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