summaryrefslogtreecommitdiff
path: root/Build/source/texk/web2c/mplibdir/mpmathbinary.w
diff options
context:
space:
mode:
authorAkira Kakuto <kakuto@fuk.kindai.ac.jp>2015-10-07 01:43:37 +0000
committerAkira Kakuto <kakuto@fuk.kindai.ac.jp>2015-10-07 01:43:37 +0000
commit4ca2692153042a0c2f7db8d681e190b867ae2faf (patch)
tree5ce3f1b1d2b21911766b40bb856a8279cec96fed /Build/source/texk/web2c/mplibdir/mpmathbinary.w
parent357d78ed0c02ea3e4c41bac6c13f2cdb7c50784c (diff)
web2c/mplibdir: Sync with the upstream.
git-svn-id: svn://tug.org/texlive/trunk@38574 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/texk/web2c/mplibdir/mpmathbinary.w')
-rw-r--r--Build/source/texk/web2c/mplibdir/mpmathbinary.w176
1 files changed, 171 insertions, 5 deletions
diff --git a/Build/source/texk/web2c/mplibdir/mpmathbinary.w b/Build/source/texk/web2c/mplibdir/mpmathbinary.w
index 38cd823499a..43edabab86c 100644
--- a/Build/source/texk/web2c/mplibdir/mpmathbinary.w
+++ b/Build/source/texk/web2c/mplibdir/mpmathbinary.w
@@ -27,6 +27,7 @@
#define MPMATHBINARY_H 1
#include "mplib.h"
#include "mpmp.h" /* internal header */
+#include <gmp.h>
#include <mpfr.h>
@<Internal library declarations@>;
#endif
@@ -61,6 +62,7 @@ static void mp_number_angle_to_scaled (mp_number *A);
static void mp_number_fraction_to_scaled (mp_number *A);
static void mp_number_scaled_to_fraction (mp_number *A);
static void mp_number_scaled_to_angle (mp_number *A);
+static void mp_binary_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig);
static void mp_binary_m_norm_rand (MP mp, mp_number *ret);
static void mp_binary_m_exp (MP mp, mp_number *ret, mp_number x_orig);
static void mp_binary_m_log (MP mp, mp_number *ret, mp_number x_orig);
@@ -396,6 +398,7 @@ void * mp_initialize_binary_math (MP mp) {
math->n_arg = mp_binary_n_arg;
math->m_log = mp_binary_m_log;
math->m_exp = mp_binary_m_exp;
+ math->m_unif_rand = mp_binary_m_unif_rand;
math->m_norm_rand = mp_binary_m_norm_rand;
math->pyth_add = mp_binary_pyth_add;
math->pyth_sub = mp_binary_pyth_sub;
@@ -870,32 +873,54 @@ The definitions below are temporarily here
@<Declarations...@>=
static void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop);
-@ Precision check is TODO
-@d too_precise(a) 0
+@ The check of the precision is based on the article "27 Bits are not enough for 8-Digit accuracy"
+@ by Bennet Goldberg which roughly says that
+@ given $p$ digits in base 10 and $q$ digits in base 2,
+@ conversion from base 10 round-trip through base 2 if and only if $10^p < $2^{q-1}$.
+@ In our case $p/\log_{10}2 + 1 < q$, or $q\geq a$
+@ where $q$ is the current precision in bits and $a=\left\lceil p/\log_{10}2 + 1\right\rceil$.
+@ Therefore if $a>q$ the number could be too precise and we emit a warning.
+@d too_precise(a) (a>precision_bits)
@c
void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop) {
int invalid = 0;
mpfr_t result;
size_t l = stop-start+1;
+ unsigned long lp, lpbit;
+ /*size_t lp = l; */
char *buf = mp_xmalloc(mp, l+1, 1);
+ char *bufp = buf;
buf[l] = '\0';
mpfr_init2(result, precision_bits);
(void)strncpy(buf,(const char *)start, l);
invalid = mpfr_set_str(result,buf, 10, ROUNDING);
//fprintf(stdout,"scan of [%s] produced %s, ", buf, mp_binnumber_tostring(result));
+ lp = (unsigned long) l;
+
+ /* strip leading - or + or 0 or .*/
+ if ( (*bufp=='-') || (*bufp=='+') || (*bufp=='0') || (*bufp=='.') ) { lp--; bufp++;}
+ /* strip also . */
+ lp = strchr(bufp,'.') ? lp-1: lp;
+ /* strip also trailing 0s */
+ bufp = buf+l-1;
+ while(*bufp == '0') {bufp--; lp--;}
+ lp = lp>0? lp: 1;
+ /* bits needed for buf */
+ lpbit = (unsigned long)ceil(lp/log10(2)+1);
free(buf);
+ bufp = NULL;
if (invalid == 0) {
set_cur_mod(result);
// fprintf(stdout,"mod=%s\n", mp_binary_number_tostring(mp,mp->cur_mod_->data.n));
- if (too_precise(l)) {
+ if (too_precise(lpbit)) {
if (mpfr_positive_p((mpfr_ptr)(internal_value (mp_warning_check).data.num)) &&
(mp->scanner_status != tex_flushing)) {
char msg[256];
const char *hlp[] = {"Continue and I'll try to cope",
- "with that big value; but it might be dangerous.",
+ "with that value; but it might be dangerous.",
"(Set warningcheck:=0 to suppress this message.)",
NULL };
- mp_snprintf (msg, 256, "Number is too large (%s)", mp_binary_number_tostring(mp,mp->cur_mod_->data.n));
+ mp_snprintf (msg, 256, "Number is too precise (%d vs. numberprecision = %f)", (unsigned int)lp,mpfr_get_d(internal_value (mp_number_precision).data.num, ROUNDING));
@.Number is too large@>;
mp_error (mp, msg, hlp, true);
}
@@ -1500,6 +1525,82 @@ void mp_binary_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_
mpfr_clear (one_eighty);
}
+@ This is the http://www-cs-faculty.stanford.edu/~uno/programs/rng.c
+with small cosmetic modifications.
+
+@c
+#define KK 100 /* the long lag */
+#define LL 37 /* the short lag */
+#define MM (1L<<30) /* the modulus */
+#define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
+/* */
+static long ran_x[KK]; /* the generator state */
+/* */
+static void ran_array(long aa[],int n) /* put n new random numbers in aa */
+ /* long aa[] destination */
+ /* int n array length (must be at least KK) */
+{
+ register int i,j;
+ for (j=0;j<KK;j++) aa[j]=ran_x[j];
+ for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
+ for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
+ for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
+}
+/* */
+/* the following routines are from exercise 3.6--15 */
+/* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */
+/* */
+#define QUALITY 1009 /* recommended quality level for high-res use */
+static long ran_arr_buf[QUALITY];
+static long ran_arr_dummy=-1, ran_arr_started=-1;
+static long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */
+/* */
+#define TT 70 /* guaranteed separation between streams */
+#define is_odd(x) ((x)&1) /* units bit of x */
+/* */
+static void ran_start(long seed) /* do this before using ran_array */
+ /* long seed selector for different streams */
+{
+ register int t,j;
+ long x[KK+KK-1]; /* the preparation buffer */
+ register long ss=(seed+2)&(MM-2);
+ for (j=0;j<KK;j++) {
+ x[j]=ss; /* bootstrap the buffer */
+ ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
+ }
+ x[1]++; /* make x[1] (and only x[1]) odd */
+ for (ss=seed&(MM-1),t=TT-1; t; ) {
+ for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
+ for (j=KK+KK-2;j>=KK;j--)
+ x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
+ x[j-KK]=mod_diff(x[j-KK],x[j]);
+ if (is_odd(ss)) { /* "multiply by z" */
+ for (j=KK;j>0;j--) x[j]=x[j-1];
+ x[0]=x[KK]; /* shift the buffer cyclically */
+ x[LL]=mod_diff(x[LL],x[KK]);
+ }
+ if (ss) ss>>=1; else t--;
+ }
+ for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
+ for (;j<KK;j++) ran_x[j-LL]=x[j];
+ for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
+ ran_arr_ptr=&ran_arr_started;
+}
+/* */
+#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
+static long ran_arr_cycle(void)
+{
+ if (ran_arr_ptr==&ran_arr_dummy)
+ ran_start(314159L); /* the user forgot to initialize */
+ ran_array(ran_arr_buf,QUALITY);
+ ran_arr_buf[KK]=-1;
+ ran_arr_ptr=ran_arr_buf+1;
+ return ran_arr_buf[0];
+}
+
+
+
+
@ To initialize the |randoms| table, we call the following routine.
@c
@@ -1522,6 +1623,9 @@ void mp_init_randoms (MP mp, int seed) {
mp_new_randoms (mp);
mp_new_randoms (mp);
mp_new_randoms (mp); /* ``warm up'' the array */
+
+ ran_start ((unsigned long)seed);
+
}
@ @c
@@ -1529,6 +1633,22 @@ void mp_binary_number_modulo (mp_number *a, mp_number b) {
mpfr_remainder (a->data.num, a->data.num, b.data.num, ROUNDING);
}
+@ To consume a random integer for the uniform generator, the program below will say `|next_unif_random|'.
+
+@c
+static void mp_next_unif_random (MP mp, mp_number *ret) {
+ mp_number rop;
+ unsigned long int op;
+ float flt_op ;
+ (void)mp;
+ mp_new_number (mp, &rop, mp_scaled_type);
+ op = (unsigned)ran_arr_next();
+ flt_op = op/(MM*1.0);
+ mpfr_set_d ((mpfr_ptr)(rop.data.num), flt_op,ROUNDING);
+ mp_number_clone (ret, rop);
+ free_number (rop);
+}
+
@ To consume a random fraction, the program below will say `|next_random|'.
@@ -1542,6 +1662,52 @@ static void mp_next_random (MP mp, mp_number *ret) {
mp_number_clone (ret, mp->randoms[mp->j_random]);
}
+@ To produce a uniform random number in the range |0<=u<x| or |0>=u>x|
+or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here.
+
+Note that the call of |take_fraction| will produce the values 0 and~|x|
+with about half the probability that it will produce any other particular
+values between 0 and~|x|, because it rounds its answers.
+
+@c
+static void mp_binary_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig) {
+ mp_number y; /* trial value */
+ mp_number x, abs_x;
+ mp_number u;
+ char *r ;mpfr_exp_t e;
+ new_fraction (y);
+ new_number (x);
+ new_number (abs_x);
+ new_number (u);
+ mp_number_clone (&x, x_orig);
+ mp_number_clone (&abs_x, x);
+ mp_binary_abs (&abs_x);
+ mp_next_unif_random(mp, &u);
+ mpfr_mul (y.data.num, abs_x.data.num, u.data.num, ROUNDING);
+ free_number (u);
+ if (mp_number_equal(y, abs_x)) {
+ mp_number_clone (ret, ((math_data *)mp->math)->zero_t);
+ } else if (mp_number_greater(x, ((math_data *)mp->math)->zero_t)) {
+ mp_number_clone (ret, y);
+ } else {
+ mp_number_clone (ret, y);
+ mp_number_negate (ret);
+ }
+ r = mpfr_get_str(NULL, // char *str,
+ &e, // mpfr_exp_t *expptr,
+ 10, // int b,
+ 0, // size_t n,
+ ret->data.num, // mpfr_t op,
+ ROUNDING // mpfr_rnd_t rnd
+ );
+ printf("\nret=%s e=%ld\n",r,e);
+ mpfr_free_str(r);
+ free_number (abs_x);
+ free_number (x);
+ free_number (y);
+}
+
+
@ Finally, a normal deviate with mean zero and unit standard deviation
can readily be obtained with the ratio method (Algorithm 3.4.1R in