diff options
Diffstat (limited to 'Build/source/texk/web2c/mplibdir/mpmath.w')
-rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmath.w | 67 |
1 files changed, 65 insertions, 2 deletions
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); +} |