/* 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); }