diff options
Diffstat (limited to 'source/luametatex/source/mp/mpc/mpmath.c')
| -rw-r--r-- | source/luametatex/source/mp/mpc/mpmath.c | 1501 |
1 files changed, 1501 insertions, 0 deletions
diff --git a/source/luametatex/source/mp/mpc/mpmath.c b/source/luametatex/source/mp/mpc/mpmath.c new file mode 100644 index 000000000..42a596dac --- /dev/null +++ b/source/luametatex/source/mp/mpc/mpmath.c @@ -0,0 +1,1501 @@ +/* This file is generated by "mtxrun --script "mtx-wtoc.lua" from the metapost cweb files. */ + + +# include "mpconfig.h" +# include "mpmath.h" +# include "mpstrings.h" + +# define coef_bound 04525252525 +# define fraction_threshold 2685 +# define half_fraction_threshold 1342 +# define scaled_threshold 8 +# define half_scaled_threshold 4 +# define near_zero_angle 26844 +# define p_over_v_threshold 0x80000 +# define equation_threshold 64 +# define unity 0x10000 +# define two (2*unity) +# define three (3*unity) +# define half_unit (unity/2) +# define three_quarter_unit (3*(unity/4)) +# define EL_GORDO 0x7fffffff +# define negative_EL_GORDO (-EL_GORDO) +# define one_third_EL_GORDO 05252525252 +# define TWEXP31 2147483648.0 +# define TWEXP28 268435456.0 +# define TWEXP16 65536.0 +# define TWEXP_16 (1.0/65536.0) +# define TWEXP_28 (1.0/268435456.0) +# define no_crossing (fraction_one + 1) +# define one_crossing fraction_one +# define zero_crossing 0 +# define fraction_half 0x08000000 +# define fraction_one 0x10000000 +# define fraction_two 0x20000000 +# define fraction_three 0x30000000 +# define fraction_four 0x40000000 +# define negate_x 1 +# define negate_y 2 +# define switch_x_and_y 4 +# define first_octant 1 +# define second_octant (first_octant + switch_x_and_y) +# define third_octant (first_octant + switch_x_and_y + negate_x) +# define fourth_octant (first_octant + negate_x) +# define fifth_octant (first_octant + negate_x + negate_y) +# define sixth_octant (first_octant + switch_x_and_y + negate_x + negate_y) +# define seventh_octant (first_octant + switch_x_and_y + negate_y) +# define eighth_octant (first_octant + negate_y) +# define forty_five_deg 0x02D00000 +# define ninety_deg 0x05A00000 +# define one_eighty_deg 0x0B400000 +# define negative_one_eighty_deg -0x0B400000 +# define three_sixty_deg 0x16800000 +# define odd(A) (abs(A)%2==1) +# define two_to_the(A) (1<<(unsigned)(A)) +# define set_cur_cmd(A) mp->cur_mod_->command = (A) +# define set_cur_mod(A) mp->cur_mod_->data.n.data.val = (A) + +static int mp_ab_vs_cd (mp_number *a, mp_number *b, mp_number *c, mp_number *d); +static void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *B); +static void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *B); +static void mp_allocate_double (MP mp, mp_number *n, double v); +static void mp_allocate_number (MP mp, mp_number *n, mp_number_type t); +static void mp_crossing_point (MP mp, mp_number *ret, mp_number *a, mp_number *b, mp_number *c); +static void mp_fraction_to_round_scaled (mp_number *x); +static void mp_free_number (MP mp, mp_number *n); +static void mp_free_scaled_math (MP mp); +static void mp_init_randoms (MP mp, int seed); +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_m_norm_rand (MP mp, mp_number *ret); +static void mp_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig); +static int mp_make_scaled (MP mp, int p, int q); +static void mp_n_arg (MP mp, mp_number *ret, mp_number *x, mp_number *y); +static void mp_n_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin); +static void mp_number_abs (mp_number *A); +static void mp_number_abs_clone (mp_number *A, mp_number *B); +static void mp_number_add (mp_number *A, mp_number *B); +static void mp_number_add_scaled (mp_number *A, int B); +static void mp_number_angle_to_scaled (mp_number *A); +static void mp_number_clone (mp_number *A, mp_number *B); +static void mp_number_divide_int (mp_number *A, int B); +static void mp_number_double (mp_number *A); +static int mp_number_equal (mp_number *A, mp_number *B); +static void mp_number_floor (mp_number *i); +static void mp_number_fraction_to_scaled (mp_number *A); +static int mp_number_greater (mp_number *A, mp_number *B); +static void mp_number_half (mp_number *A); +static int mp_number_less (mp_number *A, mp_number *B); +static void mp_number_make_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_number_make_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_number_modulo (mp_number *a, mp_number *b); +static void mp_number_multiply_int (mp_number *A, int B); +static void mp_number_negate (mp_number *A); +static void mp_number_negated_clone (mp_number *A, mp_number *B); +static int mp_number_nonequalabs (mp_number *A, mp_number *B); +static int mp_number_odd (mp_number *A); +static void mp_number_scaled_to_angle (mp_number *A); +static void mp_number_scaled_to_fraction (mp_number *A); +static void mp_number_subtract (mp_number *A, mp_number *B); +static void mp_number_swap (mp_number *A, mp_number *B); +static void mp_number_take_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_number_take_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); +static int mp_number_to_boolean (mp_number *A); +static double mp_number_to_double (mp_number *A); +static int mp_number_to_int (mp_number *A); +static int mp_number_to_scaled (mp_number *A); +static void mp_power_of (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_print_number (MP mp, mp_number *n); +static void mp_pyth_add (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_pyth_sub (MP mp, mp_number *r, mp_number *a, mp_number *b); +static int mp_round_decimals (MP mp, unsigned char *b, int k); +static int mp_round_unscaled (mp_number *x_orig); +static void mp_scaled_set_precision (MP mp); +static void mp_scan_fractional_token (MP mp, int n); +static void mp_scan_numeric_token (MP mp, int n); +static void mp_set_number_from_addition (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_number_from_boolean (mp_number *A, int B); +static void mp_set_number_from_div (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_number_from_double (mp_number *A, double B); +static void mp_set_number_from_int (mp_number *A, int B); +static void mp_set_number_from_int_div (mp_number *A, mp_number *B, int C); +static void mp_set_number_from_int_mul (mp_number *A, mp_number *B, int C); +static void mp_set_number_from_mul (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_number_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C); +static void mp_set_number_from_scaled (mp_number *A, int B); +static void mp_set_number_from_subtraction (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_number_half_from_addition (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_number_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C); +static void mp_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig); +static void mp_square_rt (MP mp, mp_number *ret, mp_number *x_orig); +static int mp_take_fraction (MP mp, int q, int f); +static int mp_take_scaled (MP mp, int q, int f); +static void mp_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t); +static void mp_wrapup_numeric_token (MP mp, int n, int f); +static char *mp_number_tostring (MP mp, mp_number *n); +static char *mp_string_scaled (MP mp, int s); +static const int mp_m_spec_log[29] = { + 0, 93032640, 38612034, 17922280, 8662214, 4261238, 2113709, 1052693, 525315, + 262400, 131136, 65552, 32772, 16385, 8192, 4096, 2048, 1024, 512, 256, 128, + 64, 32, 16, 8, 4, 2, 1, 1 +}; +static const int mp_m_spec_atan[27] = { + 0, 27855475, 14718068, 7471121, 3750058, 1876857, 938658, 469357, 234682, + 117342, 58671, 29335, 14668, 7334, 3667, 1833, 917, 458, 229, 115, 57, 29, + 14, 7, 4, 2, 1 +}; + +math_data *mp_initialize_scaled_math(MP mp) +{ + math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data)); + math->md_allocate = mp_allocate_number; + math->md_free = mp_free_number; + math->md_allocate_clone = mp_allocate_clone; + math->md_allocate_abs = mp_allocate_abs; + math->md_allocate_double = mp_allocate_double; + mp_allocate_number(mp, &math->md_precision_default, mp_scaled_type); + mp_allocate_number(mp, &math->md_precision_max, mp_scaled_type); + mp_allocate_number(mp, &math->md_precision_min, mp_scaled_type); + mp_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_inf_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_unity_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_two_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_three_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_zero_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type); + mp_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type); + mp_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type); + mp_allocate_number(mp, &math->md_one_k, mp_scaled_type); + mp_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type); + mp_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type); + mp_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type); + mp_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type); + mp_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type); + mp_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type); + mp_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type); + mp_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type); + mp_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type); + math->md_precision_default.data.val = unity * 10; + math->md_precision_max.data.val = unity * 10; + math->md_precision_min.data.val = unity * 10; + math->md_epsilon_t.data.val = 1; + math->md_inf_t.data.val = EL_GORDO; + math->md_negative_inf_t.data.val = negative_EL_GORDO; + math->md_one_third_inf_t.data.val = one_third_EL_GORDO; + math->md_warning_limit_t.data.val = fraction_one; + math->md_unity_t.data.val = unity; + math->md_two_t.data.val = two; + math->md_three_t.data.val = three; + math->md_half_unit_t.data.val = half_unit; + math->md_three_quarter_unit_t.data.val = three_quarter_unit; + math->md_arc_tol_k.data.val = (unity/4096); + math->md_fraction_one_t.data.val = fraction_one; + math->md_fraction_half_t.data.val = fraction_half; + math->md_fraction_three_t.data.val = fraction_three; + math->md_fraction_four_t.data.val = fraction_four; + math->md_three_sixty_deg_t.data.val = three_sixty_deg; + math->md_one_eighty_deg_t.data.val = one_eighty_deg; + math->md_negative_one_eighty_deg_t.data.val = negative_one_eighty_deg; + math->md_one_k.data.val = 1024; + math->md_sqrt_8_e_k.data.val = 112429; + math->md_twelve_ln_2_k.data.val = 139548960; + math->md_coef_bound_k.data.val = coef_bound; + math->md_coef_bound_minus_1.data.val = coef_bound - 1; + math->md_twelvebits_3.data.val = 1365; + math->md_twentysixbits_sqrt2_t.data.val = 94906266; + math->md_twentyeightbits_d_t.data.val = 35596755; + math->md_twentysevenbits_sqrt2_d_t.data.val = 25170707; + math->md_fraction_threshold_t.data.val = fraction_threshold; + math->md_half_fraction_threshold_t.data.val = half_fraction_threshold; + math->md_scaled_threshold_t.data.val = scaled_threshold; + math->md_half_scaled_threshold_t.data.val = half_scaled_threshold; + math->md_near_zero_angle_t.data.val = near_zero_angle; + math->md_p_over_v_threshold_t.data.val = p_over_v_threshold; + math->md_equation_threshold_t.data.val = equation_threshold; + math->md_from_int = mp_set_number_from_int; + math->md_from_boolean = mp_set_number_from_boolean; + math->md_from_scaled = mp_set_number_from_scaled; + math->md_from_double = mp_set_number_from_double; + math->md_from_addition = mp_set_number_from_addition; + math->md_half_from_addition = mp_set_number_half_from_addition; + math->md_from_subtraction = mp_set_number_from_subtraction; + math->md_half_from_subtraction = mp_set_number_half_from_subtraction; + math->md_from_oftheway = mp_set_number_from_of_the_way; + math->md_from_div = mp_set_number_from_div; + math->md_from_mul = mp_set_number_from_mul; + math->md_from_int_div = mp_set_number_from_int_div; + math->md_from_int_mul = mp_set_number_from_int_mul; + math->md_negate = mp_number_negate; + math->md_add = mp_number_add; + math->md_subtract = mp_number_subtract; + math->md_half = mp_number_half; + math->md_do_double = mp_number_double; + math->md_abs = mp_number_abs; + math->md_clone = mp_number_clone; + math->md_negated_clone = mp_number_negated_clone; + math->md_abs_clone = mp_number_abs_clone; + math->md_swap = mp_number_swap; + math->md_add_scaled = mp_number_add_scaled; + math->md_multiply_int = mp_number_multiply_int; + math->md_divide_int = mp_number_divide_int; + math->md_to_int = mp_number_to_int; + math->md_to_boolean = mp_number_to_boolean; + math->md_to_scaled = mp_number_to_scaled; + math->md_to_double = mp_number_to_double; + math->md_odd = mp_number_odd; + math->md_equal = mp_number_equal; + math->md_less = mp_number_less; + math->md_greater = mp_number_greater; + math->md_nonequalabs = mp_number_nonequalabs; + math->md_round_unscaled = mp_round_unscaled; + math->md_floor_scaled = mp_number_floor; + math->md_fraction_to_round_scaled = mp_fraction_to_round_scaled; + math->md_make_scaled = mp_number_make_scaled; + math->md_make_fraction = mp_number_make_fraction; + math->md_take_fraction = mp_number_take_fraction; + math->md_take_scaled = mp_number_take_scaled; + math->md_velocity = mp_velocity; + math->md_n_arg = mp_n_arg; + math->md_m_log = mp_m_log; + math->md_m_exp = mp_m_exp; + math->md_m_unif_rand = mp_m_unif_rand; + math->md_m_norm_rand = mp_m_norm_rand; + math->md_pyth_add = mp_pyth_add; + math->md_pyth_sub = mp_pyth_sub; + math->md_power_of = mp_power_of; + math->md_fraction_to_scaled = mp_number_fraction_to_scaled; + math->md_scaled_to_fraction = mp_number_scaled_to_fraction; + math->md_scaled_to_angle = mp_number_scaled_to_angle; + math->md_angle_to_scaled = mp_number_angle_to_scaled; + math->md_init_randoms = mp_init_randoms; + math->md_sin_cos = mp_n_sin_cos; + math->md_slow_add = mp_slow_add; + math->md_sqrt = mp_square_rt; + math->md_print = mp_print_number; + math->md_tostring = mp_number_tostring; + math->md_modulo = mp_number_modulo; + math->md_ab_vs_cd = mp_ab_vs_cd; + math->md_crossing_point = mp_crossing_point; + math->md_scan_numeric = mp_scan_numeric_token; + math->md_scan_fractional = mp_scan_fractional_token; + math->md_free_math = mp_free_scaled_math; + math->md_set_precision = mp_scaled_set_precision; + return math; +} + +void mp_scaled_set_precision (MP mp) +{ + (void) mp; +} + +void mp_free_scaled_math (MP mp) +{ + mp_free_number(mp, &(mp->math->md_epsilon_t)); + mp_free_number(mp, &(mp->math->md_inf_t)); + mp_free_number(mp, &(mp->math->md_negative_inf_t)); + mp_free_number(mp, &(mp->math->md_arc_tol_k)); + mp_free_number(mp, &(mp->math->md_three_sixty_deg_t)); + mp_free_number(mp, &(mp->math->md_one_eighty_deg_t)); + mp_free_number(mp, &(mp->math->md_negative_one_eighty_deg_t)); + mp_free_number(mp, &(mp->math->md_fraction_one_t)); + mp_free_number(mp, &(mp->math->md_fraction_half_t)); + mp_free_number(mp, &(mp->math->md_fraction_three_t)); + mp_free_number(mp, &(mp->math->md_fraction_four_t)); + mp_free_number(mp, &(mp->math->md_zero_t)); + mp_free_number(mp, &(mp->math->md_half_unit_t)); + mp_free_number(mp, &(mp->math->md_three_quarter_unit_t)); + mp_free_number(mp, &(mp->math->md_unity_t)); + mp_free_number(mp, &(mp->math->md_two_t)); + mp_free_number(mp, &(mp->math->md_three_t)); + mp_free_number(mp, &(mp->math->md_one_third_inf_t)); + mp_free_number(mp, &(mp->math->md_warning_limit_t)); + mp_free_number(mp, &(mp->math->md_one_k)); + mp_free_number(mp, &(mp->math->md_sqrt_8_e_k)); + mp_free_number(mp, &(mp->math->md_twelve_ln_2_k)); + mp_free_number(mp, &(mp->math->md_coef_bound_k)); + mp_free_number(mp, &(mp->math->md_coef_bound_minus_1)); + mp_free_number(mp, &(mp->math->md_twelvebits_3)); + mp_free_number(mp, &(mp->math->md_twentysixbits_sqrt2_t)); + mp_free_number(mp, &(mp->math->md_twentyeightbits_d_t)); + mp_free_number(mp, &(mp->math->md_twentysevenbits_sqrt2_d_t)); + mp_free_number(mp, &(mp->math->md_fraction_threshold_t)); + mp_free_number(mp, &(mp->math->md_half_fraction_threshold_t)); + mp_free_number(mp, &(mp->math->md_scaled_threshold_t)); + mp_free_number(mp, &(mp->math->md_half_scaled_threshold_t)); + mp_free_number(mp, &(mp->math->md_near_zero_angle_t)); + mp_free_number(mp, &(mp->math->md_p_over_v_threshold_t)); + mp_free_number(mp, &(mp->math->md_equation_threshold_t)); + mp_memory_free(mp->math); +} + +void mp_allocate_number (MP mp, mp_number *n, mp_number_type t) +{ + (void) mp; + n->data.val = 0; + n->type = t; +} + +void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) +{ + (void) mp; + n->type = t; + n->data.val = v->data.val; +} + +void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) +{ + (void) mp; + n->type = t; + n->data.val = abs(v->data.val); +} + +void mp_allocate_double (MP mp, mp_number *n, double v) +{ + (void) mp; + n->type = mp_scaled_type; + n->data.val = (int) (v * 65536.0); +} + +void mp_free_number (MP mp, mp_number *n) +{ + (void) mp; + n->type = mp_nan_type; +} + +void mp_set_number_from_int(mp_number *A, int B) +{ + A->data.val = B * 65536; +} + +void mp_set_number_from_boolean(mp_number *A, int B) +{ + A->data.val = B; +} + +void mp_set_number_from_scaled(mp_number *A, int B) +{ + A->data.val = B; +} + +void mp_set_number_from_double(mp_number *A, double B) +{ + A->data.val = (int) (B * 65536.0); +} + +void mp_set_number_from_addition(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.val = B->data.val + C->data.val; +} + +void mp_set_number_half_from_addition(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.val = (B->data.val + C->data.val) / 2; +} + +void mp_set_number_from_subtraction(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.val = B->data.val - C->data.val; +} + +void mp_set_number_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.val = (B->data.val - C->data.val) / 2; +} + +void mp_set_number_from_div(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.val = B->data.val / C->data.val; +} + +void mp_set_number_from_mul(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.val = B->data.val * C->data.val; +} + +void mp_set_number_from_int_div(mp_number *A, mp_number *B, int C) +{ + A->data.val = B->data.val / C; +} + +void mp_set_number_from_int_mul(mp_number *A, mp_number *B, int C) +{ + A->data.val = B->data.val * C; +} + +void mp_set_number_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C) +{ + (void) mp; + A->data.val = B->data.val - mp_take_fraction(mp, (B->data.val - C->data.val), t->data.val); +} + +void mp_number_negate(mp_number *A) +{ + A->data.val = -A->data.val; +} + +void mp_number_add(mp_number *A, mp_number *B) +{ + A->data.val = A->data.val + B->data.val; +} + +void mp_number_subtract(mp_number *A, mp_number *B) +{ + A->data.val = A->data.val - B->data.val; +} + +void mp_number_half(mp_number *A) +{ + A->data.val = A->data.val / 2; +} + +void mp_number_double(mp_number *A) +{ + A->data.val = A->data.val + A->data.val; +} + +void mp_number_add_scaled(mp_number *A, int B) +{ + A->data.val = A->data.val + B; +} + +void mp_number_multiply_int(mp_number *A, int B) +{ + A->data.val = B * A->data.val; +} + +void mp_number_divide_int(mp_number *A, int B) +{ + A->data.val = A->data.val / B; +} + +void mp_number_abs(mp_number *A) +{ + A->data.val = abs(A->data.val); +} + +void mp_number_clone(mp_number *A, mp_number *B) +{ + A->data.val = B->data.val; +} + +void mp_number_negated_clone(mp_number *A, mp_number *B) +{ + A->data.val = -B->data.val; +} + +void mp_number_abs_clone(mp_number *A, mp_number *B) +{ + A->data.val = abs(B->data.val); +} + +void mp_number_swap(mp_number *A, mp_number *B) +{ + int swap_tmp = A->data.val; + A->data.val = B->data.val; + B->data.val = swap_tmp; +} + +void mp_number_fraction_to_scaled(mp_number *A) +{ + A->type = mp_scaled_type; + A->data.val = A->data.val / 4096; +} + +void mp_number_angle_to_scaled(mp_number *A) +{ + A->type = mp_scaled_type; + if (A->data.val >= 0) { + A->data.val = (A->data.val + 8) / 16; + } else { + A->data.val = -((-A->data.val + 8) / 16); + } +} + +void mp_number_scaled_to_fraction(mp_number *A) +{ + A->type = mp_fraction_type; + A->data.val = A->data.val * 4096; +} + +void mp_number_scaled_to_angle(mp_number *A) +{ + A->type = mp_angle_type; + A->data.val = A->data.val * 16; +} + +int mp_number_to_int(mp_number *A) +{ + return A->data.val; +} + +int mp_number_to_scaled(mp_number *A) +{ + return A->data.val; +} + +int mp_number_to_boolean(mp_number *A) +{ + return A->data.val; +} + +double mp_number_to_double(mp_number *A) +{ + return A->data.val / 65536.0; +} + +int mp_number_odd(mp_number *A) +{ + return odd(mp_round_unscaled(A)); +} + +int mp_number_equal(mp_number *A, mp_number *B) { + return A->data.val == B->data.val; +} + +int mp_number_greater(mp_number *A, mp_number *B) +{ + return A->data.val > B->data.val; +} + +int mp_number_less(mp_number *A, mp_number *B) +{ + return A->data.val < B->data.val; +} + +int mp_number_nonequalabs(mp_number *A, mp_number *B) +{ + return abs(A->data.val) != abs(B->data.val); +} + +static char *mp_string_scaled (MP mp, int s) +{ + (void) mp; + static char scaled_string[32]; + int i = 0; + if (s < 0) { + scaled_string[i++] = '-'; + s = -s; + } + mp_snprintf ((scaled_string+i), 12, "%d", (int) (s / unity)); + while (*(scaled_string+i)) { + i++; + } + s = 10 * (s % unity) + 5; + if (s != 5) { + int delta = 10; + scaled_string[i++] = '.'; + do { + if (delta > unity) { + s = s + 0100000 - (delta / 2); + } + scaled_string[i++] = '0' + (s / unity); + s = 10 * (s % unity); + delta = delta * 10; + } while (s > delta); + } + scaled_string[i] = '\0'; + return scaled_string; +} + +void mp_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) +{ + int x = x_orig->data.val; + int y = y_orig->data.val; + if (x >= 0) { + if (y <= EL_GORDO - x) { + ret->data.val = x + y; + } else { + mp->arith_error = 1; + ret->data.val = EL_GORDO; + } + } else if (-y <= EL_GORDO + x) { + ret->data.val = x + y; + } else { + mp->arith_error = 1; + ret->data.val = negative_EL_GORDO; + } +} + +static int mp_make_fraction (MP mp, int p, int q) +{ + if (q == 0) { + mp_confusion (mp, "division by zero"); + return 0; + } else { + double d = TWEXP28 * (double) p / (double) q; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + mp->arith_error = 1; + return EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && (((q > 0 ? -q : q) & 077777) * (((i & 037777) << 1) - 1) & 04000) != 0) { + --i; + } + return i; + } + } else { + d -= 0.5; + if (d <= -TWEXP31) { + mp->arith_error = 1; + return -negative_EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && (((q > 0 ? q : -q) & 077777) * (((i & 037777) << 1) + 1) & 04000) != 0) { + ++i; + } + return i; + } + } + } +} + +void mp_number_make_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) +{ + ret->data.val = mp_make_fraction (mp, p->data.val, q->data.val); +} + +int mp_take_fraction (MP mp, int p, int q) +{ + double d = (double) p *(double) q *TWEXP_28; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + if (d != TWEXP31 || (((p & 077777) * (q & 077777)) & 040000) == 0) { + mp->arith_error = 1; + } + return EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && (((p & 077777) * (q & 077777)) & 040000) != 0) { + --i; + } + return i; + } + } else { + d -= 0.5; + if (d <= -TWEXP31) { + if (d != -TWEXP31 || ((-(p & 077777) * (q & 077777)) & 040000) == 0) { + mp->arith_error = 1; + } + return -negative_EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && ((-(p & 077777) * (q & 077777)) & 040000) != 0) { + ++i; + } + return i; + } + } +} + +void mp_number_take_fraction (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + ret->data.val = mp_take_fraction (mp, p_orig->data.val, q_orig->data.val); +} + +static int mp_take_scaled (MP mp, int p, int q) +{ + double d = (double) p *(double) q *TWEXP_16; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + if (d != TWEXP31 || (((p & 077777) * (q & 077777)) & 040000) == 0) { + mp->arith_error = 1; + } + return EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && (((p & 077777) * (q & 077777)) & 040000) != 0) { + --i; + } + return i; + } + } else { + d -= 0.5; + if (d <= -TWEXP31) { + if (d != -TWEXP31 || ((-(p & 077777) * (q & 077777)) & 040000) == 0) { + mp->arith_error = 1; + } + return -negative_EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && ((-(p & 077777) * (q & 077777)) & 040000) != 0) { + ++i; + } + return i; + } + } +} + +void mp_number_take_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + ret->data.val = mp_take_scaled(mp, p_orig->data.val, q_orig->data.val); +} + +int mp_make_scaled (MP mp, int p, int q) +{ + if (q == 0) { + mp_confusion (mp, "division by zero"); + return 0; + } else { + double d = TWEXP16 * (double) p / (double) q; + if ((p ^ q) >= 0) { + d += 0.5; + if (d >= TWEXP31) { + mp->arith_error = 1; + return EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && (((q > 0 ? -q : q) & 077777) * (((i & 037777) << 1) - 1) & 04000) != 0) { + --i; + } + return i; + } + } else { + d -= 0.5; + if (d <= -TWEXP31) { + mp->arith_error = 1; + return -negative_EL_GORDO; + } else { + int i = (int) d; + if (d == (double) i && (((q > 0 ? q : -q) & 077777) * (((i & 037777) << 1) + 1) & 04000) != 0) { + ++i; + } + return i; + } + } + } +} + +void mp_number_make_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + ret->data.val = mp_make_scaled(mp, p_orig->data.val, q_orig->data.val); +} + +static int mp_round_decimals (MP mp, unsigned char *b, int k) +{ + unsigned a = 0; + int l = 0; + (void) mp; + for (l = k-1; l >= 0; l-- ) { + if (l<16) { + a = (a + (unsigned) (*(b+l) - '0') * two) / 10; + } + } + return (int) (a + 1)/2; +} + +static void mp_wrapup_numeric_token (MP mp, int n, int f) +{ + if (n < 32768) { + int mod = (n * unity + f); + set_cur_mod(mod); + if (mod >= fraction_one) { + if (internal_value(mp_warning_check_internal).data.val > 0 && (mp->scanner_status != mp_tex_flushing_state)) { + char msg[256]; + mp_snprintf(msg, 256, "Number is too large (%s)", mp_string_scaled(mp,mod)); + mp_error( + mp, + msg, + "It is at least 4096. Continue and I'll try to cope with that big value;\n" + "but it might be dangerous. (Set warningcheck:=0 to suppress this message.)" + ); + } + } + } else if (mp->scanner_status != mp_tex_flushing_state) { + mp_error( + mp, + "Enormous number has been reduced", + "I can\'t handle numbers bigger than 32767.99998; so I've changed your constant\n" + "to that maximum amount." + ); + set_cur_mod(EL_GORDO); + } + set_cur_cmd(mp_numeric_command); +} + +void mp_scan_fractional_token (MP mp, int n) +{ + int f; + int k = 0; + do { + k++; + mp->cur_input.loc_field++; + } while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class); + f = mp_round_decimals(mp, (unsigned char *)(mp->buffer+mp->cur_input.loc_field-k), (int) k); + if (f == unity) { + n++; + f = 0; + } + mp_wrapup_numeric_token(mp, n, f); +} + +void mp_scan_numeric_token (MP mp, int n) +{ + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + if (n < 32768) { + n = 10 * n + mp->buffer[mp->cur_input.loc_field] - '0'; + } + mp->cur_input.loc_field++; + } + if (! (mp->buffer[mp->cur_input.loc_field] == '.' && mp->char_class[mp->buffer[mp->cur_input.loc_field + 1]] == mp_digit_class)) { + mp_wrapup_numeric_token(mp, n, 0); + } else { + mp->cur_input.loc_field++; + mp_scan_fractional_token(mp, n); + } +} + +void mp_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t) +{ + int acc, num, denom; + acc = mp_take_fraction(mp, st->data.val - (sf->data.val / 16), sf->data.val - (st->data.val / 16)); + acc = mp_take_fraction(mp, acc, ct->data.val - cf->data.val); + num = fraction_two + mp_take_fraction(mp, acc, 379625062); + denom = fraction_three + mp_take_fraction(mp, ct->data.val, 497706707) + mp_take_fraction (mp, cf->data.val, 307599661); + if (t->data.val != unity) { + num = mp_make_scaled (mp, num, t->data.val); + } + if (num / 4 >= denom) { + ret->data.val = fraction_four; + } else { + ret->data.val = mp_make_fraction(mp, num, denom); + } +} + +static int mp_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) +{ + int a = a_orig->data.val; + int b = b_orig->data.val; + int c = c_orig->data.val; + int d = d_orig->data.val; + if (a < 0) { + a = -a; + b = -b; + } + if (c < 0) { + c = -c; + d = -d; + } + if (d <= 0) { + if (b >= 0) { + if ((a == 0 || b == 0) && (c == 0 || d == 0)) { + return 0; + } else { + return 1; + } + } else if (d == 0) { + return a == 0 ? 0 : -1; + } else { + int q = a; + a = c; + c = q; + q = -b; + b = -d; + d = q; + } + } else if (b <= 0) { + if (b < 0 && a > 0) { + return -1; + } else { + return c == 0 ? 0 : -1; + } + } + while (1) { + int q = a / d; + int r = c / b; + if (q != r) { + return q > r ? 1 : -1; + } else { + q = a % d; + r = c % b; + if (r == 0) { + return q ? 1 : 0; + } else if (q == 0) { + return -1; + } else { + a = b; + b = q; + c = d; + d = r; + } + } + } +} + +static void mp_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc) +{ + int x, xx, x0, x1, x2; + int a = aa->data.val; + int b = bb->data.val; + int c = cc->data.val; + int d; + (void) mp; + if (a < 0) { + ret->data.val = zero_crossing; + return; + } else if (c >= 0) { + if (b >= 0) { + if (c > 0) { + ret->data.val = no_crossing; + } else if ((a == 0) && (b == 0)) { + ret->data.val = no_crossing; + } else { + ret->data.val = one_crossing; + } + return; + } else if (a == 0) { + ret->data.val = zero_crossing; + return; + } + } else if (a == 0) { + if (b <= 0) { + ret->data.val = zero_crossing; + return; + } + } + d = 1; + x0 = a; + x1 = a - b; + x2 = b - c; + do { + x = (x1 + x2) / 2; + if (x1 - x0 > x0) { + x2 = x; + x0 += x0; + d += d; + } else { + xx = x1 + x - x0; + if (xx > x0) { + x2 = x; + x0 += x0; + d += d; + } else { + x0 = x0 - xx; + if ((x <= x0) && (x + x2 <= x0)) { + ret->data.val = no_crossing; + return; + } else { + x1 = x; + d = d + d + 1; + } + } + } + } while (d < fraction_one); + ret->data.val = d - fraction_one; +} + +int mp_round_unscaled(mp_number *x_orig) +{ + int x = x_orig->data.val; + if (x >= 32768) { + return 1 + ((x-32768) / 65536); + } else if (x >= -32768) { + return 0; + } else { + return -(1+((-(x+1)-32768) / 65536)); + } +} + +void mp_number_floor(mp_number *i) +{ + i->data.val = i->data.val&-65536; +} + +void mp_fraction_to_round_scaled(mp_number *x_orig) +{ + int x = x_orig->data.val; + x_orig->type = mp_scaled_type; + x_orig->data.val = (x>=2048 ? 1+((x-2048) / 4096) : ( x>=-2048 ? 0 : -(1+((-(x+1)-2048) / 4096)))); +} + +void mp_square_rt (MP mp, mp_number *ret, mp_number *x_orig) +{ + int x = x_orig->data.val; + if (x <= 0) { + if (x < 0) { + char msg[256]; + mp_snprintf(msg, 256, "Square root of %s has been replaced by 0", mp_string_scaled(mp, x)); + mp_error( + mp, + msg, + "Since I don't take square roots of negative numbers, I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + } + ret->data.val = 0; + } else { + int k = 23; + int y; + int q = 2; + while (x < fraction_two) { + k--; + x = x + x + x + x; + } + if (x < fraction_four) + y = 0; + else { + x = x - fraction_four; + y = 1; + } + do { + x += x; + y += y; + if (x >= fraction_four) { + x = x - fraction_four; + y++; + }; + x += x; + y = y + y - q; + q += q; + if (x >= fraction_four) { + x = x - fraction_four; + y++; + }; + if (y > (int) q) { + y -= q; + q += 2; + } else if (y <= 0) { + q -= 2; + y += q; + }; + k--; + } while (k != 0); + ret->data.val = (int) (q/2); + } +} + +void mp_pyth_add (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + int a = abs(a_orig->data.val); + int b = abs(b_orig->data.val); + if (a < b) { + int r = b; + b = a; + a = r; + } + if (b > 0) { + int big; + if (a < fraction_two) { + big = 0; + } else { + a = a / 4; + b = b / 4; + big = 1; + } + while (1) { + int r = mp_make_fraction(mp, b, a); + r = mp_take_fraction(mp, r, r); + if (r == 0) { + break; + } else { + r = mp_make_fraction(mp, r, fraction_four + r); + a = a + mp_take_fraction(mp, a + a, r); + b = mp_take_fraction(mp, b, r); + } + } + if (big) { + if (a < fraction_two) { + a = a + a + a + a; + } else { + mp->arith_error = 1; + a = EL_GORDO; + } + } + } + ret->data.val = a; +} + +void mp_pyth_sub (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + int a = abs(a_orig->data.val); + int b = abs(b_orig->data.val); + if (a <= b) { + if (a < b) { + char msg[256]; + char *astr = mp_strdup(mp_string_scaled(mp, a)); + mp_snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, mp_string_scaled(mp, b)); + mp_memory_free(astr); + mp_error( + mp, + msg, + "Since I don't take square roots of negative numbers, I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + } + a = 0; + } else { + int big; + if (a < fraction_four) { + big = 0; + } else { + a = (int) a/2; + b = (int) b/2; + big = 1; + } + while (1) { + int r = mp_make_fraction(mp, b, a); + r = mp_take_fraction(mp, r, r); + if (r == 0) { + break; + } else { + r = mp_make_fraction(mp, r, fraction_four - r); + a = a - mp_take_fraction(mp, a + a, r); + b = mp_take_fraction(mp, b, r); + } + } + if (big) { + a *= 2; + } + } + ret->data.val = a; +} + +void mp_power_of (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + double p = pow(mp_number_to_double(a_orig), mp_number_to_double(b_orig)); + long r = lround(p * 65536.0); + if (r > 0) { + if (r >= EL_GORDO) { + mp->arith_error = 1; + r = EL_GORDO; + } + } else if (r < 0) { + if (r <= - EL_GORDO) { + mp->arith_error = 1; + r = - EL_GORDO; + } + } + ret->data.val = r; +} + +void mp_m_log (MP mp, mp_number *ret, mp_number *x_orig) +{ + int x = x_orig->data.val; + if (x <= 0) { + { + char msg[256]; + mp_snprintf(msg, 256, "Logarithm of %s has been replaced by 0", mp_string_scaled(mp, x)); + mp_error( + mp, + msg, + "Since I don't take logs of non-positive numbers, I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + ret->data.val = 0; + } + } else { + int k = 2; + int y = 1302456956 + 4 - 100; + int z = 27595 + 6553600; + while (x < fraction_four) { + x = 2*x; + y -= 93032639; + z -= 48782; + } + y = y + (z / unity); + while (x > fraction_four + 4) { + { + z = ((x - 1) / two_to_the (k)) + 1; + while (x < fraction_four + z) { + z = (z + 1)/2; + k++; + }; + y += mp_m_spec_log[k]; + x -= z; + } + } + ret->data.val = (y / 8); + } +} + +void mp_m_exp (MP mp, mp_number *ret, mp_number *x_orig) +{ + int y, z; + int x = x_orig->data.val; + if (x > 174436200) { + mp->arith_error = 1; + ret->data.val = EL_GORDO; + } else if (x < -197694359) { + ret->data.val = 0; + } else { + if (x <= 0) { + z = -8 * x; + y = 04000000; + } else { + if (x <= 127919879) { + z = 1023359037 - 8 * x; + } else { + z = 8 * (174436200 - x); + } + y = EL_GORDO; + } + { + int k = 1; + while (z > 0) { + while (z >= mp_m_spec_log[k]) { + z -= mp_m_spec_log[k]; + y = y - 1 - ((y - two_to_the(k - 1)) / two_to_the(k)); + } + k++; + } + } + if (x <= 127919879) { + ret->data.val = ((y + 8) / 16); + } else { + ret->data.val = y; + } + } +} + +void mp_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) +{ + int z; + int t; + int k; + int octant; + int x = x_orig->data.val; + int y = y_orig->data.val; + if (x >= 0) { + octant = first_octant; + } else { + x = -x; + octant = first_octant + negate_x; + } + if (y < 0) { + y = -y; + octant = octant + negate_y; + } + if (x < y) { + t = y; + y = x; + x = t; + octant = octant + switch_x_and_y; + } + if (x == 0) { + mp_error( + mp, + "angle(0,0) is taken as zero", + "The 'angle' between two identical points is undefined. I'm zeroing this one.\n" + "Proceed, with fingers crossed." + ); + ret->data.val = 0; + } else { + ret->type = mp_angle_type; + while (x >= fraction_two) { + x = x/2; + y = y/2; + } + z = 0; + if (y > 0) { + while (x < fraction_one) { + x += x; + y += y; + }; + k = 0; + do { + y += y; + k++; + if (y > x) { + z = z + mp_m_spec_atan[k]; + t = x; + x = x + (y / two_to_the(k + k)); + y = y - t; + }; + } while (k != 15); + do { + y += y; + k++; + if (y > x) { + z = z + mp_m_spec_atan[k]; + y = y - x; + }; + } while (k != 26); + } + + switch (octant) { + case first_octant: ret->data.val = z; break; + case second_octant: ret->data.val = -z + ninety_deg; break; + case third_octant: ret->data.val = z + ninety_deg; break; + case fourth_octant: ret->data.val = -z + one_eighty_deg; break; + case fifth_octant: ret->data.val = z - one_eighty_deg; break; + case sixth_octant: ret->data.val = -z - ninety_deg; break; + case seventh_octant: ret->data.val = z - ninety_deg; break; + case eighth_octant: ret->data.val = -z; break; + } + } +} + +void mp_n_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin) +{ + int k; + int q; + int x, y, t; + int z = z_orig->data.val; + mp_number x_n, y_n, ret; + mp_allocate_number(mp, &ret, mp_scaled_type); + mp_allocate_number(mp, &x_n, mp_scaled_type); + mp_allocate_number(mp, &y_n, mp_scaled_type); + while (z < 0) { + z = z + three_sixty_deg; + } + z = z % three_sixty_deg; + q = z / forty_five_deg; + z = z % forty_five_deg; + x = fraction_one; + y = x; + if (! odd(q)) { + z = forty_five_deg - z; + } + k = 1; + while (z > 0) { + if (z >= mp_m_spec_atan[k]) { + z = z - mp_m_spec_atan[k]; + t = x; + x = t + y / two_to_the(k); + y = y - t / two_to_the(k); + } + k++; + } + if (y < 0) { + y = 0; + } + switch (q) { + case 0: break; + case 1: t = x; x = y; y = t; break; + case 2: t = x; x = -y; y = t; break; + case 3: x = -x; break; + case 4: x = -x; y = -y; break; + case 5: t = x; x = -y; y = -t; break; + case 6: t = x; x = y; y = -t; break; + case 7: y = -y; break; + } + x_n.data.val = x; + y_n.data.val = y; + mp_pyth_add(mp, &ret, &x_n, &y_n); + n_cos->data.val = mp_make_fraction(mp, x, ret.data.val); + n_sin->data.val = mp_make_fraction(mp, y, ret.data.val); + mp_free_number(mp, &ret); + mp_free_number(mp, &x_n); + mp_free_number(mp, &y_n); +} + +void mp_init_randoms (MP mp, int seed) +{ + int k = 1; + int j = abs(seed); + while (j >= fraction_one) { + j = j/2; + } + for (int i = 0; i <= 54; i++) { + int jj = k; + k = j - k; + j = jj; + if (k < 0) { + k += fraction_one; + } + mp->randoms[(i * 21) % 55].data.val = j; + } + mp_new_randoms(mp); + mp_new_randoms(mp); + mp_new_randoms(mp); +} + +void mp_print_number (MP mp, mp_number *n) +{ + mp_print_e_str(mp, mp_string_scaled(mp, n->data.val)); +} + +char *mp_number_tostring (MP mp, mp_number *n) +{ + return mp_string_scaled(mp, n->data.val); +} + +void mp_number_modulo(mp_number *a, mp_number *b) +{ + a->data.val = a->data.val % b->data.val; +} + +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])); +} + +static void mp_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig) +{ + mp_number x, abs_x, u, y; + mp_allocate_number(mp, &y, mp_fraction_type); + mp_allocate_clone(mp, &x, mp_scaled_type, x_orig); + mp_allocate_abs(mp, &abs_x, mp_scaled_type, &x); + mp_allocate_number(mp, &u, mp_scaled_type); + mp_next_random(mp, &u); + mp_number_take_fraction(mp, &y, &abs_x, &u); + if (mp_number_equal(&y, &abs_x)) { + mp_number_clone(ret, &((math_data *)mp->math)->md_zero_t); + } else if (mp_number_greater(&x, &((math_data *)mp->math)->md_zero_t)) { + mp_number_clone(ret, &y); + } else { + mp_number_clone(ret, &y); + mp_number_negate(ret); + } + mp_free_number(mp, &y); + mp_free_number(mp, &abs_x); + mp_free_number(mp, &x); + mp_free_number(mp, &u); +} + +static void mp_m_norm_rand (MP mp, mp_number *ret) +{ + mp_number abs_x, u, r, la, xa; + mp_allocate_number(mp, &la, mp_scaled_type); + mp_allocate_number(mp, &xa, mp_scaled_type); + mp_allocate_number(mp, &abs_x, mp_scaled_type); + mp_allocate_number(mp, &u, mp_scaled_type); + mp_allocate_number(mp, &r, mp_scaled_type); + do { + do { + mp_number v; + mp_allocate_number(mp, &v, mp_scaled_type); + mp_next_random(mp, &v); + mp_number_subtract(&v, &((math_data *)mp->math)->md_fraction_half_t); + mp_number_take_fraction(mp, &xa, &((math_data *)mp->math)->md_sqrt_8_e_k, &v); + mp_free_number(mp, &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_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la); + } while (mp_ab_vs_cd(&((math_data *)mp->math)->md_one_k, &la, &xa, &xa) < 0); + mp_number_clone(ret, &xa); + mp_free_number(mp, &r); + mp_free_number(mp, &abs_x); + mp_free_number(mp, &la); + mp_free_number(mp, &xa); + mp_free_number(mp, &u); +} |