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