diff options
author | Denis Bitouzé <dbitouze@wanadoo.fr> | 2021-02-25 18:23:07 +0000 |
---|---|---|
committer | Denis Bitouzé <dbitouze@wanadoo.fr> | 2021-02-25 18:23:07 +0000 |
commit | c6101f91d071883b48b1b4b51e5eba0f36d9a78d (patch) | |
tree | 1bf7f5a881d7a4f5c5bf59d0b2821943dd822372 /Build/source/texk/web2c/mplibdir/mpmathdecimal.w | |
parent | 07ee7222e389b0777456b427a55c22d0e6ffd267 (diff) |
French translation for tlmgr updated
git-svn-id: svn://tug.org/texlive/trunk@57912 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Build/source/texk/web2c/mplibdir/mpmathdecimal.w')
-rw-r--r-- | Build/source/texk/web2c/mplibdir/mpmathdecimal.w | 2011 |
1 files changed, 0 insertions, 2011 deletions
diff --git a/Build/source/texk/web2c/mplibdir/mpmathdecimal.w b/Build/source/texk/web2c/mplibdir/mpmathdecimal.w deleted file mode 100644 index 5c2a8fe624e..00000000000 --- a/Build/source/texk/web2c/mplibdir/mpmathdecimal.w +++ /dev/null @@ -1,2011 +0,0 @@ -% $Id$ -% -% This file is part of MetaPost; -% the MetaPost program is in the public domain. -% See the <Show version...> code in mpost.w for more info. - -% Here is TeX material that gets inserted after \input webmac - -\font\tenlogo=logo10 % font used for the METAFONT logo -\font\logos=logosl10 -\def\MF{{\tenlogo META}\-{\tenlogo FONT}} -\def\MP{{\tenlogo META}\-{\tenlogo POST}} -\def\pct!{{\char`\%}} % percent sign in ordinary text -\def\psqrt#1{\sqrt{\mathstrut#1}} - - -\def\title{Math support functions for decNumber based math} -\pdfoutput=1 - -@ Introduction. - -@c -#include <w2c/config.h> -#include <stdio.h> -#include <stdlib.h> -#include <string.h> -#include <math.h> -#include "mpmathdecimal.h" /* internal header */ -#define ROUND(a) floor((a)+0.5) -@h - -@ @c -@<Declarations@>; - -@ @(mpmathdecimal.h@>= -#ifndef MPMATHDECIMAL_H -#define MPMATHDECIMAL_H 1 -#include "mplib.h" -#include "mpmp.h" /* internal header */ -#define DECNUMDIGITS 1000 -#include "decNumber.h" -@<Internal library declarations@>; -#endif - -@* Math initialization. - -First, here are some very important constants. - -@d E_STRING "2.7182818284590452353602874713526624977572470936999595749669676277240766303535" -@d PI_STRING "3.1415926535897932384626433832795028841971693993751058209749445923078164062862" -@d fraction_multiplier 4096 -@d angle_multiplier 16 - -@ Here are the functions that are static as they are not used elsewhere - -@<Declarations@>= -#define DEBUG 0 -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_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d); -/*|static void mp_decimal_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);|*/ -static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c); -static void mp_decimal_number_modulo (mp_number *a, mp_number b); -static void mp_decimal_print_number (MP mp, mp_number n); -static char * mp_decimal_number_tostring (MP mp, mp_number n); -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); -static void mp_decimal_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_sin); -static void mp_init_randoms (MP mp, int seed); -static void mp_number_angle_to_scaled (mp_number *A); -static void mp_number_fraction_to_scaled (mp_number *A); -static void mp_number_scaled_to_fraction (mp_number *A); -static void mp_number_scaled_to_angle (mp_number *A); -static void mp_decimal_m_unif_rand (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_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_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b); -static void mp_decimal_pyth_add (MP mp, mp_number *r, mp_number a, mp_number b); -static void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number x, mp_number y); -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_set_decimal_from_int(mp_number *A, int B); -static void mp_set_decimal_from_boolean(mp_number *A, int B); -static void mp_set_decimal_from_scaled(mp_number *A, int B); -static void mp_set_decimal_from_addition(mp_number *A, mp_number B, mp_number C); -static void mp_set_decimal_from_substraction (mp_number *A, mp_number B, mp_number C); -static void mp_set_decimal_from_div(mp_number *A, mp_number B, mp_number C); -static void mp_set_decimal_from_mul(mp_number *A, mp_number B, mp_number C); -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_of_the_way(MP mp, mp_number *A, mp_number t, mp_number B, mp_number C); -static void mp_number_negate(mp_number *A); -static void mp_number_add(mp_number *A, mp_number B); -static void mp_number_substract(mp_number *A, mp_number B); -static void mp_number_half(mp_number *A); -static void mp_number_halfp(mp_number *A); -static void mp_number_double(mp_number *A); -static void mp_number_add_scaled(mp_number *A, int B); /* also for negative B */ -static void mp_number_multiply_int(mp_number *A, int B); -static void mp_number_divide_int(mp_number *A, int B); -static void mp_decimal_abs(mp_number *A); -static void mp_number_clone(mp_number *A, mp_number B); -static void mp_number_swap(mp_number *A, mp_number *B); -static int mp_round_unscaled(mp_number x_orig); -static int mp_number_to_int(mp_number A); -static int mp_number_to_scaled(mp_number A); -static int mp_number_to_boolean(mp_number A); -static double mp_number_to_double(mp_number A); -static int mp_number_odd(mp_number A); -static int mp_number_equal(mp_number A, mp_number B); -static int mp_number_greater(mp_number A, mp_number B); -static int mp_number_less(mp_number A, mp_number B); -static int mp_number_nonequalabs(mp_number A, mp_number B); -static void mp_number_floor (mp_number *i); -static void mp_decimal_fraction_to_round_scaled (mp_number *x); -static void mp_decimal_number_make_scaled (MP mp, mp_number *r, mp_number p, mp_number q); -static void mp_decimal_number_make_fraction (MP mp, mp_number *r, mp_number p, mp_number q); -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_new_number (MP mp, mp_number *n, mp_number_type t) ; -static void mp_free_number (MP mp, mp_number *n) ; -static void mp_set_decimal_from_double(mp_number *A, double B); -static void mp_free_decimal_math (MP mp); -static void mp_decimal_set_precision (MP mp); -static void mp_check_decNumber (MP mp, decNumber *dec, decContext *context); -static int decNumber_check (decNumber *dec, decContext *context); -static char * mp_decnumber_tostring (decNumber *n); - -@ We do not want special numbers as return values for functions, so: - - -@c -int decNumber_check (decNumber *dec, decContext *context) -{ - int test = false; - if (context->status & DEC_Overflow) { - test = true; - context->status &= ~DEC_Overflow; - } - if (context->status & DEC_Underflow) { - test = true; - context->status &= ~DEC_Underflow; - } - if (context->status & DEC_Errors) { -/*|fprintf(stdout, "DEC_ERROR %x (%s)\n", context->status, decContextStatusToString(context));|*/ - test = true; - decNumberZero(dec); - } - context->status = 0; - if (decNumberIsSpecial(dec)) { - test = true; - if (decNumberIsInfinite(dec)) { - if (decNumberIsNegative(dec)) { - decNumberCopyNegate(dec, &EL_GORDO_decNumber); - } else { - decNumberCopy(dec, &EL_GORDO_decNumber); - } - } else { /* Nan */ - decNumberZero(dec); - } - } - if (decNumberIsZero(dec) && decNumberIsNegative(dec)) { - decNumberZero(dec); - } - return test; -} -void mp_check_decNumber (MP mp, decNumber *dec, decContext *context) -{ - mp->arith_error = decNumber_check (dec, context); -} - - - - -@ There are a few short decNumber functions that do not exist, but -make life easier for us: - -@d decNumberIsPositive(A) !(decNumberIsZero(A) || decNumberIsNegative(A)) - -@c -static decContext set; -static decContext limitedset; -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, &set); - return decNumberIsNegative(&comp); -} -static int decNumberGreater(decNumber *a, decNumber *b) { - decNumber comp; - decNumberCompare(&comp, a, b, &set); - return decNumberIsPositive(&comp); -} -static void decNumberFromDouble(decNumber *A, double B) { - char buf[1000]; - char *c; - snprintf(buf,1000,"%-650.325lf",B); - c = buf; - while (*c++) { - if (*c == ' ') { - *c = '\0'; - break; - } - } - decNumberFromString(A, buf, &set); -} -static double decNumberToDouble(decNumber *A) { - char *buffer = malloc(A->digits + 14); - double res = 0.0; - assert (buffer); - decNumberToString(A, buffer); - if (sscanf(buffer, "%lf", &res)) { - free(buffer); - return res; - } else { - free(buffer); - /*|mp->arith_error = 1;|*/ - return 0.0; /* whatever*/ - } -} -@ Borrowed code from libdfp: - -% x^3 x^5 x^7 -%arctan(x) = x - --- + --- - --- + ... -% 3 5 7 -$$ \arctan(x) = x - {x^3\over3} + {x^5\over5} - {x^7\over7} + \ldots$$ - - -This power series works well, if $x$ is close to zero ($|x|<0.5$). -If x is larger, the series converges too slowly, -so in order to get a smaller x, we apply the identity - -% sqrt(1+x^2) - 1 -%arctan(x) = 2*arctan --------------- -% x -$$ \arctan(x) = 2\,\arctan{{\sqrt{1+x^2}-1}\over x}$$ - -twice. The first application gives us a new $x$ with $x < 1$. -The second application gives us a new x with $x < 0.4142136$. -For that $x$, we use the power series and multiply the result by four. - - - -@c -static void decNumberAtan (decNumber *result, decNumber *x_orig, decContext *set) -{ - decNumber x, f, g, mx2, term; - int i; - decNumberCopy(&x, x_orig); - if (decNumberIsZero (&x)) { - decNumberCopy (result, &x); - return; - } - for (i=0; i<2; i++) { - decNumber y; - decNumberMultiply (&y, &x, &x, set); /* $y = x^2$ */ - decNumberAdd (&y, &y, &one, set); /* $y = y+1$*/ - decNumberSquareRoot (&y, &y, set); /* $y = sqrt(y)$ */ - decNumberSubtract (&y, &y, &one, set); /* $y = y-1$ */ - decNumberDivide (&x, &y, &x, set); /* $x = y/x$ */ - if (decNumberIsZero (&x)) { - decNumberCopy (result, &x); - return; - } - } - decNumberCopy (&f, &x); /* $f(0) = x$ */ - decNumberCopy (&g, &one); /*$ g(0) = 1$*/ - decNumberCopy (&term, &x); /*$ term = x$*/ - decNumberCopy (result, &x); /*$ sum = x $*/ - decNumberMultiply (&mx2, &x, &x, set); /*$ mx2 = x^2$*/ - decNumberMinus (&mx2, &mx2, set); /*$ mx2 = -x^2 $*/ - for (i=0; i<2*set->digits; i++) { - decNumberMultiply (&f, &f, &mx2, set); - decNumberAdd (&g, &g, &two_decNumber, set); - decNumberDivide (&term, &f, &g, set); - decNumberAdd (result, result, &term, set); - } - decNumberAdd (result, result, result, set); - decNumberAdd (result, result, result, set); - return; -} -static void decNumberAtan2 (decNumber *result, decNumber *y, decNumber *x, decContext *set) -{ - decNumber temp; - if (!decNumberIsInfinite (x) && !decNumberIsZero (y) - && !decNumberIsInfinite (y) && !decNumberIsZero (x)) { - decNumberDivide (&temp, y, x, set); - decNumberAtan (result, &temp, set); - /* decNumberAtan doesn't quite return the values in the ranges we - * want for x < 0. So we need to do some correction */ - if (decNumberIsNegative (x)) { - if (decNumberIsNegative (y)) { - decNumberSubtract(result, result, &PI_decNumber, set); - } else { - decNumberAdd(result, result, &PI_decNumber, set); - } - } - return; - } - if (decNumberIsInfinite (y) && decNumberIsInfinite (x)) { - /* If x and y are both inf, the result depends on the sign of x */ - decNumberDivide(result, &PI_decNumber, &four_decNumber, set); - if (decNumberIsNegative (x) ) { - decNumber a; - decNumberFromDouble(&a, 3.0); - decNumberMultiply(result, result, &a, set); - } - } else if (!decNumberIsZero (y) && !decNumberIsInfinite (x) ) { - /* If y is non-zero and x is non-inf, the result is +-pi/2 */ - decNumberDivide(result, &PI_decNumber, &two_decNumber, set); - } else { /* Otherwise it is +0 if x is positive, +pi if x is neg */ - if (decNumberIsNegative (x)) { - decNumberCopy(result, &PI_decNumber); - } else { - decNumberZero(result); - } - } - /* Atan2 will be negative if y<0 */ - if (decNumberIsNegative (y)) { - decNumberMinus(result, result, set); - } -} - -@ And these are the ones that {\it are} used elsewhere - -@<Internal library declarations@>= -void * mp_initialize_decimal_math (MP mp); - -@ - -@d unity 1 -@d two 2 -@d three 3 -@d four 4 -@d half_unit 0.5 -@d three_quarter_unit 0.75 -@d coef_bound ((7.0/3.0)*fraction_multiplier) /* |fraction| approximation to 7/3 */ -@d fraction_threshold 0.04096 /* a |fraction| coefficient less than this is zeroed */ -@d half_fraction_threshold (fraction_threshold/2) /* half of |fraction_threshold| */ -@d scaled_threshold 0.000122 /* a |scaled| coefficient less than this is zeroed */ -@d half_scaled_threshold (scaled_threshold/2) /* half of |scaled_threshold| */ -@d near_zero_angle (0.0256*angle_multiplier) /* an angle of about 0.0256 */ -@d p_over_v_threshold 0x80000 /* TODO */ -@d equation_threshold 0.001 -@d tfm_warn_threshold 0.0625 -@d epsilon pow(2.0,-173.0) /* almost "1E-52" */ -@d epsilonf pow(2.0,-52.0) -@d EL_GORDO "1E1000000" /* the largest value that \MP\ likes. */ -@d warning_limit "1E1000000" /* this is a large value that can just be expressed without loss of precision */ -@d DECPRECISION_DEFAULT 34 - -@<Declarations@>= -static decNumber zero; -static decNumber one; -static decNumber minusone; -static decNumber two_decNumber; -static decNumber three_decNumber; -static decNumber four_decNumber; -static decNumber fraction_multiplier_decNumber; -static decNumber angle_multiplier_decNumber; -static decNumber fraction_one_decNumber; -static decNumber fraction_one_plus_decNumber; -static decNumber PI_decNumber; -static decNumber epsilon_decNumber; -static decNumber EL_GORDO_decNumber; -static decNumber **factorials = NULL; -static int last_cached_factorial = 0; -static boolean initialized = false ; -@ @c -void * mp_initialize_decimal_math (MP mp) { - math_data *math = (math_data *)mp_xmalloc(mp,1,sizeof(math_data)); - /* various decNumber initializations */ - decContextDefault(&set, DEC_INIT_BASE); /* initialize */ - set.traps=0; /* no traps, thank you */ - decContextDefault(&limitedset, DEC_INIT_BASE); /* initialize */ - limitedset.traps=0; /* no traps, thank you */ - limitedset.emax = 999999; - limitedset.emin = -999999; - set.digits = DECPRECISION_DEFAULT; - limitedset.digits = DECPRECISION_DEFAULT; - if (!initialized) { - initialized = true ; - decNumberFromInt32(&one, 1); - decNumberFromInt32(&minusone, -1); - decNumberFromInt32(&zero, 0); - decNumberFromInt32(&two_decNumber, two); - decNumberFromInt32(&three_decNumber, three); - decNumberFromInt32(&four_decNumber, four); - decNumberFromInt32(&fraction_multiplier_decNumber, fraction_multiplier); - decNumberFromInt32(&fraction_one_decNumber, fraction_one); - decNumberFromInt32(&fraction_one_plus_decNumber, (fraction_one+1)); - decNumberFromInt32(&angle_multiplier_decNumber, angle_multiplier); - decNumberFromString(&PI_decNumber, PI_STRING, &set); - decNumberFromDouble(&epsilon_decNumber, epsilon); - decNumberFromString(&EL_GORDO_decNumber, EL_GORDO, &set); - factorials = (decNumber **)mp_xmalloc(mp,PRECALC_FACTORIALS_CACHESIZE,sizeof(decNumber *)); - factorials[0] = (decNumber *)mp_xmalloc(mp,1,sizeof(decNumber)); - decNumberCopy(factorials[0], &one); - } - - /* alloc */ - math->allocate = mp_new_number; - math->free = mp_free_number; - mp_new_number (mp, &math->precision_default, mp_scaled_type); - decNumberFromInt32(math->precision_default.data.num, DECPRECISION_DEFAULT); - mp_new_number (mp, &math->precision_max, mp_scaled_type); - decNumberFromInt32(math->precision_max.data.num, DECNUMDIGITS); - mp_new_number (mp, &math->precision_min, mp_scaled_type); - decNumberFromInt32(math->precision_min.data.num, 2); - /* here are the constants for |scaled| objects */ - mp_new_number (mp, &math->epsilon_t, mp_scaled_type); - decNumberCopy(math->epsilon_t.data.num, &epsilon_decNumber); - mp_new_number (mp, &math->inf_t, mp_scaled_type); - decNumberCopy(math->inf_t.data.num, &EL_GORDO_decNumber); - mp_new_number (mp, &math->warning_limit_t, mp_scaled_type); - decNumberFromString(math->warning_limit_t.data.num, warning_limit, &set); - mp_new_number (mp, &math->one_third_inf_t, mp_scaled_type); - decNumberDivide(math->one_third_inf_t.data.num, math->inf_t.data.num, &three_decNumber, &set); - mp_new_number (mp, &math->unity_t, mp_scaled_type); - decNumberCopy(math->unity_t.data.num, &one); - mp_new_number (mp, &math->two_t, mp_scaled_type); - decNumberFromInt32(math->two_t.data.num, two); - mp_new_number (mp, &math->three_t, mp_scaled_type); - decNumberFromInt32(math->three_t.data.num, three); - mp_new_number (mp, &math->half_unit_t, mp_scaled_type); - decNumberFromString(math->half_unit_t.data.num, "0.5", &set); - mp_new_number (mp, &math->three_quarter_unit_t, mp_scaled_type); - decNumberFromString(math->three_quarter_unit_t.data.num, "0.75", &set); - mp_new_number (mp, &math->zero_t, mp_scaled_type); - decNumberZero(math->zero_t.data.num); - /* |fractions| */ - mp_new_number (mp, &math->arc_tol_k, mp_fraction_type); - { - decNumber fourzeroninesix; - decNumberFromInt32(&fourzeroninesix, 4096); - decNumberDivide(math->arc_tol_k.data.num, &one, &fourzeroninesix, &set); - /* quit when change in arc length estimate reaches this */ - } - mp_new_number (mp, &math->fraction_one_t, mp_fraction_type); - decNumberFromInt32(math->fraction_one_t.data.num, fraction_one); - mp_new_number (mp, &math->fraction_half_t, mp_fraction_type); - decNumberFromInt32(math->fraction_half_t.data.num, fraction_half); - mp_new_number (mp, &math->fraction_three_t, mp_fraction_type); - decNumberFromInt32(math->fraction_three_t.data.num, fraction_three); - mp_new_number (mp, &math->fraction_four_t, mp_fraction_type); - decNumberFromInt32(math->fraction_four_t.data.num, fraction_four); - /* |angles| */ - mp_new_number (mp, &math->three_sixty_deg_t, mp_angle_type); - decNumberFromInt32(math->three_sixty_deg_t.data.num, 360 * angle_multiplier); - mp_new_number (mp, &math->one_eighty_deg_t, mp_angle_type); - decNumberFromInt32(math->one_eighty_deg_t.data.num, 180 * angle_multiplier); - /* various approximations */ - mp_new_number (mp, &math->one_k, mp_scaled_type); - decNumberFromDouble(math->one_k.data.num, 1.0/64); - mp_new_number (mp, &math->sqrt_8_e_k, mp_scaled_type); - { - decNumberFromDouble(math->sqrt_8_e_k.data.num, 112428.82793 / 65536.0); - /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */ - } - mp_new_number (mp, &math->twelve_ln_2_k, mp_fraction_type); - { - decNumberFromDouble(math->twelve_ln_2_k.data.num, 139548959.6165 / 65536.0); - /* $2^{24}\cdot12\ln2\approx139548959.6165$ */ - } - mp_new_number (mp, &math->coef_bound_k, mp_fraction_type); - decNumberFromDouble(math->coef_bound_k.data.num,coef_bound); - mp_new_number (mp, &math->coef_bound_minus_1, mp_fraction_type); - decNumberFromDouble(math->coef_bound_minus_1.data.num,coef_bound - 1 / 65536.0); - mp_new_number (mp, &math->twelvebits_3, mp_scaled_type); - { - decNumberFromDouble(math->twelvebits_3.data.num, 1365 / 65536.0); - /* $1365\approx 2^{12}/3$ */ - } - mp_new_number (mp, &math->twentysixbits_sqrt2_t, mp_fraction_type); - { - decNumberFromDouble(math->twentysixbits_sqrt2_t.data.num, 94906265.62 / 65536.0); - /* $2^{26}\sqrt2\approx94906265.62$ */ - } - mp_new_number (mp, &math->twentyeightbits_d_t, mp_fraction_type); - { - decNumberFromDouble(math->twentyeightbits_d_t.data.num, 35596754.69 / 65536.0); - /* $2^{28}d\approx35596754.69$ */ - } - mp_new_number (mp, &math->twentysevenbits_sqrt2_d_t, mp_fraction_type); - { - decNumberFromDouble(math->twentysevenbits_sqrt2_d_t.data.num, 25170706.63 / 65536.0); - /* $2^{27}\sqrt2\,d\approx25170706.63$ */ - } - /* thresholds */ - mp_new_number (mp, &math->fraction_threshold_t, mp_fraction_type); - decNumberFromDouble(math->fraction_threshold_t.data.num, fraction_threshold); - mp_new_number (mp, &math->half_fraction_threshold_t, mp_fraction_type); - decNumberFromDouble(math->half_fraction_threshold_t.data.num, half_fraction_threshold); - mp_new_number (mp, &math->scaled_threshold_t, mp_scaled_type); - decNumberFromDouble(math->scaled_threshold_t.data.num, scaled_threshold); - mp_new_number (mp, &math->half_scaled_threshold_t, mp_scaled_type); - decNumberFromDouble(math->half_scaled_threshold_t.data.num, half_scaled_threshold); - mp_new_number (mp, &math->near_zero_angle_t, mp_angle_type); - decNumberFromDouble(math->near_zero_angle_t.data.num, near_zero_angle); - mp_new_number (mp, &math->p_over_v_threshold_t, mp_fraction_type); - decNumberFromDouble(math->p_over_v_threshold_t.data.num, p_over_v_threshold); - mp_new_number (mp, &math->equation_threshold_t, mp_scaled_type); - decNumberFromDouble(math->equation_threshold_t.data.num, equation_threshold); - mp_new_number (mp, &math->tfm_warn_threshold_t, mp_scaled_type); - decNumberFromDouble(math->tfm_warn_threshold_t.data.num, tfm_warn_threshold); - /* functions */ - math->from_int = mp_set_decimal_from_int; - math->from_boolean = mp_set_decimal_from_boolean; - math->from_scaled = mp_set_decimal_from_scaled; - math->from_double = mp_set_decimal_from_double; - math->from_addition = mp_set_decimal_from_addition; - math->from_substraction = mp_set_decimal_from_substraction; - math->from_oftheway = mp_set_decimal_from_of_the_way; - math->from_div = mp_set_decimal_from_div; - math->from_mul = mp_set_decimal_from_mul; - math->from_int_div = mp_set_decimal_from_int_div; - math->from_int_mul = mp_set_decimal_from_int_mul; - math->negate = mp_number_negate; - math->add = mp_number_add; - math->substract = mp_number_substract; - math->half = mp_number_half; - math->halfp = mp_number_halfp; - math->do_double = mp_number_double; - math->abs = mp_decimal_abs; - math->clone = mp_number_clone; - math->swap = mp_number_swap; - math->add_scaled = mp_number_add_scaled; - math->multiply_int = mp_number_multiply_int; - math->divide_int = mp_number_divide_int; - math->to_boolean = mp_number_to_boolean; - math->to_scaled = mp_number_to_scaled; - math->to_double = mp_number_to_double; - math->to_int = mp_number_to_int; - math->odd = mp_number_odd; - math->equal = mp_number_equal; - math->less = mp_number_less; - math->greater = mp_number_greater; - math->nonequalabs = mp_number_nonequalabs; - math->round_unscaled = mp_round_unscaled; - math->floor_scaled = mp_number_floor; - math->fraction_to_round_scaled = mp_decimal_fraction_to_round_scaled; - math->make_scaled = mp_decimal_number_make_scaled; - math->make_fraction = mp_decimal_number_make_fraction; - math->take_fraction = mp_decimal_number_take_fraction; - math->take_scaled = mp_decimal_number_take_scaled; - math->velocity = mp_decimal_velocity; - math->n_arg = mp_decimal_n_arg; - math->m_log = mp_decimal_m_log; - math->m_exp = mp_decimal_m_exp; - math->m_unif_rand = mp_decimal_m_unif_rand; - math->m_norm_rand = mp_decimal_m_norm_rand; - math->pyth_add = mp_decimal_pyth_add; - math->pyth_sub = mp_decimal_pyth_sub; - math->fraction_to_scaled = mp_number_fraction_to_scaled; - math->scaled_to_fraction = mp_number_scaled_to_fraction; - math->scaled_to_angle = mp_number_scaled_to_angle; - math->angle_to_scaled = mp_number_angle_to_scaled; - math->init_randoms = mp_init_randoms; - math->sin_cos = mp_decimal_sin_cos; - math->slow_add = mp_decimal_slow_add; - math->sqrt = mp_decimal_square_rt; - math->print = mp_decimal_print_number; - math->tostring = mp_decimal_number_tostring; - math->modulo = mp_decimal_number_modulo; - math->ab_vs_cd = mp_ab_vs_cd; - math->crossing_point = mp_decimal_crossing_point; - math->scan_numeric = mp_decimal_scan_numeric_token; - math->scan_fractional = mp_decimal_scan_fractional_token; - math->free_math = mp_free_decimal_math; - math->set_precision = mp_decimal_set_precision; - return (void *)math; -} - -void mp_decimal_set_precision (MP mp) { - int i; - i = decNumberToInt32((decNumber *)internal_value (mp_number_precision).data.num, &set); - set.digits = i; - limitedset.digits = i; -} - -void mp_free_decimal_math (MP mp) { - free_number (((math_data *)mp->math)->three_sixty_deg_t); - free_number (((math_data *)mp->math)->one_eighty_deg_t); - free_number (((math_data *)mp->math)->fraction_one_t); - free_number (((math_data *)mp->math)->zero_t); - free_number (((math_data *)mp->math)->half_unit_t); - free_number (((math_data *)mp->math)->three_quarter_unit_t); - free_number (((math_data *)mp->math)->unity_t); - free_number (((math_data *)mp->math)->two_t); - free_number (((math_data *)mp->math)->three_t); - free_number (((math_data *)mp->math)->one_third_inf_t); - free_number (((math_data *)mp->math)->inf_t); - free_number (((math_data *)mp->math)->warning_limit_t); - free_number (((math_data *)mp->math)->one_k); - free_number (((math_data *)mp->math)->sqrt_8_e_k); - free_number (((math_data *)mp->math)->twelve_ln_2_k); - free_number (((math_data *)mp->math)->coef_bound_k); - free_number (((math_data *)mp->math)->coef_bound_minus_1); - free_number (((math_data *)mp->math)->fraction_threshold_t); - free_number (((math_data *)mp->math)->half_fraction_threshold_t); - free_number (((math_data *)mp->math)->scaled_threshold_t); - free_number (((math_data *)mp->math)->half_scaled_threshold_t); - free_number (((math_data *)mp->math)->near_zero_angle_t); - free_number (((math_data *)mp->math)->p_over_v_threshold_t); - free_number (((math_data *)mp->math)->equation_threshold_t); - free_number (((math_data *)mp->math)->tfm_warn_threshold_t); - /* For sake of speed, we accept this memory leak. */ - /* for (i = 0; i <= last_cached_factorial; i++) {*/ - /* free(factorials[i]);*/ - /* }*/ - /* free(factorials); */ - free(mp->math); -} - -@ Creating an destroying |mp_number| objects - -@ @c -void mp_new_number (MP mp, mp_number *n, mp_number_type t) { - (void)mp; - n->data.num = mp_xmalloc(mp,1,sizeof(decNumber)); - decNumberZero(n->data.num); - n->type = t; -} - -@ - -@c -void mp_free_number (MP mp, mp_number *n) { - (void)mp; - free(n->data.num); - n->data.num = NULL; - n->type = mp_nan_type; -} - -@ Here are the low-level functions on |mp_number| items, setters first. - -@c -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, &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, &set); -} -void mp_set_decimal_from_substraction (mp_number *A, mp_number B, mp_number C) { - decNumberSubtract(A->data.num,B.data.num,C.data.num, &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, &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, &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, &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, &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, &set); - mp_decimal_take_fraction(mp, &r1, &c, t.data.num); - decNumberSubtract(A->data.num, B.data.num, &r1, &set); - mp_check_decNumber(mp, A->data.num, &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, &set); -} -void mp_number_substract(mp_number *A, mp_number B) { - decNumberSubtract(A->data.num,A->data.num,B.data.num, &set); -} -void mp_number_half(mp_number *A) { - decNumber c; - decNumberFromInt32(&c, 2); - decNumberDivide(A->data.num,A->data.num, &c, &set); -} -void mp_number_halfp(mp_number *A) { - decNumber c; - decNumberFromInt32(&c, 2); - decNumberDivide(A->data.num,A->data.num, &c, &set); -} -void mp_number_double(mp_number *A) { - decNumber c; - decNumberFromInt32(&c, 2); - decNumberMultiply(A->data.num,A->data.num, &c, &set); -} -void mp_number_add_scaled(mp_number *A, int B) { /* also for negative B */ - decNumber b,c; - decNumberFromInt32(&c, 65536); - decNumberFromInt32(&b, B); - decNumberDivide(&b,&b, &c, &set); - decNumberAdd(A->data.num,A->data.num, &b, &set); -} -void mp_number_multiply_int(mp_number *A, int B) { - decNumber b; - decNumberFromInt32(&b, B); - decNumberMultiply(A->data.num,A->data.num, &b, &set); -} -void mp_number_divide_int(mp_number *A, int B) { - decNumber b; - decNumberFromInt32(&b, B); - decNumberDivide(A->data.num,A->data.num,&b, &set); -} -void mp_decimal_abs(mp_number *A) { - decNumberAbs(A->data.num, A->data.num, &set); -} -void mp_number_clone(mp_number *A, mp_number B) { - decNumberCopy(A->data.num, B.data.num); -} -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, &fraction_multiplier_decNumber, &set); -} -void mp_number_angle_to_scaled (mp_number *A) { - A->type = mp_scaled_type; - decNumberDivide(A->data.num, A->data.num, &angle_multiplier_decNumber, &set); -} -void mp_number_scaled_to_fraction (mp_number *A) { - A->type = mp_fraction_type; - decNumberMultiply(A->data.num, A->data.num, &fraction_multiplier_decNumber, &set); -} -void mp_number_scaled_to_angle (mp_number *A) { - A->type = mp_angle_type; - decNumberMultiply(A->data.num, A->data.num, &angle_multiplier_decNumber, &set); -} - - -@* Query functions. - -@ Convert a number to a scaled value. |decNumberToInt32| is not -able to make this conversion properly, so instead we are using -|decNumberToDouble| and a typecast. Bad! - -@c -int mp_number_to_scaled(mp_number A) { - int32_t result; - decNumber corrected; - decNumberFromInt32(&corrected, 65536); - decNumberMultiply(&corrected,&corrected,A.data.num, &set); - decNumberReduce(&corrected, &corrected, &set); - result = (int)floor(decNumberToDouble(&corrected)+0.5); - return result; -} - -@ - -@d odd(A) (abs(A)%2==1) - -@c -int mp_number_to_int(mp_number A) { - int32_t result; - set.status = 0; - result = decNumberToInt32(A.data.num, &set); - if (set.status == DEC_Invalid_operation) { - set.status = 0; - /* |mp->arith_error = 1;| */ - return 0; /* whatever */ - } else { - return result; - } -} -int mp_number_to_boolean(mp_number A) { - uint32_t result; - set.status = 0; - result = decNumberToUInt32(A.data.num, &set); - if (set.status == DEC_Invalid_operation) { - set.status = 0; - /* |mp->arith_error = 1;| */ - return mp_false_code; /* whatever */ - } else { - return result ; - } -} -double mp_number_to_double(mp_number A) { - char *buffer = malloc(((decNumber *)A.data.num)->digits + 14); - double res = 0.0; - assert (buffer); - decNumberToString(A.data.num, buffer); - if (sscanf(buffer, "%lf", &res)) { - free(buffer); - return res; - } else { - free(buffer); - /* |mp->arith_error = 1;| */ - return 0.0; /* whatever */ - } -} -int mp_number_odd(mp_number A) { - return odd(mp_number_to_int(A)); -} -int mp_number_equal(mp_number A, mp_number B) { - decNumber res; - decNumberCompare(&res,A.data.num,B.data.num, &set); - return decNumberIsZero(&res); -} -int mp_number_greater(mp_number A, mp_number B) { - decNumber res; - decNumberCompare(&res,A.data.num,B.data.num, &set); - return decNumberIsPositive(&res); -} -int mp_number_less(mp_number A, mp_number B) { - decNumber res; - decNumberCompare(&res,A.data.num,B.data.num, &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, &set); - return !decNumberIsZero(&res); -} - -@ Fixed-point arithmetic is done on {\sl scaled integers\/} that are multiples -of $2^{-16}$. In other words, a binary point is assumed to be sixteen bit -positions from the right end of a binary computer word. - -@ One of \MP's most common operations is the calculation of -$\lfloor{a+b\over2}\rfloor$, -the midpoint of two given integers |a| and~|b|. The most decent way to do -this is to write `|(a+b)/2|'; but on many machines it is more efficient -to calculate `|(a+b)>>1|'. - -Therefore the midpoint operation will always be denoted by `|half(a+b)|' -in this program. If \MP\ is being implemented with languages that permit -binary shifting, the |half| macro should be changed to make this operation -as efficient as possible. Since some systems have shift operators that can -only be trusted to work on positive numbers, there is also a macro |halfp| -that is used only when the quantity being halved is known to be positive -or zero. - -@ Here is a procedure analogous to |print_int|. The current version -is fairly stupid, and it is not round-trip safe, but this is good -enough for a beta test. - -@c -char * mp_decnumber_tostring (decNumber *n) { - decNumber corrected; - char *buffer = malloc(((decNumber *)n)->digits + 14); - assert (buffer); - decNumberCopy(&corrected,n); - decNumberTrim(&corrected); - decNumberToString(&corrected, buffer); - return buffer; -} -char * mp_decimal_number_tostring (MP mp, mp_number n) { - return mp_decnumber_tostring(n.data.num); -} - - -@ @c -void mp_decimal_print_number (MP mp, mp_number n) { - char *str = mp_decimal_number_tostring(mp, n); - mp_print (mp, str); - free (str); -} - - - - -@ Addition is not always checked to make sure that it doesn't overflow, -but in places where overflow isn't too unlikely the |slow_add| routine -is used. - -@c -void mp_decimal_slow_add (MP mp, mp_number *ret, mp_number A, mp_number B) { - decNumberAdd(ret->data.num,A.data.num,B.data.num, &set); -} - -@ The |make_fraction| routine produces the |fraction| equivalent of -|p/q|, given integers |p| and~|q|; it computes the integer -$f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are -positive. If |p| and |q| are both of the same scaled type |t|, -the ``type relation'' |make_fraction(t,t)=fraction| is valid; -and it's also possible to use the subroutine ``backwards,'' using -the relation |make_fraction(t,fraction)=t| between scaled types. - -If the result would have magnitude $2^{31}$ or more, |make_fraction| -sets |arith_error:=true|. Most of \MP's internal computations have -been designed to avoid this sort of error. - -If this subroutine were programmed in assembly language on a typical -machine, we could simply compute |(@t$2^{28}$@>*p)div q|, since a -double-precision product can often be input to a fixed-point division -instruction. But when we are restricted to int-eger arithmetic it -is necessary either to resort to multiple-precision maneuvering -or to use a simple but slow iteration. The multiple-precision technique -would be about three times faster than the code adopted here, but it -would be comparatively long and tricky, involving about sixteen -additional multiplications and divisions. - -This operation is part of \MP's ``inner loop''; indeed, it will -consume nearly 10\pct! of the running time (exclusive of input and output) -if the code below is left unchanged. A machine-dependent recoding -will therefore make \MP\ run faster. The present implementation -is highly portable, but slow; it avoids multiplication and division -except in the initial stage. System wizards should be careful to -replace it with a routine that is guaranteed to produce identical -results in all cases. -@^system dependencies@> - -As noted below, a few more routines should also be replaced by machine-dependent -code, for efficiency. But when a procedure is not part of the ``inner loop,'' -such changes aren't advisable; simplicity and robustness are -preferable to trickery, unless the cost is too high. -@^inner loop@> - -@c -void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) { - decNumberDivide(ret, p, q, &set); - mp_check_decNumber(mp, ret, &set); - decNumberMultiply(ret, ret, &fraction_multiplier_decNumber, &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); -} - -@ @<Declarations@>= -void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q); - -@ The dual of |make_fraction| is |take_fraction|, which multiplies a -given integer~|q| by a fraction~|f|. When the operands are positive, it -computes $p=\lfloor qf/2^{28}+{1\over2}\rfloor$, a symmetric function -of |q| and~|f|. - -This routine is even more ``inner loopy'' than |make_fraction|; -the present implementation consumes almost 20\pct! of \MP's computation -time during typical jobs, so a machine-language substitute is advisable. -@^inner loop@> @^system dependencies@> - -@c -void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) { - decNumberMultiply(ret, p, q, &set); - decNumberDivide(ret, ret, &fraction_multiplier_decNumber, &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); -} - -@ @<Declarations@>= -void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q); - -@ When we want to multiply something by a |scaled| quantity, we use a scheme -analogous to |take_fraction| but with a different scaling. -Given positive operands, |take_scaled| -computes the quantity $p=\lfloor qf/2^{16}+{1\over2}\rfloor$. - -Once again it is a good idea to use a machine-language replacement if -possible; otherwise |take_scaled| will use more than 2\pct! of the running time -when the Computer Modern fonts are being generated. -@^inner loop@> - -@c -void mp_decimal_number_take_scaled (MP mp, mp_number *ret, mp_number p_orig, mp_number q_orig) { - decNumberMultiply(ret->data.num, p_orig.data.num, q_orig.data.num, &set); -} - - -@ For completeness, there's also |make_scaled|, which computes a -quotient as a |scaled| number instead of as a |fraction|. -In other words, the result is $\lfloor2^{16}p/q+{1\over2}\rfloor$, if the -operands are positive. \ (This procedure is not used especially often, -so it is not part of \MP's inner loop.) - -@c -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, &set); - mp_check_decNumber(mp, ret->data.num, &set); -} - -@ -@d halfp(A) (integer)((unsigned)(A) >> 1) - -@* Scanning numbers in the input. - -The definitions below are temporarily here - -@d set_cur_cmd(A) mp->cur_mod_->type=(A) -@d set_cur_mod(A) decNumberCopy((decNumber *)(mp->cur_mod_->data.n.data.num),&A) - -@<Declarations...@>= -static void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop); - -@ -@d too_precise(a) (a == (DEC_Inexact+DEC_Rounded)) -@d too_large(a) (a & DEC_Overflow) -@c -void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop) { - decNumber result; - size_t l = stop-start+1; - char *buf = mp_xmalloc(mp, l+1, 1); - buf[l] = '\0'; - (void)strncpy(buf,(const char *)start, l); - set.status = 0; - decNumberFromString(&result,buf, &set); - free(buf); - if (set.status == 0) { - set_cur_mod(result); - } else if (mp->scanner_status != tex_flushing) { - if (too_large(set.status)) { - const char *hlp[] = {"I could not handle this number specification", - "because it is out of range.", - NULL }; - decNumber_check (&result, &set); - set_cur_mod(result); - mp_error (mp, "Enormous number has been reduced", hlp, false); - } else if (too_precise(set.status)) { - set_cur_mod(result); - if (decNumberIsPositive((decNumber *)internal_value (mp_warning_check).data.num) && - (mp->scanner_status != tex_flushing)) { - char msg[256]; - const char *hlp[] = {"Continue and I'll round the value until it fits the current numberprecision", - "(Set warningcheck:=0 to suppress this message.)", - NULL }; - mp_snprintf (msg, 256, "Number is too precise (numberprecision = %d)", set.digits); - mp_error (mp, msg, hlp, true); - } - } else { /* this also captures underflow */ - const char *hlp[] = {"I could not handle this number specification", - "Error:", - "", - NULL }; - hlp[2] = decContextStatusToString(&set); - mp_error (mp, "Erroneous number specification changed to zero", hlp, false); - decNumberZero(&result); - set_cur_mod(result); - } - } - set_cur_cmd((mp_variable_type)mp_numeric_token); -} - -@ @c -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]] == 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]] == digit_class) { - mp->cur_input.loc_field++; - } - } -} -void mp_decimal_scan_fractional_token (MP mp, int n) { /* n: scaled */ - unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; - unsigned char *stop; - while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == 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); -} - - -@ We just have to collect bytes. - -@c -void mp_decimal_scan_numeric_token (MP mp, int n) { /* n: scaled */ - unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1]; - unsigned char *stop; - while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == 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]] == 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); -} - -@ The |scaled| quantities in \MP\ programs are generally supposed to be -less than $2^{12}$ in absolute value, so \MP\ does much of its internal -arithmetic with 28~significant bits of precision. A |fraction| denotes -a scaled integer whose binary point is assumed to be 28 bit positions -from the right. - -@d fraction_half (fraction_multiplier/2) -@d fraction_one (1*fraction_multiplier) -@d fraction_two (2*fraction_multiplier) -@d fraction_three (3*fraction_multiplier) -@d fraction_four (4*fraction_multiplier) - -@ Here is a typical example of how the routines above can be used. -It computes the function -$${1\over3\tau}f(\theta,\phi)= -{\tau^{-1}\bigl(2+\sqrt2\,(\sin\theta-{1\over16}\sin\phi) - (\sin\phi-{1\over16}\sin\theta)(\cos\theta-\cos\phi)\bigr)\over -3\,\bigl(1+{1\over2}(\sqrt5-1)\cos\theta+{1\over2}(3-\sqrt5\,)\cos\phi\bigr)},$$ -where $\tau$ is a |scaled| ``tension'' parameter. This is \MP's magic -fudge factor for placing the first control point of a curve that starts -at an angle $\theta$ and ends at an angle $\phi$ from the straight path. -(Actually, if the stated quantity exceeds 4, \MP\ reduces it to~4.) - -The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$. -(It's a sum of eight terms whose absolute values can be bounded using -relations such as $\sin\theta\cos\theta\L{1\over2}$.) Thus the numerator -is positive; and since the tension $\tau$ is constrained to be at least -$3\over4$, the numerator is less than $16\over3$. The denominator is -nonnegative and at most~6. - -The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction| -arguments |st|, |ct|, |sf|, and |cf|, representing $\sin\theta$, $\cos\theta$, -$\sin\phi$, and $\cos\phi$, respectively. - -@c -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; /* registers for intermediate calculations */ - 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); /*$\sqrt{5}$*/ - decNumberSquareRoot(&sqrtfive, &sqrtfive, &set); - - - decNumberDivide(&arg1,sf.data.num, &i16, &set); /* arg1 = sf / 16*/ - decNumberSubtract(&arg1,st.data.num,&arg1, &set); /* arg1 = st - arg1*/ - decNumberDivide(&arg2,st.data.num, &i16, &set); /* arg2 = st / 16*/ - decNumberSubtract(&arg2,sf.data.num,&arg2, &set); /* arg2 = sf - arg2*/ - mp_decimal_take_fraction (mp, &acc, &arg1, &arg2); /* acc = (arg1 * arg2) / fmul*/ - - decNumberCopy(&arg1, &acc); - decNumberSubtract(&arg2, ct.data.num, cf.data.num, &set); /* arg2 = ct - cf*/ - mp_decimal_take_fraction (mp, &acc, &arg1, &arg2); /* acc = (arg1 * arg2 ) / fmul*/ - - decNumberSquareRoot(&arg1, &two_decNumber, &set); /* arg1 = $\sqrt{2}$*/ - decNumberMultiply(&arg1, &arg1, &fone, &set); /* arg1 = arg1 * fmul*/ - mp_decimal_take_fraction (mp, &r1, &acc, &arg1); /* r1 = (acc * arg1) / fmul*/ - decNumberAdd(&num, &ftwo, &r1, &set); /* num = ftwo + r1*/ - - decNumberSubtract(&arg1,&sqrtfive, &one, &set); /* arg1 = $\sqrt{5}$ - 1*/ - decNumberMultiply(&arg1,&arg1,&fhalf, &set); /* arg1 = arg1 * fmul/2*/ - decNumberMultiply(&arg1,&arg1,&three_decNumber, &set); /* arg1 = arg1 * 3*/ - - decNumberSubtract(&arg2,&three_decNumber, &sqrtfive, &set); /* arg2 = 3 - $\sqrt{5}$*/ - decNumberMultiply(&arg2,&arg2,&fhalf, &set); /* arg2 = arg2 * fmul/2*/ - decNumberMultiply(&arg2,&arg2,&three_decNumber, &set); /* arg2 = arg2 * 3*/ - mp_decimal_take_fraction (mp, &r1, ct.data.num, &arg1) ; /* r1 = (ct * arg1) / fmul*/ - mp_decimal_take_fraction (mp, &r2, cf.data.num, &arg2); /* r2 = (cf * arg2) / fmul*/ - - decNumberFromInt32(&denom, fraction_three); /* denom = 3fmul*/ - decNumberAdd(&denom, &denom, &r1, &set); /* denom = denom + r1*/ - decNumberAdd(&denom, &denom, &r2, &set); /* denom = denom + r1*/ - - decNumberCompare(&arg1, t.data.num, &one, &set); - if (!decNumberIsZero(&arg1)) { /* t != r1*/ - decNumberDivide(&num, &num, t.data.num, &set); /* num = num / t*/ - } - decNumberCopy(&r2, &num); /* r2 = num / 4*/ - decNumberDivide(&r2, &r2, &four_decNumber, &set); - if (decNumberLess(&denom,&r2)) { /* num/4 >= denom => denom < num/4*/ - decNumberFromInt32(ret->data.num,fraction_four); - } else { - mp_decimal_make_fraction (mp, ret->data.num, &num, &denom); - } -#if DEBUG - fprintf(stdout, "\n%f = velocity(%f,%f,%f,%f,%f)", mp_number_to_double(*ret), -mp_number_to_double(st),mp_number_to_double(ct), -mp_number_to_double(sf),mp_number_to_double(cf), -mp_number_to_double(t)); -#endif - mp_check_decNumber(mp, ret->data.num, &set); -} - - -@ The following somewhat different subroutine tests rigorously if $ab$ is -greater than, equal to, or less than~$cd$, -given integers $(a,b,c,d)$. In most cases a quick decision is reached. -The result is $+1$, 0, or~$-1$ in the three respective cases. - -@c -void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) { - decNumber q, r, test; /* temporary registers */ - decNumber a, b, c, d; - decNumber ab, cd; - (void)mp; - 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, &set); - decNumberMultiply (&cd, (decNumber *)c_orig.data.num, (decNumber *)d_orig.data.num, &set); - decNumberCompare(ret->data.num, &ab, &cd, &set); - mp_check_decNumber(mp, ret->data.num, &set); - if (1>0) - return; - - - @<Reduce to the case that |a,c>=0|, |b,d>0|@>; - while (1) { - decNumberDivide(&q,&a,&d, &set); - decNumberDivide(&r,&c,&b, &set); - decNumberCompare(&test,&q,&r, &set); - if (!decNumberIsZero(&test)) { - if (decNumberIsPositive(&test)) { - decNumberCopy(ret->data.num, &one); - } else { - decNumberCopy(ret->data.num, &minusone); - } - goto RETURN; - } - decNumberRemainder(&q,&a,&d, &set); - decNumberRemainder(&r,&c,&b, &set); - if (decNumberIsZero(&r)) { - if (decNumberIsZero(&q)) { - decNumberCopy(ret->data.num, &zero); - } else { - decNumberCopy(ret->data.num, &one); - } - goto RETURN; - } - if (decNumberIsZero(&q)) { - decNumberCopy(ret->data.num, &minusone); - goto RETURN; - } - decNumberCopy(&a,&b); - decNumberCopy(&b,&q); - decNumberCopy(&c,&d); - decNumberCopy(&d,&r); - } /* now |a>d>0| and |c>b>0| */ -RETURN: -#if DEBUG - fprintf(stdout, "\n%f = ab_vs_cd(%f,%f,%f,%f)", mp_number_to_double(*ret), -mp_number_to_double(a_orig),mp_number_to_double(b_orig), -mp_number_to_double(c_orig),mp_number_to_double(d_orig)); -#endif - mp_check_decNumber(mp, ret->data.num, &set); - return; -} - - -@ @<Reduce to the case that |a...@>= -if (decNumberIsNegative(&a)) { - decNumberCopyNegate(&a, &a); - decNumberCopyNegate(&b, &b); -} -if (decNumberIsNegative(&c)) { - decNumberCopyNegate(&c, &c); - decNumberCopyNegate(&d, &d); -} -if (!decNumberIsPositive(&d)) { - if (!decNumberIsNegative(&b)) { - if ((decNumberIsZero(&a) || decNumberIsZero(&b)) && (decNumberIsZero(&c) || decNumberIsZero(&d))) - decNumberCopy(ret->data.num, &zero); - else - decNumberCopy(ret->data.num, &one); - goto RETURN; - } - if (decNumberIsZero(&d)) { - if (decNumberIsZero(&a)) - decNumberCopy(ret->data.num, &zero); - else - decNumberCopy(ret->data.num, &minusone); - goto RETURN; - } - decNumberCopy(&q, &a); - decNumberCopy(&a, &c); - decNumberCopy(&c, &q); - decNumberCopyNegate(&q, &b); - decNumberCopyNegate(&b, &d); - decNumberCopy(&d, &q); -} else if (!decNumberIsPositive(&b)) { - if (decNumberIsNegative(&b) && decNumberIsPositive(&a)) { - decNumberCopy(ret->data.num, &minusone); - goto RETURN; - } - if (decNumberIsZero(&c)) - decNumberCopy(ret->data.num, &zero); - else - decNumberCopy(ret->data.num, &minusone); - goto RETURN; -} - -@ Now here's a subroutine that's handy for all sorts of path computations: -Given a quadratic polynomial $B(a,b,c;t)$, the |crossing_point| function -returns the unique |fraction| value |t| between 0 and~1 at which -$B(a,b,c;t)$ changes from positive to negative, or returns -|t=fraction_one+1| if no such value exists. If |a<0| (so that $B(a,b,c;t)$ -is already negative at |t=0|), |crossing_point| returns the value zero. - -The general bisection method is quite simple when $n=2$, hence -|crossing_point| does not take much time. At each stage in the -recursion we have a subinterval defined by |l| and~|j| such that -$B(a,b,c;2^{-l}(j+t))=B(x_0,x_1,x_2;t)$, and we want to ``zero in'' on -the subinterval where $x_0\G0$ and $\min(x_1,x_2)<0$. - -It is convenient for purposes of calculation to combine the values -of |l| and~|j| in a single variable $d=2^l+j$, because the operation -of bisection then corresponds simply to doubling $d$ and possibly -adding~1. Furthermore it proves to be convenient to modify -our previous conventions for bisection slightly, maintaining the -variables $X_0=2^lx_0$, $X_1=2^l(x_0-x_1)$, and $X_2=2^l(x_1-x_2)$. -With these variables the conditions $x_0\ge0$ and $\min(x_1,x_2)<0$ are -equivalent to $\max(X_1,X_1+X_2)>X_0\ge0$. - -The following code maintains the invariant relations -$0\L|x0|<\max(|x1|,|x1|+|x2|)$, -$\vert|x1|\vert<2^{30}$, $\vert|x2|\vert<2^{30}$; -it has been constructed in such a way that no arithmetic overflow -will occur if the inputs satisfy -$a<2^{30}$, $\vert a-b\vert<2^{30}$, and $\vert b-c\vert<2^{30}$. - -@d no_crossing { decNumberCopy(ret->data.num, &fraction_one_plus_decNumber); goto RETURN; } -@d one_crossing { decNumberCopy(ret->data.num, &fraction_one_decNumber); goto RETURN; } -@d zero_crossing { decNumberCopy(ret->data.num, &zero); goto RETURN; } - -@c -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; /* recursive counter */ - decNumber x, xx, x0, x1, x2; /* temporary registers for bisection */ - decNumber scratch, scratch2; - decNumberCopy(&a, (decNumber *)aa.data.num); - decNumberCopy(&b, (decNumber *)bb.data.num); - decNumberCopy(&c, (decNumber *)cc.data.num); - if (decNumberIsNegative(&a)) - zero_crossing; - if (!decNumberIsNegative(&c)) { - if (!decNumberIsNegative(&b)) { - if (decNumberIsPositive(&c)) { - no_crossing; - } else if (decNumberIsZero(&a) && decNumberIsZero(&b)) { - no_crossing; - } else { - one_crossing; - } - } - if (decNumberIsZero(&a)) - zero_crossing; - } else if (decNumberIsZero(&a)) { - if (!decNumberIsPositive(&b)) - zero_crossing; - } - - /* Use bisection to find the crossing point... */ - d = epsilonf; - decNumberCopy(&x0, &a); - decNumberSubtract(&x1,&a, &b, &set); - decNumberSubtract(&x2,&b, &c, &set); - /* not sure why the error correction has to be >= 1E-12 */ - decNumberFromDouble(&scratch2, 1E-12); - do { - decNumberAdd(&x, &x1, &x2, &set); - decNumberDivide(&x, &x, &two_decNumber, &set); - decNumberAdd(&x, &x, &scratch2, &set); - decNumberSubtract(&scratch, &x1, &x0, &set); - if (decNumberGreater(&scratch, &x0)) { - decNumberCopy(&x2, &x); - decNumberAdd(&x0, &x0, &x0, &set); - d += d; - } else { - decNumberAdd(&xx, &scratch, &x, &set); - if (decNumberGreater(&xx,&x0)) { - decNumberCopy(&x2,&x); - decNumberAdd(&x0, &x0, &x0, &set); - d += d; - } else { - decNumberSubtract(&x0, &x0, &xx, &set); - if (!decNumberGreater(&x,&x0)) { - decNumberAdd(&scratch, &x, &x2, &set); - if (!decNumberGreater(&scratch, &x0)) - no_crossing; - } - decNumberCopy(&x1,&x); - d = d + d + epsilonf; - } - } - } while (d < fraction_one); - decNumberFromDouble(&scratch, d); - decNumberSubtract(ret->data.num,&scratch, &fraction_one_decNumber, &set); -RETURN: -#if DEBUG - fprintf(stdout, "\n%f = crossing_point(%f,%f,%f)", mp_number_to_double(*ret), -mp_number_to_double(aa),mp_number_to_double(bb),mp_number_to_double(cc)); -#endif - mp_check_decNumber(mp, ret->data.num, &set); - return; -} - - -@ We conclude this set of elementary routines with some simple rounding -and truncation operations. - - -@ |round_unscaled| rounds a |scaled| and converts it to |int| -@c -int mp_round_unscaled(mp_number x_orig) { - double xx = mp_number_to_double(x_orig); - int x = (int)ROUND(xx); - return x; -} - -@ |number_floor| floors a number - -@c -void mp_number_floor (mp_number *i) { - int round = set.round; - set.round = DEC_ROUND_FLOOR; - decNumberToIntegralValue(i->data.num, i->data.num, &set); - set.round = round; -} - -@ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled| -@c -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, &fraction_multiplier_decNumber, &set); -} - - - -@* Algebraic and transcendental functions. -\MP\ computes all of the necessary special functions from scratch, without -relying on |real| arithmetic or system subroutines for sines, cosines, etc. - -@ - -@c -void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number x_orig) { /* return, x: scaled */ - decNumber x; - decNumberCopy(&x, x_orig.data.num); - if (!decNumberIsPositive(&x)) { - @<Handle square root of zero or negative argument@>; - } else { - decNumberSquareRoot(ret->data.num, &x, &set); - } - mp_check_decNumber(mp, ret->data.num, &set); -} - - -@ @<Handle square root of zero...@>= -{ - if (decNumberIsNegative(&x)) { - char msg[256]; - const char *hlp[] = { - "Since I don't take square roots of negative numbers,", - "I'm zeroing this one. Proceed, with fingers crossed.", - NULL }; - char *xstr = mp_decimal_number_tostring (mp, x_orig); - mp_snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr); - free(xstr); -@.Square root...replaced by 0@>; - mp_error (mp, msg, hlp, true); - } - decNumberZero(ret->data.num); - return; -} - - -@ Pythagorean addition $\psqrt{a^2+b^2}$ is implemented by a quick hack - -@c -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, &set); - decNumberMultiply(&bsq, &b, &b, &set); - decNumberAdd(&a, &asq, &bsq, &set); - decNumberSquareRoot(ret->data.num, &a, &set); - /*|if (set.status != 0) {|*/ - /*| mp->arith_error = true;|*/ - /*| decNumberCopy(ret->data.num, &EL_GORDO_decNumber);|*/ - /*|}|*/ - mp_check_decNumber(mp, ret->data.num, &set); -} - -@ Here is a similar algorithm for $\psqrt{a^2-b^2}$. Same quick hack, also. - -@c -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)) { - @<Handle erroneous |pyth_sub| and set |a:=0|@>; - } else { - decNumber asq, bsq; - decNumberMultiply(&asq, &a, &a, &set); - decNumberMultiply(&bsq, &b, &b, &set); - decNumberSubtract(&a, &asq, &bsq, &set); - decNumberSquareRoot(&a, &a, &set); - } - decNumberCopy(ret->data.num, &a); - mp_check_decNumber(mp, ret->data.num, &set); -} - - -@ @<Handle erroneous |pyth_sub| and set |a:=0|@>= -{ - if (decNumberLess(&a, &b)) { - char msg[256]; - const char *hlp[] = { - "Since I don't take square roots of negative numbers,", - "I'm zeroing this one. Proceed, with fingers crossed.", - NULL }; - 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); - free(astr); - free(bstr); -@.Pythagorean...@>; - mp_error (mp, msg, hlp, true); - } - decNumberZero(&a); -} - - -@ Here is the routine that calculates $2^8$ times the natural logarithm -of a |scaled| quantity; - -@c -void mp_decimal_m_log (MP mp, mp_number *ret, mp_number x_orig) { - if (!decNumberIsPositive((decNumber *)x_orig.data.num)) { - @<Handle non-positive logarithm@>; - } else { - decNumber twofivesix; - decNumberFromInt32(&twofivesix, 256); - decNumberLn(ret->data.num, x_orig.data.num, &limitedset); - mp_check_decNumber(mp, ret->data.num, &limitedset); - decNumberMultiply(ret->data.num, ret->data.num, &twofivesix, &set); - } - mp_check_decNumber(mp, ret->data.num, &set); -} - -@ @<Handle non-positive logarithm@>= -{ - char msg[256]; - const char *hlp[] = { - "Since I don't take logs of non-positive numbers,", - "I'm zeroing this one. Proceed, with fingers crossed.", - NULL }; - char *xstr = mp_decimal_number_tostring (mp, x_orig); - mp_snprintf (msg, 256, "Logarithm of %s has been replaced by 0", xstr); - free (xstr); -@.Logarithm...replaced by 0@>; - mp_error (mp, msg, hlp, true); - decNumberZero(ret->data.num); -} - - -@ Conversely, the exponential routine calculates $\exp(x/2^8)$, -when |x| is |scaled|. - -@c -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, &set); - limitedset.status = 0; - decNumberExp(ret->data.num, &temp, &limitedset); - if (limitedset.status & DEC_Clamped) { - if (decNumberIsPositive((decNumber *)x_orig.data.num)) { - mp->arith_error = true; - decNumberCopy(ret->data.num, &EL_GORDO_decNumber); - } else { - decNumberZero(ret->data.num); - } - } - mp_check_decNumber(mp, ret->data.num, &limitedset); - limitedset.status = 0; -} - - -@ Given integers |x| and |y|, not both zero, the |n_arg| function -returns the |angle| whose tangent points in the direction $(x,y)$. - -@c -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)) { - @<Handle undefined arg@>; - } else { - decNumber atan2val, oneeighty_angle; - ret->type = mp_angle_type; - decNumberFromInt32(&oneeighty_angle, 180 * angle_multiplier); - decNumberDivide(&oneeighty_angle, &oneeighty_angle, &PI_decNumber, &set); - checkZero(y_orig.data.num); - checkZero(x_orig.data.num); - decNumberAtan2(&atan2val, y_orig.data.num, x_orig.data.num, &set); -#if DEBUG - fprintf(stdout, "\n%g = atan2(%g,%g)", decNumberToDouble(&atan2val),mp_number_to_double(x_orig),mp_number_to_double(y_orig)); -#endif - decNumberMultiply(ret->data.num,&atan2val, &oneeighty_angle, &set); - checkZero(ret->data.num); -#if DEBUG - fprintf(stdout, "\nn_arg(%g,%g,%g)", mp_number_to_double(*ret), - mp_number_to_double(x_orig),mp_number_to_double(y_orig)); -#endif - } - mp_check_decNumber(mp, ret->data.num, &set); -} - - -@ @<Handle undefined arg@>= -{ - const char *hlp[] = { - "The `angle' between two identical points is undefined.", - "I'm zeroing this one. Proceed, with fingers crossed.", - NULL }; - mp_error (mp, "angle(0,0) is taken as zero", hlp, true); -@.angle(0,0)...zero@>; - decNumberZero(ret->data.num); -} - - -@ Conversely, the |n_sin_cos| routine takes an |angle| and produces the sine -and cosine of that angle. The results of this routine are -stored in global integer variables |n_sin| and |n_cos|. - -First, we need a decNumber function that calculates sines and cosines -using the Taylor series. This function is fairly optimized. - -@d PRECALC_FACTORIALS_CACHESIZE 50 - -@c -static void sinecosine(decNumber *theangle, decNumber *c, decNumber *s) -{ - int n, i, prec; - decNumber p, pxa, fac, cc; - decNumber n1, n2, p1; - decNumberZero(c); - decNumberZero(s); - prec = (set.digits/2); - if (prec < DECPRECISION_DEFAULT) prec = DECPRECISION_DEFAULT; - for (n=0;n<prec;n++) - { - decNumberFromInt32(&p1, n); - decNumberFromInt32(&n1, 2*n); - decNumberPower(&p, &minusone, &p1, &limitedset); - if (n==0) { - decNumberCopy(&pxa, &one); - } else { - decNumberPower(&pxa, theangle, &n1, &limitedset); - } - - if (2*n<last_cached_factorial) { - decNumberCopy(&fac,factorials[2*n]); - } else { - decNumberCopy(&fac,factorials[last_cached_factorial]); - for (i = last_cached_factorial+1; i <= 2*n; i++) { - decNumberFromInt32(&cc, i); - decNumberMultiply (&fac, &fac, &cc, &set); - if (i<PRECALC_FACTORIALS_CACHESIZE) { - factorials[i] = malloc(sizeof(decNumber)); - decNumberCopy(factorials[i],&fac); - last_cached_factorial = i; - } - } - } - - decNumberDivide (&pxa, &pxa, &fac, &set); - decNumberMultiply (&pxa, &pxa, &p, &set); - decNumberAdd (s, s, &pxa, &set); - - decNumberFromInt32(&n2, 2*n+1); - decNumberMultiply (&fac, &fac, &n2, &set); /* fac = fac * (2*n+1)*/ - decNumberPower(&pxa, theangle, &n2, &limitedset); - decNumberDivide (&pxa, &pxa, &fac, &set); - decNumberMultiply (&pxa, &pxa, &p, &set); - decNumberAdd (c, c, &pxa, &set); - /* |printf("\niteration %2d: %-42s %-42s",n,tostring(c), tostring(s));|*/ - } -} - -@ Calculate sines and cosines. -@c -void mp_decimal_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_sin) { - decNumber rad; - double tmp; - decNumber one_eighty; - tmp = mp_number_to_double(z_orig)/16.0; - -#if DEBUG - fprintf(stdout, "\nsin_cos(%f)", mp_number_to_double(z_orig)); -#endif -#if 0 - if (decNumberIsNegative(&rad)) { - while (decNumberLess(&rad,&PI_decNumber)) - decNumberAdd(&rad, &rad, &PI_decNumber, &set); - } else { - while (decNumberGreater(&rad,&PI_decNumber)) - decNumberSubtract(&rad, &rad, &PI_decNumber, &set); - } -#endif - if ((tmp == 90.0)||(tmp == -270)){ - decNumberZero(n_cos->data.num); - decNumberCopy(n_sin->data.num,&fraction_multiplier_decNumber); - } else if ((tmp == -90.0)||(tmp == 270.0)) { - decNumberZero(n_cos->data.num); - decNumberCopyNegate(n_sin->data.num,&fraction_multiplier_decNumber); - } else if ((tmp == 180.0) || (tmp == -180.0)) { - decNumberCopyNegate(n_cos->data.num,&fraction_multiplier_decNumber); - decNumberZero(n_sin->data.num); - } else { - decNumberFromInt32(&one_eighty, 180 * 16); - decNumberMultiply(&rad, z_orig.data.num, &PI_decNumber, &set); - decNumberDivide(&rad, &rad, &one_eighty, &set); - sinecosine(&rad, n_sin->data.num, n_cos->data.num); - decNumberMultiply(n_cos->data.num,n_cos->data.num,&fraction_multiplier_decNumber, &set); - decNumberMultiply(n_sin->data.num,n_sin->data.num,&fraction_multiplier_decNumber, &set); - } -#if DEBUG - fprintf(stdout, "\nsin_cos(%f,%f,%f)", decNumberToDouble(&rad), -mp_number_to_double(*n_cos), mp_number_to_double(*n_sin)); -#endif - mp_check_decNumber(mp, n_cos->data.num, &set); - mp_check_decNumber(mp, n_sin->data.num, &set); -} - -@ This is the {\tt http://www-cs-faculty.stanford.edu/~uno/programs/rng.c} -with small cosmetic modifications. - -@c -#define KK 100 /* the long lag */ -#define LL 37 /* the short lag */ -#define MM (1L<<30) /* the modulus */ -#define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */ -/* */ -static long ran_x[KK]; /* the generator state */ -/* */ -static void ran_array(long aa[],int n) /* put n new random numbers in aa */ - /* long aa[] destination */ - /* int n array length (must be at least KK) */ -{ - register int i,j; - for (j=0;j<KK;j++) aa[j]=ran_x[j]; - for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]); - for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]); - for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]); -} -/* */ -/* the following routines are from exercise 3.6--15 */ -/* after calling |ran_start|, get new randoms by, e.g., "|x=ran_arr_next()|" */ -/* */ -#define QUALITY 1009 /* recommended quality level for high-res use */ -static long ran_arr_buf[QUALITY]; -static long ran_arr_dummy=-1, ran_arr_started=-1; -static long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */ -/* */ -#define TT 70 /* guaranteed separation between streams */ -#define is_odd(x) ((x)&1) /* units bit of x */ -/* */ -static void ran_start(long seed) /* do this before using |ran_array| */ - /* |long seed| selector for different streams */ -{ - register int t,j; - long x[KK+KK-1]; /* the preparation buffer */ - register long ss=(seed+2)&(MM-2); - for (j=0;j<KK;j++) { - x[j]=ss; /* bootstrap the buffer */ - ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */ - } - x[1]++; /* make x[1] (and only x[1]) odd */ - 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; /* "square" */ - 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)) { /* "multiply by z" */ - for (j=KK;j>0;j--) x[j]=x[j-1]; - x[0]=x[KK]; /* shift the buffer cyclically */ - x[LL]=mod_diff(x[LL],x[KK]); - } - if (ss) ss>>=1; else t--; - } - for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j]; - for (;j<KK;j++) ran_x[j-LL]=x[j]; - for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */ - ran_arr_ptr=&ran_arr_started; -} -/* */ -#define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle()) -static long ran_arr_cycle(void) -{ - if (ran_arr_ptr==&ran_arr_dummy) - ran_start(314159L); /* the user forgot to initialize */ - ran_array(ran_arr_buf,QUALITY); - ran_arr_buf[KK]=-1; - ran_arr_ptr=ran_arr_buf+1; - return ran_arr_buf[0]; -} - - - -@ To initialize the |randoms| table, we call the following routine. - -@c -void mp_init_randoms (MP mp, int seed) { - int j, jj, k; /* more or less random integers */ - int i; /* index into |randoms| */ - j = abs (seed); - while (j >= fraction_one) { - j = j/2; - } - k = 1; - for (i = 0; i <= 54; i++) { - 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); /* ``warm up'' the array */ - - ran_start((unsigned long) seed); - -} - -@ @c -void mp_decimal_number_modulo (mp_number *a, mp_number b) { - decNumberRemainder(a->data.num, a->data.num, b.data.num, &set); -} - - -@ To consume a random integer for the uniform generator, the program below will say `|next_unif_random|'. - -@c -static void mp_next_unif_random (MP mp, mp_number *ret) { - decNumber a; - decNumber b; - unsigned long int op; - (void)mp; - op = (unsigned)ran_arr_next(); - decNumberFromInt32(&a, op); - decNumberFromInt32(&b, MM); - decNumberDivide (&a, &a, &b, &set); /* a = a/b */ - decNumberCopy(ret->data.num, &a); - mp_check_decNumber(mp, ret->data.num, &set); -} - - -@ To consume a random fraction, the program below will say `|next_random|'. - -@c -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]); -} - - -@ To produce a uniform random number in the range |0<=u<x| or |0>=u>x| -or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here. - -Note that the call of |take_fraction| will produce the values 0 and~|x| -with about half the probability that it will produce any other particular -values between 0 and~|x|, because it rounds its answers. - -@c -static void mp_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig) { - mp_number y; /* trial value */ - mp_number x, abs_x; - mp_number u; - new_fraction (y); - new_number (x); - new_number (abs_x); - new_number (u); - mp_number_clone (&x, x_orig); - mp_number_clone (&abs_x, x); - mp_decimal_abs (&abs_x); - mp_next_unif_random(mp, &u); - decNumberMultiply (y.data.num, abs_x.data.num, u.data.num, &set); - free_number (u); - if (mp_number_equal(y, abs_x)) { - mp_number_clone (ret, ((math_data *)mp->math)->zero_t); - } else if (mp_number_greater(x, ((math_data *)mp->math)->zero_t)) { - mp_number_clone (ret, y); - } else { - mp_number_clone (ret, y); - mp_number_negate (ret); - } - free_number (abs_x); - free_number (x); - free_number (y); -} - - - -@ Finally, a normal deviate with mean zero and unit standard deviation -can readily be obtained with the ratio method (Algorithm 3.4.1R in -{\sl The Art of Computer Programming\/}). - -@c -static void mp_decimal_m_norm_rand (MP mp, mp_number *ret) { - mp_number ab_vs_cd; - mp_number abs_x; - mp_number u; - mp_number r; - mp_number la, xa; - new_number (ab_vs_cd); - new_number (la); - new_number (xa); - new_number (abs_x); - new_number (u); - new_number (r); - - do { - do { - mp_number v; - new_number (v); - mp_next_random(mp, &v); - mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t); - mp_decimal_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v); - free_number (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_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la); - mp_ab_vs_cd (mp,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa); - } while (mp_number_less(ab_vs_cd,((math_data *)mp->math)->zero_t)); - mp_number_clone (ret, xa); - free_number (ab_vs_cd); - free_number (r); - free_number (abs_x); - free_number (la); - free_number (xa); - free_number (u); -} - - - - -@ The following subroutine could be used in |norm_rand| and tests if $ab$ is -greater than, equal to, or less than~$cd$. -The result is $+1$, 0, or~$-1$ in the three respective cases. -This is not necessary, even if it's shorter than the current |ab_vs_cd| -and looks as a native implementation. - -@c -/* -|void mp_decimal_ab_vs_cd (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig, mp_number c_orig, mp_number d_orig) {| -| decNumber a, b, c, d;| -| decNumber ab, cd;| -| (void)mp;| -|| -| 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, &set);| -| decNumberMultiply (&cd, (decNumber *)c_orig.data.num, (decNumber *)d_orig.data.num, &set);| -| decNumberCompare(ret->data.num, &ab, &cd, &set);| -| mp_check_decNumber(mp, ret->data.num, &set);| -| return;| -|| -|}| -*/ - - - |