diff options
Diffstat (limited to 'Build/source')
| -rw-r--r-- | Build/source/texk/web2c/mplibdir/ChangeLog | 5 | ||||
| -rw-r--r-- | Build/source/texk/web2c/mplibdir/mp.w | 24 | ||||
| -rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmath.w | 43 | ||||
| -rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmathbinary.w | 176 | ||||
| -rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmathdecimal.w | 143 | ||||
| -rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmathdouble.w | 132 |
6 files changed, 500 insertions, 23 deletions
diff --git a/Build/source/texk/web2c/mplibdir/ChangeLog b/Build/source/texk/web2c/mplibdir/ChangeLog index db853ec24eb..c3db5d08ce2 100644 --- a/Build/source/texk/web2c/mplibdir/ChangeLog +++ b/Build/source/texk/web2c/mplibdir/ChangeLog @@ -1,3 +1,8 @@ +2015-10-07 Akira Kakuto <kakuto@fuk.kindai.ac.jp> + + *mpmath.w, mpmathbinary.w, mpmathdecimal.w, mp.w, mpmathdouble.w: + Sync with the upstream trunk. + 2015-09-21 Akira Kakuto <kakuto@fuk.kindai.ac.jp> * mpmathdecimal.w: Sync with the upstream trunk. diff --git a/Build/source/texk/web2c/mplibdir/mp.w b/Build/source/texk/web2c/mplibdir/mp.w index d40a45fda8c..a5b6a619cc4 100644 --- a/Build/source/texk/web2c/mplibdir/mp.w +++ b/Build/source/texk/web2c/mplibdir/mp.w @@ -1,4 +1,4 @@ -% $Id: mp.w 2060 2015-05-07 11:01:39Z luigi $ +% $Id: mp.w 2070 2015-10-06 10:35:23Z luigi $ % % This file is part of MetaPost; % the MetaPost program is in the public domain. @@ -331,6 +331,7 @@ typedef struct mp_number_data mp_number; typedef void (*convert_func) (mp_number *r); typedef void (*m_log_func) (MP mp, mp_number *r, mp_number a); typedef void (*m_exp_func) (MP mp, mp_number *r, mp_number a); +typedef void (*m_unif_rand_func) (MP mp, mp_number *ret, mp_number x_orig); typedef void (*m_norm_rand_func) (MP mp, mp_number *ret); typedef void (*pyth_add_func) (MP mp, mp_number *r, mp_number a, mp_number b); typedef void (*pyth_sub_func) (MP mp, mp_number *r, mp_number a, mp_number b); @@ -475,6 +476,7 @@ typedef struct math_data { n_arg_func n_arg; m_log_func m_log; m_exp_func m_exp; + m_unif_rand_func m_unif_rand; m_norm_rand_func m_norm_rand; pyth_add_func pyth_add; pyth_sub_func pyth_sub; @@ -2566,9 +2568,12 @@ void mp_new_randoms (MP mp) { mp->j_random = 54; } -@ To consume a random fraction, the program below will say `|next_random|'. +@ To consume a random fraction, the program below will say `|next_random|'. +Now each number system has its own implementation, +true to the original as much as possibile. @c +/* Unused. static void mp_next_random (MP mp, mp_number *ret) { if ( mp->j_random==0 ) mp_new_randoms(mp); @@ -2576,7 +2581,7 @@ static void mp_next_random (MP mp, mp_number *ret) { decr(mp->j_random); 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. @@ -2584,10 +2589,15 @@ 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. +This is the original one, +that stays as reference: +As said before, now each number system has its own implementation. + @c +/*Unused. static void mp_unif_rand (MP mp, mp_number *ret, mp_number x_orig) { - mp_number y; /* trial value */ + mp_number y; // trial value mp_number x, abs_x; mp_number u; new_fraction (y); @@ -2612,7 +2622,7 @@ static void mp_unif_rand (MP mp, mp_number *ret, mp_number x_orig) { 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 @@ -10965,6 +10975,7 @@ This first set goes into the header @d n_arg(R,A,B) (((math_data *)(mp->math))->n_arg)(mp,&(R),A,B) @d m_log(R,A) (((math_data *)(mp->math))->m_log)(mp,&(R),A) @d m_exp(R,A) (((math_data *)(mp->math))->m_exp)(mp,&(R),A) +@d m_unif_rand(R,A) (((math_data *)(mp->math))->m_unif_rand)(mp,&(R),A) @d m_norm_rand(R) (((math_data *)(mp->math))->m_norm_rand)(mp,&(R)) @d velocity(R,A,B,C,D,E) (((math_data *)(mp->math))->velocity)(mp,&(R),A,B,C,D,E) @d ab_vs_cd(R,A,B,C,D) (((math_data *)(mp->math))->ab_vs_cd)(mp,&(R),A,B,C,D) @@ -24861,7 +24872,8 @@ static void mp_do_unary (MP mp, quarterword c) { { mp_number vvx; new_number (vvx); - mp_unif_rand (mp, &vvx, cur_exp_value_number ()); + /*mp_unif_rand (mp, &vvx, cur_exp_value_number ());*/ + m_unif_rand (vvx, cur_exp_value_number ()); set_cur_exp_value_number (vvx); free_number (vvx); } diff --git a/Build/source/texk/web2c/mplibdir/mpmath.w b/Build/source/texk/web2c/mplibdir/mpmath.w index b7d1948ec3e..b713dc14151 100644 --- a/Build/source/texk/web2c/mplibdir/mpmath.w +++ b/Build/source/texk/web2c/mplibdir/mpmath.w @@ -1,4 +1,4 @@ -% $Id: mpmath.w 2057 2015-03-20 01:33:59Z luigi $ +% $Id: mpmath.w 2070 2015-10-06 10:35:23Z luigi $ % % This file is part of MetaPost; % the MetaPost program is in the public domain. @@ -56,6 +56,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_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig); static void mp_m_norm_rand (MP mp, mp_number *ret); static void mp_m_exp (MP mp, mp_number *ret, mp_number x_orig); static void mp_m_log (MP mp, mp_number *ret, mp_number x_orig); @@ -254,6 +255,7 @@ void * mp_initialize_scaled_math (MP mp) { math->n_arg = mp_n_arg; math->m_log = mp_m_log; math->m_exp = mp_m_exp; + math->m_unif_rand = mp_m_unif_rand; math->m_norm_rand = mp_m_norm_rand; math->pyth_add = mp_pyth_add; math->pyth_sub = mp_pyth_sub; @@ -1836,6 +1838,45 @@ static void mp_next_random (MP mp, mp_number *ret) { } +@ 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_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; + 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_number_abs (&abs_x); + mp_next_random(mp, &u); + /*take_fraction (y, abs_x, u);*/ + mp_number_take_fraction (mp,&y, abs_x,u); + free_number (u); + if (mp_number_equal(y, abs_x)) { + /*set_number_to_zero(*ret);*/ + 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); + } + 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 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 diff --git a/Build/source/texk/web2c/mplibdir/mpmathdecimal.w b/Build/source/texk/web2c/mplibdir/mpmathdecimal.w index 2d5a7e3d0a0..d0550c51a80 100644 --- a/Build/source/texk/web2c/mplibdir/mpmathdecimal.w +++ b/Build/source/texk/web2c/mplibdir/mpmathdecimal.w @@ -68,6 +68,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_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig); static void mp_decimal_m_norm_rand (MP mp, mp_number *ret); static void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number x_orig); static void mp_decimal_m_log (MP mp, mp_number *ret, mp_number x_orig); @@ -542,6 +543,7 @@ void * mp_initialize_decimal_math (MP mp) { math->n_arg = mp_decimal_n_arg; math->m_log = mp_decimal_m_log; math->m_exp = mp_decimal_m_exp; + math->m_unif_rand = mp_decimal_m_unif_rand; math->m_norm_rand = mp_decimal_m_norm_rand; math->pyth_add = mp_decimal_pyth_add; math->pyth_sub = mp_decimal_pyth_sub; @@ -1424,8 +1426,8 @@ int mp_round_unscaled(mp_number x_orig) { @c void mp_number_floor (mp_number *i) { - int round = set.round; - set.round = DEC_ROUND_FLOOR; + int round = set.round; + set.round = DEC_ROUND_FLOOR; decNumberToIntegralValue(i->data.num, i->data.num, &set); set.round = round; } @@ -1724,6 +1726,81 @@ mp_number_to_double(*n_cos), mp_number_to_double(*n_sin)); mp_check_decNumber(mp, n_sin->data.num, &set); } +@ 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 @@ -1746,6 +1823,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 @@ -1754,14 +1834,21 @@ void mp_decimal_number_modulo (mp_number *a, mp_number b) { } +@ 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) { + decNumber a; + decNumber b; + unsigned long int op; + (void)mp; + op = (unsigned)ran_arr_next(); + decNumberFromInt32(&a, op); + decNumberFromInt32(&b, MM); + decNumberDivide (&a, &a, &b, &set); /* a = a/b */ + decNumberCopy(ret->data.num, &a); + mp_check_decNumber(mp, ret->data.num, &set); +} @ To consume a random fraction, the program below will say `|next_random|'. @@ -1776,6 +1863,43 @@ static void mp_next_random (MP mp, mp_number *ret) { } +@ 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_decimal_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; + 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_decimal_abs (&abs_x); + mp_next_unif_random(mp, &u); + decNumberMultiply (y.data.num, abs_x.data.num, u.data.num, &set); + 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); + } + 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 {\sl The Art of Computer Programming\/}). @@ -1810,7 +1934,6 @@ static void mp_decimal_m_norm_rand (MP mp, mp_number *ret) { mp_number_clone (&xa, r); mp_decimal_m_log (mp,&la, u); mp_set_decimal_from_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la); - /*mp_decimal_ab_vs_cd (mp,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa);*/ mp_ab_vs_cd (mp,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa); } while (mp_number_less(ab_vs_cd,((math_data *)mp->math)->zero_t)); mp_number_clone (ret, xa); diff --git a/Build/source/texk/web2c/mplibdir/mpmathdouble.w b/Build/source/texk/web2c/mplibdir/mpmathdouble.w index 7a6d99a4aa7..d5663b2df97 100644 --- a/Build/source/texk/web2c/mplibdir/mpmathdouble.w +++ b/Build/source/texk/web2c/mplibdir/mpmathdouble.w @@ -1,4 +1,4 @@ -% $Id: mpmathdouble.w 2057 2015-03-20 01:33:59Z luigi $ +% $Id: mpmathdouble.w 2070 2015-10-06 10:35:23Z luigi $ % % This file is part of MetaPost; % the MetaPost program is in the public domain. @@ -64,6 +64,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_double_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig); static void mp_double_m_norm_rand (MP mp, mp_number *ret); static void mp_double_m_exp (MP mp, mp_number *ret, mp_number x_orig); static void mp_double_m_log (MP mp, mp_number *ret, mp_number x_orig); @@ -262,6 +263,7 @@ void * mp_initialize_double_math (MP mp) { math->n_arg = mp_double_n_arg; math->m_log = mp_double_m_log; math->m_exp = mp_double_m_exp; + math->m_unif_rand = mp_double_m_unif_rand; math->m_norm_rand = mp_double_m_norm_rand; math->pyth_add = mp_double_pyth_add; math->pyth_sub = mp_double_pyth_sub; @@ -1215,6 +1217,81 @@ mp_number_to_double(*n_cos), mp_number_to_double(*n_sin)); #endif } +@ 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 @@ -1237,6 +1314,10 @@ 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 @@ -1255,6 +1336,18 @@ void mp_number_modulo (mp_number *a, mp_number b) { +@ 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) { + double a; + unsigned long int op; + (void)mp; + op = (unsigned)ran_arr_next(); + a = op/(MM*1.0); + ret->data.dval = a; +} + @ To consume a random fraction, the program below will say `|next_random|'. @@ -1269,6 +1362,43 @@ static void mp_next_random (MP mp, mp_number *ret) { } +@ 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_double_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; + 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_double_abs (&abs_x); + mp_next_unif_random(mp, &u); + y.data.dval = abs_x.data.dval * u.data.dval; + 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); + } + 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 {\sl The Art of Computer Programming\/}). |