summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPeter Breitenlohner <peb@mppmu.mpg.de>2015-03-20 07:58:39 +0000
committerPeter Breitenlohner <peb@mppmu.mpg.de>2015-03-20 07:58:39 +0000
commit49256ac120db7b7f06e54a931701fe8b3d2b4bdb (patch)
tree3d960eb9207936a77f4b3ef1be06dcf6c61056cd
parentcae04d40ac09c34e21a3c3efce3ddf43e7062db3 (diff)
texk/web2c/mplibdir: Sync with metapost-trunk
git-svn-id: svn://tug.org/texlive/trunk@36569 c570f23f-e606-0410-a88d-b1316a301751
-rw-r--r--Build/source/texk/web2c/mplibdir/ChangeLog8
-rw-r--r--Build/source/texk/web2c/mplibdir/lmplib.c3
-rw-r--r--Build/source/texk/web2c/mplibdir/mp.w25
-rw-r--r--Build/source/texk/web2c/mplibdir/mpmath.w67
-rw-r--r--Build/source/texk/web2c/mplibdir/mpmathbinary.w104
-rw-r--r--Build/source/texk/web2c/mplibdir/mpmathdecimal.w110
-rw-r--r--Build/source/texk/web2c/mplibdir/mpmathdouble.w96
-rw-r--r--Build/source/texk/web2c/mplibdir/mpost.w2
-rw-r--r--Build/source/texk/web2c/mplibdir/mpstrings.w2
-rw-r--r--Build/source/texk/web2c/mplibdir/mpxout.w1
-rw-r--r--Build/source/texk/web2c/mplibdir/pngout.w2
-rw-r--r--Build/source/texk/web2c/mplibdir/psout.w1
-rw-r--r--Build/source/texk/web2c/mplibdir/svgout.w2
-rw-r--r--Build/source/texk/web2c/mplibdir/tfmin.w2
14 files changed, 390 insertions, 35 deletions
diff --git a/Build/source/texk/web2c/mplibdir/ChangeLog b/Build/source/texk/web2c/mplibdir/ChangeLog
index 17480cb9de1..4cba35f7910 100644
--- a/Build/source/texk/web2c/mplibdir/ChangeLog
+++ b/Build/source/texk/web2c/mplibdir/ChangeLog
@@ -1,12 +1,18 @@
+2015-03-19 Luigi Scarso <luigi.scarso@gmail.com>
+ Fixed normal rand. generator for double/decimal/binary
+ systems.
+
2015-01-20 Luigi Scarso <luigi.scarso@gmail.com>
Added cur_sym_mod() == mp_normal_sym
to the cur_sym==0 condition mp_cur_tok
+
2015-01-19 Luigi Scarso <luigi.scarso@gmail.com>
fixed a bug in mp_recycle_value with NULL node
+
2015-01-15 Luigi Scarso <luigi.scarso@gmail.com>
- * mpmathdecimal.w: better management for static data in mpmathdecimal.w
+ mpmathdecimal.w: better management for static data in mpmathdecimal.w
(thanks to Philipp Gesang <phg@phi-gamma.net>)
2014-12-04 Luigi Scarso <luigi.scarso@gmail.com>
diff --git a/Build/source/texk/web2c/mplibdir/lmplib.c b/Build/source/texk/web2c/mplibdir/lmplib.c
index 6a3f83e8775..ab1b707b237 100644
--- a/Build/source/texk/web2c/mplibdir/lmplib.c
+++ b/Build/source/texk/web2c/mplibdir/lmplib.c
@@ -52,9 +52,6 @@
#include "mplibsvg.h"
#include "mplibpng.h"
- /*@unused@*/ static const char _svn_version[] =
- "$Id: lmplib.c 2036 2014-06-26 18:44:57Z luigi $";
-
int luaopen_mplib(lua_State * L); /* forward */
/* metatable identifiers and tests */
diff --git a/Build/source/texk/web2c/mplibdir/mp.w b/Build/source/texk/web2c/mplibdir/mp.w
index 4ba9f3586a9..9735f6f33b2 100644
--- a/Build/source/texk/web2c/mplibdir/mp.w
+++ b/Build/source/texk/web2c/mplibdir/mp.w
@@ -1,4 +1,4 @@
-% $Id: mp.w 2055 2015-01-22 15:39:12Z luigi $
+% $Id: mp.w 2057 2015-03-20 01:33:59Z 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_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);
typedef void (*n_arg_func) (MP mp, mp_number *r, mp_number a, mp_number b);
@@ -471,9 +472,10 @@ typedef struct math_data {
velocity_func velocity;
ab_vs_cd_func ab_vs_cd;
crossing_point_func crossing_point;
- n_arg_func n_arg;
- m_log_func m_log;
- m_exp_func m_exp;
+ n_arg_func n_arg;
+ m_log_func m_log;
+ m_exp_func m_exp;
+ m_norm_rand_func m_norm_rand;
pyth_add_func pyth_add;
pyth_sub_func pyth_sub;
fraction_to_round_scaled_func fraction_to_round_scaled;
@@ -2603,9 +2605,14 @@ static void mp_unif_rand (MP mp, mp_number *ret, mp_number x_orig) {
@ 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\/}).
+{\sl The Art of Computer Programming\/}). This is the original one,
+that stays as reference:
+Now each number system has its own implementation,
+true to the original as much as possibile.
+.
@c
+/* Unused.
static void mp_norm_rand (MP mp, mp_number *ret) {
mp_number ab_vs_cd;
mp_number abs_x;
@@ -2623,7 +2630,7 @@ static void mp_norm_rand (MP mp, mp_number *ret) {
mp_number v;
new_number (v);
mp_next_random(mp, &v);
- number_substract (v, fraction_half_t);
+ number_substract (v, fraction_half_t);
take_fraction (xa, sqrt_8_e_k, v);
free_number (v);
mp_next_random(mp, &u);
@@ -2644,6 +2651,7 @@ static void mp_norm_rand (MP mp, mp_number *ret) {
free_number (xa);
free_number (u);
}
+*/
@* Packed data.
@@ -10940,6 +10948,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_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)
@d crossing_point(R,A,B,C) (((math_data *)(mp->math))->crossing_point)(mp,&(R),A,B,C)
@@ -24536,7 +24545,8 @@ static void mp_do_nullary (MP mp, quarterword c) {
{
mp_number r;
new_number (r);
- mp_norm_rand (mp, &r);
+ /*mp_norm_rand (mp, &r);*/
+ m_norm_rand (r);
mp->cur_exp.type = mp_known;
set_cur_exp_value_number (r);
free_number (r);
@@ -30002,7 +30012,6 @@ Here's one of the simplest:
@ @<Declare action procedures for use by |do_statement|@>=
static void mp_do_random_seed (MP mp);
-
@ @c
void mp_do_random_seed (MP mp) {
mp_value new_expr;
diff --git a/Build/source/texk/web2c/mplibdir/mpmath.w b/Build/source/texk/web2c/mplibdir/mpmath.w
index 82d41746c72..b7d1948ec3e 100644
--- a/Build/source/texk/web2c/mplibdir/mpmath.w
+++ b/Build/source/texk/web2c/mplibdir/mpmath.w
@@ -1,4 +1,4 @@
-% $Id: mpmath.w 2037 2014-09-02 14:59:07Z luigi $
+% $Id: mpmath.w 2057 2015-03-20 01:33:59Z 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_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);
static void mp_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b);
@@ -253,6 +254,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_norm_rand = mp_m_norm_rand;
math->pyth_add = mp_pyth_add;
math->pyth_sub = mp_pyth_sub;
math->fraction_to_scaled = mp_number_fraction_to_scaled;
@@ -1816,4 +1818,65 @@ char * mp_number_tostring (MP mp, mp_number n) {
@ @c
void mp_number_modulo (mp_number *a, mp_number b) {
a->data.val = a->data.val % b.data.val;
-} \ No newline at end of file
+}
+
+
+
+
+
+@ To consume a random fraction, the program below will say `|next_random|'.
+
+@c
+static void mp_next_random (MP mp, mp_number *ret) {
+ if ( mp->j_random==0 )
+ mp_new_randoms(mp);
+ else
+ mp->j_random = mp->j_random-1;
+ mp_number_clone (ret, mp->randoms[mp->j_random]);
+}
+
+
+
+@ 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\/}).
+
+@c
+static void mp_m_norm_rand (MP mp, mp_number *ret) {
+ mp_number ab_vs_cd;
+ mp_number abs_x;
+ mp_number u;
+ mp_number r;
+ mp_number la, xa;
+ new_number (ab_vs_cd);
+ new_number (la);
+ new_number (xa);
+ new_number (abs_x);
+ new_number (u);
+ new_number (r);
+ do {
+ do {
+ mp_number v;
+ new_number (v);
+ mp_next_random(mp, &v);
+ mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t);
+ mp_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
+ free_number (v);
+ mp_next_random(mp, &u);
+ mp_number_clone (&abs_x, xa);
+ mp_number_abs (&abs_x);
+ } while (!mp_number_less(abs_x, u));
+ mp_number_make_fraction (mp, &r, xa, u);
+ mp_number_clone (&xa, r);
+ mp_m_log (mp,&la, u);
+ mp_set_number_from_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la);
+ 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);
+ free_number (ab_vs_cd);
+ free_number (r);
+ free_number (abs_x);
+ free_number (la);
+ free_number (xa);
+ free_number (u);
+}
diff --git a/Build/source/texk/web2c/mplibdir/mpmathbinary.w b/Build/source/texk/web2c/mplibdir/mpmathbinary.w
index ca223c6e353..b8b8508d8aa 100644
--- a/Build/source/texk/web2c/mplibdir/mpmathbinary.w
+++ b/Build/source/texk/web2c/mplibdir/mpmathbinary.w
@@ -47,6 +47,7 @@ First, here are some very important constants.
#define DEBUG 0
static void mp_binary_scan_fractional_token (MP mp, int n);
static void mp_binary_scan_numeric_token (MP mp, int n);
+static void mp_binary_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
static void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
static void mp_binary_crossing_point (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c);
static void mp_binary_number_modulo (mp_number *a, mp_number b);
@@ -60,6 +61,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_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);
static void mp_binary_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b);
@@ -298,7 +300,7 @@ void * mp_initialize_binary_math (MP mp) {
mpfr_set_si(math->one_eighty_deg_t.data.num, 180 * angle_multiplier, ROUNDING);
/* various approximations */
mp_new_number (mp, &math->one_k, mp_scaled_type);
- mpfr_set_si(math->one_k.data.num, 1024, ROUNDING);
+ mpfr_set_d(math->one_k.data.num, 1.0/64, ROUNDING);
mp_new_number (mp, &math->sqrt_8_e_k, mp_scaled_type);
{
mpfr_set_d(math->sqrt_8_e_k.data.num, 112428.82793 / 65536.0, ROUNDING);
@@ -394,6 +396,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_norm_rand = mp_binary_m_norm_rand;
math->pyth_add = mp_binary_pyth_add;
math->pyth_sub = mp_binary_pyth_sub;
math->fraction_to_scaled = mp_number_fraction_to_scaled;
@@ -1524,4 +1527,101 @@ void mp_init_randoms (MP mp, int seed) {
@ @c
void mp_binary_number_modulo (mp_number *a, mp_number b) {
mpfr_remainder (a->data.num, a->data.num, b.data.num, ROUNDING);
-} \ No newline at end of file
+}
+
+
+
+@ To consume a random fraction, the program below will say `|next_random|'.
+
+@c
+static void mp_next_random (MP mp, mp_number *ret) {
+ if ( mp->j_random==0 )
+ mp_new_randoms(mp);
+ else
+ mp->j_random = mp->j_random-1;
+ mp_number_clone (ret, mp->randoms[mp->j_random]);
+}
+
+
+@ 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\/}).
+
+@c
+static void mp_binary_m_norm_rand (MP mp, mp_number *ret) {
+ mp_number ab_vs_cd;
+ mp_number abs_x;
+ mp_number u;
+ mp_number r;
+ mp_number la, xa;
+ new_number (ab_vs_cd);
+ new_number (la);
+ new_number (xa);
+ new_number (abs_x);
+ new_number (u);
+ new_number (r);
+
+ do {
+ do {
+ mp_number v;
+ new_number (v);
+ mp_next_random(mp, &v);
+ mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t);
+ mp_binary_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
+ free_number (v);
+ mp_next_random(mp, &u);
+ mp_number_clone (&abs_x, xa);
+ mp_binary_abs (&abs_x);
+ } while (!mp_number_less(abs_x, u));
+ mp_binary_number_make_fraction (mp, &r, xa, u);
+ mp_number_clone (&xa, r);
+ mp_binary_m_log (mp,&la, u);
+ mp_set_binary_from_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la);
+ mp_binary_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);
+ free_number (ab_vs_cd);
+ free_number (r);
+ free_number (abs_x);
+ free_number (la);
+ free_number (xa);
+ free_number (u);
+}
+
+
+
+@ The following subroutine is used only in norm_rand and tests if $ab$ is
+greater than, equal to, or less than~$cd$.
+The result is $+1$, 0, or~$-1$ in the three respective cases.
+
+@c
+static void mp_binary_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) {
+ mpfr_t a, b, c, d;
+ mpfr_t ab, cd;
+
+ int cmp = 0;
+ (void)mp;
+ mpfr_inits2(precision_bits, a,b,c,d,ab,cd,(mpfr_ptr)0);
+ mpfr_set(a, (mpfr_ptr )a_orig.data.num, ROUNDING);
+ mpfr_set(b, (mpfr_ptr )b_orig.data.num, ROUNDING);
+ mpfr_set(c, (mpfr_ptr )c_orig.data.num, ROUNDING);
+ mpfr_set(d, (mpfr_ptr )d_orig.data.num, ROUNDING);
+
+ mpfr_mul(ab,a,b, ROUNDING);
+ mpfr_mul(cd,c,d, ROUNDING);
+
+ mpfr_set(ret->data.num, zero, ROUNDING);
+ cmp = mpfr_cmp(ab,cd);
+ if (cmp) {
+ if (cmp>0)
+ mpfr_set(ret->data.num, one, ROUNDING);
+ else
+ mpfr_set(ret->data.num, minusone, ROUNDING);
+ }
+ mp_check_mpfr_t(mp, ret->data.num);
+ mpfr_clears(a,b,c,d,ab,cd,(mpfr_ptr)0);
+ return;
+}
+
+
diff --git a/Build/source/texk/web2c/mplibdir/mpmathdecimal.w b/Build/source/texk/web2c/mplibdir/mpmathdecimal.w
index c65e2a5a420..3ae1cce3042 100644
--- a/Build/source/texk/web2c/mplibdir/mpmathdecimal.w
+++ b/Build/source/texk/web2c/mplibdir/mpmathdecimal.w
@@ -55,6 +55,7 @@ First, here are some very important constants.
static void mp_decimal_scan_fractional_token (MP mp, int n);
static void mp_decimal_scan_numeric_token (MP mp, int n);
static void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
+/*static void mp_decimal_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);*/
static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c);
static void mp_decimal_number_modulo (mp_number *a, mp_number b);
static void mp_decimal_print_number (MP mp, mp_number n);
@@ -67,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_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);
static void mp_decimal_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b);
@@ -444,7 +446,7 @@ void * mp_initialize_decimal_math (MP mp) {
decNumberFromInt32(math->one_eighty_deg_t.data.num, 180 * angle_multiplier);
/* various approximations */
mp_new_number (mp, &math->one_k, mp_scaled_type);
- decNumberFromInt32(math->one_k.data.num, 1024);
+ decNumberFromDouble(math->one_k.data.num, 1.0/64);
mp_new_number (mp, &math->sqrt_8_e_k, mp_scaled_type);
{
decNumberFromDouble(math->sqrt_8_e_k.data.num, 112428.82793 / 65536.0);
@@ -540,6 +542,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_norm_rand = mp_decimal_m_norm_rand;
math->pyth_add = mp_decimal_pyth_add;
math->pyth_sub = mp_decimal_pyth_sub;
math->fraction_to_scaled = mp_number_fraction_to_scaled;
@@ -1748,4 +1751,107 @@ void mp_init_randoms (MP mp, int seed) {
@ @c
void mp_decimal_number_modulo (mp_number *a, mp_number b) {
decNumberRemainder(a->data.num, a->data.num, b.data.num, &set);
-} \ No newline at end of file
+}
+
+
+
+
+
+
+
+
+
+
+
+
+@ To consume a random fraction, the program below will say `|next_random|'.
+
+@c
+static void mp_next_random (MP mp, mp_number *ret) {
+ if ( mp->j_random==0 )
+ mp_new_randoms(mp);
+ else
+ mp->j_random = mp->j_random-1;
+ mp_number_clone (ret, mp->randoms[mp->j_random]);
+}
+
+
+@ 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\/}).
+
+@c
+static void mp_decimal_m_norm_rand (MP mp, mp_number *ret) {
+ mp_number ab_vs_cd;
+ mp_number abs_x;
+ mp_number u;
+ mp_number r;
+ mp_number la, xa;
+ new_number (ab_vs_cd);
+ new_number (la);
+ new_number (xa);
+ new_number (abs_x);
+ new_number (u);
+ new_number (r);
+
+ do {
+ do {
+ mp_number v;
+ new_number (v);
+ mp_next_random(mp, &v);
+ mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t);
+ mp_decimal_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
+ free_number (v);
+ mp_next_random(mp, &u);
+ mp_number_clone (&abs_x, xa);
+ mp_decimal_abs (&abs_x);
+ } while (!mp_number_less(abs_x, u));
+ mp_decimal_number_make_fraction (mp, &r, xa, u);
+ 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);
+ free_number (ab_vs_cd);
+ free_number (r);
+ free_number (abs_x);
+ free_number (la);
+ free_number (xa);
+ free_number (u);
+}
+
+
+
+
+@ The following subroutine could be used in norm_rand and tests if $ab$ is
+greater than, equal to, or less than~$cd$.
+The result is $+1$, 0, or~$-1$ in the three respective cases.
+This is not necessary, even if it's shorter than the current ab_vs_cd
+and looks as a native implememtation.
+
+@c
+/*
+void mp_decimal_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) {
+ decNumber a, b, c, d;
+ decNumber ab, cd;
+ (void)mp;
+
+ decNumberCopy(&a, (decNumber *)a_orig.data.num);
+ decNumberCopy(&b, (decNumber *)b_orig.data.num);
+ decNumberCopy(&c, (decNumber *)c_orig.data.num);
+ decNumberCopy(&d, (decNumber *)d_orig.data.num);
+
+
+ decNumberMultiply (&ab, (decNumber *)a_orig.data.num, (decNumber *)b_orig.data.num, &set);
+ decNumberMultiply (&cd, (decNumber *)c_orig.data.num, (decNumber *)d_orig.data.num, &set);
+ decNumberCompare(ret->data.num, &ab, &cd, &set);
+ mp_check_decNumber(mp, ret->data.num, &set);
+ return;
+
+}
+*/
+
+
+
diff --git a/Build/source/texk/web2c/mplibdir/mpmathdouble.w b/Build/source/texk/web2c/mplibdir/mpmathdouble.w
index a6203d0e601..7a6d99a4aa7 100644
--- a/Build/source/texk/web2c/mplibdir/mpmathdouble.w
+++ b/Build/source/texk/web2c/mplibdir/mpmathdouble.w
@@ -1,4 +1,4 @@
-% $Id: mpmathdouble.w 2037 2014-09-02 14:59:07Z luigi $
+% $Id: mpmathdouble.w 2057 2015-03-20 01:33:59Z luigi $
%
% This file is part of MetaPost;
% the MetaPost program is in the public domain.
@@ -51,6 +51,7 @@ First, here are some very important constants.
static void mp_double_scan_fractional_token (MP mp, int n);
static void mp_double_scan_numeric_token (MP mp, int n);
static void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
+static void mp_double_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
static void mp_double_crossing_point (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c);
static void mp_number_modulo (mp_number *a, mp_number b);
static void mp_double_print_number (MP mp, mp_number n);
@@ -63,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_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);
static void mp_double_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b);
@@ -182,11 +184,11 @@ void * mp_initialize_double_math (MP mp) {
math->one_eighty_deg_t.data.dval = one_eighty_deg;
/* various approximations */
mp_new_number (mp, &math->one_k, mp_scaled_type);
- math->one_k.data.dval = 1024;
+ math->one_k.data.dval = 1.0/64 ;
mp_new_number (mp, &math->sqrt_8_e_k, mp_scaled_type);
- math->sqrt_8_e_k.data.dval = 112429 / 65536.0; /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
+ math->sqrt_8_e_k.data.dval = 1.71552776992141359295 ; /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
mp_new_number (mp, &math->twelve_ln_2_k, mp_fraction_type);
- math->twelve_ln_2_k.data.dval = 139548960 / 65536.0; /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
+ math->twelve_ln_2_k.data.dval = 8.31776616671934371292 *256; /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
mp_new_number (mp, &math->coef_bound_k, mp_fraction_type);
math->coef_bound_k.data.dval = coef_bound;
mp_new_number (mp, &math->coef_bound_minus_1, mp_fraction_type);
@@ -260,6 +262,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_norm_rand = mp_double_m_norm_rand;
math->pyth_add = mp_double_pyth_add;
math->pyth_sub = mp_double_pyth_sub;
math->fraction_to_scaled = mp_number_fraction_to_scaled;
@@ -1248,4 +1251,87 @@ double modulus(double left, double right) {
}
void mp_number_modulo (mp_number *a, mp_number b) {
a->data.dval = modulus (a->data.dval, b.data.dval);
-} \ No newline at end of file
+}
+
+
+
+
+
+@ To consume a random fraction, the program below will say `|next_random|'.
+
+@c
+static void mp_next_random (MP mp, mp_number *ret) {
+ if ( mp->j_random==0 )
+ mp_new_randoms(mp);
+ else
+ mp->j_random = mp->j_random-1;
+ mp_number_clone (ret, mp->randoms[mp->j_random]);
+}
+
+
+@ 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\/}).
+
+@c
+static void mp_double_m_norm_rand (MP mp, mp_number *ret) {
+ mp_number ab_vs_cd;
+ mp_number abs_x;
+ mp_number u;
+ mp_number r;
+ mp_number la, xa;
+ new_number (ab_vs_cd);
+ new_number (la);
+ new_number (xa);
+ new_number (abs_x);
+ new_number (u);
+ new_number (r);
+
+ do {
+ do {
+ mp_number v;
+ new_number (v);
+ mp_next_random(mp, &v);
+ mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t);
+ mp_double_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
+ free_number (v);
+ mp_next_random(mp, &u);
+ mp_number_clone (&abs_x, xa);
+ mp_double_abs (&abs_x);
+ } while (!mp_number_less(abs_x, u));
+ mp_double_number_make_fraction (mp, &r, xa, u);
+ mp_number_clone (&xa, r);
+ mp_double_m_log (mp,&la, u);
+ mp_set_double_from_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la);
+ mp_double_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);
+ free_number (ab_vs_cd);
+ free_number (r);
+ free_number (abs_x);
+ free_number (la);
+ free_number (xa);
+ free_number (u);
+}
+
+
+
+
+@ The following subroutine is used only in norm_rand and tests if $ab$ is
+greater than, equal to, or less than~$cd$.
+The result is $+1$, 0, or~$-1$ in the three respective cases.
+
+@c
+void mp_double_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) {
+ double ab, cd;
+ (void)mp;
+ ret->data.dval = 0 ;
+ ab = a_orig.data.dval*b_orig.data.dval;
+ cd = c_orig.data.dval*d_orig.data.dval;
+ if (ab > cd )
+ ret->data.dval = 1 ;
+ else if (ab < cd )
+ ret->data.dval = -1 ;
+ return ;
+}
+
diff --git a/Build/source/texk/web2c/mplibdir/mpost.w b/Build/source/texk/web2c/mplibdir/mpost.w
index 7c9e24b6053..cf87cceab7f 100644
--- a/Build/source/texk/web2c/mplibdir/mpost.w
+++ b/Build/source/texk/web2c/mplibdir/mpost.w
@@ -1,5 +1,3 @@
-% $Id: mpost.w 2048 2014-12-04 10:26:26Z luigi $
-%
% This file is part of MetaPost;
% the MetaPost program is in the public domain.
% See the <Show version...> code below for more info.
diff --git a/Build/source/texk/web2c/mplibdir/mpstrings.w b/Build/source/texk/web2c/mplibdir/mpstrings.w
index 0ecc4e6ff7e..ff4e5b6b44e 100644
--- a/Build/source/texk/web2c/mplibdir/mpstrings.w
+++ b/Build/source/texk/web2c/mplibdir/mpstrings.w
@@ -1,5 +1,3 @@
-% $Id: mpstrings.w 1702 2012-09-27 10:31:30Z taco $
-%
% This file is part of MetaPost;
% the MetaPost program is in the public domain.
% See the <Show version...> code in mpost.w for more info.
diff --git a/Build/source/texk/web2c/mplibdir/mpxout.w b/Build/source/texk/web2c/mplibdir/mpxout.w
index 382b907dbee..8d1dee5ed49 100644
--- a/Build/source/texk/web2c/mplibdir/mpxout.w
+++ b/Build/source/texk/web2c/mplibdir/mpxout.w
@@ -1,4 +1,3 @@
-% $Id: mpxout.w 2031 2014-06-24 13:38:58Z luigi $
% This file is part of MetaPost;
% the MetaPost program is in the public domain.
% See the <Show version...> code in mpost.w for more info.
diff --git a/Build/source/texk/web2c/mplibdir/pngout.w b/Build/source/texk/web2c/mplibdir/pngout.w
index 1ab58d6b015..7a4a0b44006 100644
--- a/Build/source/texk/web2c/mplibdir/pngout.w
+++ b/Build/source/texk/web2c/mplibdir/pngout.w
@@ -1,5 +1,3 @@
-% $Id: pngout.w 1706 2012-09-28 08:12:34Z taco $
-%
% Copyright 2012 Taco Hoekwater.
%
% This program is free software: you can redistribute it and/or modify
diff --git a/Build/source/texk/web2c/mplibdir/psout.w b/Build/source/texk/web2c/mplibdir/psout.w
index bc26e0dadcc..e11dca0cc2a 100644
--- a/Build/source/texk/web2c/mplibdir/psout.w
+++ b/Build/source/texk/web2c/mplibdir/psout.w
@@ -1,4 +1,3 @@
-% $Id: psout.w 2045 2014-11-05 11:33:13Z luigi $
% This file is part of MetaPost;
% the MetaPost program is in the public domain.
% See the <Show version...> code in mpost.w for more info.
diff --git a/Build/source/texk/web2c/mplibdir/svgout.w b/Build/source/texk/web2c/mplibdir/svgout.w
index ecc81972124..6ee2c711739 100644
--- a/Build/source/texk/web2c/mplibdir/svgout.w
+++ b/Build/source/texk/web2c/mplibdir/svgout.w
@@ -1,5 +1,3 @@
-% $Id: svgout.w 2003 2014-04-09 09:48:26Z taco $
-%
% Copyright 2008-2009 Taco Hoekwater.
%
% This program is free software: you can redistribute it and/or modify
diff --git a/Build/source/texk/web2c/mplibdir/tfmin.w b/Build/source/texk/web2c/mplibdir/tfmin.w
index 0c9e7c5385c..142824258c6 100644
--- a/Build/source/texk/web2c/mplibdir/tfmin.w
+++ b/Build/source/texk/web2c/mplibdir/tfmin.w
@@ -1,5 +1,3 @@
-% $Id: tfmin.w 1682 2011-05-30 07:24:45Z taco $
-%
% Copyright 2008-2009 Taco Hoekwater.
%
% This program is free software: you can redistribute it and/or modify