web2c/mplibdir: Sync with the upstream.

git-svn-id: svn://tug.org/texlive/trunk@38574 c570f23f-e606-0410-a88d-b1316a301751

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\/}).