source/texk/web2c/mplibdir/mpmathdecimal.w

   1 % $Id: mpmathdecimal.w 1915 2013-06-13 10:17:31Z taco $
   2 %
   3 % This file is part of MetaPost;
   4 % the MetaPost program is in the public domain.
   5 % See the <Show version...> code in mpost.w for more info.
   6
   7 % Here is TeX material that gets inserted after \input webmac
   8
   9 \font\tenlogo=logo10 % font used for the METAFONT logo
  10 \font\logos=logosl10
  11 \def\MF{{\tenlogo META}\-{\tenlogo FONT}}
  12 \def\MP{{\tenlogo META}\-{\tenlogo POST}}
  13
  14 \def\title{Math support functions for decNumber based math}
  15 \pdfoutput=1
  16
  17 @ Introduction.
  18
  19 @c
  20 #include <w2c/config.h>
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23 #include <string.h>
  24 #include <math.h>
  25 #include "mpmathdecimal.h" /* internal header */
  26 #define ROUND(a) floor((a)+0.5)
  27 @h
  28
  29 @ @c
  30 @<Declarations@>;
  31
  32 @ @(mpmathdecimal.h@>=
  33 #ifndef MPMATHDECIMAL_H
  34 #define  MPMATHDECIMAL_H 1
  35 #include "mplib.h"
  36 #include "mpmp.h" /* internal header */
  37 #define  DECNUMDIGITS 1000
  38 #include "decNumber.h"
  39 @<Internal library declarations@>;
  40 #endif
  41
  42 @* Math initialization.
  43
  44 First, here are some very important constants.
  45
  46 @d E_STRING  "2.7182818284590452353602874713526624977572470936999595749669676277240766303535"
  47 @d PI_STRING "3.1415926535897932384626433832795028841971693993751058209749445923078164062862"
  48 @d fraction_multiplier 4096
  49 @d angle_multiplier 16
  50
  51 @ Here are the functions that are static as they are not used elsewhere
  52
  53 @<Declarations@>=
  54 #define DEBUG 0
  55 static void mp_decimal_scan_fractional_token (MP mp, int n);
  56 static void mp_decimal_scan_numeric_token (MP mp, int n);
  57 static void mp_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);
  58 /*static void mp_decimal_ab_vs_cd (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c, mp_number d);*/
  59 static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number a, mp_number b, mp_number c);
  60 static void mp_decimal_number_modulo (mp_number *a, mp_number b);
  61 static void mp_decimal_print_number (MP mp, mp_number n);
  62 static char * mp_decimal_number_tostring (MP mp, mp_number n);
  63 static void mp_decimal_slow_add (MP mp, mp_number *ret, mp_number x_orig, mp_number y_orig);
  64 static void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number x_orig);
  65 static void mp_decimal_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_sin);
  66 static void mp_init_randoms (MP mp, int seed);
  67 static void mp_number_angle_to_scaled (mp_number *A);
  68 static void mp_number_fraction_to_scaled (mp_number *A);
  69 static void mp_number_scaled_to_fraction (mp_number *A);
  70 static void mp_number_scaled_to_angle (mp_number *A);
  71 static void mp_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig);
  72 static void mp_decimal_m_norm_rand (MP mp, mp_number *ret);
  73 static void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number x_orig);
  74 static void mp_decimal_m_log (MP mp, mp_number *ret, mp_number x_orig);
  75 static void mp_decimal_pyth_sub (MP mp, mp_number *r, mp_number a, mp_number b);
  76 static void mp_decimal_pyth_add (MP mp, mp_number *r, mp_number a, mp_number b);
  77 static void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number x, mp_number y);
  78 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);
  79 static void mp_set_decimal_from_int(mp_number *A, int B);
  80 static void mp_set_decimal_from_boolean(mp_number *A, int B);
  81 static void mp_set_decimal_from_scaled(mp_number *A, int B);
  82 static void mp_set_decimal_from_addition(mp_number *A, mp_number B, mp_number C);
  83 static void mp_set_decimal_from_substraction (mp_number *A, mp_number B, mp_number C);
  84 static void mp_set_decimal_from_div(mp_number *A, mp_number B, mp_number C);
  85 static void mp_set_decimal_from_mul(mp_number *A, mp_number B, mp_number C);
  86 static void mp_set_decimal_from_int_div(mp_number *A, mp_number B, int C);
  87 static void mp_set_decimal_from_int_mul(mp_number *A, mp_number B, int C);
  88 static void mp_set_decimal_from_of_the_way(MP mp, mp_number *A, mp_number t, mp_number B, mp_number C);
  89 static void mp_number_negate(mp_number *A);
  90 static void mp_number_add(mp_number *A, mp_number B);
  91 static void mp_number_substract(mp_number *A, mp_number B);
  92 static void mp_number_half(mp_number *A);
  93 static void mp_number_halfp(mp_number *A);
  94 static void mp_number_double(mp_number *A);
  95 static void mp_number_add_scaled(mp_number *A, int B); /* also for negative B */
  96 static void mp_number_multiply_int(mp_number *A, int B);
  97 static void mp_number_divide_int(mp_number *A, int B);
  98 static void mp_decimal_abs(mp_number *A);
  99 static void mp_number_clone(mp_number *A, mp_number B);
 100 static void mp_number_swap(mp_number *A, mp_number *B);
 101 static int mp_round_unscaled(mp_number x_orig);
 102 static int mp_number_to_int(mp_number A);
 103 static int mp_number_to_scaled(mp_number A);
 104 static int mp_number_to_boolean(mp_number A);
 105 static double mp_number_to_double(mp_number A);
 106 static int mp_number_odd(mp_number A);
 107 static int mp_number_equal(mp_number A, mp_number B);
 108 static int mp_number_greater(mp_number A, mp_number B);
 109 static int mp_number_less(mp_number A, mp_number B);
 110 static int mp_number_nonequalabs(mp_number A, mp_number B);
 111 static void mp_number_floor (mp_number *i);
 112 static void mp_decimal_fraction_to_round_scaled (mp_number *x);
 113 static void mp_decimal_number_make_scaled (MP mp, mp_number *r, mp_number p, mp_number q);
 114 static void mp_decimal_number_make_fraction (MP mp, mp_number *r, mp_number p, mp_number q);
 115 static void mp_decimal_number_take_fraction (MP mp, mp_number *r, mp_number p, mp_number q);
 116 static void mp_decimal_number_take_scaled (MP mp, mp_number *r, mp_number p, mp_number q);
 117 static void mp_new_number (MP mp, mp_number *n, mp_number_type t) ;
 118 static void mp_free_number (MP mp, mp_number *n) ;
 119 static void mp_set_decimal_from_double(mp_number *A, double B);
 120 static void mp_free_decimal_math (MP mp);
 121 static void mp_decimal_set_precision (MP mp);
 122 static void mp_check_decNumber (MP mp, decNumber *dec, decContext *context);
 123 static int decNumber_check (decNumber *dec, decContext *context);
 124 static char * mp_decnumber_tostring (decNumber *n);
 125
 126 @ We do not want special numbers as return values for functions, so:
 127
 128
 129 @c
 130 int decNumber_check (decNumber *dec, decContext *context)
 131 {
 132    int test = false;
 133    if (context->status & DEC_Overflow) {
 134       test = true;
 135       context->status &= ~DEC_Overflow;
 136    }
 137    if (context->status & DEC_Underflow) {
 138       test = true;
 139       context->status &= ~DEC_Underflow;
 140    }
 141    if (context->status & DEC_Errors) {
 142 //     fprintf(stdout, "DEC_ERROR %x (%s)\n", context->status, decContextStatusToString(context));
 143      test = true;
 144      decNumberZero(dec);
 145    }
 146    context->status = 0;
 147    if (decNumberIsSpecial(dec)) {
 148       test = true;
 149       if (decNumberIsInfinite(dec)) {
 150         if (decNumberIsNegative(dec)) {
 151           decNumberCopyNegate(dec, &EL_GORDO_decNumber);
 152         } else {
 153           decNumberCopy(dec, &EL_GORDO_decNumber);
 154         }
 155       } else { // Nan
 156         decNumberZero(dec);
 157       }
 158    }
 159    if (decNumberIsZero(dec) && decNumberIsNegative(dec)) {
 160      decNumberZero(dec);
 161    }
 162    return test;
 163 }
 164 void mp_check_decNumber (MP mp, decNumber *dec, decContext *context)
 165 {
 166   mp->arith_error = decNumber_check (dec, context);
 167 }
 168
 169
 170
 171
 172 @ There are a few short decNumber functions that do not exist, but
 173 make life easier for us:
 174
 175 @d decNumberIsPositive(A) !(decNumberIsZero(A) || decNumberIsNegative(A))
 176
 177 @c
 178 static decContext set;
 179 static decContext limitedset;
 180 static void checkZero (decNumber *ret) {
 181   if (decNumberIsZero(ret) && decNumberIsNegative(ret))
 182     decNumberZero(ret);
 183 }
 184 static int decNumberLess(decNumber *a, decNumber *b) {
 185   decNumber comp;
 186   decNumberCompare(&comp, a, b, &set);
 187   return decNumberIsNegative(&comp);
 188 }
 189 static int decNumberGreater(decNumber *a, decNumber *b) {
 190   decNumber comp;
 191   decNumberCompare(&comp, a, b, &set);
 192   return decNumberIsPositive(&comp);
 193 }
 194 static void decNumberFromDouble(decNumber *A, double B) {
 195   char buf[1000];
 196   char *c;
 197   snprintf(buf,1000,"%-650.325lf",B);
 198   c = buf;
 199   while (*c++) {
 200     if (*c == ' ') {
 201       *c = '\0';
 202       break;
 203     }
 204   }
 205   decNumberFromString(A, buf, &set);
 206 }
 207 static double decNumberToDouble(decNumber *A) {
 208   char *buffer = malloc(A->digits + 14);
 209   double res = 0.0;
 210   assert (buffer);
 211   decNumberToString(A, buffer);
 212   if (sscanf(buffer, "%lf", &res)) {
 213      free(buffer);
 214      return res;
 215   } else {
 216      free(buffer);
 217      //mp->arith_error = 1;
 218      return 0.0; // whatever
 219   }
 220 }
 221
 222 @ Borrowed code from libdfp:
 223
 224                 x^3   x^5   x^7
 225 arctan(x) = x - --- + --- - --- + ...
 226                  3     5     7
 227
 228 This power series works well, if x is close to zero (|x|<0.5).
 229 If x is larger, the series converges too slowly,
 230 so in order to get a smaller x, we apply the identity
 231
 232                      sqrt(1+x^2) - 1
 233 arctan(x) = 2*arctan ---------------
 234                             x
 235
 236 twice. The first application gives us a new x with x < 1.
 237 The second application gives us a new x with x < 0.4142136.
 238 For that x, we use the power series and multiply the result by four.
 239
 240 @c
 241 static void decNumberAtan (decNumber *result, decNumber *x_orig, decContext *set)
 242 {
 243   decNumber x, f, g, mx2, term;
 244   int i;
 245   decNumberCopy(&x, x_orig);
 246   if (decNumberIsZero (&x)) {
 247     decNumberCopy (result, &x);
 248     return;
 249   }
 250   for (i=0; i<2; i++) {
 251     decNumber y;
 252     decNumberMultiply (&y, &x, &x, set);     // y = x^2
 253     decNumberAdd (&y, &y, &one, set);      // y = y+1
 254     decNumberSquareRoot (&y, &y, set);     // y = sqrt(y)
 255     decNumberSubtract (&y, &y, &one, set); // y = y-1
 256     decNumberDivide (&x, &y, &x, set);       // x = y/x
 257     if (decNumberIsZero (&x)) {
 258       decNumberCopy (result, &x);
 259       return;
 260     }
 261   }
 262   decNumberCopy (&f, &x);     // f(0) = x
 263   decNumberCopy (&g, &one);  // g(0) = 1
 264   decNumberCopy (&term, &x);  // term = x
 265   decNumberCopy (result, &x); // sum  = x
 266   decNumberMultiply (&mx2, &x, &x, set); // mx2 = x^2
 267   decNumberMinus (&mx2, &mx2, set);    // mx2 = -x^2
 268   for (i=0; i<2*set->digits; i++) {
 269     decNumberMultiply (&f, &f, &mx2, set);
 270     decNumberAdd (&g, &g, &two_decNumber, set);
 271     decNumberDivide (&term, &f, &g, set);
 272     decNumberAdd (result, result, &term, set);
 273   }
 274   decNumberAdd (result, result, result, set);
 275   decNumberAdd (result, result, result, set);
 276   return;
 277 }
 278 static void decNumberAtan2 (decNumber *result, decNumber *y, decNumber *x, decContext *set)
 279 {
 280   decNumber temp;
 281   if (!decNumberIsInfinite (x) && !decNumberIsZero (y)
 282       && !decNumberIsInfinite (y) && !decNumberIsZero (x)) {
 283       decNumberDivide (&temp, y, x, set);
 284       decNumberAtan (result, &temp, set);
 285       /* decNumberAtan doesn't quite return the values in the ranges we
 286        * want for x < 0. So we need to do some correction */
 287       if (decNumberIsNegative (x)) {
 288           if (decNumberIsNegative (y)) {
 289             decNumberSubtract(result, result, &PI_decNumber, set);
 290           } else {
 291             decNumberAdd(result, result, &PI_decNumber, set);
 292           }
 293       }
 294       return;
 295   }
 296   if (decNumberIsInfinite (y) && decNumberIsInfinite (x)) {
 297       /* If x and y are both inf, the result depends on the sign of x */
 298       decNumberDivide(result, &PI_decNumber, &four_decNumber, set);
 299       if (decNumberIsNegative (x) ) {
 300         decNumber a;
 301         decNumberFromDouble(&a, 3.0);
 302         decNumberMultiply(result, result, &a, set);
 303       }
 304   } else if (!decNumberIsZero (y) && !decNumberIsInfinite (x) ) {
 305     /*  If y is non-zero and x is non-inf, the result is +-pi/2 */
 306     decNumberDivide(result, &PI_decNumber, &two_decNumber, set);
 307   } else { /*  Otherwise it is +0 if x is positive, +pi if x is neg */
 308       if (decNumberIsNegative (x)) {
 309         decNumberCopy(result, &PI_decNumber);
 310       } else {
 311         decNumberZero(result);
 312       }
 313   }
 314   /*  Atan2 will be negative if y<0 */
 315   if (decNumberIsNegative (y)) {
 316     decNumberMinus(result, result, set);
 317   }
 318 }
 319
 320 @ And these are the ones that {\it are} used elsewhere
 321
 322 @<Internal library declarations@>=
 323 void * mp_initialize_decimal_math (MP mp);
 324
 325 @
 326
 327 @d unity   1
 328 @d two 2
 329 @d three 3
 330 @d four 4
 331 @d half_unit  0.5
 332 @d three_quarter_unit 0.75
 333 @d coef_bound ((7.0/3.0)*fraction_multiplier) /* |fraction| approximation to 7/3 */
 334 @d fraction_threshold 0.04096 /* a |fraction| coefficient less than this is zeroed */
 335 @d half_fraction_threshold (fraction_threshold/2) /* half of |fraction_threshold| */
 336 @d scaled_threshold 0.000122 /* a |scaled| coefficient less than this is zeroed */
 337 @d half_scaled_threshold (scaled_threshold/2) /* half of |scaled_threshold| */
 338 @d near_zero_angle (0.0256*angle_multiplier)  /* an angle of about 0.0256 */
 339 @d p_over_v_threshold 0x80000 /* TODO */
 340 @d equation_threshold 0.001
 341 @d tfm_warn_threshold 0.0625
 342 @d epsilon "1E-52"
 343 @d epsilonf pow(2.0,-52.0)
 344 @d EL_GORDO     "1E1000000" /* the largest value that \MP\ likes. */
 345 @d warning_limit "1E1000000" /* this is a large value that can just be expressed without loss of precision */
 346 @d DECPRECISION_DEFAULT 34
 347
 348 @<Declarations@>=
 349 static decNumber zero;
 350 static decNumber one;
 351 static decNumber minusone;
 352 static decNumber two_decNumber;
 353 static decNumber three_decNumber;
 354 static decNumber four_decNumber;
 355 static decNumber fraction_multiplier_decNumber;
 356 static decNumber angle_multiplier_decNumber;
 357 static decNumber fraction_one_decNumber;
 358 static decNumber fraction_one_plus_decNumber;
 359 static decNumber PI_decNumber;
 360 static decNumber epsilon_decNumber;
 361 static decNumber EL_GORDO_decNumber;
 362 static decNumber **factorials = NULL;
 363 static int last_cached_factorial = 0;
 364 static boolean initialized = false ;
 365 @ @c
 366 void * mp_initialize_decimal_math (MP mp) {
 367   math_data *math = (math_data *)mp_xmalloc(mp,1,sizeof(math_data));
 368   // various decNumber initializations
 369   decContextDefault(&set, DEC_INIT_BASE); // initialize
 370   set.traps=0;                     // no traps, thank you
 371   decContextDefault(&limitedset, DEC_INIT_BASE); // initialize
 372   limitedset.traps=0;                     // no traps, thank you
 373   limitedset.emax = 999999;
 374   limitedset.emin = -999999;
 375   set.digits = DECPRECISION_DEFAULT;
 376   limitedset.digits = DECPRECISION_DEFAULT;
 377   if (!initialized) {
 378     initialized = true ;
 379     decNumberFromInt32(&one, 1);
 380     decNumberFromInt32(&minusone, -1);
 381     decNumberFromInt32(&zero, 0);
 382     decNumberFromInt32(&two_decNumber, two);
 383     decNumberFromInt32(&three_decNumber, three);
 384     decNumberFromInt32(&four_decNumber, four);
 385     decNumberFromInt32(&fraction_multiplier_decNumber, fraction_multiplier);
 386     decNumberFromInt32(&fraction_one_decNumber, fraction_one);
 387     decNumberFromInt32(&fraction_one_plus_decNumber, (fraction_one+1));
 388     decNumberFromInt32(&angle_multiplier_decNumber, angle_multiplier);
 389     decNumberFromString(&PI_decNumber, PI_STRING, &set);
 390     decNumberFromString(&epsilon_decNumber, epsilon, &set);
 391     decNumberFromString(&EL_GORDO_decNumber, EL_GORDO, &set);
 392     factorials = (decNumber **)mp_xmalloc(mp,PRECALC_FACTORIALS_CACHESIZE,sizeof(decNumber *));
 393     factorials[0] = (decNumber *)mp_xmalloc(mp,1,sizeof(decNumber));
 394     decNumberCopy(factorials[0], &one);
 395   }
 396
 397   /* alloc */
 398   math->allocate = mp_new_number;
 399   math->free = mp_free_number;
 400   mp_new_number (mp, &math->precision_default, mp_scaled_type);
 401   decNumberFromInt32(math->precision_default.data.num, DECPRECISION_DEFAULT);
 402   mp_new_number (mp, &math->precision_max, mp_scaled_type);
 403   decNumberFromInt32(math->precision_max.data.num, DECNUMDIGITS);
 404   mp_new_number (mp, &math->precision_min, mp_scaled_type);
 405   decNumberFromInt32(math->precision_min.data.num, 1);
 406   /* here are the constants for |scaled| objects */
 407   mp_new_number (mp, &math->epsilon_t, mp_scaled_type);
 408   decNumberCopy(math->epsilon_t.data.num, &epsilon_decNumber);
 409   mp_new_number (mp, &math->inf_t, mp_scaled_type);
 410   decNumberCopy(math->inf_t.data.num, &EL_GORDO_decNumber);
 411   mp_new_number (mp, &math->warning_limit_t, mp_scaled_type);
 412   decNumberFromString(math->warning_limit_t.data.num, warning_limit, &set);
 413   mp_new_number (mp, &math->one_third_inf_t, mp_scaled_type);
 414   decNumberDivide(math->one_third_inf_t.data.num, math->inf_t.data.num, &three_decNumber, &set);
 415   mp_new_number (mp, &math->unity_t, mp_scaled_type);
 416   decNumberCopy(math->unity_t.data.num, &one);
 417   mp_new_number (mp, &math->two_t, mp_scaled_type);
 418   decNumberFromInt32(math->two_t.data.num, two);
 419   mp_new_number (mp, &math->three_t, mp_scaled_type);
 420   decNumberFromInt32(math->three_t.data.num, three);
 421   mp_new_number (mp, &math->half_unit_t, mp_scaled_type);
 422   decNumberFromString(math->half_unit_t.data.num, "0.5", &set);
 423   mp_new_number (mp, &math->three_quarter_unit_t, mp_scaled_type);
 424   decNumberFromString(math->three_quarter_unit_t.data.num, "0.75", &set);
 425   mp_new_number (mp, &math->zero_t, mp_scaled_type);
 426   decNumberZero(math->zero_t.data.num);
 427   /* |fractions| */
 428   mp_new_number (mp, &math->arc_tol_k, mp_fraction_type);
 429   {
 430      decNumber fourzeroninesix;
 431      decNumberFromInt32(&fourzeroninesix, 4096);
 432      decNumberDivide(math->arc_tol_k.data.num, &one, &fourzeroninesix, &set);
 433      /* quit when change in arc length estimate reaches this */
 434   }
 435   mp_new_number (mp, &math->fraction_one_t, mp_fraction_type);
 436   decNumberFromInt32(math->fraction_one_t.data.num, fraction_one);
 437   mp_new_number (mp, &math->fraction_half_t, mp_fraction_type);
 438   decNumberFromInt32(math->fraction_half_t.data.num, fraction_half);
 439   mp_new_number (mp, &math->fraction_three_t, mp_fraction_type);
 440   decNumberFromInt32(math->fraction_three_t.data.num, fraction_three);
 441   mp_new_number (mp, &math->fraction_four_t, mp_fraction_type);
 442   decNumberFromInt32(math->fraction_four_t.data.num, fraction_four);
 443   /* |angles| */
 444   mp_new_number (mp, &math->three_sixty_deg_t, mp_angle_type);
 445   decNumberFromInt32(math->three_sixty_deg_t.data.num, 360  * angle_multiplier);
 446   mp_new_number (mp, &math->one_eighty_deg_t, mp_angle_type);
 447   decNumberFromInt32(math->one_eighty_deg_t.data.num, 180 * angle_multiplier);
 448   /* various approximations */
 449   mp_new_number (mp, &math->one_k, mp_scaled_type);
 450   decNumberFromDouble(math->one_k.data.num, 1.0/64);
 451   mp_new_number (mp, &math->sqrt_8_e_k, mp_scaled_type);
 452   {
 453     decNumberFromDouble(math->sqrt_8_e_k.data.num, 112428.82793 / 65536.0);
 454     /* $2^{16}\sqrt{8/e}\approx 112428.82793$ */
 455   }
 456   mp_new_number (mp, &math->twelve_ln_2_k, mp_fraction_type);
 457   {
 458     decNumberFromDouble(math->twelve_ln_2_k.data.num, 139548959.6165 / 65536.0);
 459     /* $2^{24}\cdot12\ln2\approx139548959.6165$ */
 460   }
 461   mp_new_number (mp, &math->coef_bound_k, mp_fraction_type);
 462   decNumberFromDouble(math->coef_bound_k.data.num,coef_bound);
 463   mp_new_number (mp, &math->coef_bound_minus_1, mp_fraction_type);
 464   decNumberFromDouble(math->coef_bound_minus_1.data.num,coef_bound - 1 / 65536.0);
 465   mp_new_number (mp, &math->twelvebits_3, mp_scaled_type);
 466   {
 467     decNumberFromDouble(math->twelvebits_3.data.num, 1365 / 65536.0);
 468     /* $1365\approx 2^{12}/3$ */
 469   }
 470   mp_new_number (mp, &math->twentysixbits_sqrt2_t, mp_fraction_type);
 471   {
 472     decNumberFromDouble(math->twentysixbits_sqrt2_t.data.num, 94906265.62 / 65536.0);
 473     /* $2^{26}\sqrt2\approx94906265.62$ */
 474   }
 475   mp_new_number (mp, &math->twentyeightbits_d_t, mp_fraction_type);
 476   {
 477     decNumberFromDouble(math->twentyeightbits_d_t.data.num, 35596754.69  / 65536.0);
 478     /* $2^{28}d\approx35596754.69$ */
 479   }
 480   mp_new_number (mp, &math->twentysevenbits_sqrt2_d_t, mp_fraction_type);
 481   {
 482     decNumberFromDouble(math->twentysevenbits_sqrt2_d_t.data.num, 25170706.63  / 65536.0);
 483     /* $2^{27}\sqrt2\,d\approx25170706.63$ */
 484   }
 485   /* thresholds */
 486   mp_new_number (mp, &math->fraction_threshold_t, mp_fraction_type);
 487   decNumberFromDouble(math->fraction_threshold_t.data.num, fraction_threshold);
 488   mp_new_number (mp, &math->half_fraction_threshold_t, mp_fraction_type);
 489   decNumberFromDouble(math->half_fraction_threshold_t.data.num, half_fraction_threshold);
 490   mp_new_number (mp, &math->scaled_threshold_t, mp_scaled_type);
 491   decNumberFromDouble(math->scaled_threshold_t.data.num, scaled_threshold);
 492   mp_new_number (mp, &math->half_scaled_threshold_t, mp_scaled_type);
 493   decNumberFromDouble(math->half_scaled_threshold_t.data.num, half_scaled_threshold);
 494   mp_new_number (mp, &math->near_zero_angle_t, mp_angle_type);
 495   decNumberFromDouble(math->near_zero_angle_t.data.num, near_zero_angle);
 496   mp_new_number (mp, &math->p_over_v_threshold_t, mp_fraction_type);
 497   decNumberFromDouble(math->p_over_v_threshold_t.data.num, p_over_v_threshold);
 498   mp_new_number (mp, &math->equation_threshold_t, mp_scaled_type);
 499   decNumberFromDouble(math->equation_threshold_t.data.num, equation_threshold);
 500   mp_new_number (mp, &math->tfm_warn_threshold_t, mp_scaled_type);
 501   decNumberFromDouble(math->tfm_warn_threshold_t.data.num, tfm_warn_threshold);
 502   /* functions */
 503   math->from_int = mp_set_decimal_from_int;
 504   math->from_boolean = mp_set_decimal_from_boolean;
 505   math->from_scaled = mp_set_decimal_from_scaled;
 506   math->from_double = mp_set_decimal_from_double;
 507   math->from_addition  = mp_set_decimal_from_addition;
 508   math->from_substraction  = mp_set_decimal_from_substraction;
 509   math->from_oftheway  = mp_set_decimal_from_of_the_way;
 510   math->from_div  = mp_set_decimal_from_div;
 511   math->from_mul  = mp_set_decimal_from_mul;
 512   math->from_int_div  = mp_set_decimal_from_int_div;
 513   math->from_int_mul  = mp_set_decimal_from_int_mul;
 514   math->negate = mp_number_negate;
 515   math->add  = mp_number_add;
 516   math->substract = mp_number_substract;
 517   math->half = mp_number_half;
 518   math->halfp = mp_number_halfp;
 519   math->do_double = mp_number_double;
 520   math->abs = mp_decimal_abs;
 521   math->clone = mp_number_clone;
 522   math->swap = mp_number_swap;
 523   math->add_scaled = mp_number_add_scaled;
 524   math->multiply_int = mp_number_multiply_int;
 525   math->divide_int = mp_number_divide_int;
 526   math->to_boolean = mp_number_to_boolean;
 527   math->to_scaled = mp_number_to_scaled;
 528   math->to_double = mp_number_to_double;
 529   math->to_int = mp_number_to_int;
 530   math->odd = mp_number_odd;
 531   math->equal = mp_number_equal;
 532   math->less = mp_number_less;
 533   math->greater = mp_number_greater;
 534   math->nonequalabs = mp_number_nonequalabs;
 535   math->round_unscaled = mp_round_unscaled;
 536   math->floor_scaled = mp_number_floor;
 537   math->fraction_to_round_scaled = mp_decimal_fraction_to_round_scaled;
 538   math->make_scaled = mp_decimal_number_make_scaled;
 539   math->make_fraction = mp_decimal_number_make_fraction;
 540   math->take_fraction = mp_decimal_number_take_fraction;
 541   math->take_scaled = mp_decimal_number_take_scaled;
 542   math->velocity = mp_decimal_velocity;
 543   math->n_arg = mp_decimal_n_arg;
 544   math->m_log = mp_decimal_m_log;
 545   math->m_exp = mp_decimal_m_exp;
 546   math->m_unif_rand = mp_decimal_m_unif_rand;
 547   math->m_norm_rand = mp_decimal_m_norm_rand;
 548   math->pyth_add = mp_decimal_pyth_add;
 549   math->pyth_sub = mp_decimal_pyth_sub;
 550   math->fraction_to_scaled = mp_number_fraction_to_scaled;
 551   math->scaled_to_fraction = mp_number_scaled_to_fraction;
 552   math->scaled_to_angle = mp_number_scaled_to_angle;
 553   math->angle_to_scaled = mp_number_angle_to_scaled;
 554   math->init_randoms = mp_init_randoms;
 555   math->sin_cos = mp_decimal_sin_cos;
 556   math->slow_add = mp_decimal_slow_add;
 557   math->sqrt = mp_decimal_square_rt;
 558   math->print = mp_decimal_print_number;
 559   math->tostring = mp_decimal_number_tostring;
 560   math->modulo = mp_decimal_number_modulo;
 561   math->ab_vs_cd = mp_ab_vs_cd;
 562   math->crossing_point = mp_decimal_crossing_point;
 563   math->scan_numeric = mp_decimal_scan_numeric_token;
 564   math->scan_fractional = mp_decimal_scan_fractional_token;
 565   math->free_math = mp_free_decimal_math;
 566   math->set_precision = mp_decimal_set_precision;
 567   return (void *)math;
 568 }
 569
 570 void mp_decimal_set_precision (MP mp) {
 571   int i;
 572   i = decNumberToInt32((decNumber *)internal_value (mp_number_precision).data.num, &set);
 573   set.digits = i;
 574   limitedset.digits = i;
 575 }
 576
 577 void mp_free_decimal_math (MP mp) {
 578   int i;
 579   free_number (((math_data *)mp->math)->three_sixty_deg_t);
 580   free_number (((math_data *)mp->math)->one_eighty_deg_t);
 581   free_number (((math_data *)mp->math)->fraction_one_t);
 582   free_number (((math_data *)mp->math)->zero_t);
 583   free_number (((math_data *)mp->math)->half_unit_t);
 584   free_number (((math_data *)mp->math)->three_quarter_unit_t);
 585   free_number (((math_data *)mp->math)->unity_t);
 586   free_number (((math_data *)mp->math)->two_t);
 587   free_number (((math_data *)mp->math)->three_t);
 588   free_number (((math_data *)mp->math)->one_third_inf_t);
 589   free_number (((math_data *)mp->math)->inf_t);
 590   free_number (((math_data *)mp->math)->warning_limit_t);
 591   free_number (((math_data *)mp->math)->one_k);
 592   free_number (((math_data *)mp->math)->sqrt_8_e_k);
 593   free_number (((math_data *)mp->math)->twelve_ln_2_k);
 594   free_number (((math_data *)mp->math)->coef_bound_k);
 595   free_number (((math_data *)mp->math)->coef_bound_minus_1);
 596   free_number (((math_data *)mp->math)->fraction_threshold_t);
 597   free_number (((math_data *)mp->math)->half_fraction_threshold_t);
 598   free_number (((math_data *)mp->math)->scaled_threshold_t);
 599   free_number (((math_data *)mp->math)->half_scaled_threshold_t);
 600   free_number (((math_data *)mp->math)->near_zero_angle_t);
 601   free_number (((math_data *)mp->math)->p_over_v_threshold_t);
 602   free_number (((math_data *)mp->math)->equation_threshold_t);
 603   free_number (((math_data *)mp->math)->tfm_warn_threshold_t);
 604   for (i = 0; i <= last_cached_factorial; i++) {
 605     free(factorials[i]);
 606   }
 607   free(factorials);
 608   free(mp->math);
 609 }
 610
 611 @ Creating an destroying |mp_number| objects
 612
 613 @ @c
 614 void mp_new_number (MP mp, mp_number *n, mp_number_type t) {
 615   (void)mp;
 616   n->data.num = mp_xmalloc(mp,1,sizeof(decNumber));
 617   decNumberZero(n->data.num);
 618   n->type = t;
 619 }
 620
 621 @
 622
 623 @c
 624 void mp_free_number (MP mp, mp_number *n) {
 625   (void)mp;
 626   free(n->data.num);
 627   n->data.num = NULL;
 628   n->type = mp_nan_type;
 629 }
 630
 631 @ Here are the low-level functions on |mp_number| items, setters first.
 632
 633 @c
 634 void mp_set_decimal_from_int(mp_number *A, int B) {
 635   decNumberFromInt32(A->data.num,B);
 636 }
 637 void mp_set_decimal_from_boolean(mp_number *A, int B) {
 638   decNumberFromInt32(A->data.num,B);
 639 }
 640 void mp_set_decimal_from_scaled(mp_number *A, int B) {
 641   decNumber c;
 642   decNumberFromInt32(&c, 65536);
 643   decNumberFromInt32(A->data.num,B);
 644   decNumberDivide(A->data.num,A->data.num,&c, &set);
 645 }
 646 void mp_set_decimal_from_double(mp_number *A, double B) {
 647   decNumberFromDouble(A->data.num, B);
 648 }
 649 void mp_set_decimal_from_addition(mp_number *A, mp_number B, mp_number C) {
 650   decNumberAdd(A->data.num,B.data.num,C.data.num, &set);
 651 }
 652 void mp_set_decimal_from_substraction (mp_number *A, mp_number B, mp_number C) {
 653   decNumberSubtract(A->data.num,B.data.num,C.data.num, &set);
 654 }
 655 void mp_set_decimal_from_div(mp_number *A, mp_number B, mp_number C) {
 656   decNumberDivide(A->data.num,B.data.num,C.data.num, &set);
 657 }
 658 void mp_set_decimal_from_mul(mp_number *A, mp_number B, mp_number C) {
 659   decNumberMultiply(A->data.num,B.data.num,C.data.num, &set);
 660 }
 661 void mp_set_decimal_from_int_div(mp_number *A, mp_number B, int C) {
 662   decNumber c;
 663   decNumberFromInt32(&c, C);
 664   decNumberDivide(A->data.num,B.data.num,&c, &set);
 665 }
 666 void mp_set_decimal_from_int_mul(mp_number *A, mp_number B, int C) {
 667   decNumber c;
 668   decNumberFromInt32(&c, C);
 669   decNumberMultiply(A->data.num,B.data.num,&c, &set);
 670 }
 671 void mp_set_decimal_from_of_the_way(MP mp, mp_number *A, mp_number t, mp_number B, mp_number C) {
 672   decNumber c;
 673   decNumber r1;
 674   decNumberSubtract(&c,B.data.num, C.data.num, &set);
 675   mp_decimal_take_fraction(mp, &r1, &c, t.data.num);
 676   decNumberSubtract(A->data.num, B.data.num, &r1, &set);
 677   mp_check_decNumber(mp, A->data.num, &set);
 678 }
 679 void mp_number_negate(mp_number *A) {
 680   decNumberCopyNegate(A->data.num, A->data.num);
 681   checkZero(A->data.num);
 682 }
 683 void mp_number_add(mp_number *A, mp_number B) {
 684   decNumberAdd(A->data.num,A->data.num,B.data.num, &set);
 685 }
 686 void mp_number_substract(mp_number *A, mp_number B) {
 687   decNumberSubtract(A->data.num,A->data.num,B.data.num, &set);
 688 }
 689 void mp_number_half(mp_number *A) {
 690   decNumber c;
 691   decNumberFromInt32(&c, 2);
 692   decNumberDivide(A->data.num,A->data.num, &c, &set);
 693 }
 694 void mp_number_halfp(mp_number *A) {
 695   decNumber c;
 696   decNumberFromInt32(&c, 2);
 697   decNumberDivide(A->data.num,A->data.num, &c, &set);
 698 }
 699 void mp_number_double(mp_number *A) {
 700   decNumber c;
 701   decNumberFromInt32(&c, 2);
 702   decNumberMultiply(A->data.num,A->data.num, &c, &set);
 703 }
 704 void mp_number_add_scaled(mp_number *A, int B) { /* also for negative B */
 705   decNumber b,c;
 706   decNumberFromInt32(&c, 65536);
 707   decNumberFromInt32(&b, B);
 708   decNumberDivide(&b,&b, &c, &set);
 709   decNumberAdd(A->data.num,A->data.num, &b, &set);
 710 }
 711 void mp_number_multiply_int(mp_number *A, int B) {
 712   decNumber b;
 713   decNumberFromInt32(&b, B);
 714   decNumberMultiply(A->data.num,A->data.num, &b, &set);
 715 }
 716 void mp_number_divide_int(mp_number *A, int B) {
 717   decNumber b;
 718   decNumberFromInt32(&b, B);
 719   decNumberDivide(A->data.num,A->data.num,&b, &set);
 720 }
 721 void mp_decimal_abs(mp_number *A) {
 722   decNumberAbs(A->data.num, A->data.num, &set);
 723 }
 724 void mp_number_clone(mp_number *A, mp_number B) {
 725   decNumberCopy(A->data.num, B.data.num);
 726 }
 727 void mp_number_swap(mp_number *A, mp_number *B) {
 728   decNumber swap_tmp;
 729   decNumberCopy(&swap_tmp,   A->data.num);
 730   decNumberCopy(A->data.num, B->data.num);
 731   decNumberCopy(B->data.num, &swap_tmp);
 732 }
 733 void mp_number_fraction_to_scaled (mp_number *A) {
 734     A->type = mp_scaled_type;
 735     decNumberDivide(A->data.num, A->data.num, &fraction_multiplier_decNumber, &set);
 736 }
 737 void mp_number_angle_to_scaled (mp_number *A) {
 738     A->type = mp_scaled_type;
 739     decNumberDivide(A->data.num, A->data.num, &angle_multiplier_decNumber, &set);
 740 }
 741 void mp_number_scaled_to_fraction (mp_number *A) {
 742     A->type = mp_fraction_type;
 743     decNumberMultiply(A->data.num, A->data.num, &fraction_multiplier_decNumber, &set);
 744 }
 745 void mp_number_scaled_to_angle (mp_number *A) {
 746     A->type = mp_angle_type;
 747     decNumberMultiply(A->data.num, A->data.num, &angle_multiplier_decNumber, &set);
 748 }
 749
 750
 751 @* Query functions
 752
 753 @ Convert a number to a scaled value. |decNumberToInt32| is not
 754 able to make this conversion properly, so instead we are using
 755 |decNumberToDouble| and a typecast. Bad!
 756
 757 @c
 758 int mp_number_to_scaled(mp_number A) {
 759   int32_t result;
 760   decNumber corrected;
 761   decNumberFromInt32(&corrected, 65536);
 762   decNumberMultiply(&corrected,&corrected,A.data.num, &set);
 763   decNumberReduce(&corrected, &corrected, &set);
 764   result = (int)floor(decNumberToDouble(&corrected)+0.5);
 765   return result;
 766 }
 767
 768 @
 769
 770 @d odd(A)   (abs(A)%2==1)
 771
 772 @c
 773 int mp_number_to_int(mp_number A) {
 774   int32_t result;
 775   set.status = 0;
 776   result = decNumberToInt32(A.data.num, &set);
 777   if (set.status == DEC_Invalid_operation) {
 778      set.status = 0;
 779      // mp->arith_error = 1;
 780      return 0; // whatever
 781   } else {
 782      return result;
 783   }
 784 }
 785 int mp_number_to_boolean(mp_number A) {
 786   uint32_t result;
 787   set.status = 0;
 788   result = decNumberToUInt32(A.data.num, &set);
 789   if (set.status == DEC_Invalid_operation) {
 790      set.status = 0;
 791      // mp->arith_error = 1;
 792      return mp_false_code; // whatever
 793   } else {
 794      return result ;
 795   }
 796 }
 797 double mp_number_to_double(mp_number A) {
 798   char *buffer = malloc(((decNumber *)A.data.num)->digits + 14);
 799   double res = 0.0;
 800   assert (buffer);
 801   decNumberToString(A.data.num, buffer);
 802   if (sscanf(buffer, "%lf", &res)) {
 803      free(buffer);
 804      return res;
 805   } else {
 806      free(buffer);
 807      //mp->arith_error = 1;
 808      return 0.0; // whatever
 809   }
 810 }
 811 int mp_number_odd(mp_number A) {
 812   return odd(mp_number_to_int(A));
 813 }
 814 int mp_number_equal(mp_number A, mp_number B) {
 815   decNumber res;
 816   decNumberCompare(&res,A.data.num,B.data.num, &set);
 817   return decNumberIsZero(&res);
 818 }
 819 int mp_number_greater(mp_number A, mp_number B) {
 820   decNumber res;
 821   decNumberCompare(&res,A.data.num,B.data.num, &set);
 822   return decNumberIsPositive(&res);
 823 }
 824 int mp_number_less(mp_number A, mp_number B) {
 825   decNumber res;
 826   decNumberCompare(&res,A.data.num,B.data.num, &set);
 827   return decNumberIsNegative(&res);
 828 }
 829 int mp_number_nonequalabs(mp_number A, mp_number B) {
 830   decNumber res, a, b;
 831   decNumberCopyAbs(&a, A.data.num);
 832   decNumberCopyAbs(&b, B.data.num);
 833   decNumberCompare(&res, &a, &b, &set);
 834   return !decNumberIsZero(&res);
 835 }
 836
 837 @ Fixed-point arithmetic is done on {\sl scaled integers\/} that are multiples
 838 of $2^{-16}$. In other words, a binary point is assumed to be sixteen bit
 839 positions from the right end of a binary computer word.
 840
 841 @ One of \MP's most common operations is the calculation of
 842 $\lfloor{a+b\over2}\rfloor$,
 843 the midpoint of two given integers |a| and~|b|. The most decent way to do
 844 this is to write `|(a+b)/2|'; but on many machines it is more efficient
 845 to calculate `|(a+b)>>1|'.
 846
 847 Therefore the midpoint operation will always be denoted by `|half(a+b)|'
 848 in this program. If \MP\ is being implemented with languages that permit
 849 binary shifting, the |half| macro should be changed to make this operation
 850 as efficient as possible.  Since some systems have shift operators that can
 851 only be trusted to work on positive numbers, there is also a macro |halfp|
 852 that is used only when the quantity being halved is known to be positive
 853 or zero.
 854
 855 @ Here is a procedure analogous to |print_int|.  The current version
 856 is fairly stupid, and it is not round-trip safe, but this is good
 857 enough for a beta test.
 858
 859 @c
 860 char * mp_decnumber_tostring (decNumber *n) {
 861   decNumber corrected;
 862   char *buffer = malloc(((decNumber *)n)->digits + 14);
 863   assert (buffer);
 864   decNumberCopy(&corrected,n);
 865   decNumberTrim(&corrected);
 866   decNumberToString(&corrected, buffer);
 867   return buffer;
 868 }
 869 char * mp_decimal_number_tostring (MP mp, mp_number n) {
 870   return mp_decnumber_tostring(n.data.num);
 871 }
 872
 873
 874 @ @c
 875 void mp_decimal_print_number (MP mp, mp_number n) {
 876   char *str = mp_decimal_number_tostring(mp, n);
 877   mp_print (mp, str);
 878   free (str);
 879 }
 880
 881
 882
 883
 884 @ Addition is not always checked to make sure that it doesn't overflow,
 885 but in places where overflow isn't too unlikely the |slow_add| routine
 886 is used.
 887
 888 @c
 889 void mp_decimal_slow_add (MP mp, mp_number *ret, mp_number A, mp_number B) {
 890   decNumberAdd(ret->data.num,A.data.num,B.data.num, &set);
 891 }
 892
 893 @ The |make_fraction| routine produces the |fraction| equivalent of
 894 |p/q|, given integers |p| and~|q|; it computes the integer
 895 $f=\lfloor2^{28}p/q+{1\over2}\rfloor$, when $p$ and $q$ are
 896 positive. If |p| and |q| are both of the same scaled type |t|,
 897 the ``type relation'' |make_fraction(t,t)=fraction| is valid;
 898 and it's also possible to use the subroutine ``backwards,'' using
 899 the relation |make_fraction(t,fraction)=t| between scaled types.
 900
 901 If the result would have magnitude $2^{31}$ or more, |make_fraction|
 902 sets |arith_error:=true|. Most of \MP's internal computations have
 903 been designed to avoid this sort of error.
 904
 905 If this subroutine were programmed in assembly language on a typical
 906 machine, we could simply compute |(@t$2^{28}$@>*p)div q|, since a
 907 double-precision product can often be input to a fixed-point division
 908 instruction. But when we are restricted to int-eger arithmetic it
 909 is necessary either to resort to multiple-precision maneuvering
 910 or to use a simple but slow iteration. The multiple-precision technique
 911 would be about three times faster than the code adopted here, but it
 912 would be comparatively long and tricky, involving about sixteen
 913 additional multiplications and divisions.
 914
 915 This operation is part of \MP's ``inner loop''; indeed, it will
 916 consume nearly 10\pct! of the running time (exclusive of input and output)
 917 if the code below is left unchanged. A machine-dependent recoding
 918 will therefore make \MP\ run faster. The present implementation
 919 is highly portable, but slow; it avoids multiplication and division
 920 except in the initial stage. System wizards should be careful to
 921 replace it with a routine that is guaranteed to produce identical
 922 results in all cases.
 923 @^system dependencies@>
 924
 925 As noted below, a few more routines should also be replaced by machine-dependent
 926 code, for efficiency. But when a procedure is not part of the ``inner loop,''
 927 such changes aren't advisable; simplicity and robustness are
 928 preferable to trickery, unless the cost is too high.
 929 @^inner loop@>
 930
 931 @c
 932 void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) {
 933   decNumberDivide(ret, p, q, &set);
 934   mp_check_decNumber(mp, ret, &set);
 935   decNumberMultiply(ret, ret, &fraction_multiplier_decNumber, &set);
 936 }
 937 void mp_decimal_number_make_fraction (MP mp, mp_number *ret, mp_number p, mp_number q) {
 938   mp_decimal_make_fraction (mp, ret->data.num, p.data.num, q.data.num);
 939 }
 940
 941 @ @<Declarations@>=
 942 void mp_decimal_make_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q);
 943
 944 @ The dual of |make_fraction| is |take_fraction|, which multiplies a
 945 given integer~|q| by a fraction~|f|. When the operands are positive, it
 946 computes $p=\lfloor qf/2^{28}+{1\over2}\rfloor$, a symmetric function
 947 of |q| and~|f|.
 948
 949 This routine is even more ``inner loopy'' than |make_fraction|;
 950 the present implementation consumes almost 20\pct! of \MP's computation
 951 time during typical jobs, so a machine-language substitute is advisable.
 952 @^inner loop@> @^system dependencies@>
 953
 954 @c
 955 void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q) {
 956   decNumberMultiply(ret, p, q, &set);
 957   decNumberDivide(ret, ret, &fraction_multiplier_decNumber, &set);
 958 }
 959 void mp_decimal_number_take_fraction (MP mp, mp_number *ret, mp_number p, mp_number q) {
 960   mp_decimal_take_fraction (mp, ret->data.num, p.data.num, q.data.num);
 961 }
 962
 963 @ @<Declarations@>=
 964 void mp_decimal_take_fraction (MP mp, decNumber *ret, decNumber *p, decNumber *q);
 965
 966 @ When we want to multiply something by a |scaled| quantity, we use a scheme
 967 analogous to |take_fraction| but with a different scaling.
 968 Given positive operands, |take_scaled|
 969 computes the quantity $p=\lfloor qf/2^{16}+{1\over2}\rfloor$.
 970
 971 Once again it is a good idea to use a machine-language replacement if
 972 possible; otherwise |take_scaled| will use more than 2\pct! of the running time
 973 when the Computer Modern fonts are being generated.
 974 @^inner loop@>
 975
 976 @c
 977 void mp_decimal_number_take_scaled (MP mp, mp_number *ret, mp_number p_orig, mp_number q_orig) {
 978   decNumberMultiply(ret->data.num, p_orig.data.num, q_orig.data.num, &set);
 979 }
 980
 981
 982 @ For completeness, there's also |make_scaled|, which computes a
 983 quotient as a |scaled| number instead of as a |fraction|.
 984 In other words, the result is $\lfloor2^{16}p/q+{1\over2}\rfloor$, if the
 985 operands are positive. \ (This procedure is not used especially often,
 986 so it is not part of \MP's inner loop.)
 987
 988 @c
 989 void mp_decimal_number_make_scaled (MP mp, mp_number *ret, mp_number p_orig, mp_number q_orig) {
 990   decNumberDivide(ret->data.num, p_orig.data.num, q_orig.data.num, &set);
 991   mp_check_decNumber(mp, ret->data.num, &set);
 992 }
 993
 994 @
 995 @d halfp(A) (integer)((unsigned)(A) >> 1)
 996
 997 @* Scanning numbers in the input
 998
 999 The definitions below are temporarily here
1000
1001 @d set_cur_cmd(A) mp->cur_mod_->type=(A)
1002 @d set_cur_mod(A) decNumberCopy((decNumber *)(mp->cur_mod_->data.n.data.num),&A)
1003
1004 @<Declarations...@>=
1005 static void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop);
1006
1007 @
1008 @d too_precise(a) (a == (DEC_Inexact+DEC_Rounded))
1009 @d too_large(a) (a & DEC_Overflow)
1010 @c
1011 void mp_wrapup_numeric_token(MP mp, unsigned char *start, unsigned char *stop) {
1012   decNumber result;
1013   size_t l = stop-start+1;
1014   char *buf = mp_xmalloc(mp, l+1, 1);
1015   buf[l] = '\0';
1016   (void)strncpy(buf,(const char *)start, l);
1017   set.status = 0;
1018   decNumberFromString(&result,buf, &set);
1019   free(buf);
1020   if (set.status == 0) {
1021     set_cur_mod(result);
1022   } else if (mp->scanner_status != tex_flushing) {
1023     if (too_large(set.status)) {
1024        const char *hlp[] = {"I could not handle this number specification",
1025                            "because it is out of range.",
1026                             NULL };
1027        decNumber_check (&result, &set);
1028        set_cur_mod(result);
1029        mp_error (mp, "Enormous number has been reduced", hlp, false);
1030     } else if (too_precise(set.status)) {
1031        set_cur_mod(result);
1032        if (decNumberIsPositive((decNumber *)internal_value (mp_warning_check).data.num) &&
1033           (mp->scanner_status != tex_flushing)) {
1034         char msg[256];
1035         const char *hlp[] = {"Continue and I'll round the value until it fits the current numberprecision",
1036                "(Set warningcheck:=0 to suppress this message.)",
1037                NULL };
1038         mp_snprintf (msg, 256, "Number is too precise (numberprecision = %d)", set.digits);
1039         mp_error (mp, msg, hlp, true);
1040       }
1041     } else { // this also captures underflow
1042       const char *hlp[] = {"I could not handle this number specification",
1043                            "Error:",
1044                            "",
1045                             NULL };
1046       hlp[2] = decContextStatusToString(&set);
1047       mp_error (mp, "Erroneous number specification changed to zero", hlp, false);
1048       decNumberZero(&result);
1049       set_cur_mod(result);
1050     }
1051   }
1052   set_cur_cmd((mp_variable_type)mp_numeric_token);
1053 }
1054
1055 @ @c
1056 static void find_exponent (MP mp)  {
1057   if (mp->buffer[mp->cur_input.loc_field] == 'e' ||
1058       mp->buffer[mp->cur_input.loc_field] == 'E') {
1059      mp->cur_input.loc_field++;
1060      if (!(mp->buffer[mp->cur_input.loc_field] == '+' ||
1061         mp->buffer[mp->cur_input.loc_field] == '-' ||
1062         mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class)) {
1063        mp->cur_input.loc_field--;
1064        return;
1065      }
1066      if (mp->buffer[mp->cur_input.loc_field] == '+' ||
1067         mp->buffer[mp->cur_input.loc_field] == '-') {
1068         mp->cur_input.loc_field++;
1069      }
1070      while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
1071        mp->cur_input.loc_field++;
1072      }
1073   }
1074 }
1075 void mp_decimal_scan_fractional_token (MP mp, int n) { /* n: scaled */
1076   unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
1077   unsigned char *stop;
1078   while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
1079      mp->cur_input.loc_field++;
1080   }
1081   find_exponent(mp);
1082   stop = &mp->buffer[mp->cur_input.loc_field-1];
1083   mp_wrapup_numeric_token (mp, start, stop);
1084 }
1085
1086
1087 @ We just have to collect bytes.
1088
1089 @c
1090 void mp_decimal_scan_numeric_token (MP mp, int n) { /* n: scaled */
1091   unsigned char *start = &mp->buffer[mp->cur_input.loc_field -1];
1092   unsigned char *stop;
1093   while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
1094      mp->cur_input.loc_field++;
1095   }
1096   if (mp->buffer[mp->cur_input.loc_field] == '.' &&
1097       mp->buffer[mp->cur_input.loc_field+1] != '.') {
1098      mp->cur_input.loc_field++;
1099      while (mp->char_class[mp->buffer[mp->cur_input.loc_field]] == digit_class) {
1100        mp->cur_input.loc_field++;
1101      }
1102   }
1103   find_exponent(mp);
1104   stop = &mp->buffer[mp->cur_input.loc_field-1];
1105   mp_wrapup_numeric_token (mp, start, stop);
1106 }
1107
1108 @ The |scaled| quantities in \MP\ programs are generally supposed to be
1109 less than $2^{12}$ in absolute value, so \MP\ does much of its internal
1110 arithmetic with 28~significant bits of precision. A |fraction| denotes
1111 a scaled integer whose binary point is assumed to be 28 bit positions
1112 from the right.
1113
1114 @d fraction_half (fraction_multiplier/2)
1115 @d fraction_one (1*fraction_multiplier)
1116 @d fraction_two (2*fraction_multiplier)
1117 @d fraction_three (3*fraction_multiplier)
1118 @d fraction_four (4*fraction_multiplier)
1119
1120 @ Here is a typical example of how the routines above can be used.
1121 It computes the function
1122 $${1\over3\tau}f(\theta,\phi)=
1123 {\tau^{-1}\bigl(2+\sqrt2\,(\sin\theta-{1\over16}\sin\phi)
1124  (\sin\phi-{1\over16}\sin\theta)(\cos\theta-\cos\phi)\bigr)\over
1125 3\,\bigl(1+{1\over2}(\sqrt5-1)\cos\theta+{1\over2}(3-\sqrt5\,)\cos\phi\bigr)},$$
1126 where $\tau$ is a |scaled| ``tension'' parameter. This is \MP's magic
1127 fudge factor for placing the first control point of a curve that starts
1128 at an angle $\theta$ and ends at an angle $\phi$ from the straight path.
1129 (Actually, if the stated quantity exceeds 4, \MP\ reduces it to~4.)
1130
1131 The trigonometric quantity to be multiplied by $\sqrt2$ is less than $\sqrt2$.
1132 (It's a sum of eight terms whose absolute values can be bounded using
1133 relations such as $\sin\theta\cos\theta\L{1\over2}$.) Thus the numerator
1134 is positive; and since the tension $\tau$ is constrained to be at least
1135 $3\over4$, the numerator is less than $16\over3$. The denominator is
1136 nonnegative and at most~6.
1137
1138 The angles $\theta$ and $\phi$ are given implicitly in terms of |fraction|
1139 arguments |st|, |ct|, |sf|, and |cf|, representing $\sin\theta$, $\cos\theta$,
1140 $\sin\phi$, and $\cos\phi$, respectively.
1141
1142 @c
1143 void mp_decimal_velocity (MP mp, mp_number *ret, mp_number st, mp_number ct, mp_number sf,
1144                   mp_number cf, mp_number t) {
1145   decNumber acc, num, denom;      /* registers for intermediate calculations */
1146   decNumber r1, r2;
1147   decNumber arg1, arg2;
1148   decNumber i16, fone, fhalf, ftwo, sqrtfive;
1149   decNumberFromInt32(&i16, 16);
1150   decNumberFromInt32(&fone, fraction_one);
1151   decNumberFromInt32(&fhalf, fraction_half);
1152   decNumberFromInt32(&ftwo, fraction_two);
1153   decNumberFromInt32(&sqrtfive, 5);                     // sqrt(5)
1154   decNumberSquareRoot(&sqrtfive, &sqrtfive, &set);
1155
1156
1157   decNumberDivide(&arg1,sf.data.num, &i16, &set); // arg1 = sf / 16
1158   decNumberSubtract(&arg1,st.data.num,&arg1, &set); // arg1 = st - arg1
1159   decNumberDivide(&arg2,st.data.num, &i16, &set); // arg2 = st / 16
1160   decNumberSubtract(&arg2,sf.data.num,&arg2, &set); // arg2 = sf - arg2
1161   mp_decimal_take_fraction (mp, &acc, &arg1, &arg2); // acc = (arg1 * arg2) / fmul
1162
1163   decNumberCopy(&arg1, &acc);
1164   decNumberSubtract(&arg2, ct.data.num, cf.data.num, &set); // arg2 = ct - cf
1165   mp_decimal_take_fraction (mp, &acc, &arg1, &arg2); // acc = (arg1 * arg2 ) / fmul
1166
1167   decNumberSquareRoot(&arg1, &two_decNumber, &set); // arg1 = sqrt(2)
1168   decNumberMultiply(&arg1, &arg1, &fone, &set);     // arg1 = arg1 * fmul
1169   mp_decimal_take_fraction (mp, &r1, &acc, &arg1);  // r1 = (acc * arg1) / fmul
1170   decNumberAdd(&num, &ftwo, &r1, &set);             // num = ftwo + r1
1171
1172   decNumberSubtract(&arg1,&sqrtfive, &one, &set);   // arg1 = sqrt(5) - 1
1173   decNumberMultiply(&arg1,&arg1,&fhalf, &set);      // arg1 = arg1 * fmul/2
1174   decNumberMultiply(&arg1,&arg1,&three_decNumber, &set); // arg1 = arg1 * 3
1175
1176   decNumberSubtract(&arg2,&three_decNumber, &sqrtfive, &set); // arg2 = 3 - sqrt(5)
1177   decNumberMultiply(&arg2,&arg2,&fhalf, &set);            // arg2 = arg2 * fmul/2
1178   decNumberMultiply(&arg2,&arg2,&three_decNumber, &set);  // arg2 = arg2 * 3
1179   mp_decimal_take_fraction (mp, &r1, ct.data.num, &arg1) ; // r1 = (ct * arg1) / fmul
1180   mp_decimal_take_fraction (mp, &r2, cf.data.num, &arg2);  // r2 = (cf * arg2) / fmul
1181
1182   decNumberFromInt32(&denom, fraction_three);  // denom = 3fmul
1183   decNumberAdd(&denom, &denom, &r1, &set);     // denom = denom + r1
1184   decNumberAdd(&denom, &denom, &r2, &set);     // denom = denom + r1
1185
1186   decNumberCompare(&arg1, t.data.num, &one, &set);
1187   if (!decNumberIsZero(&arg1)) {                 // t != r1
1188     decNumberDivide(&num, &num, t.data.num, &set); // num = num / t
1189   }
1190   decNumberCopy(&r2, &num);                        // r2 = num / 4
1191   decNumberDivide(&r2, &r2, &four_decNumber, &set);
1192   if (decNumberLess(&denom,&r2)) { // num/4 >= denom => denom < num/4
1193     decNumberFromInt32(ret->data.num,fraction_four);
1194   } else {
1195     mp_decimal_make_fraction (mp, ret->data.num, &num, &denom);
1196   }
1197 #if DEBUG
1198   fprintf(stdout, "\n%f = velocity(%f,%f,%f,%f,%f)", mp_number_to_double(*ret),
1199 mp_number_to_double(st),mp_number_to_double(ct),
1200 mp_number_to_double(sf),mp_number_to_double(cf),
1201 mp_number_to_double(t));
1202 #endif
1203   mp_check_decNumber(mp, ret->data.num, &set);
1204 }
1205
1206
1207 @ The following somewhat different subroutine tests rigorously if $ab$ is
1208 greater than, equal to, or less than~$cd$,
1209 given integers $(a,b,c,d)$. In most cases a quick decision is reached.
1210 The result is $+1$, 0, or~$-1$ in the three respective cases.
1211
1212 @c
1213 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) {
1214   decNumber q, r, test; /* temporary registers */
1215   decNumber a, b, c, d;
1216   (void)mp;
1217   decNumberCopy(&a, (decNumber *)a_orig.data.num);
1218   decNumberCopy(&b, (decNumber *)b_orig.data.num);
1219   decNumberCopy(&c, (decNumber *)c_orig.data.num);
1220   decNumberCopy(&d, (decNumber *)d_orig.data.num);
1221   @<Reduce to the case that |a,c>=0|, |b,d>0|@>;
1222   while (1) {
1223     decNumberDivide(&q,&a,&d, &set);
1224     decNumberDivide(&r,&c,&b, &set);
1225     decNumberCompare(&test,&q,&r, &set);
1226     if (!decNumberIsZero(&test)) {
1227       if (decNumberIsPositive(&test)) {
1228          decNumberCopy(ret->data.num, &one);
1229       } else {
1230          decNumberCopy(ret->data.num, &minusone);
1231       }
1232       goto RETURN;
1233     }
1234     decNumberRemainder(&q,&a,&d, &set);
1235     decNumberRemainder(&r,&c,&b, &set);
1236     if (decNumberIsZero(&r)) {
1237       if (decNumberIsZero(&q)) {
1238          decNumberCopy(ret->data.num, &zero);
1239       } else {
1240          decNumberCopy(ret->data.num, &one);
1241       }
1242       goto RETURN;
1243     }
1244     if (decNumberIsZero(&q)) {
1245       decNumberCopy(ret->data.num, &minusone);
1246       goto RETURN;
1247     }
1248     decNumberCopy(&a,&b);
1249     decNumberCopy(&b,&q);
1250     decNumberCopy(&c,&d);
1251     decNumberCopy(&d,&r);
1252   }                             /* now |a>d>0| and |c>b>0| */
1253 RETURN:
1254 #if DEBUG
1255   fprintf(stdout, "\n%f = ab_vs_cd(%f,%f,%f,%f)", mp_number_to_double(*ret),
1256 mp_number_to_double(a_orig),mp_number_to_double(b_orig),
1257 mp_number_to_double(c_orig),mp_number_to_double(d_orig));
1258 #endif
1259   mp_check_decNumber(mp, ret->data.num, &set);
1260   return;
1261 }
1262
1263
1264 @ @<Reduce to the case that |a...@>=
1265 if (decNumberIsNegative(&a)) {
1266   decNumberCopyNegate(&a, &a);
1267   decNumberCopyNegate(&b, &b);
1268 }
1269 if (decNumberIsNegative(&c)) {
1270   decNumberCopyNegate(&c, &c);
1271   decNumberCopyNegate(&d, &d);
1272 }
1273 if (!decNumberIsPositive(&d)) {
1274   if (!decNumberIsNegative(&b)) {
1275     if ((decNumberIsZero(&a) || decNumberIsZero(&b)) && (decNumberIsZero(&c) || decNumberIsZero(&d)))
1276          decNumberCopy(ret->data.num, &zero);
1277     else
1278          decNumberCopy(ret->data.num, &one);
1279     goto RETURN;
1280   }
1281   if (decNumberIsZero(&d)) {
1282     if (decNumberIsZero(&a))
1283          decNumberCopy(ret->data.num, &zero);
1284     else
1285          decNumberCopy(ret->data.num, &minusone);
1286     goto RETURN;
1287   }
1288   decNumberCopy(&q, &a);
1289   decNumberCopy(&a, &c);
1290   decNumberCopy(&c, &q);
1291   decNumberCopyNegate(&q, &b);
1292   decNumberCopyNegate(&b, &d);
1293   decNumberCopy(&d, &q);
1294 } else if (!decNumberIsPositive(&b)) {
1295   if (decNumberIsNegative(&b) && decNumberIsPositive(&a)) {
1296     decNumberCopy(ret->data.num, &minusone);
1297     goto RETURN;
1298   }
1299   if (decNumberIsZero(&c))
1300     decNumberCopy(ret->data.num, &zero);
1301   else
1302     decNumberCopy(ret->data.num, &minusone);
1303   goto RETURN;
1304 }
1305
1306 @ Now here's a subroutine that's handy for all sorts of path computations:
1307 Given a quadratic polynomial $B(a,b,c;t)$, the |crossing_point| function
1308 returns the unique |fraction| value |t| between 0 and~1 at which
1309 $B(a,b,c;t)$ changes from positive to negative, or returns
1310 |t=fraction_one+1| if no such value exists. If |a<0| (so that $B(a,b,c;t)$
1311 is already negative at |t=0|), |crossing_point| returns the value zero.
1312
1313 The general bisection method is quite simple when $n=2$, hence
1314 |crossing_point| does not take much time. At each stage in the
1315 recursion we have a subinterval defined by |l| and~|j| such that
1316 $B(a,b,c;2^{-l}(j+t))=B(x_0,x_1,x_2;t)$, and we want to ``zero in'' on
1317 the subinterval where $x_0\G0$ and $\min(x_1,x_2)<0$.
1318
1319 It is convenient for purposes of calculation to combine the values
1320 of |l| and~|j| in a single variable $d=2^l+j$, because the operation
1321 of bisection then corresponds simply to doubling $d$ and possibly
1322 adding~1. Furthermore it proves to be convenient to modify
1323 our previous conventions for bisection slightly, maintaining the
1324 variables $X_0=2^lx_0$, $X_1=2^l(x_0-x_1)$, and $X_2=2^l(x_1-x_2)$.
1325 With these variables the conditions $x_0\ge0$ and $\min(x_1,x_2)<0$ are
1326 equivalent to $\max(X_1,X_1+X_2)>X_0\ge0$.
1327
1328 The following code maintains the invariant relations
1329 $0\L|x0|<\max(|x1|,|x1|+|x2|)$,
1330 $\vert|x1|\vert<2^{30}$, $\vert|x2|\vert<2^{30}$;
1331 it has been constructed in such a way that no arithmetic overflow
1332 will occur if the inputs satisfy
1333 $a<2^{30}$, $\vert a-b\vert<2^{30}$, and $\vert b-c\vert<2^{30}$.
1334
1335 @d no_crossing   { decNumberCopy(ret->data.num, &fraction_one_plus_decNumber); goto RETURN; }
1336 @d one_crossing  { decNumberCopy(ret->data.num, &fraction_one_decNumber); goto RETURN; }
1337 @d zero_crossing { decNumberCopy(ret->data.num, &zero); goto RETURN; }
1338
1339 @c
1340 static void mp_decimal_crossing_point (MP mp, mp_number *ret, mp_number aa, mp_number bb, mp_number cc) {
1341   decNumber a,b,c;
1342   double d;    /* recursive counter */
1343   decNumber x, xx, x0, x1, x2;    /* temporary registers for bisection */
1344   decNumber scratch, scratch2;
1345   decNumberCopy(&a, (decNumber *)aa.data.num);
1346   decNumberCopy(&b, (decNumber *)bb.data.num);
1347   decNumberCopy(&c, (decNumber *)cc.data.num);
1348   if (decNumberIsNegative(&a))
1349     zero_crossing;
1350   if (!decNumberIsNegative(&c)) {
1351     if (!decNumberIsNegative(&b)) {
1352       if (decNumberIsPositive(&c)) {
1353         no_crossing;
1354       } else if (decNumberIsZero(&a) && decNumberIsZero(&b)) {
1355         no_crossing;
1356       } else {
1357         one_crossing;
1358       }
1359     }
1360     if (decNumberIsZero(&a))
1361       zero_crossing;
1362   } else if (decNumberIsZero(&a)) {
1363     if (!decNumberIsPositive(&b))
1364       zero_crossing;
1365   }
1366
1367   /* Use bisection to find the crossing point... */
1368   d = epsilonf;
1369   decNumberCopy(&x0, &a);
1370   decNumberSubtract(&x1,&a, &b, &set);
1371   decNumberSubtract(&x2,&b, &c, &set);
1372   /* not sure why the error correction has to be >= 1E-12 */
1373   decNumberFromDouble(&scratch2, 1E-12);
1374   do {
1375     decNumberAdd(&x, &x1, &x2, &set);
1376     decNumberDivide(&x, &x, &two_decNumber, &set);
1377     decNumberAdd(&x, &x, &scratch2, &set);
1378     decNumberSubtract(&scratch, &x1, &x0, &set);
1379     if (decNumberGreater(&scratch, &x0)) {
1380       decNumberCopy(&x2, &x);
1381       decNumberAdd(&x0, &x0, &x0, &set);
1382       d += d;
1383     } else {
1384       decNumberAdd(&xx, &scratch, &x, &set);
1385       if (decNumberGreater(&xx,&x0)) {
1386         decNumberCopy(&x2,&x);
1387         decNumberAdd(&x0, &x0, &x0, &set);
1388         d += d;
1389       } else {
1390         decNumberSubtract(&x0, &x0, &xx, &set);
1391         if (!decNumberGreater(&x,&x0)) {
1392           decNumberAdd(&scratch, &x, &x2, &set);
1393           if (!decNumberGreater(&scratch, &x0))
1394             no_crossing;
1395         }
1396         decNumberCopy(&x1,&x);
1397         d = d + d + epsilonf;
1398       }
1399     }
1400   } while (d < fraction_one);
1401   decNumberFromDouble(&scratch, d);
1402   decNumberSubtract(ret->data.num,&scratch, &fraction_one_decNumber, &set);
1403 RETURN:
1404 #if DEBUG
1405   fprintf(stdout, "\n%f = crossing_point(%f,%f,%f)", mp_number_to_double(*ret),
1406 mp_number_to_double(aa),mp_number_to_double(bb),mp_number_to_double(cc));
1407 #endif
1408   mp_check_decNumber(mp, ret->data.num, &set);
1409   return;
1410 }
1411
1412
1413 @ We conclude this set of elementary routines with some simple rounding
1414 and truncation operations.
1415
1416
1417 @ |round_unscaled| rounds a |scaled| and converts it to |int|
1418 @c
1419 int mp_round_unscaled(mp_number x_orig) {
1420   double xx = mp_number_to_double(x_orig);
1421   int x = (int)ROUND(xx);
1422   return x;
1423 }
1424
1425 @ |number_floor| floors a number
1426
1427 @c
1428 void mp_number_floor (mp_number *i) {
1429   int round = set.round;
1430   set.round = DEC_ROUND_FLOOR;
1431   decNumberToIntegralValue(i->data.num, i->data.num, &set);
1432   set.round = round;
1433 }
1434
1435 @ |fraction_to_scaled| rounds a |fraction| and converts it to |scaled|
1436 @c
1437 void mp_decimal_fraction_to_round_scaled (mp_number *x_orig) {
1438   x_orig->type = mp_scaled_type;
1439   decNumberDivide(x_orig->data.num, x_orig->data.num, &fraction_multiplier_decNumber, &set);
1440 }
1441
1442
1443
1444 @* Algebraic and transcendental functions.
1445 \MP\ computes all of the necessary special functions from scratch, without
1446 relying on |real| arithmetic or system subroutines for sines, cosines, etc.
1447
1448 @
1449
1450 @c
1451 void mp_decimal_square_rt (MP mp, mp_number *ret, mp_number x_orig) { /* return, x: scaled */
1452   decNumber x;
1453   decNumberCopy(&x, x_orig.data.num);
1454   if (!decNumberIsPositive(&x)) {
1455     @<Handle square root of zero or negative argument@>;
1456   } else {
1457     decNumberSquareRoot(ret->data.num, &x, &set);
1458   }
1459   mp_check_decNumber(mp, ret->data.num, &set);
1460 }
1461
1462
1463 @ @<Handle square root of zero...@>=
1464 {
1465   if (decNumberIsNegative(&x)) {
1466     char msg[256];
1467     const char *hlp[] = {
1468            "Since I don't take square roots of negative numbers,",
1469            "I'm zeroing this one. Proceed, with fingers crossed.",
1470            NULL };
1471     char *xstr = mp_decimal_number_tostring (mp, x_orig);
1472     mp_snprintf(msg, 256, "Square root of %s has been replaced by 0", xstr);
1473     free(xstr);
1474 @.Square root...replaced by 0@>;
1475     mp_error (mp, msg, hlp, true);
1476   }
1477   decNumberZero(ret->data.num);
1478   return;
1479 }
1480
1481
1482 @ Pythagorean addition $\psqrt{a^2+b^2}$ is implemented by a quick hack
1483
1484 @c
1485 void mp_decimal_pyth_add (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig) {
1486   decNumber a, b;
1487   decNumber asq, bsq;
1488   decNumberCopyAbs(&a, a_orig.data.num);
1489   decNumberCopyAbs(&b, b_orig.data.num);
1490   decNumberMultiply(&asq, &a, &a, &set);
1491   decNumberMultiply(&bsq, &b, &b, &set);
1492   decNumberAdd(&a, &asq, &bsq, &set);
1493   decNumberSquareRoot(ret->data.num, &a, &set);
1494   //if (set.status != 0) {
1495   //  mp->arith_error = true;
1496   //  decNumberCopy(ret->data.num, &EL_GORDO_decNumber);
1497   //}
1498   mp_check_decNumber(mp, ret->data.num, &set);
1499 }
1500
1501 @ Here is a similar algorithm for $\psqrt{a^2-b^2}$. Same quick hack, also.
1502
1503 @c
1504 void mp_decimal_pyth_sub (MP mp, mp_number *ret, mp_number a_orig, mp_number b_orig) {
1505   decNumber a, b;
1506   decNumberCopyAbs(&a, a_orig.data.num);
1507   decNumberCopyAbs(&b, b_orig.data.num);
1508   if (!decNumberGreater(&a,&b)) {
1509     @<Handle erroneous |pyth_sub| and set |a:=0|@>;
1510   } else {
1511     decNumber asq, bsq;
1512     decNumberMultiply(&asq, &a, &a, &set);
1513     decNumberMultiply(&bsq, &b, &b, &set);
1514     decNumberSubtract(&a, &asq, &bsq, &set);
1515     decNumberSquareRoot(&a, &a, &set);
1516   }
1517   decNumberCopy(ret->data.num, &a);
1518   mp_check_decNumber(mp, ret->data.num, &set);
1519 }
1520
1521
1522 @ @<Handle erroneous |pyth_sub| and set |a:=0|@>=
1523 {
1524   if (decNumberLess(&a, &b)) {
1525     char msg[256];
1526     const char *hlp[] = {
1527          "Since I don't take square roots of negative numbers,",
1528          "I'm zeroing this one. Proceed, with fingers crossed.",
1529          NULL };
1530     char *astr = mp_decimal_number_tostring (mp, a_orig);
1531     char *bstr = mp_decimal_number_tostring (mp, b_orig);
1532     mp_snprintf (msg, 256, "Pythagorean subtraction %s+-+%s has been replaced by 0", astr, bstr);
1533     free(astr);
1534     free(bstr);
1535 @.Pythagorean...@>;
1536     mp_error (mp, msg, hlp, true);
1537   }
1538   decNumberZero(&a);
1539 }
1540
1541
1542 @ Here is the routine that calculates $2^8$ times the natural logarithm
1543 of a |scaled| quantity;
1544
1545 @c
1546 void mp_decimal_m_log (MP mp, mp_number *ret, mp_number x_orig) {
1547   if (!decNumberIsPositive((decNumber *)x_orig.data.num)) {
1548     @<Handle non-positive logarithm@>;
1549   } else {
1550     decNumber twofivesix;
1551     decNumberFromInt32(&twofivesix, 256);
1552     decNumberLn(ret->data.num, x_orig.data.num, &limitedset);
1553     mp_check_decNumber(mp, ret->data.num, &limitedset);
1554     decNumberMultiply(ret->data.num, ret->data.num, &twofivesix, &set);
1555   }
1556   mp_check_decNumber(mp, ret->data.num, &set);
1557 }
1558
1559 @ @<Handle non-positive logarithm@>=
1560 {
1561   char msg[256];
1562   const char *hlp[] = {
1563          "Since I don't take logs of non-positive numbers,",
1564          "I'm zeroing this one. Proceed, with fingers crossed.",
1565           NULL };
1566   char *xstr = mp_decimal_number_tostring (mp, x_orig);
1567   mp_snprintf (msg, 256, "Logarithm of %s has been replaced by 0", xstr);
1568   free (xstr);
1569 @.Logarithm...replaced by 0@>;
1570   mp_error (mp, msg, hlp, true);
1571   decNumberZero(ret->data.num);
1572 }
1573
1574
1575 @ Conversely, the exponential routine calculates $\exp(x/2^8)$,
1576 when |x| is |scaled|.
1577
1578 @c
1579 void mp_decimal_m_exp (MP mp, mp_number *ret, mp_number x_orig) {
1580   decNumber temp, twofivesix;
1581   decNumberFromInt32(&twofivesix, 256);
1582   decNumberDivide(&temp, x_orig.data.num, &twofivesix, &set);
1583   limitedset.status = 0;
1584   decNumberExp(ret->data.num, &temp, &limitedset);
1585   if (limitedset.status & DEC_Clamped) {
1586     if (decNumberIsPositive((decNumber *)x_orig.data.num)) {
1587       mp->arith_error = true;
1588       decNumberCopy(ret->data.num, &EL_GORDO_decNumber);
1589     } else {
1590       decNumberZero(ret->data.num);
1591     }
1592   }
1593   mp_check_decNumber(mp, ret->data.num, &limitedset);
1594   limitedset.status = 0;
1595 }
1596
1597
1598 @ Given integers |x| and |y|, not both zero, the |n_arg| function
1599 returns the |angle| whose tangent points in the direction $(x,y)$.
1600
1601 @c
1602 void mp_decimal_n_arg (MP mp, mp_number *ret, mp_number x_orig, mp_number y_orig) {
1603   if (decNumberIsZero((decNumber *)x_orig.data.num) && decNumberIsZero((decNumber *)y_orig.data.num)) {
1604     @<Handle undefined arg@>;
1605   } else {
1606     decNumber atan2val, oneeighty_angle;
1607     ret->type = mp_angle_type;
1608     decNumberFromInt32(&oneeighty_angle, 180 * angle_multiplier);
1609     decNumberDivide(&oneeighty_angle, &oneeighty_angle, &PI_decNumber, &set);
1610     checkZero(y_orig.data.num);
1611     checkZero(x_orig.data.num);
1612     decNumberAtan2(&atan2val, y_orig.data.num, x_orig.data.num, &set);
1613 #if DEBUG
1614     fprintf(stdout, "\n%g = atan2(%g,%g)", decNumberToDouble(&atan2val),mp_number_to_double(x_orig),mp_number_to_double(y_orig));
1615 #endif
1616     decNumberMultiply(ret->data.num,&atan2val, &oneeighty_angle, &set);
1617     checkZero(ret->data.num);
1618 #if DEBUG
1619     fprintf(stdout, "\nn_arg(%g,%g,%g)", mp_number_to_double(*ret),
1620     mp_number_to_double(x_orig),mp_number_to_double(y_orig));
1621 #endif
1622   }
1623   mp_check_decNumber(mp, ret->data.num, &set);
1624 }
1625
1626
1627 @ @<Handle undefined arg@>=
1628 {
1629   const char *hlp[] = {
1630          "The `angle' between two identical points is undefined.",
1631          "I'm zeroing this one. Proceed, with fingers crossed.",
1632          NULL };
1633   mp_error (mp, "angle(0,0) is taken as zero", hlp, true);
1634 @.angle(0,0)...zero@>;
1635   decNumberZero(ret->data.num);
1636 }
1637
1638
1639 @ Conversely, the |n_sin_cos| routine takes an |angle| and produces the sine
1640 and cosine of that angle. The results of this routine are
1641 stored in global integer variables |n_sin| and |n_cos|.
1642
1643 First, we need a decNumber function that calculates sines and cosines
1644 using the Taylor series. This function is fairly optimized.
1645
1646 @d PRECALC_FACTORIALS_CACHESIZE 50
1647
1648 @c
1649 static void sinecosine(decNumber *theangle, decNumber *c, decNumber *s)
1650 {
1651     int n, i, prec;
1652     decNumber p, pxa, fac, cc;
1653     decNumber n1, n2, p1;
1654     decNumberZero(c);
1655     decNumberZero(s);
1656     prec = (set.digits/2);
1657     if (prec < DECPRECISION_DEFAULT) prec = DECPRECISION_DEFAULT;
1658     for (n=0;n<prec;n++)
1659     {
1660         decNumberFromInt32(&p1, n);
1661         decNumberFromInt32(&n1, 2*n);
1662         decNumberPower(&p,  &minusone, &p1, &limitedset);
1663         if (n==0) {
1664           decNumberCopy(&pxa, &one);
1665         } else {
1666           decNumberPower(&pxa, theangle, &n1, &limitedset);
1667         }
1668
1669         if (2*n<last_cached_factorial) {
1670           decNumberCopy(&fac,factorials[2*n]);
1671         } else {
1672           decNumberCopy(&fac,factorials[last_cached_factorial]);
1673           for (i = last_cached_factorial+1; i <= 2*n; i++) {
1674             decNumberFromInt32(&cc, i);
1675             decNumberMultiply (&fac, &fac, &cc, &set);
1676             if (i<PRECALC_FACTORIALS_CACHESIZE) {
1677               factorials[i] = malloc(sizeof(decNumber));
1678               decNumberCopy(factorials[i],&fac);
1679               last_cached_factorial = i;
1680             }
1681           }
1682         }
1683
1684         decNumberDivide   (&pxa, &pxa, &fac, &set);
1685         decNumberMultiply (&pxa, &pxa, &p,   &set);
1686         decNumberAdd      (s,    s,    &pxa, &set);
1687
1688         decNumberFromInt32(&n2, 2*n+1);
1689         decNumberMultiply (&fac, &fac, &n2,  &set); // fac = fac * (2*n+1)
1690         decNumberPower(&pxa, theangle, &n2,  &limitedset);
1691         decNumberDivide   (&pxa, &pxa, &fac, &set);
1692         decNumberMultiply (&pxa, &pxa, &p,   &set);
1693         decNumberAdd      (c,    c,    &pxa, &set);
1694         // printf("\niteration %2d: %-42s %-42s",n,tostring(c), tostring(s));
1695     }
1696 }
1697
1698 @ Calculate sines and cosines.
1699 @c
1700 void mp_decimal_sin_cos (MP mp, mp_number z_orig, mp_number *n_cos, mp_number *n_sin) {
1701   decNumber rad;
1702   decNumber one_eighty;
1703   decNumberFromInt32(&one_eighty, 180 * 16);
1704 #if DEBUG
1705   fprintf(stdout, "\nsin_cos(%f)", mp_number_to_double(z_orig));
1706 #endif
1707   decNumberMultiply(&rad, z_orig.data.num, &PI_decNumber, &set);
1708   decNumberDivide(&rad, &rad, &one_eighty, &set);
1709 #if 0
1710   if (decNumberIsNegative(&rad)) {
1711     while (decNumberLess(&rad,&PI_decNumber))
1712       decNumberAdd(&rad, &rad, &PI_decNumber, &set);
1713   } else {
1714     while (decNumberGreater(&rad,&PI_decNumber))
1715       decNumberSubtract(&rad, &rad, &PI_decNumber, &set);
1716   }
1717 #endif
1718   sinecosine(&rad, n_sin->data.num, n_cos->data.num);
1719   decNumberMultiply(n_cos->data.num,n_cos->data.num,&fraction_multiplier_decNumber, &set);
1720   decNumberMultiply(n_sin->data.num,n_sin->data.num,&fraction_multiplier_decNumber, &set);
1721 #if DEBUG
1722   fprintf(stdout, "\nsin_cos(%f,%f,%f)", decNumberToDouble(&rad),
1723 mp_number_to_double(*n_cos), mp_number_to_double(*n_sin));
1724 #endif
1725    mp_check_decNumber(mp, n_cos->data.num, &set);
1726    mp_check_decNumber(mp, n_sin->data.num, &set);
1727 }
1728
1729 @ This is the http://www-cs-faculty.stanford.edu/~uno/programs/rng.c
1730 with  small cosmetic modifications.
1731
1732 @c
1733 #define KK 100                     /* the long lag  */
1734 #define LL  37                     /* the short lag */
1735 #define MM (1L<<30)                /* the modulus   */
1736 #define mod_diff(x,y) (((x)-(y))&(MM-1)) /* subtraction mod MM */
1737 /* */
1738 static long ran_x[KK];                    /* the generator state */
1739 /* */
1740 static void ran_array(long aa[],int n) /* put n new random numbers in aa */
1741   /* long aa[]    destination */
1742   /* int n       array length (must be at least KK) */
1743 {
1744   register int i,j;
1745   for (j=0;j<KK;j++) aa[j]=ran_x[j];
1746   for (;j<n;j++) aa[j]=mod_diff(aa[j-KK],aa[j-LL]);
1747   for (i=0;i<LL;i++,j++) ran_x[i]=mod_diff(aa[j-KK],aa[j-LL]);
1748   for (;i<KK;i++,j++) ran_x[i]=mod_diff(aa[j-KK],ran_x[i-LL]);
1749 }
1750 /* */
1751 /* the following routines are from exercise 3.6--15 */
1752 /* after calling ran_start, get new randoms by, e.g., "x=ran_arr_next()" */
1753 /* */
1754 #define QUALITY 1009 /* recommended quality level for high-res use */
1755 static long ran_arr_buf[QUALITY];
1756 static long ran_arr_dummy=-1, ran_arr_started=-1;
1757 static long *ran_arr_ptr=&ran_arr_dummy; /* the next random number, or -1 */
1758 /* */
1759 #define TT  70   /* guaranteed separation between streams */
1760 #define is_odd(x)  ((x)&1)          /* units bit of x */
1761 /* */
1762 static void ran_start(long seed) /* do this before using ran_array */
1763   /* long seed             selector for different streams */
1764 {
1765   register int t,j;
1766   long x[KK+KK-1];              /* the preparation buffer */
1767   register long ss=(seed+2)&(MM-2);
1768   for (j=0;j<KK;j++) {
1769     x[j]=ss;                      /* bootstrap the buffer */
1770     ss<<=1; if (ss>=MM) ss-=MM-2; /* cyclic shift 29 bits */
1771   }
1772   x[1]++;              /* make x[1] (and only x[1]) odd */
1773   for (ss=seed&(MM-1),t=TT-1; t; ) {
1774     for (j=KK-1;j>0;j--) x[j+j]=x[j], x[j+j-1]=0; /* "square" */
1775     for (j=KK+KK-2;j>=KK;j--)
1776       x[j-(KK-LL)]=mod_diff(x[j-(KK-LL)],x[j]),
1777       x[j-KK]=mod_diff(x[j-KK],x[j]);
1778     if (is_odd(ss)) {              /* "multiply by z" */
1779       for (j=KK;j>0;j--)  x[j]=x[j-1];
1780       x[0]=x[KK];            /* shift the buffer cyclically */
1781       x[LL]=mod_diff(x[LL],x[KK]);
1782     }
1783     if (ss) ss>>=1; else t--;
1784   }
1785   for (j=0;j<LL;j++) ran_x[j+KK-LL]=x[j];
1786   for (;j<KK;j++) ran_x[j-LL]=x[j];
1787   for (j=0;j<10;j++) ran_array(x,KK+KK-1); /* warm things up */
1788   ran_arr_ptr=&ran_arr_started;
1789 }
1790 /* */
1791 #define ran_arr_next() (*ran_arr_ptr>=0? *ran_arr_ptr++: ran_arr_cycle())
1792 static long ran_arr_cycle(void)
1793 {
1794   if (ran_arr_ptr==&ran_arr_dummy)
1795     ran_start(314159L); /* the user forgot to initialize */
1796   ran_array(ran_arr_buf,QUALITY);
1797   ran_arr_buf[KK]=-1;
1798   ran_arr_ptr=ran_arr_buf+1;
1799   return ran_arr_buf[0];
1800 }
1801
1802
1803
1804 @ To initialize the |randoms| table, we call the following routine.
1805
1806 @c
1807 void mp_init_randoms (MP mp, int seed) {
1808   int j, jj, k;    /* more or less random integers */
1809   int i;        /* index into |randoms| */
1810   j =  abs (seed);
1811   while (j >= fraction_one) {
1812     j = j/2;
1813   }
1814   k = 1;
1815   for (i = 0; i <= 54; i++) {
1816     jj = k;
1817     k = j - k;
1818     j = jj;
1819     if (k<0)
1820       k += fraction_one;
1821     decNumberFromInt32(mp->randoms[(i * 21) % 55].data.num, j);
1822   }
1823   mp_new_randoms (mp);
1824   mp_new_randoms (mp);
1825   mp_new_randoms (mp);          /* ``warm up'' the array */
1826
1827   ran_start((unsigned long) seed);
1828
1829 }
1830
1831 @ @c
1832 void mp_decimal_number_modulo (mp_number *a, mp_number b) {
1833    decNumberRemainder(a->data.num, a->data.num, b.data.num, &set);
1834 }
1835
1836
1837 @ To consume a random  integer for the uniform generator, the program below will say `|next_unif_random|'.
1838
1839 @c
1840 static void mp_next_unif_random (MP mp, mp_number *ret) {
1841   decNumber a;
1842   decNumber b;
1843   unsigned long int op;
1844   (void)mp;
1845   op = (unsigned)ran_arr_next();
1846   decNumberFromInt32(&a, op);
1847   decNumberFromInt32(&b, MM);
1848   decNumberDivide (&a, &a, &b, &set); /* a = a/b */
1849   decNumberCopy(ret->data.num, &a);
1850   mp_check_decNumber(mp, ret->data.num, &set);
1851 }
1852
1853
1854 @ To consume a random fraction, the program below will say `|next_random|'.
1855
1856 @c
1857 static void mp_next_random (MP mp, mp_number *ret) {
1858   if ( mp->j_random==0 )
1859     mp_new_randoms(mp);
1860   else
1861     mp->j_random = mp->j_random-1;
1862   mp_number_clone (ret, mp->randoms[mp->j_random]);
1863 }
1864
1865
1866 @ To produce a uniform random number in the range |0<=u<x| or |0>=u>x|
1867 or |0=u=x|, given a |scaled| value~|x|, we proceed as shown here.
1868
1869 Note that the call of |take_fraction| will produce the values 0 and~|x|
1870 with about half the probability that it will produce any other particular
1871 values between 0 and~|x|, because it rounds its answers.
1872
1873 @c
1874 static void mp_decimal_m_unif_rand (MP mp, mp_number *ret, mp_number x_orig) {
1875   mp_number y;     /* trial value */
1876   mp_number x, abs_x;
1877   mp_number u;
1878   new_fraction (y);
1879   new_number (x);
1880   new_number (abs_x);
1881   new_number (u);
1882   mp_number_clone (&x, x_orig);
1883   mp_number_clone (&abs_x, x);
1884   mp_decimal_abs (&abs_x);
1885   mp_next_unif_random(mp, &u);
1886   decNumberMultiply (y.data.num, abs_x.data.num, u.data.num, &set);
1887   free_number (u);
1888   if (mp_number_equal(y, abs_x)) {
1889     mp_number_clone (ret, ((math_data *)mp->math)->zero_t);
1890   } else if (mp_number_greater(x, ((math_data *)mp->math)->zero_t)) {
1891     mp_number_clone (ret, y);
1892   } else {
1893     mp_number_clone (ret, y);
1894     mp_number_negate (ret);
1895   }
1896   free_number (abs_x);
1897   free_number (x);
1898   free_number (y);
1899 }
1900
1901
1902
1903 @ Finally, a normal deviate with mean zero and unit standard deviation
1904 can readily be obtained with the ratio method (Algorithm 3.4.1R in
1905 {\sl The Art of Computer Programming\/}).
1906
1907 @c
1908 static void mp_decimal_m_norm_rand (MP mp, mp_number *ret) {
1909   mp_number ab_vs_cd;
1910   mp_number abs_x;
1911   mp_number u;
1912   mp_number r;
1913   mp_number la, xa;
1914   new_number (ab_vs_cd);
1915   new_number (la);
1916   new_number (xa);
1917   new_number (abs_x);
1918   new_number (u);
1919   new_number (r);
1920
1921   do {
1922     do {
1923       mp_number v;
1924       new_number (v);
1925       mp_next_random(mp, &v);
1926       mp_number_substract (&v, ((math_data *)mp->math)->fraction_half_t);
1927       mp_decimal_number_take_fraction (mp,&xa, ((math_data *)mp->math)->sqrt_8_e_k, v);
1928       free_number (v);
1929       mp_next_random(mp, &u);
1930       mp_number_clone (&abs_x, xa);
1931       mp_decimal_abs (&abs_x);
1932     } while (!mp_number_less(abs_x, u));
1933     mp_decimal_number_make_fraction (mp, &r, xa, u);
1934     mp_number_clone (&xa, r);
1935     mp_decimal_m_log (mp,&la, u);
1936     mp_set_decimal_from_substraction(&la, ((math_data *)mp->math)->twelve_ln_2_k, la);
1937     mp_ab_vs_cd (mp,&ab_vs_cd, ((math_data *)mp->math)->one_k, la, xa, xa);
1938   } while (mp_number_less(ab_vs_cd,((math_data *)mp->math)->zero_t));
1939   mp_number_clone (ret, xa);
1940   free_number (ab_vs_cd);
1941   free_number (r);
1942   free_number (abs_x);
1943   free_number (la);
1944   free_number (xa);
1945   free_number (u);
1946 }
1947
1948
1949
1950
1951 @ The following subroutine could be used  in norm_rand and tests  if $ab$ is
1952 greater than, equal to, or less than~$cd$.
1953 The result is $+1$, 0, or~$-1$ in the three respective cases.
1954 This is not necessary, even if it's shorter than the current ab_vs_cd
1955 and looks as a native implememtation.
1956
1957 @c
1958 /*
1959 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) {
1960   decNumber a, b, c, d;
1961   decNumber ab, cd;
1962   (void)mp;
1963
1964   decNumberCopy(&a, (decNumber *)a_orig.data.num);
1965   decNumberCopy(&b, (decNumber *)b_orig.data.num);
1966   decNumberCopy(&c, (decNumber *)c_orig.data.num);
1967   decNumberCopy(&d, (decNumber *)d_orig.data.num);
1968
1969
1970   decNumberMultiply (&ab, (decNumber *)a_orig.data.num, (decNumber *)b_orig.data.num, &set);
1971   decNumberMultiply (&cd, (decNumber *)c_orig.data.num, (decNumber *)d_orig.data.num, &set);
1972   decNumberCompare(ret->data.num, &ab, &cd, &set);
1973   mp_check_decNumber(mp, ret->data.num, &set);
1974   return;
1975
1976 }
1977 */
1978
1979
1980