/* This file is generated by "mtxrun --script "mtx-wtoc.lua" from the metapost cweb files. */ # include "mpconfig.h" # include "mpmathdecimal.h" # define DECNUMDIGITS 1000 # include "decNumber.h" # define E_STRING "2.7182818284590452353602874713526624977572470936999595749669676277240766303535" # define PI_STRING "3.1415926535897932384626433832795028841971693993751058209749445923078164062862" # define fraction_multiplier 4096 # define angle_multiplier 16 # define unity 1 # define two 2 # define three 3 # define four 4 # define half_unit 0.5 # define three_quarter_unit 0.75 # 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 epsilon pow(2.0,-173.0) # define epsilonf pow(2.0,-52.0) # define EL_GORDO "1E1000000" # define negative_EL_GORDO "-1E1000000" # define warning_limit "1E1000000" # define DECPRECISION_DEFAULT 34 # define FACTORIALS_CACHESIZE 50 # define too_precise(a) (a == (DEC_Inexact+DEC_Rounded)) # define too_large(a) (a & DEC_Overflow) # define fraction_half (fraction_multiplier/2) # define fraction_one (1*fraction_multiplier) # define fraction_two (2*fraction_multiplier) # define fraction_three (3*fraction_multiplier) # define fraction_four (4*fraction_multiplier) # define no_crossing mp_decimal_data.fraction_one_plus_decNumber # define one_crossing mp_decimal_data.fraction_one_decNumber # define zero_crossing mp_decimal_data.zero # define odd(A) (abs(A) % 2 == 1) # define set_cur_cmd(A) mp->cur_mod_->type = (A) # define set_cur_mod(A) decNumberCopy((decNumber *) (mp->cur_mod_->data.n.data.num), A) # define decNumberIsPositive(A) (! (decNumberIsZero(A) || decNumberIsNegative(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 void mp_decnumber_check (MP mp, decNumber *dec, decContext *context); static void mp_decimal_abs (mp_number *A); static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number *a, mp_number *b, mp_number *c); static void mp_decimal_fraction_to_round_scaled (mp_number *x); static void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number *x_orig); static void mp_decimal_m_log (MP mp, mp_number *ret, mp_number *x_orig); static void mp_decimal_m_norm_rand (MP mp, mp_number *ret); static void mp_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number *x_orig); void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q); static void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number *x, mp_number *y); static void mp_decimal_number_make_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); static void mp_decimal_number_make_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); static void mp_decimal_number_modulo (mp_number *a, mp_number *b); static void mp_decimal_number_take_fraction (MP mp, mp_number *r, mp_number *p, mp_number *q); static void mp_decimal_number_take_scaled (MP mp, mp_number *r, mp_number *p, mp_number *q); static void mp_decimal_power_of (MP mp, mp_number *r, mp_number *a, mp_number *b); static void mp_decimal_print_number (MP mp, mp_number *n); static void mp_decimal_pyth_add (MP mp, mp_number *r, mp_number *a, mp_number *b); static void mp_decimal_pyth_sub (MP mp, mp_number *r, mp_number *a, mp_number *b); static void mp_decimal_scan_fractional_token (MP mp, int n); static void mp_decimal_scan_numeric_token (MP mp, int n); static void mp_decimal_set_precision (MP mp); static void mp_decimal_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin); static void mp_decimal_slow_add (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig); static void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number *x_orig); void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q); static void mp_decimal_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_decimal_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_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_decimal_from_addition (mp_number *A, mp_number *B, mp_number *C); static void mp_set_decimal_from_boolean (mp_number *A, int B); static void mp_set_decimal_from_div (mp_number *A, mp_number *B, mp_number *C); static void mp_set_decimal_from_double (mp_number *A, double B); static void mp_set_decimal_from_int (mp_number *A, int B); static void mp_set_decimal_from_int_div (mp_number *A, mp_number *B, int C); static void mp_set_decimal_from_int_mul (mp_number *A, mp_number *B, int C); static void mp_set_decimal_from_mul (mp_number *A, mp_number *B, mp_number *C); static void mp_set_decimal_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C); static void mp_set_decimal_from_scaled (mp_number *A, int B); static void mp_set_decimal_from_subtraction (mp_number *A, mp_number *B, mp_number *C); static void mp_set_decimal_half_from_addition (mp_number *A, mp_number *B, mp_number *C); static void mp_set_decimal_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_decimal_number_tostring (MP mp, mp_number *n); static char *mp_decnumber_tostring (decNumber *n); typedef struct mp_decimal_info { decContext set; decContext limitedset; decNumber zero; decNumber one; decNumber minusone; decNumber two_decNumber; decNumber three_decNumber; decNumber four_decNumber; decNumber fraction_multiplier_decNumber; decNumber angle_multiplier_decNumber; decNumber fraction_one_decNumber; decNumber fraction_one_plus_decNumber; decNumber PI_decNumber; decNumber epsilon_decNumber; decNumber EL_GORDO_decNumber; decNumber negative_EL_GORDO_decNumber; decNumber **factorials; int last_cached_factorial; int initialized; } mp_decimal_info; mp_decimal_info mp_decimal_data = { .factorials = NULL, .last_cached_factorial = 0, .initialized = 0, }; void mp_decnumber_check(MP mp, decNumber *dec, decContext *context) { int test = 0; (void) mp; if (context->status & DEC_Overflow) { test = 1; context->status &= ~DEC_Overflow; } if (context->status & DEC_Underflow) { test = 1; context->status &= ~DEC_Underflow; } if (context->status & DEC_Errors) { test = 1; decNumberZero(dec); } context->status = 0; if (decNumberIsSpecial(dec)) { test = 1; if (decNumberIsInfinite(dec)) { if (decNumberIsNegative(dec)) { decNumberCopyNegate(dec, &mp_decimal_data.EL_GORDO_decNumber); } else { decNumberCopy(dec, &mp_decimal_data.EL_GORDO_decNumber); } } else { decNumberZero(dec); } } if (decNumberIsZero(dec) && decNumberIsNegative(dec)) { decNumberZero(dec); } mp->arith_error = test; } static void checkZero(decNumber *ret) { if (decNumberIsZero(ret) && decNumberIsNegative(ret)) { decNumberZero(ret); } } static int decNumberLess(decNumber *a, decNumber *b) { decNumber comp; decNumberCompare(&comp, a, b, &mp_decimal_data.set); return decNumberIsNegative(&comp); } static int decNumberGreater(decNumber *a, decNumber *b) { decNumber comp; decNumberCompare(&comp, a, b, &mp_decimal_data.set); return decNumberIsPositive(&comp); } static void decNumberFromDouble(decNumber *A, double B) { char buffer[1000]; char *c = buffer; snprintf(buffer, 1000, "%-650.325lf", B); while (*c++) { if (*c == ' ') { *c = '\0'; break; } } decNumberFromString(A, buffer, &mp_decimal_data.set); } static double decNumberToDouble(decNumber *A) { char *buffer = mp_memory_allocate(A->digits + 14); double res = 0.0; decNumberToString(A, buffer); if (sscanf(buffer, "%lf", &res)) { mp_memory_free(buffer); return res; } else { mp_memory_free(buffer); return 0.0; } } static void decNumberAtan(decNumber *result, decNumber *x_orig, decContext *localset) { decNumber x; decNumberCopy(&x, x_orig); if (decNumberIsZero(&x)) { decNumberCopy(result, &x); } else { decNumber f, g, mx2, term; for (int i = 0; i<2; i++) { decNumber y; decNumberMultiply(&y, &x, &x, localset); decNumberAdd(&y, &y, &mp_decimal_data.one, localset); decNumberSquareRoot(&y, &y, localset); decNumberSubtract(&y, &y, &mp_decimal_data.one, localset); decNumberDivide(&x, &y, &x, localset); if (decNumberIsZero(&x)) { decNumberCopy(result, &x); return; } } decNumberCopy(&f, &x); decNumberCopy(&g, &mp_decimal_data.one); decNumberCopy(&term, &x); decNumberCopy(result, &x); decNumberMultiply(&mx2, &x, &x, localset); decNumberMinus (&mx2, &mx2, localset); for (int i = 0; i < 2 * localset->digits; i++) { decNumberMultiply(&f, &f, &mx2, localset); decNumberAdd(&g, &g, &mp_decimal_data.two_decNumber, localset); decNumberDivide(&term, &f, &g, localset); decNumberAdd(result, result, &term, localset); } decNumberAdd(result, result, result, localset); decNumberAdd(result, result, result, localset); } } static void decNumberAtan2(decNumber *result, decNumber *y, decNumber *x, decContext *localset) { if (! decNumberIsInfinite(x) && ! decNumberIsZero(y) && ! decNumberIsInfinite(y) && ! decNumberIsZero(x)) { decNumber temp; decNumberDivide(&temp, y, x, localset); decNumberAtan(result, &temp, localset); if (decNumberIsNegative(x)) { if (decNumberIsNegative(y)) { decNumberSubtract(result, result, &mp_decimal_data.PI_decNumber, localset); } else { decNumberAdd(result, result, &mp_decimal_data.PI_decNumber, localset); } } } else { if (decNumberIsInfinite(y) && decNumberIsInfinite(x)) { decNumberDivide(result, &mp_decimal_data.PI_decNumber, &mp_decimal_data.four_decNumber, localset); if (decNumberIsNegative(x) ) { decNumber a; decNumberFromDouble(&a, 3.0); decNumberMultiply(result, result, &a, localset); } } else if (!decNumberIsZero(y) && !decNumberIsInfinite(x) ) { decNumberDivide(result, &mp_decimal_data.PI_decNumber, &mp_decimal_data.two_decNumber, localset); } else { if (decNumberIsNegative(x)) { decNumberCopy(result, &mp_decimal_data.PI_decNumber); } else { decNumberZero(result); } } if (decNumberIsNegative(y)) { decNumberMinus(result, result, localset); } } } math_data *mp_initialize_decimal_math (MP mp) { math_data *math = (math_data *) mp_memory_allocate(sizeof(math_data)); decContextDefault(&mp_decimal_data.set, DEC_INIT_BASE); mp_decimal_data.set.traps = 0; decContextDefault(&mp_decimal_data.limitedset, DEC_INIT_BASE); mp_decimal_data.limitedset.traps = 0; mp_decimal_data.limitedset.emax = 999999; mp_decimal_data.limitedset.emin = -999999; mp_decimal_data.set.digits = DECPRECISION_DEFAULT; mp_decimal_data.limitedset.digits = DECPRECISION_DEFAULT; if (! mp_decimal_data.initialized) { mp_decimal_data.initialized = 1 ; decNumberFromInt32(&mp_decimal_data.one, 1); decNumberFromInt32(&mp_decimal_data.minusone, -1); decNumberFromInt32(&mp_decimal_data.zero, 0); decNumberFromInt32(&mp_decimal_data.two_decNumber, two); decNumberFromInt32(&mp_decimal_data.three_decNumber, three); decNumberFromInt32(&mp_decimal_data.four_decNumber, four); decNumberFromInt32(&mp_decimal_data.fraction_multiplier_decNumber, fraction_multiplier); decNumberFromInt32(&mp_decimal_data.fraction_one_decNumber, fraction_one); decNumberFromInt32(&mp_decimal_data.fraction_one_plus_decNumber, (fraction_one+1)); decNumberFromInt32(&mp_decimal_data.angle_multiplier_decNumber, angle_multiplier); decNumberFromString(&mp_decimal_data.PI_decNumber, PI_STRING, &mp_decimal_data.set); decNumberFromDouble(&mp_decimal_data.epsilon_decNumber, epsilon); decNumberFromString(&mp_decimal_data.EL_GORDO_decNumber, EL_GORDO, &mp_decimal_data.set); decNumberFromString(&mp_decimal_data.negative_EL_GORDO_decNumber, negative_EL_GORDO, &mp_decimal_data.set); mp_decimal_data.factorials = (decNumber **) mp_memory_allocate(FACTORIALS_CACHESIZE * sizeof(decNumber *)); mp_decimal_data.factorials[0] = (decNumber *) mp_memory_allocate(sizeof(decNumber)); decNumberCopy(mp_decimal_data.factorials[0], &mp_decimal_data.one); } 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); decNumberFromInt32( math->md_precision_default.data.num, DECPRECISION_DEFAULT); mp_allocate_number(mp, &math->md_precision_max, mp_scaled_type); decNumberFromInt32( math->md_precision_max.data.num, DECNUMDIGITS); mp_allocate_number(mp, &math->md_precision_min, mp_scaled_type); decNumberFromInt32( math->md_precision_min.data.num, 2); mp_allocate_number(mp, &math->md_epsilon_t, mp_scaled_type); decNumberCopy( math->md_epsilon_t.data.num, &mp_decimal_data.epsilon_decNumber); mp_allocate_number(mp, &math->md_inf_t, mp_scaled_type); decNumberCopy( math->md_inf_t.data.num, &mp_decimal_data.EL_GORDO_decNumber); mp_allocate_number(mp, &math->md_negative_inf_t, mp_scaled_type); decNumberCopy( math->md_negative_inf_t.data.num, &mp_decimal_data.negative_EL_GORDO_decNumber); mp_allocate_number(mp, &math->md_warning_limit_t, mp_scaled_type); decNumberFromString( math->md_warning_limit_t.data.num, warning_limit, &mp_decimal_data.set); mp_allocate_number(mp, &math->md_one_third_inf_t, mp_scaled_type); decNumberDivide( math->md_one_third_inf_t.data.num, math->md_inf_t.data.num, &mp_decimal_data.three_decNumber, &mp_decimal_data.set); mp_allocate_number(mp, &math->md_unity_t, mp_scaled_type); decNumberCopy( math->md_unity_t.data.num, &mp_decimal_data.one); mp_allocate_number(mp, &math->md_two_t, mp_scaled_type); decNumberFromInt32( math->md_two_t.data.num, two); mp_allocate_number(mp, &math->md_three_t, mp_scaled_type); decNumberFromInt32( math->md_three_t.data.num, three); mp_allocate_number(mp, &math->md_half_unit_t, mp_scaled_type); decNumberFromString( math->md_half_unit_t.data.num, "0.5", &mp_decimal_data.set); mp_allocate_number(mp, &math->md_three_quarter_unit_t, mp_scaled_type); decNumberFromString( math->md_three_quarter_unit_t.data.num, "0.75", &mp_decimal_data.set); mp_allocate_number(mp, &math->md_zero_t, mp_scaled_type); decNumberZero( math->md_zero_t.data.num); { decNumber fourzeroninesix; decNumberFromInt32(&fourzeroninesix, 4096); mp_allocate_number(mp, &math->md_arc_tol_k, mp_fraction_type); decNumberDivide( math->md_arc_tol_k.data.num, &mp_decimal_data.one, &fourzeroninesix, &mp_decimal_data.set); } mp_allocate_number(mp, &math->md_fraction_one_t, mp_fraction_type); decNumberFromInt32( math->md_fraction_one_t.data.num, fraction_one); mp_allocate_number(mp, &math->md_fraction_half_t, mp_fraction_type); decNumberFromInt32( math->md_fraction_half_t.data.num, fraction_half); mp_allocate_number(mp, &math->md_fraction_three_t, mp_fraction_type); decNumberFromInt32( math->md_fraction_three_t.data.num, fraction_three); mp_allocate_number(mp, &math->md_fraction_four_t, mp_fraction_type); decNumberFromInt32( math->md_fraction_four_t.data.num, fraction_four); mp_allocate_number(mp, &math->md_three_sixty_deg_t, mp_angle_type); decNumberFromInt32( math->md_three_sixty_deg_t.data.num, 360 * angle_multiplier); mp_allocate_number(mp, &math->md_one_eighty_deg_t, mp_angle_type); decNumberFromInt32( math->md_one_eighty_deg_t.data.num, 180 * angle_multiplier); mp_allocate_number(mp, &math->md_negative_one_eighty_deg_t, mp_angle_type); decNumberFromInt32( math->md_negative_one_eighty_deg_t.data.num, -180 * angle_multiplier); mp_allocate_number(mp, &math->md_one_k, mp_scaled_type); decNumberFromDouble( math->md_one_k.data.num, 1.0/64); mp_allocate_number(mp, &math->md_sqrt_8_e_k, mp_scaled_type); decNumberFromDouble( math->md_sqrt_8_e_k.data.num, 112428.82793 / 65536.0); mp_allocate_number(mp, &math->md_twelve_ln_2_k, mp_fraction_type); decNumberFromDouble( math->md_twelve_ln_2_k.data.num, 139548959.6165 / 65536.0); mp_allocate_number(mp, &math->md_coef_bound_k, mp_fraction_type); decNumberFromDouble( math->md_coef_bound_k.data.num,coef_bound); mp_allocate_number(mp, &math->md_coef_bound_minus_1, mp_fraction_type); decNumberFromDouble( math->md_coef_bound_minus_1.data.num,coef_bound - 1 / 65536.0); mp_allocate_number(mp, &math->md_twelvebits_3, mp_scaled_type); decNumberFromDouble( math->md_twelvebits_3.data.num, 1365 / 65536.0); mp_allocate_number(mp, &math->md_twentysixbits_sqrt2_t, mp_fraction_type); decNumberFromDouble( math->md_twentysixbits_sqrt2_t.data.num, 94906265.62 / 65536.0); mp_allocate_number(mp, &math->md_twentyeightbits_d_t, mp_fraction_type); decNumberFromDouble( math->md_twentyeightbits_d_t.data.num, 35596754.69 / 65536.0); mp_allocate_number(mp, &math->md_twentysevenbits_sqrt2_d_t, mp_fraction_type); decNumberFromDouble( math->md_twentysevenbits_sqrt2_d_t.data.num, 25170706.63 / 65536.0); mp_allocate_number(mp, &math->md_fraction_threshold_t, mp_fraction_type); decNumberFromDouble( math->md_fraction_threshold_t.data.num, fraction_threshold); mp_allocate_number(mp, &math->md_half_fraction_threshold_t, mp_fraction_type); decNumberFromDouble( math->md_half_fraction_threshold_t.data.num, half_fraction_threshold); mp_allocate_number(mp, &math->md_scaled_threshold_t, mp_scaled_type); decNumberFromDouble( math->md_scaled_threshold_t.data.num, scaled_threshold); mp_allocate_number(mp, &math->md_half_scaled_threshold_t, mp_scaled_type); decNumberFromDouble( math->md_half_scaled_threshold_t.data.num, half_scaled_threshold); mp_allocate_number(mp, &math->md_near_zero_angle_t, mp_angle_type); decNumberFromDouble( math->md_near_zero_angle_t.data.num, near_zero_angle); mp_allocate_number(mp, &math->md_p_over_v_threshold_t, mp_fraction_type); decNumberFromDouble( math->md_p_over_v_threshold_t.data.num, p_over_v_threshold); mp_allocate_number(mp, &math->md_equation_threshold_t, mp_scaled_type); decNumberFromDouble( math->md_equation_threshold_t.data.num, equation_threshold); math->md_from_int = mp_set_decimal_from_int; math->md_from_boolean = mp_set_decimal_from_boolean; math->md_from_scaled = mp_set_decimal_from_scaled; math->md_from_double = mp_set_decimal_from_double; math->md_from_addition = mp_set_decimal_from_addition; math->md_half_from_addition = mp_set_decimal_half_from_addition; math->md_from_subtraction = mp_set_decimal_from_subtraction; math->md_half_from_subtraction = mp_set_decimal_half_from_subtraction; math->md_from_oftheway = mp_set_decimal_from_of_the_way; math->md_from_div = mp_set_decimal_from_div; math->md_from_mul = mp_set_decimal_from_mul; math->md_from_int_div = mp_set_decimal_from_int_div; math->md_from_int_mul = mp_set_decimal_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_decimal_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_decimal_fraction_to_round_scaled; math->md_make_scaled = mp_decimal_number_make_scaled; math->md_make_fraction = mp_decimal_number_make_fraction; math->md_take_fraction = mp_decimal_number_take_fraction; math->md_take_scaled = mp_decimal_number_take_scaled; math->md_velocity = mp_decimal_velocity; math->md_n_arg = mp_decimal_n_arg; math->md_m_log = mp_decimal_m_log; math->md_m_exp = mp_decimal_m_exp; math->md_m_unif_rand = mp_decimal_m_unif_rand; math->md_m_norm_rand = mp_decimal_m_norm_rand; math->md_pyth_add = mp_decimal_pyth_add; math->md_pyth_sub = mp_decimal_pyth_sub; math->md_power_of = mp_decimal_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_decimal_sin_cos; math->md_slow_add = mp_decimal_slow_add; math->md_sqrt = mp_decimal_square_rt; math->md_print = mp_decimal_print_number; math->md_tostring = mp_decimal_number_tostring; math->md_modulo = mp_decimal_number_modulo; math->md_ab_vs_cd = mp_ab_vs_cd; math->md_crossing_point = mp_decimal_crossing_point; math->md_scan_numeric = mp_decimal_scan_numeric_token; math->md_scan_fractional = mp_decimal_scan_fractional_token; math->md_free_math = mp_free_decimal_math; math->md_set_precision = mp_decimal_set_precision; return math; } void mp_decimal_set_precision (MP mp) { int i = decNumberToInt32((decNumber *) internal_value(mp_number_precision_internal).data.num, &mp_decimal_data.set); mp_decimal_data.set.digits = i; mp_decimal_data.limitedset.digits = i; } void mp_free_decimal_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.num = mp_memory_allocate(sizeof(decNumber)); n->type = t; decNumberZero(n->data.num); } void mp_allocate_clone (MP mp, mp_number *n, mp_number_type t, mp_number *v) { (void) mp; n->data.num = mp_memory_allocate(sizeof(decNumber)); n->type = t; decNumberZero(n->data.num); decNumberCopy(n->data.num, v->data.num); } void mp_allocate_abs (MP mp, mp_number *n, mp_number_type t, mp_number *v) { (void) mp; n->data.num = mp_memory_allocate(sizeof(decNumber)); n->type = t; decNumberZero(n->data.num); decNumberAbs(n->data.num, v->data.num, &mp_decimal_data.set); } void mp_allocate_double (MP mp, mp_number *n, double v) { (void) mp; n->data.num = mp_memory_allocate(sizeof(decNumber)); n->type = mp_scaled_type; decNumberZero(n->data.num); decNumberFromDouble(n->data.num, v); } void mp_free_number (MP mp, mp_number *n) { (void) mp; if (n->data.num) { mp_memory_free(n->data.num); n->data.num = NULL; n->type = mp_nan_type; } } void mp_set_decimal_from_int(mp_number *A, int B) { decNumberFromInt32(A->data.num, B); } void mp_set_decimal_from_boolean(mp_number *A, int B) { decNumberFromInt32(A->data.num, B); } void mp_set_decimal_from_scaled(mp_number *A, int B) { decNumber c; decNumberFromInt32(&c, 65536); decNumberFromInt32(A->data.num,B); decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set); } void mp_set_decimal_from_double(mp_number *A, double B) { decNumberFromDouble(A->data.num, B); } void mp_set_decimal_from_addition(mp_number *A, mp_number *B, mp_number *C) { decNumberAdd(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set); } void mp_set_decimal_half_from_addition(mp_number *A, mp_number *B, mp_number *C) { decNumber c; decNumberAdd(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set); decNumberFromInt32(&c, 2); decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set); } void mp_set_decimal_from_subtraction(mp_number *A, mp_number *B, mp_number *C) { decNumberSubtract(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set); } void mp_set_decimal_half_from_subtraction(mp_number *A, mp_number *B, mp_number *C) { decNumber c; decNumberSubtract(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set); decNumberFromInt32(&c, 2); decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set); } void mp_set_decimal_from_div(mp_number *A, mp_number *B, mp_number *C) { decNumberDivide(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set); } void mp_set_decimal_from_mul(mp_number *A, mp_number *B, mp_number *C) { decNumberMultiply(A->data.num, B->data.num, C->data.num, &mp_decimal_data.set); } void mp_set_decimal_from_int_div(mp_number *A, mp_number *B, int C) { decNumber c; decNumberFromInt32(&c, C); decNumberDivide(A->data.num, B->data.num, &c, &mp_decimal_data.set); } void mp_set_decimal_from_int_mul(mp_number *A, mp_number *B, int C) { decNumber c; decNumberFromInt32(&c, C); decNumberMultiply(A->data.num, B->data.num, &c, &mp_decimal_data.set); } void mp_set_decimal_from_of_the_way (MP mp, mp_number *A, mp_number *t, mp_number *B, mp_number *C) { decNumber c; decNumber r1; decNumberSubtract(&c, B->data.num, C->data.num, &mp_decimal_data.set); mp_decimal_take_fraction(mp, &r1, &c, t->data.num); decNumberSubtract(A->data.num, B->data.num, &r1, &mp_decimal_data.set); mp_decnumber_check(mp, A->data.num, &mp_decimal_data.set); } void mp_number_negate(mp_number *A) { decNumberCopyNegate(A->data.num, A->data.num); checkZero(A->data.num); } void mp_number_add(mp_number *A, mp_number *B) { decNumberAdd(A->data.num, A->data.num, B->data.num, &mp_decimal_data.set); } void mp_number_subtract(mp_number *A, mp_number *B) { decNumberSubtract(A->data.num, A->data.num, B->data.num, &mp_decimal_data.set); } void mp_number_half(mp_number *A) { decNumber c; decNumberFromInt32(&c, 2); decNumberDivide(A->data.num, A->data.num, &c, &mp_decimal_data.set); } void mp_number_double(mp_number *A) { decNumber c; decNumberFromInt32(&c, 2); decNumberMultiply(A->data.num, A->data.num, &c, &mp_decimal_data.set); } void mp_number_add_scaled(mp_number *A, int B) { decNumber b, c; decNumberFromInt32(&c, 65536); decNumberFromInt32(&b, B); decNumberDivide(&b, &b, &c, &mp_decimal_data.set); decNumberAdd(A->data.num, A->data.num, &b, &mp_decimal_data.set); } void mp_number_multiply_int(mp_number *A, int B) { decNumber b; decNumberFromInt32(&b, B); decNumberMultiply(A->data.num, A->data.num, &b, &mp_decimal_data.set); } void mp_number_divide_int(mp_number *A, int B) { decNumber b; decNumberFromInt32(&b, B); decNumberDivide(A->data.num, A->data.num, &b, &mp_decimal_data.set); } void mp_decimal_abs(mp_number *A) { decNumberAbs(A->data.num, A->data.num, &mp_decimal_data.set); } void mp_number_clone(mp_number *A, mp_number *B) { decNumberCopy(A->data.num, B->data.num); } void mp_number_negated_clone(mp_number *A, mp_number *B) { decNumberCopyNegate(A->data.num, B->data.num); checkZero(A->data.num); } void mp_number_abs_clone(mp_number *A, mp_number *B) { decNumberAbs(A->data.num, B->data.num, &mp_decimal_data.set); } void mp_number_swap(mp_number *A, mp_number *B) { decNumber swap_tmp; decNumberCopy(&swap_tmp, A->data.num); decNumberCopy(A->data.num, B->data.num); decNumberCopy(B->data.num, &swap_tmp); } void mp_number_fraction_to_scaled(mp_number *A) { A->type = mp_scaled_type; decNumberDivide(A->data.num, A->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } void mp_number_angle_to_scaled(mp_number *A) { A->type = mp_scaled_type; decNumberDivide(A->data.num, A->data.num, &mp_decimal_data.angle_multiplier_decNumber, &mp_decimal_data.set); } void mp_number_scaled_to_fraction(mp_number *A) { A->type = mp_fraction_type; decNumberMultiply(A->data.num, A->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } void mp_number_scaled_to_angle(mp_number *A) { A->type = mp_angle_type; decNumberMultiply(A->data.num, A->data.num, &mp_decimal_data.angle_multiplier_decNumber, &mp_decimal_data.set); } int mp_number_to_scaled(mp_number *A) { int result; decNumber corrected; decNumberFromInt32(&corrected, 65536); decNumberMultiply(&corrected, &corrected, A->data.num, &mp_decimal_data.set); decNumberReduce(&corrected, &corrected, &mp_decimal_data.set); result = lround(decNumberToDouble(&corrected)); return result; } int mp_number_to_int(mp_number *A) { int result; mp_decimal_data.set.status = 0; result = decNumberToInt32(A->data.num, &mp_decimal_data.set); if (mp_decimal_data.set.status == DEC_Invalid_operation) { mp_decimal_data.set.status = 0; return 0; } else { return result; } } int mp_number_to_boolean(mp_number *A) { unsigned int result; mp_decimal_data.set.status = 0; result = decNumberToUInt32(A->data.num, &mp_decimal_data.set); if (mp_decimal_data.set.status == DEC_Invalid_operation) { mp_decimal_data.set.status = 0; return mp_false_operation; } else { return result ; } } double mp_number_to_double(mp_number *A) { char *buffer = mp_memory_allocate((size_t) ((decNumber *) A->data.num)->digits + 14); double res = 0.0; decNumberToString(A->data.num, buffer); if (sscanf(buffer, "%lf", &res)) { mp_memory_free(buffer); return res; } else { mp_memory_free(buffer); return 0.0; } } int mp_number_odd(mp_number *A) { decNumber r1, r2; decNumberAbs(&r1, A->data.num, &mp_decimal_data.set); decNumberRemainder(&r2, &r1, &mp_decimal_data.two_decNumber, &mp_decimal_data.set); decNumberCompare(&r1, &r2, &mp_decimal_data.one, &mp_decimal_data.set); return decNumberIsZero(&r1); } int mp_number_equal(mp_number *A, mp_number *B) { decNumber res; decNumberCompare(&res, A->data.num, B->data.num, &mp_decimal_data.set); return decNumberIsZero(&res); } int mp_number_greater(mp_number *A, mp_number *B) { decNumber res; decNumberCompare(&res, A->data.num, B->data.num, &mp_decimal_data.set); return decNumberIsPositive(&res); } int mp_number_less(mp_number *A, mp_number *B) { decNumber res; decNumberCompare(&res, A->data.num, B->data.num, &mp_decimal_data.set); return decNumberIsNegative(&res); } int mp_number_nonequalabs(mp_number *A, mp_number *B) { decNumber res, a, b; decNumberCopyAbs(&a, A->data.num); decNumberCopyAbs(&b, B->data.num); decNumberCompare(&res, &a, &b, &mp_decimal_data.set); return ! decNumberIsZero(&res); } char *mp_decnumber_tostring(decNumber *n) { decNumber corrected; char *buffer = mp_memory_allocate((size_t) ((decNumber *) n)->digits + 14); decNumberCopy(&corrected, n); decNumberTrim(&corrected); decNumberToString(&corrected, buffer); return buffer; } char *mp_decimal_number_tostring (MP mp, mp_number *n) { (void) mp; return mp_decnumber_tostring(n->data.num); } void mp_decimal_print_number (MP mp, mp_number *n) { char *str = mp_decnumber_tostring(n->data.num); mp_print_e_str(mp, str); mp_memory_free(str); } void mp_decimal_slow_add (MP mp, mp_number *ret, mp_number *A, mp_number *B) { (void) mp; decNumberAdd(ret->data.num, A->data.num, B->data.num, &mp_decimal_data.set); } void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) { decNumberDivide(ret, p, q, &mp_decimal_data.set); mp_decnumber_check(mp, ret, &mp_decimal_data.set); decNumberMultiply(ret, ret, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } void mp_decimal_number_make_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) { mp_decimal_make_fraction(mp, ret->data.num, p->data.num, q->data.num); } void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) { (void) mp; decNumberMultiply(ret, p, q, &mp_decimal_data.set); decNumberDivide(ret, ret, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } void mp_decimal_number_take_fraction (MP mp, mp_number *ret, mp_number *p, mp_number *q) { mp_decimal_take_fraction(mp, ret->data.num, p->data.num, q->data.num); } void mp_decimal_number_take_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) { (void) mp; decNumberMultiply(ret->data.num, p_orig->data.num, q_orig->data.num, &mp_decimal_data.set); } void mp_decimal_number_make_scaled (MP mp, mp_number *ret, mp_number *p_orig, mp_number *q_orig) { decNumberDivide(ret->data.num, p_orig->data.num, q_orig->data.num, &mp_decimal_data.set); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } void mp_wrapup_numeric_token (MP mp, unsigned char *start, unsigned char *stop) { decNumber result; size_t l = stop-start+1; char *buf = mp_memory_allocate(l + 1); buf[l] = '\0'; (void) strncpy(buf, (const char *) start, l); mp_decimal_data.set.status = 0; decNumberFromString(&result,buf, &mp_decimal_data.set); mp_memory_free(buf); if (mp_decimal_data.set.status == 0) { set_cur_mod(&result); } else if (mp->scanner_status != mp_tex_flushing_state) { if (too_large(mp_decimal_data.set.status)) { mp_decnumber_check(mp, &result, &mp_decimal_data.set); set_cur_mod(&result); mp_error( mp, "Enormous number has been reduced", "I could not handle this number specification because it is out of range." ); } else if (too_precise(mp_decimal_data.set.status)) { set_cur_mod(&result); if (decNumberIsPositive((decNumber *) internal_value(mp_warning_check_internal).data.num) && (mp->scanner_status != mp_tex_flushing_state)) { char msg[256]; mp_snprintf (msg, 256, "Number is too precise (numberprecision = %d)", mp_decimal_data.set.digits); mp_error( mp, msg, "Continue and I'll round the value until it fits the current numberprecision\n" "(Set warningcheck:=0 to suppress this message.)" ); } } else { mp_error( mp, "Erroneous number specification changed to zero", "I could not handle this number specification" ); decNumberZero(&result); set_cur_mod(&result); } } set_cur_cmd((mp_variable_type) mp_numeric_command); } static void 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_decimal_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++; } find_exponent(mp); stop = &mp->buffer[mp->cur_input.loc_field-1]; mp_wrapup_numeric_token(mp, start, stop); } void mp_decimal_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++; } } find_exponent(mp); stop = &mp->buffer[mp->cur_input.loc_field-1]; mp_wrapup_numeric_token(mp, start, stop); } void mp_decimal_velocity (MP mp, mp_number *ret, mp_number *st, mp_number *ct, mp_number *sf, mp_number *cf, mp_number *t) { decNumber acc, num, denom; decNumber r1, r2; decNumber arg1, arg2; decNumber i16, fone, fhalf, ftwo, sqrtfive; decNumberFromInt32(&i16, 16); decNumberFromInt32(&fone, fraction_one); decNumberFromInt32(&fhalf, fraction_half); decNumberFromInt32(&ftwo, fraction_two); decNumberFromInt32(&sqrtfive, 5); decNumberSquareRoot(&sqrtfive, &sqrtfive, &mp_decimal_data.set); decNumberDivide(&arg1, sf->data.num, &i16, &mp_decimal_data.set); decNumberSubtract(&arg1, st->data.num,&arg1, &mp_decimal_data.set); decNumberDivide(&arg2, st->data.num, &i16, &mp_decimal_data.set); decNumberSubtract(&arg2, sf->data.num,&arg2, &mp_decimal_data.set); mp_decimal_take_fraction(mp, &acc, &arg1, &arg2); decNumberCopy(&arg1, &acc); decNumberSubtract(&arg2, ct->data.num, cf->data.num, &mp_decimal_data.set); mp_decimal_take_fraction(mp, &acc, &arg1, &arg2); decNumberSquareRoot(&arg1, &mp_decimal_data.two_decNumber, &mp_decimal_data.set); decNumberMultiply(&arg1, &arg1, &fone, &mp_decimal_data.set); mp_decimal_take_fraction(mp, &r1, &acc, &arg1); decNumberAdd(&num, &ftwo, &r1, &mp_decimal_data.set); decNumberSubtract(&arg1,&sqrtfive, &mp_decimal_data.one, &mp_decimal_data.set); decNumberMultiply(&arg1,&arg1,&fhalf, &mp_decimal_data.set); decNumberMultiply(&arg1,&arg1,&mp_decimal_data.three_decNumber, &mp_decimal_data.set); decNumberSubtract(&arg2,&mp_decimal_data.three_decNumber, &sqrtfive, &mp_decimal_data.set); decNumberMultiply(&arg2,&arg2, &fhalf, &mp_decimal_data.set); decNumberMultiply(&arg2,&arg2, &mp_decimal_data.three_decNumber, &mp_decimal_data.set); mp_decimal_take_fraction(mp, &r1, ct->data.num, &arg1) ; mp_decimal_take_fraction(mp, &r2, cf->data.num, &arg2); decNumberFromInt32(&denom, fraction_three); decNumberAdd(&denom, &denom, &r1, &mp_decimal_data.set); decNumberAdd(&denom, &denom, &r2, &mp_decimal_data.set); decNumberCompare(&arg1, t->data.num, &mp_decimal_data.one, &mp_decimal_data.set); if (! decNumberIsZero(&arg1)) { decNumberDivide(&num, &num, t->data.num, &mp_decimal_data.set); } decNumberCopy(&r2, &num); decNumberDivide(&r2, &r2, &mp_decimal_data.four_decNumber, &mp_decimal_data.set); if (decNumberLess(&denom, &r2)) { decNumberFromInt32(ret->data.num, fraction_four); } else { mp_decimal_make_fraction(mp, ret->data.num, &num, &denom); } mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } int mp_ab_vs_cd (mp_number *a_orig, mp_number *b_orig, mp_number *c_orig, mp_number *d_orig) { decNumber a, b, c, d; decNumber ab, cd; decNumberCopy(&a, (decNumber *) a_orig->data.num); decNumberCopy(&b, (decNumber *) b_orig->data.num); decNumberCopy(&c, (decNumber *) c_orig->data.num); decNumberCopy(&d, (decNumber *) d_orig->data.num); decNumberMultiply(&ab, (decNumber *) a_orig->data.num, (decNumber *)b_orig->data.num, &mp_decimal_data.set); decNumberMultiply(&cd, (decNumber *) c_orig->data.num, (decNumber *)d_orig->data.num, &mp_decimal_data.set); if (decNumberLess(&ab, &cd)) { return -1; } else if (decNumberGreater(&ab, &cd)) { return 1; } else { return 0; } } static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number *aa, mp_number *bb, mp_number *cc) { decNumber a, b, c; double d; decNumber x, xx, x0, x1, x2; decNumber scratch, scratch2; decNumberCopy(&a, (decNumber *) aa->data.num); decNumberCopy(&b, (decNumber *) bb->data.num); decNumberCopy(&c, (decNumber *) cc->data.num); if (decNumberIsNegative(&a)) { decNumberCopy(ret->data.num, &zero_crossing); goto RETURN; } if (! decNumberIsNegative(&c)) { if (! decNumberIsNegative(&b)) { if (decNumberIsPositive(&c)) { decNumberCopy(ret->data.num, &no_crossing); } else if (decNumberIsZero(&a) && decNumberIsZero(&b)) { decNumberCopy(ret->data.num, &no_crossing); } else { decNumberCopy(ret->data.num, &one_crossing); } goto RETURN; } if (decNumberIsZero(&a)) { decNumberCopy(ret->data.num, &zero_crossing); goto RETURN; } } else if (decNumberIsZero(&a) && ! decNumberIsPositive(&b)) { decNumberCopy(ret->data.num, &zero_crossing); goto RETURN; } d = epsilonf; decNumberCopy(&x0, &a); decNumberSubtract(&x1, &a, &b, &mp_decimal_data.set); decNumberSubtract(&x2, &b, &c, &mp_decimal_data.set); decNumberFromDouble(&scratch2, 1E-12); do { decNumberAdd(&x, &x1, &x2, &mp_decimal_data.set); decNumberDivide(&x, &x, &mp_decimal_data.two_decNumber, &mp_decimal_data.set); decNumberAdd(&x, &x, &scratch2, &mp_decimal_data.set); decNumberSubtract(&scratch, &x1, &x0, &mp_decimal_data.set); if (decNumberGreater(&scratch, &x0)) { decNumberCopy(&x2, &x); decNumberAdd(&x0, &x0, &x0, &mp_decimal_data.set); d += d; } else { decNumberAdd(&xx, &scratch, &x, &mp_decimal_data.set); if (decNumberGreater(&xx,&x0)) { decNumberCopy(&x2,&x); decNumberAdd(&x0, &x0, &x0, &mp_decimal_data.set); d += d; } else { decNumberSubtract(&x0, &x0, &xx, &mp_decimal_data.set); if (! decNumberGreater(&x,&x0)) { decNumberAdd(&scratch, &x, &x2, &mp_decimal_data.set); if (! decNumberGreater(&scratch, &x0)) { decNumberCopy(ret->data.num, &no_crossing); goto RETURN; } } decNumberCopy(&x1,&x); d = d + d + epsilonf; } } } while (d < fraction_one); decNumberFromDouble(&scratch, d); decNumberSubtract(ret->data.num,&scratch, &mp_decimal_data.fraction_one_decNumber, &mp_decimal_data.set); RETURN: mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } int mp_round_unscaled(mp_number *x_orig) { return (int) lround(mp_number_to_double(x_orig)); } void mp_number_floor(mp_number *i) { int round = mp_decimal_data.set.round; mp_decimal_data.set.round = DEC_ROUND_FLOOR; decNumberToIntegralValue(i->data.num, i->data.num, &mp_decimal_data.set); mp_decimal_data.set.round = round; } void mp_decimal_fraction_to_round_scaled(mp_number *x_orig) { x_orig->type = mp_scaled_type; decNumberDivide(x_orig->data.num, x_orig->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number *x_orig) { decNumber x; decNumberCopy(&x, x_orig->data.num); if (! decNumberIsPositive(&x)) { if (decNumberIsNegative(&x)) { char msg[256]; char *xstr = mp_decimal_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." ); } decNumberZero(ret->data.num); } else { decNumberSquareRoot(ret->data.num, &x, &mp_decimal_data.set); } mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } void mp_decimal_pyth_add (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) { decNumber a, b; decNumber asq, bsq; decNumberCopyAbs(&a, a_orig->data.num); decNumberCopyAbs(&b, b_orig->data.num); decNumberMultiply(&asq, &a, &a, &mp_decimal_data.set); decNumberMultiply(&bsq, &b, &b, &mp_decimal_data.set); decNumberAdd(&a, &asq, &bsq, &mp_decimal_data.set); decNumberSquareRoot(ret->data.num, &a, &mp_decimal_data.set); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } void mp_decimal_pyth_sub (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) { decNumber a, b; decNumberCopyAbs(&a, a_orig->data.num); decNumberCopyAbs(&b, b_orig->data.num); if (! decNumberGreater(&a, &b)) { if (decNumberLess(&a, &b)) { char msg[256]; char *astr = mp_decimal_number_tostring(mp, a_orig); char *bstr = mp_decimal_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, I'm zeroing this one.\n" "Proceed, with fingers crossed." ); } decNumberZero(&a); } else { decNumber asq, bsq; decNumberMultiply(&asq, &a, &a, &mp_decimal_data.set); decNumberMultiply(&bsq, &b, &b, &mp_decimal_data.set); decNumberSubtract(&a, &asq, &bsq, &mp_decimal_data.set); decNumberSquareRoot(&a, &a, &mp_decimal_data.set); } decNumberCopy(ret->data.num, &a); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } void mp_decimal_power_of (MP mp, mp_number *ret, mp_number *a_orig, mp_number *b_orig) { decNumberPower(ret->data.num, a_orig->data.num, b_orig->data.num, &mp_decimal_data.set); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } void mp_decimal_m_log (MP mp, mp_number *ret, mp_number *x_orig) { if (! decNumberIsPositive((decNumber *) x_orig->data.num)) { char msg[256]; char *xstr = mp_decimal_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." ); decNumberZero(ret->data.num); } else { decNumber twofivesix; decNumberFromInt32(&twofivesix, 256); decNumberLn(ret->data.num, x_orig->data.num, &mp_decimal_data.limitedset); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.limitedset); decNumberMultiply(ret->data.num, ret->data.num, &twofivesix, &mp_decimal_data.set); } mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number *x_orig) { decNumber temp, twofivesix; decNumberFromInt32(&twofivesix, 256); decNumberDivide(&temp, x_orig->data.num, &twofivesix, &mp_decimal_data.set); mp_decimal_data.limitedset.status = 0; decNumberExp(ret->data.num, &temp, &mp_decimal_data.limitedset); if (mp_decimal_data.limitedset.status & DEC_Clamped) { if (decNumberIsPositive((decNumber *) x_orig->data.num)) { mp->arith_error = 1; decNumberCopy(ret->data.num, &mp_decimal_data.EL_GORDO_decNumber); } else { decNumberZero(ret->data.num); } } mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.limitedset); mp_decimal_data.limitedset.status = 0; } void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number *x_orig, mp_number *y_orig) { if (decNumberIsZero((decNumber *) x_orig->data.num) && decNumberIsZero((decNumber *) y_orig->data.num)) { 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." ); decNumberZero(ret->data.num); } else { decNumber atan2val, oneeighty_angle; ret->type = mp_angle_type; decNumberFromInt32(&oneeighty_angle, 180 * angle_multiplier); decNumberDivide(&oneeighty_angle, &oneeighty_angle, &mp_decimal_data.PI_decNumber, &mp_decimal_data.set); checkZero(y_orig->data.num); checkZero(x_orig->data.num); decNumberAtan2(&atan2val, y_orig->data.num, x_orig->data.num, &mp_decimal_data.set); decNumberMultiply(ret->data.num,&atan2val, &oneeighty_angle, &mp_decimal_data.set); checkZero(ret->data.num); } mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } static void sinecosine(decNumber *theangle, decNumber *c, decNumber *s) { int prec = mp_decimal_data.set.digits/2; decNumber p, pxa, fac, cc; decNumber n1, n2, p1; decNumberZero(c); decNumberZero(s); if (prec < DECPRECISION_DEFAULT) { prec = DECPRECISION_DEFAULT; } for (int n = 0; n < prec; n++) { decNumberFromInt32(&p1, n); decNumberFromInt32(&n1, 2*n); decNumberPower(&p, &mp_decimal_data.minusone, &p1, &mp_decimal_data.limitedset); if (n == 0) { decNumberCopy(&pxa, &mp_decimal_data.one); } else { decNumberPower(&pxa, theangle, &n1, &mp_decimal_data.limitedset); } if (2*n < mp_decimal_data.last_cached_factorial) { decNumberCopy(&fac,mp_decimal_data.factorials[2*n]); } else { decNumberCopy(&fac,mp_decimal_data.factorials[mp_decimal_data.last_cached_factorial]); for (int i = mp_decimal_data.last_cached_factorial+1; i <= 2*n; i++) { decNumberFromInt32(&cc, i); decNumberMultiply (&fac, &fac, &cc, &mp_decimal_data.set); if (i < FACTORIALS_CACHESIZE) { mp_decimal_data.factorials[i] = mp_memory_allocate(sizeof(decNumber)); decNumberCopy(mp_decimal_data.factorials[i], &fac); mp_decimal_data.last_cached_factorial = i; } } } decNumberDivide(&pxa, &pxa, &fac, &mp_decimal_data.set); decNumberMultiply(&pxa, &pxa, &p, &mp_decimal_data.set); decNumberAdd(s, s, &pxa, &mp_decimal_data.set); decNumberFromInt32(&n2, 2*n+1); decNumberMultiply(&fac, &fac, &n2, &mp_decimal_data.set); decNumberPower(&pxa, theangle, &n2, &mp_decimal_data.limitedset); decNumberDivide(&pxa, &pxa, &fac, &mp_decimal_data.set); decNumberMultiply(&pxa, &pxa, &p, &mp_decimal_data.set); decNumberAdd(c, c, &pxa, &mp_decimal_data.set); } } void mp_decimal_sin_cos (MP mp, mp_number *z_orig, mp_number *n_cos, mp_number *n_sin) { decNumber rad; decNumber one_eighty; double tmp = mp_number_to_double(z_orig)/16.0; if ((tmp == 90.0)||(tmp == -270)){ decNumberZero(n_cos->data.num); decNumberCopy(n_sin->data.num, &mp_decimal_data.fraction_multiplier_decNumber); } else if ((tmp == -90.0)||(tmp == 270.0)) { decNumberZero(n_cos->data.num); decNumberCopyNegate(n_sin->data.num, &mp_decimal_data.fraction_multiplier_decNumber); } else if ((tmp == 180.0) || (tmp == -180.0)) { decNumberCopyNegate(n_cos->data.num, &mp_decimal_data.fraction_multiplier_decNumber); decNumberZero(n_sin->data.num); } else { decNumberFromInt32(&one_eighty, 180 * 16); decNumberMultiply(&rad, z_orig->data.num, &mp_decimal_data.PI_decNumber, &mp_decimal_data.set); decNumberDivide(&rad, &rad, &one_eighty, &mp_decimal_data.set); sinecosine(&rad, n_sin->data.num, n_cos->data.num); decNumberMultiply(n_cos->data.num, n_cos->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); decNumberMultiply(n_sin->data.num, n_sin->data.num, &mp_decimal_data.fraction_multiplier_decNumber, &mp_decimal_data.set); } mp_decnumber_check(mp, n_cos->data.num, &mp_decimal_data.set); mp_decnumber_check(mp, n_sin->data.num, &mp_decimal_data.set); } # define KK 100 # define LL 37 # define MM (1L<<30) # define mod_diff(x,y) (((x)-(y))&(MM-1)) # define QUALITY 1009 # define TT 70 # define is_odd(x) ((x)&1) typedef struct mp_decimal_random_info { long x[KK]; long buf[QUALITY]; long dummy; long started; long *ptr; } mp_decimal_random_info; static mp_decimal_random_info mp_decimal_random_data = { .dummy = -1, .started = -1, .ptr = &mp_decimal_random_data.dummy }; static void ran_array(long aa[],int n) { int i, j; for (j = 0; j < KK;j++) { aa[j] = mp_decimal_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_decimal_random_data.x[i] = mod_diff(aa[j - KK], aa[j - LL]); } for (;i < KK; i++, j++) { mp_decimal_random_data.x[i] = mod_diff(aa[j - KK], mp_decimal_random_data.x[i - LL]); } } static void 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_decimal_random_data.x[j + KK -LL] = x[j]; } for (; j < KK; j++) { mp_decimal_random_data.x[j - LL] = x[j]; } for (j = 0; j < 10; j++) { ran_array(x, KK + KK - 1); } mp_decimal_random_data.ptr = &mp_decimal_random_data.started; } static long ran_arr_cycle(void) { if (mp_decimal_random_data.ptr == &mp_decimal_random_data.dummy) { ran_start(314159L); } ran_array(mp_decimal_random_data.buf, QUALITY); mp_decimal_random_data.buf[KK] = -1; mp_decimal_random_data.ptr = mp_decimal_random_data.buf + 1; return mp_decimal_random_data.buf[0]; } 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; } decNumberFromInt32(mp->randoms[(i * 21) % 55].data.num, j); } mp_new_randoms(mp); mp_new_randoms(mp); mp_new_randoms(mp); ran_start((unsigned long) seed); } void mp_decimal_number_modulo(mp_number *a, mp_number *b) { decNumberRemainder(a->data.num, a->data.num, b->data.num, &mp_decimal_data.set); } static void mp_next_unif_random (MP mp, mp_number *ret) { decNumber a; decNumber b; unsigned long int op = (unsigned) (*mp_decimal_random_data.ptr>=0? *mp_decimal_random_data.ptr++: ran_arr_cycle()); (void) mp; decNumberFromInt32(&a, op); decNumberFromInt32(&b, MM); decNumberDivide(&a, &a, &b, &mp_decimal_data.set); decNumberCopy(ret->data.num, &a); mp_decnumber_check(mp, ret->data.num, &mp_decimal_data.set); } 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_decimal_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); decNumberMultiply(y.data.num, abs_x.data.num, u.data.num, &mp_decimal_data.set); 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, &x); mp_free_number(mp, &abs_x); mp_free_number(mp, &y); mp_free_number(mp, &u); } static void mp_decimal_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_decimal_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_decimal_abs(&abs_x); } while (! mp_number_less(&abs_x, &u)); mp_decimal_number_make_fraction(mp, &r, &xa, &u); mp_number_clone(&xa, &r); mp_decimal_m_log(mp, &la, &u); mp_set_decimal_from_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); }