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