diff options
Diffstat (limited to 'source/luametatex/source/mp/mpc/mpmathdouble.c')
| -rw-r--r-- | source/luametatex/source/mp/mpc/mpmathdouble.c | 1160 |
1 files changed, 1160 insertions, 0 deletions
diff --git a/source/luametatex/source/mp/mpc/mpmathdouble.c b/source/luametatex/source/mp/mpc/mpmathdouble.c new file mode 100644 index 000000000..99be94727 --- /dev/null +++ b/source/luametatex/source/mp/mpc/mpmathdouble.c @@ -0,0 +1,1160 @@ +/* This file is generated by "mtxrun --script "mtx-wtoc.lua" from the metapost cweb files. */ + + +# include "mpconfig.h" +# include "mpmathdouble.h" + + +# define PI 3.1415926535897932384626433832795028841971 +# define fraction_multiplier 4096.0 +# define angle_multiplier 16.0 +# define coef_bound ((7.0/3.0)*fraction_multiplier) +# define fraction_threshold 0.04096 +# define half_fraction_threshold (fraction_threshold/2) +# define scaled_threshold 0.000122 +# define half_scaled_threshold (scaled_threshold/2) +# define near_zero_angle (0.0256*angle_multiplier) +# define p_over_v_threshold 0x80000 +# define equation_threshold 0.001 +# define warning_limit pow(2.0,52.0) +# define epsilon pow(2.0,-52.0) +# define unity 1.0 +# define two 2.0 +# define three 3.0 +# define half_unit 0.5 +# define three_quarter_unit 0.75 +# define EL_GORDO (DBL_MAX/2.0-1.0) +# define negative_EL_GORDO (-EL_GORDO) +# define one_third_EL_GORDO (EL_GORDO/3.0) +# define fraction_half (0.5*fraction_multiplier) +# define fraction_one (1.0*fraction_multiplier) +# define fraction_two (2.0*fraction_multiplier) +# define fraction_three (3.0*fraction_multiplier) +# define fraction_four (4.0*fraction_multiplier) +# define no_crossing (fraction_one + 1) +# define one_crossing fraction_one +# define zero_crossing 0 +# define one_eighty_deg (180.0*angle_multiplier) +# define negative_one_eighty_deg (-180.0*angle_multiplier) +# define three_sixty_deg (360.0*angle_multiplier) +# 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.dval = (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 *v); +static void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v); +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 int mp_double_ab_vs_cd (mp_number *a, mp_number *b, mp_number *c, mp_number *d); +static void mp_double_abs (mp_number *A); +static void mp_double_crossing_point (MP mp, mp_number *ret, mp_number *a, mp_number *b, mp_number *c); +static void mp_double_fraction_to_round_scaled (mp_number *x); +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_m_norm_rand (MP mp, mp_number *ret); +static void mp_double_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig); +static void mp_double_n_arg (MP mp, mp_number *ret, mp_number *x, mp_number *y); +static void mp_double_number_make_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_double_number_make_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_double_number_take_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_double_number_take_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); +static void mp_double_power_of (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_double_print_number (MP mp, mp_number *n); +static void mp_double_pyth_add (MP mp, mp_number *r, mp_number *a, mp_number *b); +static void mp_double_pyth_sub (MP mp, mp_number *r, mp_number *a, mp_number *b); +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_double_set_precision (MP mp); +static void mp_double_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin); +static void mp_double_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig); +static void mp_double_square_rt (MP mp, mp_number *ret, mp_number *x_orig); +static void mp_double_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t); +static void mp_free_double_math (MP mp); +static void mp_free_number (MP mp, mp_number *n); +static void mp_init_randoms (MP mp, int seed); +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_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 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 int mp_round_unscaled (mp_number *x_orig); +static void mp_set_double_from_addition (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_double_from_boolean (mp_number *A, int B); +static void mp_set_double_from_div (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_double_from_double (mp_number *A, double B); +static void mp_set_double_from_int (mp_number *A, int B); +static void mp_set_double_from_int_div (mp_number *A, mp_number *B, int C); +static void mp_set_double_from_int_mul (mp_number *A, mp_number *B, int C); +static void mp_set_double_from_mul (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_double_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C); +static void mp_set_double_from_scaled (mp_number *A, int B); +static void mp_set_double_from_subtraction (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_double_half_from_addition (mp_number *A, mp_number *B, mp_number *C); +static void mp_set_double_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C); +static void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop); +static char *mp_double_number_tostring (MP mp, mp_number *n); +inline double mp_double_make_fraction (double p, double q) { return (p / q) * fraction_multiplier; } +inline double mp_double_take_fraction (double p, double q) { return (p * q) / fraction_multiplier; } +inline double mp_double_make_scaled (double p, double q) { return p / q; } + +math_data *mp_initialize_double_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.dval = 16 * unity; + math->md_precision_max.data.dval = 16 * unity; + math->md_precision_min.data.dval = 16 * unity; + math->md_epsilon_t.data.dval = epsilon; + math->md_inf_t.data.dval = EL_GORDO; + math->md_negative_inf_t.data.dval = negative_EL_GORDO; + math->md_one_third_inf_t.data.dval = one_third_EL_GORDO; + math->md_warning_limit_t.data.dval = warning_limit; + math->md_unity_t.data.dval = unity; + math->md_two_t.data.dval = two; + math->md_three_t.data.dval = three; + math->md_half_unit_t.data.dval = half_unit; + math->md_three_quarter_unit_t.data.dval = three_quarter_unit; + math->md_arc_tol_k.data.dval = (unity/4096); + math->md_fraction_one_t.data.dval = fraction_one; + math->md_fraction_half_t.data.dval = fraction_half; + math->md_fraction_three_t.data.dval = fraction_three; + math->md_fraction_four_t.data.dval = fraction_four; + math->md_three_sixty_deg_t.data.dval = three_sixty_deg; + math->md_one_eighty_deg_t.data.dval = one_eighty_deg; + math->md_negative_one_eighty_deg_t.data.dval = negative_one_eighty_deg; + math->md_one_k.data.dval = 1.0/64 ; + math->md_sqrt_8_e_k.data.dval = 1.71552776992141359295; + math->md_twelve_ln_2_k.data.dval = 8.31776616671934371292 *256; + math->md_coef_bound_k.data.dval = coef_bound; + math->md_coef_bound_minus_1.data.dval = coef_bound - 1/65536.0; + math->md_twelvebits_3.data.dval = 1365 / 65536.0; + math->md_twentysixbits_sqrt2_t.data.dval = 94906266 / 65536.0; + math->md_twentyeightbits_d_t.data.dval = 35596755 / 65536.0; + math->md_twentysevenbits_sqrt2_d_t.data.dval = 25170707 / 65536.0; + math->md_fraction_threshold_t.data.dval = fraction_threshold; + math->md_half_fraction_threshold_t.data.dval = half_fraction_threshold; + math->md_scaled_threshold_t.data.dval = scaled_threshold; + math->md_half_scaled_threshold_t.data.dval = half_scaled_threshold; + math->md_near_zero_angle_t.data.dval = near_zero_angle; + math->md_p_over_v_threshold_t.data.dval = p_over_v_threshold; + math->md_equation_threshold_t.data.dval = equation_threshold; + math->md_from_int = mp_set_double_from_int; + math->md_from_boolean = mp_set_double_from_boolean; + math->md_from_scaled = mp_set_double_from_scaled; + math->md_from_double = mp_set_double_from_double; + math->md_from_addition = mp_set_double_from_addition; + math->md_half_from_addition = mp_set_double_half_from_addition; + math->md_from_subtraction = mp_set_double_from_subtraction; + math->md_half_from_subtraction = mp_set_double_half_from_subtraction; + math->md_from_oftheway = mp_set_double_from_of_the_way; + math->md_from_div = mp_set_double_from_div; + math->md_from_mul = mp_set_double_from_mul; + math->md_from_int_div = mp_set_double_from_int_div; + math->md_from_int_mul = mp_set_double_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_double_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_boolean = mp_number_to_boolean; + math->md_to_scaled = mp_number_to_scaled; + math->md_to_double = mp_number_to_double; + math->md_to_int = mp_number_to_int; + 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_double_fraction_to_round_scaled; + math->md_make_scaled = mp_double_number_make_scaled; + math->md_make_fraction = mp_double_number_make_fraction; + math->md_take_fraction = mp_double_number_take_fraction; + math->md_take_scaled = mp_double_number_take_scaled; + math->md_velocity = mp_double_velocity; + math->md_n_arg = mp_double_n_arg; + math->md_m_log = mp_double_m_log; + math->md_m_exp = mp_double_m_exp; + math->md_m_unif_rand = mp_double_m_unif_rand; + math->md_m_norm_rand = mp_double_m_norm_rand; + math->md_pyth_add = mp_double_pyth_add; + math->md_pyth_sub = mp_double_pyth_sub; + math->md_power_of = mp_double_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_double_sin_cos; + math->md_slow_add = mp_double_slow_add; + math->md_sqrt = mp_double_square_rt; + math->md_print = mp_double_print_number; + math->md_tostring = mp_double_number_tostring; + math->md_modulo = mp_number_modulo; + math->md_ab_vs_cd = mp_ab_vs_cd; + math->md_crossing_point = mp_double_crossing_point; + math->md_scan_numeric = mp_double_scan_numeric_token; + math->md_scan_fractional = mp_double_scan_fractional_token; + math->md_free_math = mp_free_double_math; + math->md_set_precision = mp_double_set_precision; + return math; +} + +void mp_double_set_precision (MP mp) +{ + (void) mp; +} + +void mp_free_double_math (MP mp) +{ + 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_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_inf_t)); + mp_free_number(mp, &(mp->math->md_negative_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_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.dval = 0.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.dval = v->data.dval; +} + +void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) +{ + (void) mp; + n->type = t; + n->data.dval = fabs(v->data.dval); +} + +void mp_allocate_double (MP mp, mp_number *n, double v) +{ + (void) mp; + n->type = mp_scaled_type; + n->data.dval = v; +} + +void mp_free_number (MP mp, mp_number *n) +{ + (void) mp; + n->type = mp_nan_type; +} + +void mp_set_double_from_int(mp_number *A, int B) +{ + A->data.dval = B; +} + +void mp_set_double_from_boolean(mp_number *A, int B) +{ + A->data.dval = B; +} + +void mp_set_double_from_scaled(mp_number *A, int B) +{ + A->data.dval = B / 65536.0; +} + +void mp_set_double_from_double(mp_number *A, double B) +{ + A->data.dval = B; +} + +void mp_set_double_from_addition(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.dval = B->data.dval + C->data.dval; +} + +void mp_set_double_half_from_addition(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.dval = (B->data.dval + C->data.dval) / 2.0; +} + +void mp_set_double_from_subtraction(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.dval = B->data.dval - C->data.dval; +} + +void mp_set_double_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.dval = (B->data.dval - C->data.dval) / 2.0; +} + +void mp_set_double_from_div(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.dval = B->data.dval / C->data.dval; +} + +void mp_set_double_from_mul(mp_number *A, mp_number *B, mp_number *C) +{ + A->data.dval = B->data.dval * C->data.dval; +} + +void mp_set_double_from_int_div(mp_number *A, mp_number *B, int C) +{ + A->data.dval = B->data.dval / C; +} + +void mp_set_double_from_int_mul(mp_number *A, mp_number *B, int C) +{ + A->data.dval = B->data.dval * C; +} + +void mp_set_double_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C) +{ + (void) mp; + A->data.dval = B->data.dval - mp_double_take_fraction(B->data.dval - C->data.dval, t->data.dval); +} + +void mp_number_negate(mp_number *A) +{ + A->data.dval = -A->data.dval; + if (A->data.dval == -0.0) { + A->data.dval = 0.0; + } +} + +void mp_number_add(mp_number *A, mp_number *B) +{ + A->data.dval = A->data.dval + B->data.dval; +} + +void mp_number_subtract(mp_number *A, mp_number *B) +{ + A->data.dval = A->data.dval - B->data.dval; +} + +void mp_number_half(mp_number *A) +{ + A->data.dval = A->data.dval / 2.0; +} + +void mp_number_double(mp_number *A) +{ + A->data.dval = A->data.dval * 2.0; +} + +void mp_number_add_scaled(mp_number *A, int B) +{ + A->data.dval = A->data.dval + (B / 65536.0); +} + +void mp_number_multiply_int(mp_number *A, int B) +{ + A->data.dval = (double)(A->data.dval * B); +} + +void mp_number_divide_int(mp_number *A, int B) +{ + A->data.dval = A->data.dval / (double)B; +} + +void mp_double_abs(mp_number *A) +{ + A->data.dval = fabs(A->data.dval); +} + +void mp_number_clone(mp_number *A, mp_number *B) +{ + A->data.dval = B->data.dval; +} + +void mp_number_negated_clone(mp_number *A, mp_number *B) +{ + A->data.dval = -B->data.dval; + if (A->data.dval == -0.0) { + A->data.dval = 0.0; + } +} + +void mp_number_abs_clone(mp_number *A, mp_number *B) +{ + A->data.dval = fabs(B->data.dval); +} + +void mp_number_swap(mp_number *A, mp_number *B) +{ + double swap_tmp = A->data.dval; + A->data.dval = B->data.dval; + B->data.dval = swap_tmp; +} + +void mp_number_fraction_to_scaled(mp_number *A) +{ + A->type = mp_scaled_type; + A->data.dval = A->data.dval / fraction_multiplier; +} + +void mp_number_angle_to_scaled(mp_number *A) +{ + A->type = mp_scaled_type; + A->data.dval = A->data.dval / angle_multiplier; +} + +void mp_number_scaled_to_fraction(mp_number *A) +{ + A->type = mp_fraction_type; + A->data.dval = A->data.dval * fraction_multiplier; +} + +void mp_number_scaled_to_angle(mp_number *A) +{ + A->type = mp_angle_type; + A->data.dval = A->data.dval * angle_multiplier; +} + +int mp_number_to_scaled(mp_number *A) +{ + return (int) lround(A->data.dval * 65536.0); +} + +int mp_number_to_int(mp_number *A) +{ + return (int) (A->data.dval); +} + +int mp_number_to_boolean(mp_number *A) +{ + return (int) (A->data.dval); +} + +double mp_number_to_double(mp_number *A) +{ + return A->data.dval; +} + +int mp_number_odd(mp_number *A) +{ + return odd((int) lround(A->data.dval)); +} + +int mp_number_equal(mp_number *A, mp_number *B) +{ + return A->data.dval == B->data.dval; +} + +int mp_number_greater(mp_number *A, mp_number *B) +{ + return A->data.dval > B->data.dval; +} + +int mp_number_less(mp_number *A, mp_number *B) +{ + return A->data.dval < B->data.dval; +} + +int mp_number_nonequalabs(mp_number *A, mp_number *B) +{ + return fabs(A->data.dval) != fabs(B->data.dval); +} + +char *mp_double_number_tostring (MP mp, mp_number *n) +{ + static char set[64]; + int l = 0; + char *ret = mp_memory_allocate(64); + (void) mp; + snprintf(set, 64, "%.17g", n->data.dval); + while (set[l] == ' ') { + l++; + } + strcpy(ret, set+l); + return ret; +} + +void mp_double_print_number (MP mp, mp_number *n) +{ + char *str = mp_double_number_tostring(mp, n); + mp_print_e_str(mp, str); + mp_memory_free(str); +} + +void mp_double_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) +{ + double x = x_orig->data.dval; + double y = y_orig->data.dval; + if (x >= 0.0) { + if (y <= EL_GORDO - x) { + ret->data.dval = x + y; + } else { + mp->arith_error = 1; + ret->data.dval = EL_GORDO; + } + } else if (-y <= EL_GORDO + x) { + ret->data.dval = x + y; + } else { + mp->arith_error = 1; + ret->data.dval = negative_EL_GORDO; + } +} + +void mp_double_number_make_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) { + (void) mp; + ret->data.dval = mp_double_make_fraction(p->data.dval, q->data.dval); +} + +void mp_double_number_take_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) { + (void) mp; + ret->data.dval = mp_double_take_fraction(p->data.dval, q->data.dval); +} + +void mp_double_number_take_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + (void) mp; + ret->data.dval = p_orig->data.dval * q_orig->data.dval; +} + +void mp_double_number_make_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) +{ + (void) mp; + ret->data.dval = p_orig->data.dval / q_orig->data.dval; +} + +void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop) +{ + double result; + char *end = (char *) stop; + errno = 0; + result = strtod((char *) start, &end); + if (errno == 0) { + set_cur_mod(result); + if (result >= warning_limit) { + if (internal_value(mp_warning_check_internal).data.dval > 0 && (mp->scanner_status != mp_tex_flushing_state)) { + char msg[256]; + mp_snprintf(msg, 256, "Number is too large (%g)", result); + mp_error( + mp, + msg, + "Continue and I'll try to cope with that big value; 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 could not handle this number specification probably because it is out of" + "range." + ); + set_cur_mod(EL_GORDO); + } + set_cur_cmd(mp_numeric_command); +} + +static void mp_double_aux_find_exponent (MP mp) +{ + if (mp->buffer[mp->cur_input.loc_field] == 'e' || mp->buffer[mp->cur_input.loc_field] == 'E') { + mp->cur_input.loc_field++; + if (!(mp->buffer[mp->cur_input.loc_field] == '+' + || mp->buffer[mp->cur_input.loc_field] == '-' + || mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class)) { + mp->cur_input.loc_field--; + return; + } + if (mp->buffer[mp->cur_input.loc_field] == '+' + || mp->buffer[mp->cur_input.loc_field] == '-') { + mp->cur_input.loc_field++; + } + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + } +} + +void mp_double_scan_fractional_token (MP mp, int n) +{ + unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; + unsigned char *stop; + (void) n; + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + mp_double_aux_find_exponent(mp); + stop = &mp->buffer[mp->cur_input.loc_field-1]; + mp_wrapup_numeric_token(mp, start, stop); +} + +void mp_double_scan_numeric_token (MP mp, int n) +{ + unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; + unsigned char *stop; + (void) n; + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + if (mp->buffer[mp->cur_input.loc_field] == '.' && mp->buffer[mp->cur_input.loc_field+1] != '.') { + mp->cur_input.loc_field++; + while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == mp_digit_class) { + mp->cur_input.loc_field++; + } + } + mp_double_aux_find_exponent(mp); + stop = &mp->buffer[mp->cur_input.loc_field-1]; + mp_wrapup_numeric_token(mp, start, stop); +} + +void mp_double_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t) +{ + double acc, num, denom; + (void) mp; + acc = mp_double_take_fraction(st->data.dval - (sf->data.dval / 16.0), sf->data.dval - (st->data.dval / 16.0)); + acc = mp_double_take_fraction(acc, ct->data.dval - cf->data.dval); + num = fraction_two + mp_double_take_fraction(acc, sqrt(2)*fraction_one); + denom = fraction_three + + mp_double_take_fraction(ct->data.dval, 3*fraction_half*(sqrt(5.0)-1.0)) + + mp_double_take_fraction(cf->data.dval, 3*fraction_half*(3.0-sqrt(5.0))); + if (t->data.dval != unity) { + num = mp_double_make_scaled(num, t->data.dval); + } + if (num / 4 >= denom) { + ret->data.dval = fraction_four; + } else { + ret->data.dval = mp_double_make_fraction(num, denom); + } +} + +int mp_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) +{ + return mp_double_ab_vs_cd(a_orig, b_orig, c_orig, d_orig); +} + +static void mp_double_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc) +{ + double d; + double xx, x0, x1, x2; + double a = aa->data.dval; + double b = bb->data.dval; + double c = cc->data.dval; + (void) mp; + if (a < 0.0) { + ret->data.dval = zero_crossing; + return; + } + if (c >= 0.0) { + if (b >= 0.0) { + if (c > 0.0) { + ret->data.dval = no_crossing; + } else if ((a == 0.0) && (b == 0.0)) { + ret->data.dval = no_crossing; + } else { + ret->data.dval = one_crossing; + } + return; + } + if (a == 0.0) { + ret->data.dval = zero_crossing; + return; + } + } else if ((a == 0.0) && (b <= 0.0)) { + ret->data.dval = zero_crossing; + return; + } + d = epsilon; + x0 = a; + x1 = a - b; + x2 = b - c; + do { + double x = (x1 + x2) / 2 + 1E-12; + 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.dval = no_crossing; + return; + } + x1 = x; + d = d + d + epsilon; + } + } + } while (d < fraction_one); + ret->data.dval = (d - fraction_one); +} + +int mp_round_unscaled(mp_number *x_orig) +{ + return (int) lround(x_orig->data.dval); +} + +void mp_number_floor(mp_number *i) +{ + i->data.dval = floor(i->data.dval); +} + +void mp_double_fraction_to_round_scaled(mp_number *x_orig) +{ + double x = x_orig->data.dval; + x_orig->type = mp_scaled_type; + x_orig->data.dval = x/fraction_multiplier; +} + +void mp_double_square_rt (MP mp, mp_number *ret, mp_number *x_orig) +{ + double x = x_orig->data.dval; + if (x > 0) { + ret->data.dval = sqrt(x); + } else { + if (x < 0) { + char msg[256]; + char *xstr = mp_double_number_tostring(mp, x_orig); + mp_snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr); + mp_memory_free(xstr); + 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.dval = 0; + } +} + +void mp_double_pyth_add (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + double a = fabs(a_orig->data.dval); + double b = fabs(b_orig->data.dval); + errno = 0; + ret->data.dval = sqrt(a*a + b*b); + if (errno) { + mp->arith_error = 1; + ret->data.dval = EL_GORDO; + } +} + +void mp_double_pyth_sub (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + double a = fabs(a_orig->data.dval); + double b = fabs(b_orig->data.dval); + if (a > b) { + a = sqrt(a*a - b*b); + } else { + if (a < b) { + char msg[256]; + char *astr = mp_double_number_tostring(mp, a_orig); + char *bstr = mp_double_number_tostring(mp, b_orig); + mp_snprintf(msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr); + mp_memory_free(astr); + mp_memory_free(bstr); + mp_error( + mp, + msg, + "Since I don't take square roots of negative numbers, Im zeroing this one.\n" + "Proceed, with fingers crossed." + ); + } + a = 0; + } + ret->data.dval = a; +} + +void mp_double_power_of (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) +{ + errno = 0; + ret->data.dval = pow(a_orig->data.dval, b_orig->data.dval); + if (errno) { + mp->arith_error = 1; + ret->data.dval = EL_GORDO; + } +} + +void mp_double_m_log (MP mp, mp_number *ret, mp_number *x_orig) +{ + if (x_orig->data.dval > 0) { + ret->data.dval = log(x_orig->data.dval)*256.0; + } else { + char msg[256]; + char *xstr = mp_double_number_tostring(mp, x_orig); + mp_snprintf(msg, 256, "Logarithm of %s has been replaced by 0", xstr); + mp_memory_free(xstr); + 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.dval = 0; + } +} + +void mp_double_m_exp (MP mp, mp_number *ret, mp_number *x_orig) +{ + errno = 0; + ret->data.dval = exp(x_orig->data.dval/256.0); + if (errno) { + if (x_orig->data.dval > 0) { + mp->arith_error = 1; + ret->data.dval = EL_GORDO; + } else { + ret->data.dval = 0; + } + } +} + +void mp_double_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) +{ + if (x_orig->data.dval == 0.0 && y_orig->data.dval == 0.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.dval = 0; + } else { + ret->type = mp_angle_type; + ret->data.dval = atan2(y_orig->data.dval, x_orig->data.dval) * (180.0 / PI) * angle_multiplier; + if (ret->data.dval == -0.0) + ret->data.dval = 0.0; + } +} + +void mp_double_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin) +{ + double rad = (z_orig->data.dval / angle_multiplier); + (void) mp; + if ((rad == 90.0) || (rad == -270)){ + n_cos->data.dval = 0.0; + n_sin->data.dval = fraction_multiplier; + } else if ((rad == -90.0) || (rad == 270.0)) { + n_cos->data.dval = 0.0; + n_sin->data.dval = -fraction_multiplier; + } else if ((rad == 180.0) || (rad == -180.0)) { + n_cos->data.dval = -fraction_multiplier; + n_sin->data.dval = 0.0; + } else { + rad = rad * PI/180.0; + n_cos->data.dval = cos(rad) * fraction_multiplier; + n_sin->data.dval = sin(rad) * fraction_multiplier; + } +} + +# define KK 100 +# define LL 37 +# define MM (1L<<30) +# define mod_diff(x,y) (((x)-(y))&(MM-1)) +# define TT 70 +# define is_odd(x) ((x)&1) +# define QUALITY 1009 + + +typedef struct mp_double_random_info { + long x[KK]; + long buf[QUALITY]; + long dummy; + long started; + long *ptr; +} mp_double_random_info; + +mp_double_random_info mp_double_random_data = { + .dummy = -1, + .started = -1, + .ptr = &mp_double_random_data.dummy +}; + +static void mp_double_aux_ran_array(long aa[], int n) +{ + int i, j; + for (j = 0; j < KK; j++) { + aa[j] = mp_double_random_data.x[j]; + } + for (; j < n; j++) { + aa[j] = mod_diff(aa[j - KK], aa[j - LL]); + } + for (i = 0; i < LL; i++, j++) { + mp_double_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]); + } + for (; i < KK; i++, j++) { + mp_double_random_data.x[i] = mod_diff(aa[j - KK], mp_double_random_data.x[i - LL]); + } +} + + +static void mp_double_aux_ran_start(long seed) +{ + int t, j; + long x[KK + KK - 1]; + long ss = (seed+2) & (MM - 2); + for (j = 0; j < KK; j++) { + x[j] = ss; + ss <<= 1; + if (ss >= MM) { + ss -= MM - 2; + } + } + x[1]++; + for (ss = seed & (MM - 1), t = TT - 1; t;) { + for (j = KK - 1; j > 0; j--) { + x[j + j] = x[j]; + x[j + j - 1] = 0; + } + for (j = KK + KK - 2; j >= KK; j--) { + x[j - (KK -LL)] = mod_diff(x[j - (KK - LL)], x[j]); + x[j - KK] = mod_diff(x[j - KK], x[j]); + } + if (is_odd(ss)) { + for (j = KK; j>0; j--) { + x[j] = x[j-1]; + } + x[0] = x[KK]; + x[LL] = mod_diff(x[LL], x[KK]); + } + if (ss) { + ss >>= 1; + } else { + t--; + } + } + for (j = 0; j < LL; j++) { + mp_double_random_data.x[j + KK - LL] = x[j]; + } + for (;j < KK; j++) { + mp_double_random_data.x[j - LL] = x[j]; + } + for (j = 0; j < 10; j++) { + mp_double_aux_ran_array(x, KK + KK - 1); + } + mp_double_random_data.ptr = &mp_double_random_data.started; +} + +# define mp_double_aux_ran_arr_next() (*mp_double_random_data.ptr>=0? *mp_double_random_data.ptr++: mp_double_aux_ran_arr_cycle()) + +static long mp_double_aux_ran_arr_cycle(void) +{ + if (mp_double_random_data.ptr == &mp_double_random_data.dummy) { + mp_double_aux_ran_start(314159L); + } + mp_double_aux_ran_array(mp_double_random_data.buf, QUALITY); + mp_double_random_data.buf[KK] = -1; + mp_double_random_data.ptr = mp_double_random_data.buf + 1; + return mp_double_random_data.buf[0]; +} + +void mp_init_randoms (MP mp, int seed) +{ + int k = 1; + int j = abs(seed); + int f = (int) fraction_one; + while (j >= f) { + j = j/2; + } + for (int i = 0; i <= 54; i++) { + int jj = k; + k = j - k; + j = jj; + if (k < 0) { + k += f; + } + mp->randoms[(i * 21) % 55].data.dval = j; + } + mp_new_randoms(mp); + mp_new_randoms(mp); + mp_new_randoms(mp); + mp_double_aux_ran_start((unsigned long) seed); +} + +void mp_number_modulo(mp_number *a, mp_number *b) +{ + double tmp; + a->data.dval = modf((double) a->data.dval / (double) b->data.dval, &tmp) * (double) b->data.dval; +} + +static void mp_next_unif_random (MP mp, mp_number *ret) +{ + unsigned long int op = (unsigned) mp_double_aux_ran_arr_next(); + double a = op / (MM * 1.0); + (void) mp; + ret->data.dval = a; +} + +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_double_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_unif_random(mp, &u); + y.data.dval = abs_x.data.dval * u.data.dval; + mp_free_number(mp, &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_negated_clone(ret, &y); + } + mp_free_number(mp, &abs_x); + mp_free_number(mp, &x); + mp_free_number(mp, &y); +} + +static void mp_double_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_double_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_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_subtraction(&la, &((math_data *)mp->math)->md_twelve_ln_2_k, &la); + } while (mp_double_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); +} + +int mp_double_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) +{ + double ab = a_orig->data.dval * b_orig->data.dval; + double cd = c_orig->data.dval * d_orig->data.dval; + if (ab > cd) { + return 1; + } else if (ab < cd) { + return -1; + } else { + return 0; + } +} |