gcalctool/mp.c

   1
   2 /*  $Header$
   3  *
   4  *  Copyright (c) 1987-2006 Sun Microsystems, Inc. All Rights Reserved.
   5  *
   6  *  This program is free software; you can redistribute it and/or modify
   7  *  it under the terms of the GNU General Public License as published by
   8  *  the Free Software Foundation; either version 2, or (at your option)
   9  *  any later version.
  10  *
  11  *  This program is distributed in the hope that it will be useful, but
  12  *  WITHOUT ANY WARRANTY; without even the implied warranty of
  13  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  14  *  General Public License for more details.
  15  *
  16  *  You should have received a copy of the GNU General Public License
  17  *  along with this program; if not, write to the Free Software
  18  *  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  19  *  02111-1307, USA.
  20  */
  21
  22 /*  This maths library is based on the MP multi-precision floating-point
  23  *  arithmetic package originally written in FORTRAN by Richard Brent,
  24  *  Computer Centre, Australian National University in the 1970's.
  25  *
  26  *  It has been converted from FORTRAN into C using the freely available
  27  *  f2c translator, available via netlib on research.att.com.
  28  *
  29  *  The subsequently converted C code has then been tidied up, mainly to
  30  *  remove any dependencies on the libI77 and libF77 support libraries.
  31  *
  32  *  FOR A GENERAL DESCRIPTION OF THE PHILOSOPHY AND DESIGN OF MP,
  33  *  SEE - R. P. BRENT, A FORTRAN MULTIPLE-PRECISION ARITHMETIC
  34  *  PACKAGE, ACM TRANS. MATH. SOFTWARE 4 (MARCH 1978), 57-70.
  35  *  SOME ADDITIONAL DETAILS ARE GIVEN IN THE SAME ISSUE, 71-81.
  36  *  FOR DETAILS OF THE IMPLEMENTATION, CALLING SEQUENCES ETC. SEE
  37  *  THE MP USERS GUIDE.
  38  */
  39
  40 #include <stdio.h>
  41 #include "calctool.h"
  42 #include "extern.h"
  43
  44 static struct {
  45     int b, t, m, mxr, r[MP_SIZE];
  46 } MP;
  47
  48 /* Table of constant values */
  49
  50 static int c__0   = 0;
  51 static int c__1   = 1;
  52 static int c__4   = 4;
  53 static int c__2   = 2;
  54 static int c__6   = 6;
  55 static int c__5   = 5;
  56 static int c__12  = 12;
  57 static int c_n2   = -2;
  58 static int c__10  = 10;
  59 static int c__32  = 32;
  60 static int c__3   = 3;
  61 static int c__8   = 8;
  62 static int c__14  = 14;
  63 static int c_n1   = -1;
  64 static int c__239 = 239;
  65 static int c__7   = 7;
  66 static int c__16  = 16;
  67
  68 static double mppow_di(double *, int *);
  69 static double mppow_ri(float *, int *);
  70
  71 static int mpcmpi(int *, int *);
  72 static int mpcmpr(int *, float *);
  73 static int mpcomp(int *, int *);
  74 static int pow_ii(int *, int *);
  75
  76 static void mpadd2(int *, int *, int *, int *, int *);
  77 static void mpadd3(int *, int *, int *, int *, int *);
  78 static void mpaddq(int *, int *, int *, int *);
  79 static void mpart1(int *, int *);
  80 static void mpchk(int *, int *);
  81 static void mpcmr(int *, float *);
  82 static void mpcqm(int *, int *, int *);
  83 static void mpcrm(float *, int *);
  84 static void mpexp1(int *, int *);
  85 static void mpext(int *, int *, int *);
  86 static void mpgcd(int *, int *);
  87 static void mplns(int *, int *);
  88 static void mpmaxr(int *);
  89 static void mpmlp(int *, int *, int *, int *);
  90 static void mpmul2(int *, int *, int *, int *);
  91 static void mpmulq(int *, int *, int *, int *);
  92 static void mpnzr(int *, int *, int *, int *);
  93 static void mpovfl(int *);
  94 static void mprec(int *, int *);
  95 static void mproot(int *, int *, int *);
  96 static void mpsin1(int *, int *, int *);
  97 static void mpunfl(int *);
  98
  99
 100 void
 101 mpabs(int *x, int *y)
 102 {
 103
 104 /* SETS Y = ABS(X) FOR MP NUMBERS X AND Y */
 105
 106     --y;                      /* Parameter adjustments */
 107     --x;
 108
 109     mpstr(&x[1], &y[1]);
 110     y[1] = C_abs(y[1]);
 111 }
 112
 113
 114 void
 115 mpadd(int *x, int *y, int *z)
 116 {
 117
 118 /*  ADDS X AND Y, FORMING RESULT IN Z, WHERE X, Y AND Z ARE MP
 119  *  NUMBERS.  FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING.
 120  */
 121
 122     --z;                    /* Parameter adjustments */
 123     --y;
 124     --x;
 125
 126     mpadd2(&x[1], &y[1], &z[1], &y[1], &c__0);
 127 }
 128
 129
 130 static void
 131 mpadd2(int *x, int *y, int *z, int *y1, int *trunc)
 132 {
 133     int i__1, i__2;
 134
 135     static int j, s;
 136     static int ed, re, rs, med;
 137
 138 /*  CALLED BY MPADD, MPSUB ETC.
 139  *  X, Y AND Z ARE MP NUMBERS, Y1 AND TRUNC ARE INTEGERS.
 140  *  TO FORCE CALL BY REFERENCE RATHER THAN VALUE/RESULT, Y1 IS
 141  *  DECLARED AS AN ARRAY, BUT ONLY Y1(1) IS EVER USED.
 142  *  SETS Z = X + Y1(1)*ABS(Y), WHERE Y1(1) = +- Y(1).
 143  *  IF TRUNC.EQ.0 R*-ROUNDING IS USED, OTHERWISE TRUNCATION.
 144  *  R*-ROUNDING IS DEFINED IN KUKI AND CODI, COMM. ACM
 145  *  16(1973), 223.  (SEE ALSO BRENT, IEEE TC-22(1973), 601.)
 146  *  CHECK FOR X OR Y ZERO
 147  */
 148
 149     --y1;                      /* Parameter adjustments */
 150     --z;
 151     --y;
 152     --x;
 153
 154     if (x[1] != 0) goto L20;
 155
 156 /* X = 0 OR NEGLIGIBLE, SO RESULT = +-Y */
 157
 158 L10:
 159     mpstr(&y[1], &z[1]);
 160     z[1] = y1[1];
 161     return;
 162
 163 L20:
 164     if (y1[1] != 0) goto L40;
 165
 166 /* Y = 0 OR NEGLIGIBLE, SO RESULT = X */
 167
 168 L30:
 169     mpstr(&x[1], &z[1]);
 170     return;
 171
 172 /* COMPARE SIGNS */
 173
 174 L40:
 175     s = x[1] * y1[1];
 176     if (C_abs(s) <= 1) goto L60;
 177
 178     mpchk(&c__1, &c__4);
 179     if (v->MPerrors) {
 180         FPRINTF(stderr, "*** SIGN NOT 0, +1 OR -1 IN MPADD2 CALL.\nPOSSIBLE OVERWRITING PROBLEM ***\n");
 181     }
 182
 183     mperr();
 184     z[1] = 0;
 185     return;
 186
 187 /* COMPARE EXPONENTS */
 188
 189 L60:
 190     ed = x[2] - y[2];
 191     med = C_abs(ed);
 192     if (ed < 0)  goto L90;
 193     else if (ed == 0) goto L70;
 194     else goto L120;
 195
 196 /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
 197
 198 L70:
 199     if (s > 0) goto L100;
 200
 201     i__1 = MP.t;
 202     for (j = 1; j <= i__1; ++j) {
 203         if ((i__2 = x[j + 2] - y[j + 2]) < 0) goto L100;
 204         else if (i__2 == 0) continue;
 205         else goto L130;
 206     }
 207
 208 /* RESULT IS ZERO */
 209
 210     z[1] = 0;
 211     return;
 212
 213 /* HERE EXPONENT(Y) .GE. EXPONENT(X) */
 214
 215 L90:
 216     if (med > MP.t) goto L10;
 217
 218 L100:
 219     rs = y1[1];
 220     re = y[2];
 221     mpadd3(&x[1], &y[1], &s, &med, &re);
 222
 223 /* NORMALIZE, ROUND OR TRUNCATE, AND RETURN */
 224
 225 L110:
 226     mpnzr(&rs, &re, &z[1], trunc);
 227     return;
 228
 229 /* ABS(X) .GT. ABS(Y) */
 230
 231 L120:
 232     if (med > MP.t) goto L30;
 233
 234 L130:
 235     rs = x[1];
 236     re = x[2];
 237     mpadd3(&y[1], &x[1], &s, &med, &re);
 238     goto L110;
 239 }
 240
 241
 242 static void
 243 mpadd3(int *x, int *y, int *s, int *med, int *re)
 244 {
 245     int i__1;
 246
 247     static int c, i, j, i2, i2p, ted;
 248
 249 /* CALLED BY MPADD2, DOES INNER LOOPS OF ADDITION */
 250
 251     --y;                 /* Parameter adjustments */
 252     --x;
 253
 254     ted = MP.t + *med;
 255     i2 = MP.t + 4;
 256     i = i2;
 257     c = 0;
 258
 259 /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
 260
 261 L10:
 262     if (i <= ted) goto L20;
 263
 264     MP.r[i - 1] = 0;
 265     --i;
 266     goto L10;
 267
 268 L20:
 269     if (*s < 0) goto L130;
 270
 271 /* HERE DO ADDITION, EXPONENT(Y) .GE. EXPONENT(X) */
 272
 273     if (i <= MP.t) goto L40;
 274
 275 L30:
 276     j = i - *med;
 277     MP.r[i - 1] = x[j + 2];
 278     --i;
 279     if (i > MP.t) goto L30;
 280
 281 L40:
 282     if (i <= *med) goto L60;
 283
 284     j = i - *med;
 285     c = y[i + 2] + x[j + 2] + c;
 286     if (c < MP.b) goto L50;
 287
 288 /* CARRY GENERATED HERE */
 289
 290     MP.r[i - 1] = c - MP.b;
 291     c = 1;
 292     --i;
 293     goto L40;
 294
 295 /* NO CARRY GENERATED HERE */
 296
 297 L50:
 298     MP.r[i - 1] = c;
 299     c = 0;
 300     --i;
 301     goto L40;
 302
 303 L60:
 304     if (i <= 0) goto L90;
 305
 306     c = y[i + 2] + c;
 307     if (c < MP.b) goto L70;
 308
 309     MP.r[i - 1] = 0;
 310     c = 1;
 311     --i;
 312     goto L60;
 313
 314 L70:
 315     MP.r[i - 1] = c;
 316     --i;
 317
 318 /* NO CARRY POSSIBLE HERE */
 319
 320 L80:
 321     if (i <= 0) return;
 322
 323     MP.r[i - 1] = y[i + 2];
 324     --i;
 325     goto L80;
 326
 327 L90:
 328     if (c == 0) return;
 329
 330 /* MUST SHIFT RIGHT HERE AS CARRY OFF END */
 331
 332     i2p = i2 + 1;
 333     i__1 = i2;
 334     for (j = 2; j <= i__1; ++j) {
 335         i = i2p - j;
 336         MP.r[i] = MP.r[i - 1];
 337     }
 338     MP.r[0] = 1;
 339     ++(*re);
 340     return;
 341
 342 /* HERE DO SUBTRACTION, ABS(Y) .GT. ABS(X) */
 343
 344 L110:
 345     j = i - *med;
 346     MP.r[i - 1] = c - x[j + 2];
 347     c = 0;
 348     if (MP.r[i - 1] >= 0) goto L120;
 349
 350 /* BORROW GENERATED HERE */
 351
 352     c = -1;
 353     MP.r[i - 1] += MP.b;
 354
 355 L120:
 356     --i;
 357
 358 L130:
 359     if (i > MP.t) goto L110;
 360
 361 L140:
 362     if (i <= *med) goto L160;
 363
 364     j = i - *med;
 365     c = y[i + 2] + c - x[j + 2];
 366     if (c >= 0) goto L150;
 367
 368 /* BORROW GENERATED HERE */
 369
 370     MP.r[i - 1] = c + MP.b;
 371     c = -1;
 372     --i;
 373     goto L140;
 374
 375 /* NO BORROW GENERATED HERE */
 376
 377 L150:
 378     MP.r[i - 1] = c;
 379     c = 0;
 380     --i;
 381     goto L140;
 382
 383 L160:
 384     if (i <= 0) return;
 385
 386     c = y[i + 2] + c;
 387     if (c >= 0) goto L70;
 388
 389     MP.r[i - 1] = c + MP.b;
 390     c = -1;
 391     --i;
 392     goto L160;
 393 }
 394
 395
 396 void
 397 mpaddi(int *x, int *iy, int *z)
 398 {
 399
 400 /*  ADDS MULTIPLE-PRECISION X TO INTEGER IY
 401  *  GIVING MULTIPLE-PRECISION Z.
 402  *  DIMENSION OF R IN CALLING PROGRAM MUST BE
 403  *  AT LEAST 2T+6 (BUT Z(1) MAY BE R(T+5)).
 404  *  DIMENSION R(6) BECAUSE RALPH COMPILER ON UNIVAC 1100 COMPUTERS
 405  *  OBJECTS TO DIMENSION R(1).
 406  *  CHECK LEGALITY OF B, T, M AND MXR
 407  */
 408
 409     --z;               /* Parameter adjustments */
 410     --x;
 411
 412     mpchk(&c__2, &c__6);
 413     mpcim(iy, &MP.r[MP.t + 4]);
 414     mpadd(&x[1], &MP.r[MP.t + 4], &z[1]);
 415 }
 416
 417
 418 static void
 419 mpaddq(int *x, int *i, int *j, int *y)
 420 {
 421
 422 /*  ADDS THE RATIONAL NUMBER I/J TO MP NUMBER X, MP RESULT IN Y
 423  *  DIMENSION OF R MUST BE AT LEAST 2T+6
 424  *  CHECK LEGALITY OF B, T, M AND MXR
 425  */
 426
 427     --y;               /* Parameter adjustments */
 428     --x;
 429
 430     mpchk(&c__2, &c__6);
 431     mpcqm(i, j, &MP.r[MP.t + 4]);
 432     mpadd(&x[1], &MP.r[MP.t + 4], &y[1]);
 433 }
 434
 435
 436 static void
 437 mpart1(int *n, int *y)
 438 {
 439     int i__1, i__2;
 440
 441     static int i, b2, i2, id, ts;
 442
 443 /*  COMPUTES MP Y = ARCTAN(1/N), ASSUMING INTEGER N .GT. 1.
 444  *  USES SERIES ARCTAN(X) = X - X**3/3 + X**5/5 - ...
 445  *  DIMENSION OF R IN CALLING PROGRAM MUST BE
 446  *  AT LEAST 2T+6
 447  *  CHECK LEGALITY OF B, T, M AND MXR
 448  */
 449
 450     --y;                  /* Parameter adjustments */
 451
 452     mpchk(&c__2, &c__6);
 453     if (*n > 1) goto L20;
 454
 455     if (v->MPerrors) {
 456         FPRINTF(stderr, "*** N .LE. 1 IN CALL TO MPART1 ***\n");
 457     }
 458
 459     mperr();
 460     y[1] = 0;
 461     return;
 462
 463 L20:
 464     i2 = MP.t + 5;
 465     ts = MP.t;
 466
 467 /* SET SUM TO X = 1/N */
 468
 469     mpcqm(&c__1, n, &y[1]);
 470
 471 /* SET ADDITIVE TERM TO X */
 472
 473     mpstr(&y[1], &MP.r[i2 - 1]);
 474     i = 1;
 475     id = 0;
 476
 477 /* ASSUME AT LEAST 16-BIT WORD. */
 478
 479     b2 = max(MP.b, 64);
 480     if (*n < b2) id = b2 * 7 * b2 / (*n * *n);
 481
 482 /* MAIN LOOP.  FIRST REDUCE T IF POSSIBLE */
 483
 484 L30:
 485     MP.t = ts + 2 + MP.r[i2] - y[2];
 486     if (MP.t < 2) goto L60;
 487
 488     MP.t = min(MP.t,ts);
 489
 490 /*  IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
 491  *  FOLLOWING HAS TO BE PERFORMED IN SEVERAL STEPS.
 492  */
 493
 494     if (i >= id) goto L40;
 495
 496     i__1 = -i;
 497     i__2 = (i + 2) * *n * *n;
 498     mpmulq(&MP.r[i2 - 1], &i__1, &i__2, &MP.r[i2 - 1]);
 499     goto L50;
 500
 501 L40:
 502     i__1 = -i;
 503     i__2 = i + 2;
 504     mpmulq(&MP.r[i2 - 1], &i__1, &i__2, &MP.r[i2 - 1]);
 505     mpdivi(&MP.r[i2 - 1], n, &MP.r[i2 - 1]);
 506     mpdivi(&MP.r[i2 - 1], n, &MP.r[i2 - 1]);
 507
 508 L50:
 509     i += 2;
 510
 511 /* RESTORE T */
 512
 513     MP.t = ts;
 514
 515 /* ADD TO SUM, USING MPADD2 (FASTER THAN MPADD) */
 516
 517     mpadd2(&MP.r[i2 - 1], &y[1], &y[1], &y[1], &c__0);
 518     if (MP.r[i2 - 1] != 0) goto L30;
 519
 520 L60:
 521     MP.t = ts;
 522 }
 523
 524
 525 void
 526 mpasin(int *x, int *y)
 527 {
 528     int i__1;
 529
 530     static int i2, i3;
 531
 532 /*  RETURNS Y = ARCSIN(X), ASSUMING ABS(X) .LE. 1,
 533  *  FOR MP NUMBERS X AND Y.
 534  *  Y IS IN THE RANGE -PI/2 TO +PI/2.
 535  *  METHOD IS TO USE MPATAN, SO TIME IS O(M(T)T).
 536  *  DIMENSION OF R MUST BE AT LEAST 5T+12
 537  *  CHECK LEGALITY OF B, T, M AND MXR
 538  */
 539
 540     --y;                 /* Parameter adjustments */
 541     --x;
 542
 543     mpchk(&c__5, &c__12);
 544     i3 = (MP.t << 2) + 11;
 545     if (x[1] == 0) goto L30;
 546
 547     if (x[2] <= 0) goto L40;
 548
 549 /* HERE ABS(X) .GE. 1.  SEE IF X = +-1 */
 550
 551     mpcim(&x[1], &MP.r[i3 - 1]);
 552     if (mpcomp(&x[1], &MP.r[i3 - 1]) != 0) goto L10;
 553
 554 /* X = +-1 SO RETURN +-PI/2 */
 555
 556     mppi(&y[1]);
 557     i__1 = MP.r[i3 - 1] << 1;
 558     mpdivi(&y[1], &i__1, &y[1]);
 559     return;
 560
 561 L10:
 562     if (v->MPerrors) {
 563         FPRINTF(stderr, "*** ABS(X) .GT. 1 IN CALL TO MPASIN ***\n");
 564     }
 565
 566     mperr();
 567
 568 L30:
 569     y[1] = 0;
 570     return;
 571
 572 /* HERE ABS(X) .LT. 1 SO USE ARCTAN(X/SQRT(1 - X**2)) */
 573
 574 L40:
 575     i2 = i3 - (MP.t + 2);
 576     mpcim(&c__1, &MP.r[i2 - 1]);
 577     mpstr(&MP.r[i2 - 1], &MP.r[i3 - 1]);
 578     mpsub(&MP.r[i2 - 1], &x[1], &MP.r[i2 - 1]);
 579     mpadd(&MP.r[i3 - 1], &x[1], &MP.r[i3 - 1]);
 580     mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
 581     mproot(&MP.r[i3 - 1], &c_n2, &MP.r[i3 - 1]);
 582     mpmul(&x[1], &MP.r[i3 - 1], &y[1]);
 583     mpatan(&y[1], &y[1]);
 584 }
 585
 586
 587 void
 588 mpatan(int *x, int *y)
 589 {
 590     int i__1, i__2;
 591     float r__1;
 592
 593     static int i, q, i2, i3, ie, ts;
 594     static float rx, ry;
 595
 596 /*  RETURNS Y = ARCTAN(X) FOR MP X AND Y, USING AN O(T.M(T)) METHOD
 597  *  WHICH COULD EASILY BE MODIFIED TO AN O(SQRT(T)M(T))
 598  *  METHOD (AS IN MPEXP1). Y IS IN THE RANGE -PI/2 TO +PI/2.
 599  *  FOR AN ASYMPTOTICALLY FASTER METHOD, SEE - FAST MULTIPLE-
 600  *  PRECISION EVALUATION OF ELEMENTARY FUNCTIONS
 601  *  (BY R. P. BRENT), J. ACM 23 (1976), 242-251,
 602  *  AND THE COMMENTS IN MPPIGL.
 603  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 5T+12
 604  *  CHECK LEGALITY OF B, T, M AND MXR
 605  */
 606
 607     --y;                 /* Parameter adjustments */
 608     --x;
 609
 610     mpchk(&c__5, &c__12);
 611     i2 = MP.t * 3 + 9;
 612     i3 = i2 + MP.t + 2;
 613     if (x[1] != 0) {
 614         goto L10;
 615     }
 616     y[1] = 0;
 617     return;
 618
 619 L10:
 620     mpstr(&x[1], &MP.r[i3 - 1]);
 621     ie = C_abs(x[2]);
 622     if (ie <= 2) mpcmr(&x[1], &rx);
 623
 624     q = 1;
 625
 626 /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
 627
 628 L20:
 629     if (MP.r[i3] < 0) goto L30;
 630
 631     if (MP.r[i3] == 0 && (MP.r[i3 + 1] + 1) << 1 <= MP.b)
 632     goto L30;
 633
 634     q <<= 1;
 635     mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], &y[1]);
 636     mpaddi(&y[1], &c__1, &y[1]);
 637     mpsqrt(&y[1], &y[1]);
 638     mpaddi(&y[1], &c__1, &y[1]);
 639     mpdiv(&MP.r[i3 - 1], &y[1], &MP.r[i3 - 1]);
 640     goto L20;
 641
 642 /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
 643
 644 L30:
 645     mpstr(&MP.r[i3 - 1], &y[1]);
 646     mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]);
 647     i = 1;
 648     ts = MP.t;
 649
 650 /* SERIES LOOP.  REDUCE T IF POSSIBLE. */
 651
 652 L40:
 653     MP.t = ts + 2 + MP.r[i3];
 654     if (MP.t <= 2) goto L50;
 655
 656     MP.t = min(MP.t,ts);
 657     mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i3 - 1]);
 658     i__1 = -i;
 659     i__2 = i + 2;
 660     mpmulq(&MP.r[i3 - 1], &i__1, &i__2, &MP.r[i3 - 1]);
 661     i += 2;
 662     MP.t = ts;
 663     mpadd(&y[1], &MP.r[i3 - 1], &y[1]);
 664     if (MP.r[i3 - 1] != 0) goto L40;
 665
 666 /* RESTORE T, CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
 667
 668 L50:
 669     MP.t = ts;
 670     mpmuli(&y[1], &q, &y[1]);
 671
 672 /*  CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
 673  *  OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
 674  */
 675
 676     if (ie > 2) return;
 677
 678     mpcmr(&y[1], &ry);
 679     if ((r__1 = ry - atan(rx), dabs(r__1)) < dabs(ry) * (float).01)
 680         return;
 681
 682 /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
 683
 684     if (v->MPerrors) {
 685         FPRINTF(stderr, "*** ERROR OCCURRED IN MPATAN, RESULT INCORRECT ***\n");
 686     }
 687
 688     mperr();
 689 }
 690
 691
 692 void
 693 mpcdm(double *dx, int *z)
 694 {
 695     int i__1;
 696
 697     static int i, k, i2, ib, ie, re, tp, rs;
 698     static double db, dj;
 699
 700 /*  CONVERTS DOUBLE-PRECISION NUMBER DX TO MULTIPLE-PRECISION Z.
 701  *  SOME NUMBERS WILL NOT CONVERT EXACTLY ON MACHINES
 702  *  WITH BASE OTHER THAN TWO, FOUR OR SIXTEEN.
 703  *  THIS ROUTINE IS NOT CALLED BY ANY OTHER ROUTINE IN MP,
 704  *  SO MAY BE OMITTED IF DOUBLE-PRECISION IS NOT AVAILABLE.
 705  *  CHECK LEGALITY OF B, T, M AND MXR
 706  */
 707
 708     --z;              /* Parameter adjustments */
 709
 710     mpchk(&c__1, &c__4);
 711     i2 = MP.t + 4;
 712
 713 /* CHECK SIGN */
 714
 715     if (*dx < 0.) goto L20;
 716     else if (*dx == 0) goto L10;
 717     else goto L30;
 718
 719 /* IF DX = 0D0 RETURN 0 */
 720
 721 L10:
 722     z[1] = 0;
 723     return;
 724
 725 /* DX .LT. 0D0 */
 726
 727 L20:
 728     rs = -1;
 729     dj = -(*dx);
 730     goto L40;
 731
 732 /* DX .GT. 0D0 */
 733
 734 L30:
 735     rs = 1;
 736     dj = *dx;
 737
 738 L40:
 739     ie = 0;
 740
 741 L50:
 742     if (dj < 1.) goto L60;
 743
 744 /* INCREASE IE AND DIVIDE DJ BY 16. */
 745
 746     ++ie;
 747     dj *= .0625;
 748     goto L50;
 749
 750 L60:
 751     if (dj >= .0625) goto L70;
 752
 753     --ie;
 754     dj *= 16.;
 755     goto L60;
 756
 757 /*  NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
 758  *  SET EXPONENT TO 0
 759  */
 760
 761 L70:
 762     re = 0;
 763
 764 /* DB = DFLOAT(B) IS NOT ANSI STANDARD SO USE FLOAT AND DBLE */
 765
 766     db = (double) ((float) MP.b);
 767
 768 /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
 769
 770     i__1 = i2;
 771     for (i = 1; i <= i__1; ++i) {
 772         dj = db * dj;
 773         MP.r[i - 1] = (int) dj;
 774         dj -= (double) ((float) MP.r[i - 1]);
 775     }
 776
 777 /* NORMALIZE RESULT */
 778
 779     mpnzr(&rs, &re, &z[1], &c__0);
 780
 781 /* Computing MAX */
 782
 783     i__1 = MP.b * 7 * MP.b;
 784     ib = max(i__1, 32767) / 16;
 785     tp = 1;
 786
 787 /* NOW MULTIPLY BY 16**IE */
 788
 789     if (ie < 0)  goto L90;
 790     else if (ie == 0) goto L130;
 791     else goto L110;
 792
 793 L90:
 794     k = -ie;
 795     i__1 = k;
 796     for (i = 1; i <= i__1; ++i) {
 797         tp <<= 4;
 798         if (tp <= ib && tp != MP.b && i < k) continue;
 799         mpdivi(&z[1], &tp, &z[1]);
 800         tp = 1;
 801     }
 802     return;
 803
 804 L110:
 805     i__1 = ie;
 806     for (i = 1; i <= i__1; ++i) {
 807         tp <<= 4;
 808         if (tp <= ib && tp != MP.b && i < ie) continue;
 809         mpmuli(&z[1], &tp, &z[1]);
 810         tp = 1;
 811     }
 812
 813 L130:
 814     return;
 815 }
 816
 817
 818 static void
 819 mpchk(int *i, int *j)
 820 {
 821     static int ib, mx;
 822
 823 /*  CHECKS LEGALITY OF B, T, M AND MXR.
 824  *  THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT
 825  *  MXR .GE. (I*T + J)
 826  */
 827
 828 /* CHECK LEGALITY OF B, T AND M */
 829
 830     if (MP.b > 1) goto L40;
 831
 832     if (v->MPerrors) {
 833         FPRINTF(stderr, "*** B = %d ILLEGAL IN CALL TO MPCHK.\nPERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***\n", MP.b);
 834     }
 835     mperr();
 836
 837 L40:
 838     if (MP.t > 1) goto L60;
 839
 840     if (v->MPerrors) {
 841         FPRINTF(stderr, "*** T = %d ILLEGAL IN CALL TO MPCHK.\nPERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***\n", MP.t);
 842     }
 843     mperr();
 844
 845 L60:
 846     if (MP.m > MP.t) goto L80;
 847
 848     if (v->MPerrors) {
 849         FPRINTF(stderr, "*** M .LE. T IN CALL TO MPCHK.\nPERHAPS NOT SET BEFORE CALL TO AN MP ROUTINE ***\n");
 850     }
 851     mperr();
 852
 853 /*  8*B*B-1 SHOULD BE REPRESENTABLE, IF NOT WILL OVERFLOW
 854  *  AND MAY BECOME NEGATIVE, SO CHECK FOR THIS
 855  */
 856
 857 L80:
 858     ib = (MP.b << 2) * MP.b - 1;
 859     if (ib > 0 && (ib << 1) + 1 > 0) goto L100;
 860
 861     if (v->MPerrors) {
 862         FPRINTF(stderr, "*** B TOO LARGE IN CALL TO MPCHK ***\n");
 863     }
 864
 865     mperr();
 866
 867 /* CHECK THAT SPACE IN COMMON IS SUFFICIENT */
 868
 869 L100:
 870     mx = *i * MP.t + *j;
 871     if (MP.mxr >= mx) return;
 872
 873 /* HERE COMMON IS TOO SMALL, SO GIVE ERROR MESSAGE. */
 874
 875     if (v->MPerrors) {
 876         FPRINTF(stderr,
 877           "*** MXR TOO SMALL OR NOT SET TO DIM(R) BEFORE CALL TO AN MP ROUTINE ***\n");
 878         FPRINTF(stderr,
 879           "*** MXR SHOULD BE AT LEAST %d*T + %d = %d  ***\n*** ACTUALLY MXR = %d, AND T = %d  ***\n",
 880           *i, *j, mx, MP.mxr, MP.t);
 881     }
 882
 883     mperr();
 884 }
 885
 886
 887 void
 888 mpcim(int *ix, int *z)
 889 {
 890     int i__1;
 891
 892     static int i, n;
 893
 894 /*  CONVERTS INTEGER IX TO MULTIPLE-PRECISION Z.
 895  *  CHECK LEGALITY OF B, T, M AND MXR
 896  */
 897
 898     --z;            /* Parameter adjustments */
 899
 900     mpchk(&c__1, &c__4);
 901     n = *ix;
 902     if (n < 0)  goto L20;
 903     else if (n == 0) goto L10;
 904     else goto L30;
 905
 906 L10:
 907     z[1] = 0;
 908     return;
 909
 910 L20:
 911     n = -n;
 912     z[1] = -1;
 913     goto L40;
 914
 915 L30:
 916     z[1] = 1;
 917
 918 /* SET EXPONENT TO T */
 919
 920 L40:
 921     z[2] = MP.t;
 922
 923 /* CLEAR FRACTION */
 924
 925     i__1 = MP.t;
 926     for (i = 2; i <= i__1; ++i) z[i + 1] = 0;
 927
 928 /* INSERT N */
 929
 930     z[MP.t + 2] = n;
 931
 932 /* NORMALIZE BY CALLING MPMUL2 */
 933
 934     mpmul2(&z[1], &c__1, &z[1], &c__1);
 935 }
 936
 937
 938 void
 939 mpcmd(int *x, double *dz)
 940 {
 941     int i__1;
 942     double d__1;
 943
 944     static int i, tm;
 945     static double db, dz2;
 946
 947 /*  CONVERTS MULTIPLE-PRECISION X TO DOUBLE-PRECISION DZ.
 948  *  ASSUMES X IS IN ALLOWABLE RANGE FOR DOUBLE-PRECISION
 949  *  NUMBERS.   THERE IS SOME LOSS OF ACCURACY IF THE
 950  *  EXPONENT IS LARGE.
 951  *  CHECK LEGALITY OF B, T, M AND MXR
 952  */
 953
 954     --x;         /* Parameter adjustments */
 955
 956     mpchk(&c__1, &c__4);
 957     *dz = 0.;
 958     if (x[1] == 0) return;
 959
 960 /* DB = DFLOAT(B) IS NOT ANSI STANDARD, SO USE FLOAT AND DBLE */
 961
 962     db = (double) ((float) MP.b);
 963     i__1 = MP.t;
 964     for (i = 1; i <= i__1; ++i) {
 965         *dz = db * *dz + (double) ((float) x[i + 2]);
 966         tm = i;
 967
 968 /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
 969
 970         dz2 = *dz + 1.;
 971
 972 /*  TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20,
 973  *  FOR EXAMPLE ON CYBER 76.
 974  */
 975         if (dz2 - *dz <= 0.) goto L20;
 976     }
 977
 978 /* NOW ALLOW FOR EXPONENT */
 979
 980 L20:
 981     i__1 = x[2] - tm;
 982     *dz *= mppow_di(&db, &i__1);
 983
 984 /* CHECK REASONABLENESS OF RESULT. */
 985
 986     if (*dz <= 0.) goto L30;
 987
 988 /* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */
 989
 990     if ((d__1 = (double) ((float) x[2]) - (log(*dz) / log((double)
 991                 ((float) MP.b)) + .5), C_abs(d__1)) > .6) {
 992         goto L30;
 993     }
 994
 995     if (x[1] < 0) *dz = -(*dz);
 996     return;
 997
 998 /*  FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
 999  *  TRY USING MPCMDE INSTEAD.
1000  */
1001
1002 L30:
1003     if (v->MPerrors) {
1004         FPRINTF(stderr, "*** FLOATING-POINT OVER/UNDER-FLOW IN MPCMD ***\n");
1005     }
1006
1007     mperr();
1008 }
1009
1010
1011 void
1012 mpcmf(int *x, int *y)
1013 {
1014     int i__1;
1015
1016     static int i;
1017     static int x2, il, ip, xs;
1018
1019 /*  FOR MP X AND Y, RETURNS FRACTIONAL PART OF X IN Y,
1020  *  I.E., Y = X - INTEGER PART OF X (TRUNCATED TOWARDS 0).
1021  */
1022
1023     --y;                   /* Parameter adjustments */
1024     --x;
1025
1026     if (x[1] != 0) goto L20;
1027
1028 /* RETURN 0 IF X = 0 */
1029
1030 L10:
1031     y[1] = 0;
1032     return;
1033
1034 L20:
1035     x2 = x[2];
1036
1037 /* RETURN 0 IF EXPONENT SO LARGE THAT NO FRACTIONAL PART */
1038
1039     if (x2 >= MP.t) goto L10;
1040
1041 /* IF EXPONENT NOT POSITIVE CAN RETURN X */
1042
1043     if (x2 > 0) goto L30;
1044
1045     mpstr(&x[1], &y[1]);
1046     return;
1047
1048 /* CLEAR ACCUMULATOR */
1049
1050 L30:
1051     i__1 = x2;
1052     for (i = 1; i <= i__1; ++i) MP.r[i - 1] = 0;
1053
1054     il = x2 + 1;
1055
1056 /* MOVE FRACTIONAL PART OF X TO ACCUMULATOR */
1057
1058     i__1 = MP.t;
1059     for (i = il; i <= i__1; ++i) MP.r[i - 1] = x[i + 2];
1060
1061     for (i = 1; i <= 4; ++i) {
1062         ip = i + MP.t;
1063         MP.r[ip - 1] = 0;
1064     }
1065     xs = x[1];
1066
1067 /* NORMALIZE RESULT AND RETURN */
1068
1069     mpnzr(&xs, &x2, &y[1], &c__1);
1070 }
1071
1072
1073 void
1074 mpcmi(int *x, int *iz)
1075 {
1076     int i__1;
1077
1078     static int i, j, k, j1, x2, kx, xs, izs;
1079
1080 /*  CONVERTS MULTIPLE-PRECISION X TO INTEGER IZ,
1081  *  ASSUMING THAT X NOT TOO LARGE (ELSE USE MPCMIM).
1082  *  X IS TRUNCATED TOWARDS ZERO.
1083  *  IF INT(X)IS TOO LARGE TO BE REPRESENTED AS A SINGLE-
1084  *  PRECISION INTEGER, IZ IS RETURNED AS ZERO.  THE USER
1085  *  MAY CHECK FOR THIS POSSIBILITY BY TESTING IF
1086  *  ((X(1).NE.0).AND.(X(2).GT.0).AND.(IZ.EQ.0)) IS TRUE ON
1087  *  RETURN FROM MPCMI.
1088  */
1089
1090     --x;             /* Parameter adjustments */
1091
1092     xs = x[1];
1093     *iz = 0;
1094     if (xs == 0) return;
1095
1096     if (x[2] <= 0) return;
1097
1098     x2 = x[2];
1099     i__1 = x2;
1100     for (i = 1; i <= i__1; ++i) {
1101         izs = *iz;
1102         *iz = MP.b * *iz;
1103         if (i <= MP.t) *iz += x[i + 2];
1104
1105 /* CHECK FOR SIGNS OF INTEGER OVERFLOW */
1106
1107         if (*iz <= 0 || *iz <= izs) goto L30;
1108     }
1109
1110 /*  CHECK THAT RESULT IS CORRECT (AN UNDETECTED OVERFLOW MAY
1111  *  HAVE OCCURRED).
1112  */
1113
1114     j = *iz;
1115     i__1 = x2;
1116     for (i = 1; i <= i__1; ++i) {
1117         j1 = j / MP.b;
1118         k = x2 + 1 - i;
1119         kx = 0;
1120         if (k <= MP.t) kx = x[k + 2];
1121         if (kx != j - MP.b * j1) goto L30;
1122         j = j1;
1123     }
1124     if (j != 0) goto L30;
1125
1126 /* RESULT CORRECT SO RESTORE SIGN AND RETURN */
1127
1128     *iz = xs * *iz;
1129     return;
1130
1131 /*  HERE OVERFLOW OCCURRED (OR X WAS UNNORMALIZED), SO
1132  *  RETURN ZERO.
1133  */
1134
1135 L30:
1136     *iz = 0;
1137 }
1138
1139
1140 void
1141 mpcmim(int *x, int *y)
1142 {
1143     int tmp[MP_SIZE];     /* Temporary store for the number. */
1144     int accuracy;         /* Temporary story for the accuracy. */
1145     int pointed;          /* Whether a decimal point has been given. */
1146     int toclear;          /* Indicates if display should be cleared. */
1147     char disp[MAXLINE];   /* Setup a string to store what would be displayed */
1148     enum num_type dtype;  /* Setup a temp display type variable */
1149
1150     int i__1;
1151
1152     static int i, il, ll;
1153
1154 /* RETURNS Y = INTEGER PART OF X (TRUNCATED TOWARDS 0), FOR MP X AND Y.
1155  * USE IF Y TOO LARGE TO BE REPRESENTABLE AS A SINGLE-PRECISION INTEGER.
1156  * (ELSE COULD USE MPCMI).
1157  * CHECK LEGALITY OF B, T, M AND MXR
1158  */
1159
1160     --y;                 /* Parameter adjustments */
1161     --x;
1162
1163     mpchk(&c__1, &c__4);
1164     mpstr(&x[1], &y[1]);
1165     if (y[1] == 0) {
1166         return;
1167     }
1168
1169     il = y[2] + 1;
1170     ll = il;
1171
1172 /* IF EXPONENT LARGE ENOUGH RETURN Y = X */
1173
1174     if (il > MP.t) {
1175         return;
1176     }
1177
1178 /* IF EXPONENT SMALL ENOUGH RETURN ZERO */
1179
1180     if (il > 1) {
1181         goto L10;
1182     }
1183
1184     y[1] = 0;
1185     return;
1186
1187 /* SET FRACTION TO ZERO */
1188
1189 L10:
1190     i__1 = MP.t;
1191     for (i = il; i <= i__1; ++i) {
1192         y[i + 2] = 0;
1193     }
1194
1195 /* Fix for Sun bugtraq bug #4006391 - problem with Int function for numbers
1196  * like 4800 when internal representation in something like 4799.999999999.
1197  */
1198
1199     accuracy = v->accuracy;
1200     dtype = v->dtype;
1201     pointed = v->pointed;
1202     toclear = v->toclear;
1203
1204     v->dtype = FIX;
1205     v->accuracy = MAX_DIGITS;
1206     mpcmf(&x[1], tmp);
1207     STRCPY(disp, make_number(tmp, v->base, FALSE));
1208
1209     if (disp[0] == '1') {
1210         y[ll+1]++;
1211     }
1212
1213     v->accuracy = accuracy;
1214     v->dtype = dtype;
1215     v->pointed = pointed;
1216     v->toclear = toclear;
1217 }
1218
1219
1220 static int
1221 mpcmpi(int *x, int *i)
1222 {
1223     int ret_val;
1224
1225 /*  COMPARES MP NUMBER X WITH INTEGER I, RETURNING
1226  *      +1 IF X .GT. I,
1227  *       0 IF X .EQ. I,
1228  *      -1 IF X .LT. I
1229  *  DIMENSION OF R IN COMMON AT LEAST 2T+6
1230  *  CHECK LEGALITY OF B, T, M AND MXR
1231  */
1232
1233     --x;              /* Parameter adjustments */
1234
1235     mpchk(&c__2, &c__6);
1236
1237 /* CONVERT I TO MULTIPLE-PRECISION AND COMPARE */
1238
1239     mpcim(i, &MP.r[MP.t + 4]);
1240     ret_val = mpcomp(&x[1], &MP.r[MP.t + 4]);
1241     return(ret_val);
1242 }
1243
1244
1245 static int
1246 mpcmpr(int *x, float *ri)
1247 {
1248     int ret_val;
1249
1250 /*  COMPARES MP NUMBER X WITH REAL NUMBER RI, RETURNING
1251  *      +1 IF X .GT. RI,
1252  *       0 IF X .EQ. RI,
1253  *      -1 IF X .LT. RI
1254  *  DIMENSION OF R IN COMMON AT LEAST 2T+6
1255  *  CHECK LEGALITY OF B, T, M AND MXR
1256  */
1257
1258     --x;              /* Parameter adjustments */
1259
1260     mpchk(&c__2, &c__6);
1261
1262 /* CONVERT RI TO MULTIPLE-PRECISION AND COMPARE */
1263
1264     mpcrm(ri, &MP.r[MP.t + 4]);
1265     ret_val = mpcomp(&x[1], &MP.r[MP.t + 4]);
1266     return(ret_val);
1267 }
1268
1269
1270 static void
1271 mpcmr(int *x, float *rz)
1272 {
1273     int i__1;
1274     float r__1;
1275
1276     static int i, tm;
1277     static float rb, rz2;
1278
1279 /*  CONVERTS MULTIPLE-PRECISION X TO SINGLE-PRECISION RZ.
1280  *  ASSUMES X IN ALLOWABLE RANGE.  THERE IS SOME LOSS OF
1281  *  ACCURACY IF THE EXPONENT IS LARGE.
1282  *  CHECK LEGALITY OF B, T, M AND MXR
1283  */
1284
1285     --x;               /* Parameter adjustments */
1286
1287     mpchk(&c__1, &c__4);
1288     *rz = (float) 0.0;
1289     if (x[1] == 0) return;
1290
1291     rb = (float) MP.b;
1292     i__1 = MP.t;
1293     for (i = 1; i <= i__1; ++i) {
1294         *rz = rb * *rz + (float) x[i + 2];
1295         tm = i;
1296
1297 /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
1298
1299         rz2 = *rz + (float) 1.0;
1300         if (rz2 <= *rz) goto L20;
1301     }
1302
1303 /* NOW ALLOW FOR EXPONENT */
1304
1305 L20:
1306     i__1 = x[2] - tm;
1307     *rz *= mppow_ri(&rb, &i__1);
1308
1309 /* CHECK REASONABLENESS OF RESULT */
1310
1311     if (*rz <= (float)0.) goto L30;
1312
1313 /* LHS SHOULD BE .LE. 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
1314
1315     if ((r__1 = (float) x[2] - (log(*rz) / log((float) MP.b) + (float).5),
1316                  dabs(r__1)) > (float).6) {
1317         goto L30;
1318     }
1319
1320     if (x[1] < 0) *rz = -(double)(*rz);
1321     return;
1322
1323 /*  FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
1324  *  TRY USING MPCMRE INSTEAD.
1325  */
1326
1327 L30:
1328     if (v->MPerrors) {
1329         FPRINTF(stderr, "*** FLOATING-POINT OVER/UNDER-FLOW IN MPCMR ***\n");
1330     }
1331
1332     mperr();
1333 }
1334
1335
1336 static int
1337 mpcomp(int *x, int *y)
1338 {
1339     int ret_val, i__1, i__2;
1340
1341     static int i, t2;
1342
1343 /*  COMPARES THE MULTIPLE-PRECISION NUMBERS X AND Y,
1344  *  RETURNING +1 IF X .GT. Y,
1345  *            -1 IF X .LT. Y,
1346  *  AND        0 IF X .EQ. Y.
1347  */
1348
1349     --y;             /* Parameter adjustments */
1350     --x;
1351
1352     if ((i__1 = x[1] - y[1]) < 0) goto L10;
1353     else if (i__1 == 0) goto L30;
1354     else goto L20;
1355
1356 /* X .LT. Y */
1357
1358 L10:
1359     ret_val = -1;
1360     return(ret_val);
1361
1362 /* X .GT. Y */
1363
1364 L20:
1365     ret_val = 1;
1366     return(ret_val);
1367
1368 /* SIGN(X) = SIGN(Y), SEE IF ZERO */
1369
1370 L30:
1371     if (x[1] != 0) goto L40;
1372
1373 /* X = Y = 0 */
1374
1375     ret_val = 0;
1376     return(ret_val);
1377
1378 /* HAVE TO COMPARE EXPONENTS AND FRACTIONS */
1379
1380 L40:
1381     t2 = MP.t + 2;
1382     i__1 = t2;
1383     for (i = 2; i <= i__1; ++i) {
1384         if ((i__2 = x[i] - y[i]) < 0) goto L60;
1385         else if (i__2 == 0) continue;
1386         else goto L70;
1387     }
1388
1389 /* NUMBERS EQUAL */
1390
1391     ret_val = 0;
1392     return(ret_val);
1393
1394 /* ABS(X) .LT. ABS(Y) */
1395
1396 L60:
1397     ret_val = -x[1];
1398     return(ret_val);
1399
1400 /* ABS(X) .GT. ABS(Y) */
1401
1402 L70:
1403     ret_val = x[1];
1404     return(ret_val);
1405 }
1406
1407
1408 void
1409 mpcos(int *x, int *y)
1410 {
1411     static int i2;
1412
1413 /*  RETURNS Y = COS(X) FOR MP X AND Y, USING MPSIN AND MPSIN1.
1414  *  DIMENSION OF R IN COMMON AT LEAST 5T+12.
1415  */
1416
1417     --y;              /* Parameter adjustments */
1418     --x;
1419
1420     if (x[1] != 0) goto L10;
1421
1422 /* COS(0) = 1 */
1423
1424     mpcim(&c__1, &y[1]);
1425     return;
1426
1427 /* CHECK LEGALITY OF B, T, M AND MXR */
1428
1429 L10:
1430     mpchk(&c__5, &c__12);
1431     i2 = MP.t * 3 + 12;
1432
1433 /* SEE IF ABS(X) .LE. 1 */
1434
1435     mpabs(&x[1], &y[1]);
1436     if (mpcmpi(&y[1], &c__1) <= 0) goto L20;
1437
1438 /*  HERE ABS(X) .GT. 1 SO USE COS(X) = SIN(PI/2 - ABS(X)),
1439  *  COMPUTING PI/2 WITH ONE GUARD DIGIT.
1440  */
1441
1442     ++MP.t;
1443     mppi(&MP.r[i2 - 1]);
1444     mpdivi(&MP.r[i2 - 1], &c__2, &MP.r[i2 - 1]);
1445     --MP.t;
1446     mpsub(&MP.r[i2 - 1], &y[1], &y[1]);
1447     mpsin(&y[1], &y[1]);
1448     return;
1449
1450 /* HERE ABS(X) .LE. 1 SO USE POWER SERIES */
1451
1452 L20:
1453     mpsin1(&y[1], &y[1], &c__0);
1454 }
1455
1456
1457 void
1458 mpcosh(int *x, int *y)
1459 {
1460     static int i2;
1461
1462 /*  RETURNS Y = COSH(X) FOR MP NUMBERS X AND Y, X NOT TOO LARGE.
1463  *  USES MPEXP, DIMENSION OF R IN COMMON AT LEAST 5T+12
1464  */
1465
1466     --y;               /* Parameter adjustments */
1467     --x;
1468
1469     if (x[1] != 0) goto L10;
1470
1471 /* COSH(0) = 1 */
1472
1473     mpcim(&c__1, &y[1]);
1474     return;
1475
1476 /* CHECK LEGALITY OF B, T, M AND MXR */
1477
1478 L10:
1479     mpchk(&c__5, &c__12);
1480     i2 = (MP.t << 2) + 11;
1481     mpabs(&x[1], &MP.r[i2 - 1]);
1482
1483 /*  IF ABS(X) TOO LARGE MPEXP WILL PRINT ERROR MESSAGE
1484  *  INCREASE M TO AVOID OVERFLOW WHEN COSH(X) REPRESENTABLE
1485  */
1486
1487     MP.m += 2;
1488     mpexp(&MP.r[i2 - 1], &MP.r[i2 - 1]);
1489     mprec(&MP.r[i2 - 1], &y[1]);
1490     mpadd(&MP.r[i2 - 1], &y[1], &y[1]);
1491
1492 /*  RESTORE M.  IF RESULT OVERFLOWS OR UNDERFLOWS, MPDIVI WILL
1493  *  ACT ACCORDINGLY.
1494  */
1495
1496     MP.m += -2;
1497     mpdivi(&y[1], &c__2, &y[1]);
1498 }
1499
1500
1501 static void
1502 mpcqm(int *i, int *j, int *q)
1503 {
1504     static int i1, j1;
1505
1506 /* CONVERTS THE RATIONAL NUMBER I/J TO MULTIPLE PRECISION Q. */
1507
1508     --q;               /* Parameter adjustments */
1509
1510     i1 = *i;
1511     j1 = *j;
1512     mpgcd(&i1, &j1);
1513     if (j1 < 0)  goto L30;
1514     else if (j1 == 0) goto L10;
1515     else goto L40;
1516
1517 L10:
1518     if (v->MPerrors) {
1519         FPRINTF(stderr, "*** J = 0 IN CALL TO MPCQM ***\n");
1520     }
1521
1522     mperr();
1523     q[1] = 0;
1524     return;
1525
1526 L30:
1527     i1 = -i1;
1528     j1 = -j1;
1529
1530 L40:
1531     mpcim(&i1, &q[1]);
1532     if (j1 != 1) mpdivi(&q[1], &j1, &q[1]);
1533 }
1534
1535
1536 static void
1537 mpcrm(float *rx, int *z)
1538 {
1539     int i__1;
1540
1541     static int i, k, i2, ib, ie, re, tp, rs;
1542     static float rb, rj;
1543
1544 /*  CONVERTS SINGLE-PRECISION NUMBER RX TO MULTIPLE-PRECISION Z.
1545  *  SOME NUMBERS WILL NOT CONVERT EXACTLY ON MACHINES
1546  *  WITH BASE OTHER THAN TWO, FOUR OR SIXTEEN.
1547  *  CHECK LEGALITY OF B, T, M AND MXR
1548  */
1549
1550     --z;                   /* Parameter adjustments */
1551
1552     mpchk(&c__1, &c__4);
1553     i2 = MP.t + 4;
1554
1555 /* CHECK SIGN */
1556
1557     if (*rx < (float) 0.0) goto L20;
1558     else if (*rx == 0) goto L10;
1559     else goto L30;
1560
1561 /* IF RX = 0E0 RETURN 0 */
1562
1563 L10:
1564     z[1] = 0;
1565     return;
1566
1567 /* RX .LT. 0E0 */
1568
1569 L20:
1570     rs = -1;
1571     rj = -(double)(*rx);
1572     goto L40;
1573
1574 /* RX .GT. 0E0 */
1575
1576 L30:
1577     rs = 1;
1578     rj = *rx;
1579
1580 L40:
1581     ie = 0;
1582
1583 L50:
1584     if (rj < (float)1.0) goto L60;
1585
1586 /* INCREASE IE AND DIVIDE RJ BY 16. */
1587
1588     ++ie;
1589     rj *= (float) 0.0625;
1590     goto L50;
1591
1592 L60:
1593     if (rj >= (float).0625) goto L70;
1594
1595     --ie;
1596     rj *= (float)16.0;
1597     goto L60;
1598
1599 /*  NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
1600  *  SET EXPONENT TO 0
1601  */
1602
1603 L70:
1604     re = 0;
1605     rb = (float) MP.b;
1606
1607 /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
1608
1609     i__1 = i2;
1610     for (i = 1; i <= i__1; ++i) {
1611         rj = rb * rj;
1612         MP.r[i - 1] = (int) rj;
1613         rj -= (float) MP.r[i - 1];
1614     }
1615
1616 /* NORMALIZE RESULT */
1617
1618     mpnzr(&rs, &re, &z[1], &c__0);
1619
1620 /* Computing MAX */
1621
1622     i__1 = MP.b * 7 * MP.b;
1623     ib = max(i__1, 32767) / 16;
1624     tp = 1;
1625
1626 /* NOW MULTIPLY BY 16**IE */
1627
1628     if (ie < 0)  goto L90;
1629     else if (ie == 0) goto L130;
1630     else goto L110;
1631
1632 L90:
1633     k = -ie;
1634     i__1 = k;
1635     for (i = 1; i <= i__1; ++i) {
1636         tp <<= 4;
1637         if (tp <= ib && tp != MP.b && i < k) continue;
1638         mpdivi(&z[1], &tp, &z[1]);
1639         tp = 1;
1640     }
1641     return;
1642
1643 L110:
1644     i__1 = ie;
1645     for (i = 1; i <= i__1; ++i) {
1646         tp <<= 4;
1647         if (tp <= ib && tp != MP.b && i < ie) continue;
1648         mpmuli(&z[1], &tp, &z[1]);
1649         tp = 1;
1650     }
1651
1652 L130:
1653     return;
1654 }
1655
1656
1657 void
1658 mpdiv(int *x, int *y, int *z)
1659 {
1660     static int i, i2, ie, iz3;
1661
1662 /*  SETS Z = X/Y, FOR MP X, Y AND Z.
1663  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 4T+10
1664  *  (BUT Z(1) MAY BE R(3T+9)).
1665  *  CHECK LEGALITY OF B, T, M AND MXR
1666  */
1667
1668     --z;             /* Parameter adjustments */
1669     --y;
1670     --x;
1671
1672     mpchk(&c__4, &c__10);
1673
1674 /* CHECK FOR DIVISION BY ZERO */
1675
1676     if (y[1] != 0) goto L20;
1677
1678     if (v->MPerrors) {
1679         FPRINTF(stderr, "*** ATTEMPTED DIVISION BY ZERO IN CALL TO MPDIV ***\n");
1680     }
1681
1682     mperr();
1683     z[1] = 0;
1684     return;
1685
1686 /* SPACE USED BY MPREC IS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
1687
1688 L20:
1689     i2 = MP.t * 3 + 9;
1690
1691 /* CHECK FOR X = 0 */
1692
1693     if (x[1] != 0) goto L30;
1694
1695     z[1] = 0;
1696     return;
1697
1698 /* INCREASE M TO AVOID OVERFLOW IN MPREC */
1699
1700 L30:
1701     MP.m += 2;
1702
1703 /* FORM RECIPROCAL OF Y */
1704
1705     mprec(&y[1], &MP.r[i2 - 1]);
1706
1707 /* SET EXPONENT OF R(I2) TO ZERO TO AVOID OVERFLOW IN MPMUL */
1708
1709     ie = MP.r[i2];
1710     MP.r[i2] = 0;
1711     i = MP.r[i2 + 1];
1712
1713 /* MULTIPLY BY X */
1714
1715     mpmul(&x[1], &MP.r[i2 - 1], &z[1]);
1716     iz3 = z[3];
1717     mpext(&i, &iz3, &z[1]);
1718
1719 /* RESTORE M, CORRECT EXPONENT AND RETURN */
1720
1721     MP.m += -2;
1722     z[2] += ie;
1723     if (z[2] >= -MP.m) goto L40;
1724
1725 /* UNDERFLOW HERE */
1726
1727     mpunfl(&z[1]);
1728     return;
1729
1730 L40:
1731     if (z[2] <= MP.m) return;
1732
1733 /* OVERFLOW HERE */
1734
1735     if (v->MPerrors) {
1736         FPRINTF(stderr, "*** OVERFLOW OCCURRED IN MPDIV ***\n");
1737     }
1738
1739     mpovfl(&z[1]);
1740 }
1741
1742
1743 void
1744 mpdivi(int *x, int *iy, int *z)
1745 {
1746     int i__1, i__2;
1747
1748     static int c, i, j, k, b2, c2, i2, j1, j2, r1;
1749     static int j11, kh, re, iq, ir, rs, iqj;
1750
1751 /*  DIVIDES MP X BY THE SINGLE-PRECISION INTEGER IY GIVING MP Z.
1752  *  THIS IS MUCH FASTER THAN DIVISION BY AN MP NUMBER.
1753  */
1754
1755     --z;                  /* Parameter adjustments */
1756     --x;
1757
1758     rs = x[1];
1759     j = *iy;
1760     if (j < 0)  goto L30;
1761     else if (j == 0) goto L10;
1762     else             goto L40;
1763
1764 L10:
1765     if (v->MPerrors) {
1766         FPRINTF(stderr, "*** ATTEMPTED DIVISION BY ZERO IN CALL TO MPDIVI ***\n");
1767     }
1768
1769     goto L230;
1770
1771 L30:
1772     j = -j;
1773     rs = -rs;
1774
1775 L40:
1776     re = x[2];
1777
1778 /* CHECK FOR ZERO DIVIDEND */
1779
1780     if (rs == 0) goto L120;
1781
1782 /* CHECK FOR DIVISION BY B */
1783
1784     if (j != MP.b) goto L50;
1785     mpstr(&x[1], &z[1]);
1786     if (re <= -MP.m) goto L240;
1787     z[1] = rs;
1788     z[2] = re - 1;
1789     return;
1790
1791 /* CHECK FOR DIVISION BY 1 OR -1 */
1792
1793 L50:
1794     if (j != 1) goto L60;
1795     mpstr(&x[1], &z[1]);
1796     z[1] = rs;
1797     return;
1798
1799 L60:
1800     c = 0;
1801     i2 = MP.t + 4;
1802     i = 0;
1803
1804 /*  IF J*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
1805  *  LONG DIVISION.  ASSUME AT LEAST 16-BIT WORD.
1806  */
1807
1808 /* Computing MAX */
1809
1810     i__1 = MP.b << 3, i__2 = 32767 / MP.b;
1811     b2 = max(i__1,i__2);
1812     if (j >= b2) goto L130;
1813
1814 /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
1815
1816 L70:
1817     ++i;
1818     c = MP.b * c;
1819     if (i <= MP.t) c += x[i + 2];
1820     r1 = c / j;
1821     if (r1 < 0)  goto L210;
1822     else if (r1 == 0) goto L70;
1823     else goto L80;
1824
1825 /* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
1826
1827 L80:
1828     re = re + 1 - i;
1829     MP.r[0] = r1;
1830     c = MP.b * (c - j * r1);
1831     kh = 2;
1832     if (i >= MP.t) goto L100;
1833     kh = MP.t + 1 - i;
1834     i__1 = kh;
1835     for (k = 2; k <= i__1; ++k) {
1836         ++i;
1837         c += x[i + 2];
1838         MP.r[k - 1] = c / j;
1839         c = MP.b * (c - j * MP.r[k - 1]);
1840     }
1841     if (c < 0) goto L210;
1842     ++kh;
1843
1844 L100:
1845     i__1 = i2;
1846     for (k = kh; k <= i__1; ++k) {
1847         MP.r[k - 1] = c / j;
1848         c = MP.b * (c - j * MP.r[k - 1]);
1849     }
1850     if (c < 0) goto L210;
1851
1852 /* NORMALIZE AND ROUND RESULT */
1853
1854 L120:
1855     mpnzr(&rs, &re, &z[1], &c__0);
1856     return;
1857
1858 /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
1859
1860 L130:
1861     c2 = 0;
1862     j1 = j / MP.b;
1863     j2 = j - j1 * MP.b;
1864     j11 = j1 + 1;
1865
1866 /* LOOK FOR FIRST NONZERO DIGIT */
1867
1868 L140:
1869     ++i;
1870     c = MP.b * c + c2;
1871     c2 = 0;
1872     if (i <= MP.t) c2 = x[i + 2];
1873     if ((i__1 = c - j1) < 0) goto L140;
1874     else if (i__1 == 0) goto L150;
1875     else goto L160;
1876
1877 L150:
1878     if (c2 < j2) goto L140;
1879
1880 /* COMPUTE T+4 QUOTIENT DIGITS */
1881
1882 L160:
1883     re = re + 1 - i;
1884     k = 1;
1885     goto L180;
1886
1887 /* MAIN LOOP FOR LARGE ABS(IY) CASE */
1888
1889 L170:
1890     ++k;
1891     if (k > i2) goto L120;
1892     ++i;
1893
1894 /* GET APPROXIMATE QUOTIENT FIRST */
1895
1896 L180:
1897     ir = c / j11;
1898
1899 /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */
1900
1901     iq = c - ir * j1;
1902     if (iq < b2) goto L190;
1903
1904 /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */
1905
1906     ++ir;
1907     iq -= j1;
1908
1909 L190:
1910     iq = iq * MP.b - ir * j2;
1911     if (iq >= 0) goto L200;
1912
1913 /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
1914
1915     --ir;
1916     iq += j;
1917
1918 L200:
1919     if (i <= MP.t) iq += x[i + 2];
1920     iqj = iq / j;
1921
1922 /* R(K) = QUOTIENT, C = REMAINDER */
1923
1924     MP.r[k - 1] = iqj + ir;
1925     c = iq - j * iqj;
1926     if (c >= 0) goto L170;
1927
1928 /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
1929
1930 L210:
1931     mpchk(&c__1, &c__4);
1932
1933     if (v->MPerrors) {
1934         FPRINTF(stderr, "*** INTEGER OVERFLOW IN MPDIVI, B TOO LARGE ***\n");
1935     }
1936
1937 L230:
1938     mperr();
1939     z[1] = 0;
1940     return;
1941
1942 /* UNDERFLOW HERE */
1943
1944 L240:
1945     mpunfl(&z[1]);
1946 }
1947
1948
1949 int
1950 mpeq(int *x, int *y)
1951 {
1952     int ret_val;
1953
1954 /* RETURNS LOGICAL VALUE OF (X .EQ. Y) FOR MP X AND Y. */
1955
1956     --y;               /* Parameter adjustments */
1957     --x;
1958
1959     ret_val = mpcomp(&x[1], &y[1]) == 0;
1960
1961     return(ret_val);
1962 }
1963
1964
1965 void
1966 mperr()
1967 {
1968
1969 /*  THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
1970  *  AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
1971  */
1972
1973     doerr(_("Error"));
1974 }
1975
1976
1977 void
1978 mpexp(int *x, int *y)
1979 {
1980     int i__1, i__2;
1981     float r__1;
1982
1983     static int i, i2, i3, ie, ix, ts, xs, tss;
1984     static float rx, ry, rlb;
1985
1986 /*  RETURNS Y = EXP(X) FOR MP X AND Y.
1987  *  EXP OF INTEGER AND FRACTIONAL PARTS OF X ARE COMPUTED
1988  *  SEPARATELY.  SEE ALSO COMMENTS IN MPEXP1.
1989  *  TIME IS O(SQRT(T)M(T)).
1990  *  DIMENSION OF R MUST BE AT LEAST 4T+10 IN CALLING PROGRAM
1991  *  CHECK LEGALITY OF B, T, M AND MXR
1992  */
1993
1994     --y;                 /* Parameter adjustments */
1995     --x;
1996
1997     mpchk(&c__4, &c__10);
1998     i2 = (MP.t << 1) + 7;
1999     i3 = i2 + MP.t + 2;
2000
2001 /* CHECK FOR X = 0 */
2002
2003     if (x[1] != 0) goto L10;
2004     mpcim(&c__1, &y[1]);
2005     return;
2006
2007 /* CHECK IF ABS(X) .LT. 1 */
2008
2009 L10:
2010     if (x[2] > 0) goto L20;
2011
2012 /* USE MPEXP1 HERE */
2013
2014     mpexp1(&x[1], &y[1]);
2015     mpaddi(&y[1], &c__1, &y[1]);
2016     return;
2017
2018 /*  SEE IF ABS(X) SO LARGE THAT EXP(X) WILL CERTAINLY OVERFLOW
2019  *  OR UNDERFLOW.  1.01 IS TO ALLOW FOR ERRORS IN ALOG.
2020  */
2021
2022 L20:
2023     rlb = log((float) MP.b) * (float)1.01;
2024     r__1 = -(double)((float) (MP.m + 1)) * rlb;
2025     if (mpcmpr(&x[1], &r__1) >= 0) goto L40;
2026
2027 /* UNDERFLOW SO CALL MPUNFL AND RETURN */
2028
2029 L30:
2030     mpunfl(&y[1]);
2031     return;
2032
2033 L40:
2034     r__1 = (float) MP.m * rlb;
2035     if (mpcmpr(&x[1], &r__1) <= 0) goto L70;
2036
2037 /* OVERFLOW HERE */
2038
2039 L50:
2040     if (v->MPerrors) {
2041         FPRINTF(stderr, "*** OVERFLOW IN SUBROUTINE MPEXP ***\n");
2042     }
2043
2044     mpovfl(&y[1]);
2045     return;
2046
2047 /* NOW SAFE TO CONVERT X TO REAL */
2048
2049 L70:
2050     mpcmr(&x[1], &rx);
2051
2052 /* SAVE SIGN AND WORK WITH ABS(X) */
2053
2054     xs = x[1];
2055     mpabs(&x[1], &MP.r[i3 - 1]);
2056
2057 /*  IF ABS(X) .GT. M POSSIBLE THAT INT(X) OVERFLOWS,
2058  *  SO DIVIDE BY 32.
2059  */
2060
2061     if (dabs(rx) > (float) MP.m) {
2062         mpdivi(&MP.r[i3 - 1], &c__32, &MP.r[i3 - 1]);
2063     }
2064
2065 /* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
2066
2067     mpcmi(&MP.r[i3 - 1], &ix);
2068     mpcmf(&MP.r[i3 - 1], &MP.r[i3 - 1]);
2069
2070 /* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */
2071
2072     MP.r[i3 - 1] = xs * MP.r[i3 - 1];
2073     mpexp1(&MP.r[i3 - 1], &y[1]);
2074     mpaddi(&y[1], &c__1, &y[1]);
2075
2076 /*  COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
2077  *  (BUT ONLY ONE EXTRA DIGIT IF T .LT. 4)
2078  */
2079
2080     tss = MP.t;
2081     ts = MP.t + 2;
2082     if (MP.t < 4) ts = MP.t + 1;
2083     MP.t = ts;
2084     i2 = MP.t + 5;
2085     i3 = i2 + MP.t + 2;
2086     MP.r[i3 - 1] = 0;
2087     mpcim(&xs, &MP.r[i2 - 1]);
2088     i = 1;
2089
2090 /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
2091
2092 L80:
2093
2094 /* Computing MIN */
2095
2096     i__1 = ts, i__2 = ts + 2 + MP.r[i2];
2097     MP.t = min(i__1,i__2);
2098     if (MP.t <= 2) goto L90;
2099     ++i;
2100     i__1 = i * xs;
2101     mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]);
2102     MP.t = ts;
2103     mpadd2(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i3 - 1], &
2104             MP.r[i2 - 1], &c__0);
2105     if (MP.r[i2 - 1] != 0) goto L80;
2106
2107 /* RAISE E OR 1/E TO POWER IX */
2108
2109 L90:
2110     MP.t = ts;
2111     if (xs > 0) {
2112         mpaddi(&MP.r[i3 - 1], &c__2, &MP.r[i3 - 1]);
2113     }
2114     mppwr(&MP.r[i3 - 1], &ix, &MP.r[i3 - 1]);
2115
2116 /* RESTORE T NOW */
2117
2118     MP.t = tss;
2119
2120 /* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
2121
2122     mpmul(&y[1], &MP.r[i3 - 1], &y[1]);
2123
2124 /* MUST CORRECT RESULT IF DIVIDED BY 32 ABOVE. */
2125
2126     if (dabs(rx) <= (float) MP.m || y[1] == 0) goto L110;
2127
2128     for (i = 1; i <= 5; ++i) {
2129
2130 /* SAVE EXPONENT TO AVOID OVERFLOW IN MPMUL */
2131
2132         ie = y[2];
2133         y[2] = 0;
2134         mpmul(&y[1], &y[1], &y[1]);
2135         y[2] += ie << 1;
2136
2137 /* CHECK FOR UNDERFLOW AND OVERFLOW */
2138
2139         if (y[2] < -MP.m) goto L30;
2140         if (y[2] > MP.m)  goto L50;
2141     }
2142
2143 /*  CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
2144  *  (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
2145  */
2146
2147 L110:
2148     if (dabs(rx) > (float)10.0) return;
2149
2150     mpcmr(&y[1], &ry);
2151     if ((r__1 = ry - exp(rx), dabs(r__1)) < ry * (float)0.01) return;
2152
2153 /*  THE FOLLOWING MESSAGE MAY INDICATE THAT
2154  *  B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE
2155  *  RESULT UNDERFLOWED.
2156  */
2157
2158     if (v->MPerrors) {
2159         FPRINTF(stderr, "*** ERROR OCCURRED IN MPEXP, RESULT INCORRECT ***\n");
2160     }
2161
2162     mperr();
2163 }
2164
2165
2166 static void
2167 mpexp1(int *x, int *y)
2168 {
2169     int i__1;
2170
2171     static int i, q, i2, i3, ib, ic, ts;
2172     static float rlb;
2173
2174 /*  ASSUMES THAT X AND Y ARE MP NUMBERS,  -1 .LT. X .LT. 1.
2175  *  RETURNS Y = EXP(X) - 1 USING AN O(SQRT(T).M(T)) ALGORITHM
2176  *  DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
2177  *  PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
2178  *  SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
2179  *  ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
2180  *  UNLESS T IS VERY LARGE. SEE COMMENTS TO MPATAN AND MPPIGL.
2181  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 3T+8
2182  *  CHECK LEGALITY OF B, T, M AND MXR
2183  */
2184
2185     --y;                 /* Parameter adjustments */
2186     --x;
2187
2188     mpchk(&c__3, &c__8);
2189     i2 = MP.t + 5;
2190     i3 = i2 + MP.t + 2;
2191
2192 /* CHECK FOR X = 0 */
2193
2194     if (x[1] != 0) goto L20;
2195
2196 L10:
2197     y[1] = 0;
2198     return;
2199
2200 /* CHECK THAT ABS(X) .LT. 1 */
2201
2202 L20:
2203     if (x[2] <= 0) goto L40;
2204
2205     if (v->MPerrors) {
2206         FPRINTF(stderr, "*** ABS(X) NOT LESS THAN 1 IN CALL TO MPEXP1 ***\n");
2207     }
2208
2209     mperr();
2210     goto L10;
2211
2212 L40:
2213     mpstr(&x[1], &MP.r[i2 - 1]);
2214     rlb = log((float) MP.b);
2215
2216 /* COMPUTE APPROXIMATELY OPTIMAL Q (AND DIVIDE X BY 2**Q) */
2217
2218     q = (int) (sqrt((float) MP.t * (float).48 * rlb) + (float) x[2] *
2219               (float)1.44 * rlb);
2220
2221 /* HALVE Q TIMES */
2222
2223     if (q <= 0) goto L60;
2224     ib = MP.b << 2;
2225     ic = 1;
2226     i__1 = q;
2227     for (i = 1; i <= i__1; ++i) {
2228         ic <<= 1;
2229         if (ic < ib && ic != MP.b && i < q) continue;
2230         mpdivi(&MP.r[i2 - 1], &ic, &MP.r[i2 - 1]);
2231         ic = 1;
2232     }
2233
2234 L60:
2235     if (MP.r[i2 - 1] == 0) goto L10;
2236     mpstr(&MP.r[i2 - 1], &y[1]);
2237     mpstr(&MP.r[i2 - 1], &MP.r[i3 - 1]);
2238     i = 1;
2239     ts = MP.t;
2240
2241 /* SUM SERIES, REDUCING T WHERE POSSIBLE */
2242
2243 L70:
2244     MP.t = ts + 2 + MP.r[i3] - y[2];
2245     if (MP.t <= 2) goto L80;
2246
2247     MP.t = min(MP.t,ts);
2248     mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
2249     ++i;
2250     mpdivi(&MP.r[i3 - 1], &i, &MP.r[i3 - 1]);
2251     MP.t = ts;
2252     mpadd2(&MP.r[i3 - 1], &y[1], &y[1], &y[1], &c__0);
2253     if (MP.r[i3 - 1] != 0) goto L70;
2254
2255 L80:
2256     MP.t = ts;
2257     if (q <= 0) return;
2258
2259 /* APPLY (X+1)**2 - 1 = X(2 + X) FOR Q ITERATIONS */
2260
2261     i__1 = q;
2262     for (i = 1; i <= i__1; ++i) {
2263         mpaddi(&y[1], &c__2, &MP.r[i2 - 1]);
2264         mpmul(&MP.r[i2 - 1], &y[1], &y[1]);
2265     }
2266 }
2267
2268
2269 static void
2270 mpext(int *i, int *j, int *x)
2271 {
2272     static int q, s;
2273
2274 /*  ROUTINE CALLED BY MPDIV AND MPSQRT TO ENSURE THAT
2275  *  RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
2276  *  CAN BE.  X IS AN MP NUMBER, I AND J ARE INTEGERS.
2277  */
2278
2279     --x;               /* Parameter adjustments */
2280
2281     if (x[1] == 0 || MP.t <= 2 || *i == 0) return;
2282
2283 /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
2284
2285     q = (*j + 1) / *i + 1;
2286     s = MP.b * x[MP.t + 1] + x[MP.t + 2];
2287     if (s > q) goto L10;
2288
2289 /* SET LAST TWO DIGITS TO ZERO */
2290
2291     x[MP.t + 1] = 0;
2292     x[MP.t + 2] = 0;
2293     return;
2294
2295 L10:
2296     if (s + q < MP.b * MP.b) return;
2297
2298 /* ROUND UP HERE */
2299
2300     x[MP.t + 1] = MP.b - 1;
2301     x[MP.t + 2] = MP.b;
2302
2303 /* NORMALIZE X (LAST DIGIT B IS OK IN MPMULI) */
2304
2305     mpmuli(&x[1], &c__1, &x[1]);
2306 }
2307
2308
2309 static void
2310 mpgcd(int *k, int *l)
2311 {
2312     static int i, j, is, js;
2313
2314 /*  RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
2315  *  GREATEST COMMON DIVISOR OF K AND L.
2316  *  SAVE INPUT PARAMETERS IN LOCAL VARIABLES
2317  */
2318
2319     i = *k;
2320     j = *l;
2321     is = C_abs(i);
2322     js = C_abs(j);
2323     if (js == 0) goto L30;
2324
2325 /* EUCLIDEAN ALGORITHM LOOP */
2326
2327 L10:
2328     is %= js;
2329     if (is == 0) goto L20;
2330     js %= is;
2331     if (js != 0) goto L10;
2332     js = is;
2333
2334 /* HERE JS IS THE GCD OF I AND J */
2335
2336 L20:
2337     *k = i / js;
2338     *l = j / js;
2339     return;
2340
2341 /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */
2342
2343 L30:
2344     *k = 1;
2345     if (is == 0) *k = 0;
2346     *l = 0;
2347 }
2348
2349
2350 int
2351 mpge(int *x, int *y)
2352 {
2353     int ret_val;
2354
2355 /* RETURNS LOGICAL VALUE OF (X .GE. Y) FOR MP X AND Y. */
2356
2357     --y;               /* Parameter adjustments */
2358     --x;
2359
2360     ret_val = mpcomp(&x[1], &y[1]) >= 0;
2361
2362     return(ret_val);
2363 }
2364
2365
2366 int
2367 mpgt(int *x, int *y)
2368 {
2369     int ret_val;
2370
2371 /* RETURNS LOGICAL VALUE OF (X .GT. Y) FOR MP X AND Y. */
2372
2373     --y;             /* Parameter adjustments */
2374     --x;
2375
2376     ret_val = mpcomp(&x[1], &y[1]) > 0;
2377
2378     return(ret_val);
2379 }
2380
2381
2382 int
2383 mple(int *x, int *y)
2384 {
2385     int ret_val;
2386
2387 /* RETURNS LOGICAL VALUE OF (X .LE. Y) FOR MP X AND Y. */
2388
2389     --y;               /* Parameter adjustments */
2390     --x;
2391
2392     ret_val = mpcomp(&x[1], &y[1]) <= 0;
2393
2394     return(ret_val);
2395 }
2396
2397
2398 void
2399 mpln(int *x, int *y)
2400 {
2401     float r__1;
2402
2403     static int e, k, i2, i3;
2404     static float rx, rlx;
2405
2406 /*  RETURNS Y = LN(X), FOR MP X AND Y, USING MPLNS.
2407  *  RESTRICTION - INTEGER PART OF LN(X) MUST BE REPRESENTABLE
2408  *  AS A SINGLE-PRECISION INTEGER.  TIME IS O(SQRT(T).M(T)).
2409  *  FOR SMALL INTEGER X, MPLNI IS FASTER.
2410  *  ASYMPTOTICALLY FASTER METHODS EXIST (EG THE GAUSS-SALAMIN
2411  *  METHOD, SEE MPLNGS), BUT ARE NOT USEFUL UNLESS T IS LARGE.
2412  *  SEE COMMENTS TO MPATAN, MPEXP1 AND MPPIGL.
2413  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 6T+14.
2414  *  CHECK LEGALITY OF B, T, M AND MXR
2415  */
2416
2417     --y;                 /* Parameter adjustments */
2418     --x;
2419
2420     mpchk(&c__6, &c__14);
2421     i2 = (MP.t << 2) + 11;
2422     i3 = i2 + MP.t + 2;
2423
2424 /* CHECK THAT X IS POSITIVE */
2425     if (x[1] > 0) goto L20;
2426
2427     if (v->MPerrors) {
2428         FPRINTF(stderr, "*** X NONPOSITIVE IN CALL TO MPLN ***\n");
2429     }
2430
2431     mperr();
2432     y[1] = 0;
2433     return;
2434
2435 /* MOVE X TO LOCAL STORAGE */
2436
2437 L20:
2438     mpstr(&x[1], &MP.r[i2 - 1]);
2439     y[1] = 0;
2440     k = 0;
2441
2442 /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
2443
2444 L30:
2445     mpaddi(&MP.r[i2 - 1], &c_n1, &MP.r[i3 - 1]);
2446
2447 /* IF POSSIBLE GO TO CALL MPLNS */
2448
2449     if (MP.r[i3 - 1] == 0 || MP.r[i3] + 1 <= 0) goto L50;
2450
2451 /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
2452
2453     e = MP.r[i2];
2454     MP.r[i2] = 0;
2455     mpcmr(&MP.r[i2 - 1], &rx);
2456
2457 /* RESTORE EXPONENT AND COMPUTE SINGLE-PRECISION LOG */
2458
2459     MP.r[i2] = e;
2460     rlx = log(rx) + (float) e * log((float) MP.b);
2461     r__1 = -(double)rlx;
2462     mpcrm(&r__1, &MP.r[i3 - 1]);
2463
2464 /* UPDATE Y AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
2465
2466     mpsub(&y[1], &MP.r[i3 - 1], &y[1]);
2467     mpexp(&MP.r[i3 - 1], &MP.r[i3 - 1]);
2468
2469 /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
2470
2471     mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]);
2472
2473 /* MAKE SURE NOT LOOPING INDEFINITELY */
2474     ++k;
2475     if (k < 10) goto L30;
2476
2477     if (v->MPerrors) {
2478         FPRINTF(stderr, "*** ERROR IN MPLN, ITERATION NOT CONVERGING ***\n");
2479     }
2480
2481     mperr();
2482     return;
2483
2484 /* COMPUTE FINAL CORRECTION ACCURATELY USING MPLNS */
2485
2486 L50:
2487     mplns(&MP.r[i3 - 1], &MP.r[i3 - 1]);
2488     mpadd(&y[1], &MP.r[i3 - 1], &y[1]);
2489 }
2490
2491
2492 static void
2493 mplns(int *x, int *y)
2494 {
2495     int i__1, i__2;
2496
2497     static int i2, i3, i4, ts, it0, ts2, ts3;
2498
2499 /*  RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
2500  *  CONDITION ABS(X) .LT. 1/B, ERROR OTHERWISE.
2501  *  USES NEWTONS METHOD TO SOLVE THE EQUATION
2502  *  EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
2503  *  (HERE EXP1(Y) = EXP(Y) - 1 IS COMPUTED USING MPEXP1).
2504  *  TIME IS O(SQRT(T).M(T)) AS FOR MPEXP1, SPACE = 5T+12.
2505  *  CHECK LEGALITY OF B, T, M AND MXR
2506  */
2507
2508     --y;                 /* Parameter adjustments */
2509     --x;
2510
2511     mpchk(&c__5, &c__12);
2512     i2 = (MP.t << 1) + 7;
2513     i3 = i2 + MP.t + 2;
2514     i4 = i3 + MP.t + 2;
2515
2516 /* CHECK FOR X = 0 EXACTLY */
2517
2518     if (x[1] != 0) goto L10;
2519     y[1] = 0;
2520     return;
2521
2522 /* CHECK THAT ABS(X) .LT. 1/B */
2523
2524 L10:
2525     if (x[2] + 1 <= 0) goto L30;
2526
2527     if (v->MPerrors) {
2528         FPRINTF(stderr, "*** ABS(X) .GE. 1/B IN CALL TO MPLNS ***\n");
2529     }
2530
2531     mperr();
2532     y[1] = 0;
2533     return;
2534
2535 /* SAVE T AND GET STARTING APPROXIMATION TO -LN(1+X) */
2536
2537 L30:
2538     ts = MP.t;
2539     mpstr(&x[1], &MP.r[i3 - 1]);
2540     mpdivi(&x[1], &c__4, &MP.r[i2 - 1]);
2541     mpaddq(&MP.r[i2 - 1], &c_n1, &c__3, &MP.r[i2 - 1]);
2542     mpmul(&x[1], &MP.r[i2 - 1], &MP.r[i2 - 1]);
2543     mpaddq(&MP.r[i2 - 1], &c__1, &c__2, &MP.r[i2 - 1]);
2544     mpmul(&x[1], &MP.r[i2 - 1], &MP.r[i2 - 1]);
2545     mpaddi(&MP.r[i2 - 1], &c_n1, &MP.r[i2 - 1]);
2546     mpmul(&x[1], &MP.r[i2 - 1], &y[1]);
2547
2548 /* START NEWTON ITERATION USING SMALL T, LATER INCREASE */
2549
2550 /* Computing MAX */
2551
2552     i__1 = 5, i__2 = 13 - (MP.b << 1);
2553     MP.t = max(i__1,i__2);
2554     if (MP.t > ts) goto L80;
2555
2556     it0 = (MP.t + 5) / 2;
2557
2558 L40:
2559     mpexp1(&y[1], &MP.r[i4 - 1]);
2560     mpmul(&MP.r[i3 - 1], &MP.r[i4 - 1], &MP.r[i2 - 1]);
2561     mpadd(&MP.r[i4 - 1], &MP.r[i2 - 1], &MP.r[i4 - 1]);
2562     mpadd(&MP.r[i3 - 1], &MP.r[i4 - 1], &MP.r[i4 - 1]);
2563     mpsub(&y[1], &MP.r[i4 - 1], &y[1]);
2564     if (MP.t >= ts) goto L60;
2565
2566 /*  FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
2567  *  BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
2568  *  WE CAN ALMOST DOUBLE T EACH TIME.
2569  */
2570
2571     ts3 = MP.t;
2572     MP.t = ts;
2573
2574 L50:
2575     ts2 = MP.t;
2576     MP.t = (MP.t + it0) / 2;
2577     if (MP.t > ts3) goto L50;
2578
2579     MP.t = ts2;
2580     goto L40;
2581
2582 /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
2583
2584 L60:
2585     if (MP.r[i4 - 1] == 0 || MP.r[i4] << 1 <= it0 - MP.t) {
2586         goto L80;
2587     }
2588
2589     if (v->MPerrors) {
2590         FPRINTF(stderr, "*** ERROR OCCURRED IN MPLNS.\nNEWTON ITERATION NOT CONVERGING PROPERLY ***\n");
2591     }
2592
2593     mperr();
2594
2595 /* REVERSE SIGN OF Y AND RETURN */
2596
2597 L80:
2598     y[1] = -y[1];
2599     MP.t = ts;
2600 }
2601
2602
2603 int
2604 mplt(int *x, int *y)
2605 {
2606     int ret_val;
2607
2608 /* RETURNS LOGICAL VALUE OF (X .LT. Y) FOR MP X AND Y. */
2609
2610     --y;               /* Parameter adjustments */
2611     --x;
2612
2613     ret_val = mpcomp(&x[1], &y[1]) < 0;
2614
2615     return(ret_val);
2616 }
2617
2618
2619 static void
2620 mpmaxr(int *x)
2621 {
2622     int i__1;
2623
2624     static int i, it;
2625
2626 /*  SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
2627  *  CHECK LEGALITY OF B, T, M AND MXR
2628  */
2629
2630     --x;                /* Parameter adjustments */
2631
2632     mpchk(&c__1, &c__4);
2633     it = MP.b - 1;
2634
2635 /* SET FRACTION DIGITS TO B-1 */
2636
2637     i__1 = MP.t;
2638     for (i = 1; i <= i__1; ++i) x[i + 2] = it;
2639
2640 /* SET SIGN AND EXPONENT */
2641
2642     x[1] = 1;
2643     x[2] = MP.m;
2644 }
2645
2646
2647 static void
2648 mpmlp(int *u, int *v, int *w, int *j)
2649 {
2650     int i__1;
2651
2652     static int i;
2653
2654 /*  PERFORMS INNER MULTIPLICATION LOOP FOR MPMUL
2655  *  NOTE THAT CARRIES ARE NOT PROPAGATED IN INNER LOOP,
2656  *  WHICH SAVES TIME AT THE EXPENSE OF SPACE.
2657  */
2658
2659     --v;                 /* Parameter adjustments */
2660     --u;
2661
2662     i__1 = *j;
2663     for (i = 1; i <= i__1; ++i) u[i] += *w * v[i];
2664 }
2665
2666
2667 void
2668 mpmul(int *x, int *y, int *z)
2669 {
2670     int i__1, i__2, i__3, i__4;
2671
2672     static int c, i, j, i2, j1, re, ri, xi, rs, i2p;
2673
2674 /*  MULTIPLIES X AND Y, RETURNING RESULT IN Z, FOR MP X, Y AND Z.
2675  *  THE SIMPLE O(T**2) ALGORITHM IS USED, WITH
2676  *  FOUR GUARD DIGITS AND R*-ROUNDING.
2677  *  ADVANTAGE IS TAKEN OF ZERO DIGITS IN X, BUT NOT IN Y.
2678  *  ASYMPTOTICALLY FASTER ALGORITHMS ARE KNOWN (SEE KNUTH,
2679  *  VOL. 2), BUT ARE DIFFICULT TO IMPLEMENT IN FORTRAN IN AN
2680  *  EFFICIENT AND MACHINE-INDEPENDENT MANNER.
2681  *  IN COMMENTS TO OTHER MP ROUTINES, M(T) IS THE TIME
2682  *  TO PERFORM T-DIGIT MP MULTIPLICATION.   THUS
2683  *  M(T) = O(T**2) WITH THE PRESENT VERSION OF MPMUL,
2684  *  BUT M(T) = O(T.LOG(T).LOG(LOG(T))) IS THEORETICALLY POSSIBLE.
2685  *  CHECK LEGALITY OF B, T, M AND MXR
2686  */
2687
2688     --z;                 /* Parameter adjustments */
2689     --y;
2690     --x;
2691
2692     mpchk(&c__1, &c__4);
2693     i2 = MP.t + 4;
2694     i2p = i2 + 1;
2695
2696 /* FORM SIGN OF PRODUCT */
2697
2698     rs = x[1] * y[1];
2699     if (rs != 0) goto L10;
2700
2701 /* SET RESULT TO ZERO */
2702
2703     z[1] = 0;
2704     return;
2705
2706 /* FORM EXPONENT OF PRODUCT */
2707
2708 L10:
2709     re = x[2] + y[2];
2710
2711 /* CLEAR ACCUMULATOR */
2712
2713     i__1 = i2;
2714     for (i = 1; i <= i__1; ++i) MP.r[i - 1] = 0;
2715
2716 /* PERFORM MULTIPLICATION */
2717
2718     c = 8;
2719     i__1 = MP.t;
2720     for (i = 1; i <= i__1; ++i) {
2721         xi = x[i + 2];
2722
2723 /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
2724
2725         if (xi == 0) continue;
2726
2727 /* Computing MIN */
2728
2729         i__3 = MP.t, i__4 = i2 - i;
2730         i__2 = min(i__3,i__4);
2731         mpmlp(&MP.r[i], &y[3], &xi, &i__2);
2732         --c;
2733         if (c > 0) continue;
2734
2735 /* CHECK FOR LEGAL BASE B DIGIT */
2736
2737         if (xi < 0 || xi >= MP.b) goto L90;
2738
2739 /*  PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
2740  *  FASTER THAN DOING IT EVERY TIME.
2741  */
2742
2743         i__2 = i2;
2744         for (j = 1; j <= i__2; ++j) {
2745             j1 = i2p - j;
2746             ri = MP.r[j1 - 1] + c;
2747             if (ri < 0) goto L70;
2748             c = ri / MP.b;
2749             MP.r[j1 - 1] = ri - MP.b * c;
2750         }
2751         if (c != 0) goto L90;
2752         c = 8;
2753     }
2754
2755     if (c == 8) goto L60;
2756     if (xi < 0 || xi >= MP.b) goto L90;
2757     c = 0;
2758     i__1 = i2;
2759     for (j = 1; j <= i__1; ++j) {
2760         j1 = i2p - j;
2761         ri = MP.r[j1 - 1] + c;
2762         if (ri < 0) goto L70;
2763         c = ri / MP.b;
2764         MP.r[j1 - 1] = ri - MP.b * c;
2765     }
2766     if (c != 0) goto L90;
2767
2768 /* NORMALIZE AND ROUND RESULT */
2769
2770 L60:
2771     mpnzr(&rs, &re, &z[1], &c__0);
2772     return;
2773
2774 L70:
2775     if (v->MPerrors) {
2776         FPRINTF(stderr, "*** INTEGER OVERFLOW IN MPMUL, B TOO LARGE ***\n");
2777     }
2778
2779     goto L110;
2780
2781 L90:
2782     if (v->MPerrors) {
2783         FPRINTF(stderr, "*** ILLEGAL BASE B DIGIT IN CALL TO MPMUL.\nPOSSIBLE OVERWRITING PROBLEM ***\n");
2784     }
2785
2786 L110:
2787     mperr();
2788     z[1] = 0;
2789 }
2790
2791
2792 static void
2793 mpmul2(int *x, int *iy, int *z, int *trunc)
2794 {
2795     int i__1, i__2;
2796
2797     static int c, i, j, c1, c2, j1, j2;
2798     static int t1, t3, t4, ij, re, ri, is, ix, rs;
2799
2800 /*  MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z.
2801  *  MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER
2802  *  EVEN IF SOME DIGITS ARE GREATER THAN B-1.
2803  *  RESULT IS ROUNDED IF TRUNC.EQ.0, OTHERWISE TRUNCATED.
2804  */
2805
2806     --z;                /* Parameter adjustments */
2807     --x;
2808
2809     rs = x[1];
2810     if (rs == 0) goto L10;
2811     j = *iy;
2812     if (j < 0)  goto L20;
2813     else if (j == 0) goto L10;
2814     else goto L50;
2815
2816 /* RESULT ZERO */
2817
2818 L10:
2819     z[1] = 0;
2820     return;
2821
2822 L20:
2823     j = -j;
2824     rs = -rs;
2825
2826 /* CHECK FOR MULTIPLICATION BY B */
2827
2828     if (j != MP.b)   goto L50;
2829     if (x[2] < MP.m) goto L40;
2830
2831     mpchk(&c__1, &c__4);
2832
2833     if (v->MPerrors) {
2834         FPRINTF(stderr, "*** OVERFLOW OCCURRED IN MPMUL2 ***\n");
2835     }
2836
2837     mpovfl(&z[1]);
2838     return;
2839
2840 L40:
2841     mpstr(&x[1], &z[1]);
2842     z[1] = rs;
2843     z[2] = x[2] + 1;
2844     return;
2845
2846 /* SET EXPONENT TO EXPONENT(X) + 4 */
2847
2848 L50:
2849     re = x[2] + 4;
2850
2851 /* FORM PRODUCT IN ACCUMULATOR */
2852
2853     c = 0;
2854     t1 = MP.t + 1;
2855     t3 = MP.t + 3;
2856     t4 = MP.t + 4;
2857
2858 /*  IF J*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
2859  *  DOUBLE-PRECISION MULTIPLICATION.
2860  */
2861
2862 /* Computing MAX */
2863
2864     i__1 = MP.b << 3, i__2 = 32767 / MP.b;
2865     if (j >= max(i__1,i__2)) goto L110;
2866
2867     i__1 = MP.t;
2868     for (ij = 1; ij <= i__1; ++ij) {
2869         i = t1 - ij;
2870         ri = j * x[i + 2] + c;
2871         c = ri / MP.b;
2872         MP.r[i + 3] = ri - MP.b * c;
2873     }
2874
2875 /* CHECK FOR INTEGER OVERFLOW */
2876
2877     if (ri < 0) goto L130;
2878
2879 /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */
2880
2881     for (ij = 1; ij <= 4; ++ij) {
2882         i = 5 - ij;
2883         ri = c;
2884         c = ri / MP.b;
2885         MP.r[i - 1] = ri - MP.b * c;
2886     }
2887     if (c == 0) goto L100;
2888
2889 /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
2890
2891 L80:
2892     i__1 = t3;
2893     for (ij = 1; ij <= i__1; ++ij) {
2894         i = t4 - ij;
2895         MP.r[i] = MP.r[i - 1];
2896     }
2897     ri = c;
2898     c = ri / MP.b;
2899     MP.r[0] = ri - MP.b * c;
2900     ++re;
2901     if (c < 0)  goto L130;
2902     else if (c == 0) goto L100;
2903     else goto L80;
2904
2905 /* NORMALIZE AND ROUND OR TRUNCATE RESULT */
2906
2907 L100:
2908     mpnzr(&rs, &re, &z[1], trunc);
2909     return;
2910
2911 /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
2912
2913 L110:
2914     j1 = j / MP.b;
2915     j2 = j - j1 * MP.b;
2916
2917 /* FORM PRODUCT */
2918
2919     i__1 = t4;
2920     for (ij = 1; ij <= i__1; ++ij) {
2921         c1 = c / MP.b;
2922         c2 = c - MP.b * c1;
2923         i = t1 - ij;
2924         ix = 0;
2925         if (i > 0) ix = x[i + 2];
2926         ri = j2 * ix + c2;
2927         is = ri / MP.b;
2928         c = j1 * ix + c1 + is;
2929         MP.r[i + 3] = ri - MP.b * is;
2930     }
2931
2932     if (c < 0)  goto L130;
2933     else if (c == 0) goto L100;
2934     else goto L80;
2935
2936 /* CAN ONLY GET HERE IF INTEGER OVERFLOW OCCURRED */
2937
2938 L130:
2939     mpchk(&c__1, &c__4);
2940
2941     if (v->MPerrors) {
2942         FPRINTF(stderr, "*** INTEGER OVERFLOW IN MPMUL2, B TOO LARGE ***\n");
2943     }
2944
2945     mperr();
2946     goto L10;
2947 }
2948
2949
2950 void
2951 mpmuli(int *x, int *iy, int *z)
2952 {
2953
2954 /*  MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z.
2955  *  THIS IS FASTER THAN USING MPMUL.  RESULT IS ROUNDED.
2956  *  MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER
2957  *  EVEN IF THE LAST DIGIT IS B.
2958  */
2959
2960     --z;             /* Parameter adjustments */
2961     --x;
2962
2963     mpmul2(&x[1], iy, &z[1], &c__0);
2964 }
2965
2966
2967 static void
2968 mpmulq(int *x, int *i, int *j, int *y)
2969 {
2970     int i__1;
2971
2972     static int is, js;
2973
2974 /* MULTIPLIES MP X BY I/J, GIVING MP Y */
2975
2976     --y;                /* Parameter adjustments */
2977     --x;
2978
2979     if (*j != 0) goto L20;
2980     mpchk(&c__1, &c__4);
2981
2982     if (v->MPerrors) {
2983         FPRINTF(stderr, "*** ATTEMPTED DIVISION BY ZERO IN MPMULQ ***\n");
2984     }
2985
2986     mperr();
2987     goto L30;
2988
2989 L20:
2990     if (*i != 0) goto L40;
2991
2992 L30:
2993     y[1] = 0;
2994     return;
2995
2996 /* REDUCE TO LOWEST TERMS */
2997
2998 L40:
2999     is = *i;
3000     js = *j;
3001     mpgcd(&is, &js);
3002     if (C_abs(is) == 1) goto L50;
3003
3004     mpdivi(&x[1], &js, &y[1]);
3005     mpmul2(&y[1], &is, &y[1], &c__0);
3006     return;
3007
3008 /* HERE IS = +-1 */
3009
3010 L50:
3011     i__1 = is * js;
3012     mpdivi(&x[1], &i__1, &y[1]);
3013 }
3014
3015
3016 void
3017 mpneg(int *x, int *y)
3018 {
3019
3020 /* SETS Y = -X FOR MP NUMBERS X AND Y */
3021
3022     --y;             /* Parameter adjustments */
3023     --x;
3024
3025     mpstr(&x[1], &y[1]);
3026     y[1] = -y[1];
3027 }
3028
3029
3030 static void
3031 mpnzr(int *rs, int *re, int *z, int *trunc)
3032 {
3033     int i__1;
3034
3035     static int i, j, k, b2, i2, is, it, i2m, i2p;
3036
3037 /*  ASSUMES LONG (I.E. (T+4)-DIGIT) FRACTION IN
3038  *  R, SIGN = RS, EXPONENT = RE.  NORMALIZES,
3039  *  AND RETURNS MP RESULT IN Z.  INTEGER ARGUMENTS RS AND RE
3040  *  ARE NOT PRESERVED. R*-ROUNDING IS USED IF TRUNC.EQ.0
3041  */
3042
3043     --z;                 /* Parameter adjustments */
3044
3045     i2 = MP.t + 4;
3046     if (*rs != 0) goto L20;
3047
3048 /* STORE ZERO IN Z */
3049
3050 L10:
3051     z[1] = 0;
3052     return;
3053
3054 /* CHECK THAT SIGN = +-1 */
3055
3056 L20:
3057     if (C_abs(*rs) <= 1) goto L40;
3058
3059     if (v->MPerrors) {
3060         FPRINTF(stderr, "*** SIGN NOT 0, +1 OR -1 IN CALL TO MPNZR.\nPOSSIBLE OVERWRITING PROBLEM ***\n");
3061     }
3062
3063     mperr();
3064     goto L10;
3065
3066 /* LOOK FOR FIRST NONZERO DIGIT */
3067
3068 L40:
3069     i__1 = i2;
3070     for (i = 1; i <= i__1; ++i) {
3071         is = i - 1;
3072         if (MP.r[i - 1] > 0) goto L60;
3073     }
3074
3075 /* FRACTION ZERO */
3076
3077     goto L10;
3078
3079 L60:
3080     if (is == 0) goto L90;
3081
3082 /* NORMALIZE */
3083
3084     *re -= is;
3085     i2m = i2 - is;
3086     i__1 = i2m;
3087     for (j = 1; j <= i__1; ++j) {
3088         k = j + is;
3089         MP.r[j - 1] = MP.r[k - 1];
3090     }
3091     i2p = i2m + 1;
3092     i__1 = i2;
3093     for (j = i2p; j <= i__1; ++j) MP.r[j - 1] = 0;
3094
3095 /* CHECK TO SEE IF TRUNCATION IS DESIRED */
3096
3097 L90:
3098     if (*trunc != 0) goto L150;
3099
3100 /*  SEE IF ROUNDING NECESSARY
3101  *  TREAT EVEN AND ODD BASES DIFFERENTLY
3102  */
3103
3104     b2 = MP.b / 2;
3105     if (b2 << 1 != MP.b) goto L130;
3106
3107 /*  B EVEN.  ROUND IF R(T+1).GE.B2 UNLESS R(T) ODD AND ALL ZEROS
3108  *  AFTER R(T+2).
3109  */
3110
3111     if ((i__1 = MP.r[MP.t] - b2) < 0) goto L150;
3112     else if (i__1 == 0) goto L100;
3113     else goto L110;
3114
3115 L100:
3116     if (MP.r[MP.t - 1] % 2 == 0) goto L110;
3117
3118     if (MP.r[MP.t + 1] + MP.r[MP.t + 2] + MP.r[MP.t + 3] == 0) {
3119         goto L150;
3120     }
3121
3122 /* ROUND */
3123
3124 L110:
3125     i__1 = MP.t;
3126     for (j = 1; j <= i__1; ++j) {
3127         i = MP.t + 1 - j;
3128         ++MP.r[i - 1];
3129         if (MP.r[i - 1] < MP.b) goto L150;
3130         MP.r[i - 1] = 0;
3131     }
3132
3133 /* EXCEPTIONAL CASE, ROUNDED UP TO .10000... */
3134
3135     ++(*re);
3136     MP.r[0] = 1;
3137     goto L150;
3138
3139 /* ODD BASE, ROUND IF R(T+1)... .GT. 1/2 */
3140
3141 L130:
3142     for (i = 1; i <= 4; ++i) {
3143       it = MP.t + i;
3144       if ((i__1 = MP.r[it - 1] - b2) < 0) goto L150;
3145       else if (i__1 == 0) continue;
3146       else goto L110;
3147     }
3148
3149 /* CHECK FOR OVERFLOW */
3150
3151 L150:
3152     if (*re <= MP.m) goto L170;
3153
3154     if (v->MPerrors) {
3155         FPRINTF(stderr, "*** OVERFLOW OCCURRED IN MPNZR ***\n");
3156     }
3157
3158     mpovfl(&z[1]);
3159     return;
3160
3161 /* CHECK FOR UNDERFLOW */
3162
3163 L170:
3164     if (*re < -MP.m) goto L190;
3165
3166 /* STORE RESULT IN Z */
3167
3168     z[1] = *rs;
3169     z[2] = *re;
3170     i__1 = MP.t;
3171     for (i = 1; i <= i__1; ++i) z[i + 2] = MP.r[i - 1];
3172     return;
3173
3174 /* UNDERFLOW HERE */
3175
3176 L190:
3177     mpunfl(&z[1]);
3178 }
3179
3180
3181 static void
3182 mpovfl(int *x)
3183 {
3184
3185 /*  CALLED ON MULTIPLE-PRECISION OVERFLOW, IE WHEN THE
3186  *  EXPONENT OF MP NUMBER X WOULD EXCEED M.
3187  *  AT PRESENT EXECUTION IS TERMINATED WITH AN ERROR MESSAGE
3188  *  AFTER CALLING MPMAXR(X), BUT IT WOULD BE POSSIBLE TO RETURN,
3189  *  POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION AFTER
3190  *  A PRESET NUMBER OF OVERFLOWS.  ACTION COULD EASILY BE DETERMINED
3191  *  BY A FLAG IN LABELLED COMMON.
3192  *  M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC.
3193  */
3194
3195     --x;                 /* Parameter adjustments */
3196
3197     mpchk(&c__1, &c__4);
3198
3199 /* SET X TO LARGEST POSSIBLE POSITIVE NUMBER */
3200
3201     mpmaxr(&x[1]);
3202
3203     if (v->MPerrors) {
3204         FPRINTF(stderr, "*** CALL TO MPOVFL, MP OVERFLOW OCCURRED ***\n");
3205     }
3206
3207 /* TERMINATE EXECUTION BY CALLING MPERR */
3208
3209     mperr();
3210 }
3211
3212
3213 void
3214 mppi(int *x)
3215 {
3216     float r__1;
3217
3218     static int i2;
3219     static float rx;
3220
3221 /*  SETS MP X = PI TO THE AVAILABLE PRECISION.
3222  *  USES PI/4 = 4.ARCTAN(1/5) - ARCTAN(1/239).
3223  *  TIME IS O(T**2).
3224  *  DIMENSION OF R MUST BE AT LEAST 3T+8
3225  *  CHECK LEGALITY OF B, T, M AND MXR
3226  */
3227
3228     --x;                /* Parameter adjustments */
3229
3230     mpchk(&c__3, &c__8);
3231
3232 /* ALLOW SPACE FOR MPART1 */
3233
3234     i2 = (MP.t << 1) + 7;
3235     mpart1(&c__5, &MP.r[i2 - 1]);
3236     mpmuli(&MP.r[i2 - 1], &c__4, &MP.r[i2 - 1]);
3237     mpart1(&c__239, &x[1]);
3238     mpsub(&MP.r[i2 - 1], &x[1], &x[1]);
3239     mpmuli(&x[1], &c__4, &x[1]);
3240
3241 /* RETURN IF ERROR IS LESS THAN 0.01 */
3242
3243     mpcmr(&x[1], &rx);
3244     if ((r__1 = rx - (float)3.1416, dabs(r__1)) < (float) 0.01) return;
3245
3246 /* FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL */
3247
3248     if (v->MPerrors) {
3249         FPRINTF(stderr, "*** ERROR OCCURRED IN MPPI, RESULT INCORRECT ***\n");
3250     }
3251
3252     mperr();
3253 }
3254
3255
3256 void
3257 mppwr(int *x, int *n, int *y)
3258 {
3259     static int i2, n2, ns;
3260
3261 /*  RETURNS Y = X**N, FOR MP X AND Y, INTEGER N, WITH 0**0 = 1.
3262  *  R MUST BE DIMENSIONED AT LEAST 4T+10 IN CALLING PROGRAM
3263  *  (2T+6 IS ENOUGH IF N NONNEGATIVE).
3264  */
3265
3266     --y;            /* Parameter adjustments */
3267     --x;
3268
3269     i2 = MP.t + 5;
3270     n2 = *n;
3271     if (n2 < 0)  goto L20;
3272     else if (n2 == 0) goto L10;
3273     else              goto L40;
3274
3275 /* N = 0, RETURN Y = 1. */
3276
3277 L10:
3278     mpcim(&c__1, &y[1]);
3279     return;
3280
3281 /* N .LT. 0 */
3282
3283 L20:
3284     mpchk(&c__4, &c__10);
3285     n2 = -n2;
3286     if (x[1] != 0) goto L60;
3287
3288     if (v->MPerrors) {
3289         FPRINTF(stderr, "*** ATTEMPT TO RAISE ZERO TO NEGATIVE POWER IN CALL TO SUBROUTINE MPPWR ***\n");
3290     }
3291
3292     mperr();
3293     goto L50;
3294
3295 /* N .GT. 0 */
3296
3297 L40:
3298     mpchk(&c__2, &c__6);
3299     if (x[1] != 0) goto L60;
3300
3301 /* X = 0, N .GT. 0, SO Y = 0 */
3302
3303 L50:
3304     y[1] = 0;
3305     return;
3306
3307 /* MOVE X */
3308
3309 L60:
3310     mpstr(&x[1], &y[1]);
3311
3312 /* IF N .LT. 0 FORM RECIPROCAL */
3313
3314     if (*n < 0) mprec(&y[1], &y[1]);
3315     mpstr(&y[1], &MP.r[i2 - 1]);
3316
3317 /* SET PRODUCT TERM TO ONE */
3318
3319     mpcim(&c__1, &y[1]);
3320
3321 /* MAIN LOOP, LOOK AT BITS OF N2 FROM RIGHT */
3322
3323 L70:
3324     ns = n2;
3325     n2 /= 2;
3326     if (n2 << 1 != ns) mpmul(&y[1], &MP.r[i2 - 1], &y[1]);
3327     if (n2 <= 0) return;
3328
3329     mpmul(&MP.r[i2 - 1], &MP.r[i2 - 1], &MP.r[i2 - 1]);
3330     goto L70;
3331 }
3332
3333
3334 void
3335 mppwr2(int *x, int *y, int *z)
3336 {
3337     static int i2;
3338
3339 /*  RETURNS Z = X**Y FOR MP NUMBERS X, Y AND Z, WHERE X IS
3340  *  POSITIVE (X .EQ. 0 ALLOWED IF Y .GT. 0).  SLOWER THAN
3341  *  MPPWR AND MPQPWR, SO USE THEM IF POSSIBLE.
3342  *  DIMENSION OF R IN COMMON AT LEAST 7T+16
3343  *  CHECK LEGALITY OF B, T, M AND MXR
3344  */
3345
3346     --z;            /* Parameter adjustments */
3347     --y;
3348     --x;
3349
3350     mpchk(&c__7, &c__16);
3351     if (x[1] < 0)  goto L10;
3352     else if (x[1] == 0) goto L30;
3353     else goto L70;
3354
3355 L10:
3356     show_error(_("Negative X and non-integer Y not supported"));
3357     goto L50;
3358
3359 /* HERE X IS ZERO, RETURN ZERO IF Y POSITIVE, OTHERWISE ERROR */
3360
3361 L30:
3362     if (y[1] > 0) goto L60;
3363
3364     if (v->MPerrors) {
3365         FPRINTF(stderr, "*** X ZERO AND Y NONPOSITIVE IN CALL TO MPPWR2 ***\n");
3366     }
3367
3368 L50:
3369     mperr();
3370
3371 /* RETURN ZERO HERE */
3372
3373 L60:
3374     z[1] = 0;
3375     return;
3376
3377 /*  USUAL CASE HERE, X POSITIVE
3378  *  USE MPLN AND MPEXP TO COMPUTE POWER
3379  */
3380
3381 L70:
3382     i2 = MP.t * 6 + 15;
3383     mpln(&x[1], &MP.r[i2 - 1]);
3384     mpmul(&y[1], &MP.r[i2 - 1], &z[1]);
3385
3386 /* IF X**Y IS TOO LARGE, MPEXP WILL PRINT ERROR MESSAGE */
3387
3388     mpexp(&z[1], &z[1]);
3389 }
3390
3391
3392 static void
3393 mprec(int *x, int *y)
3394 {
3395
3396 /* Initialized data */
3397
3398     static int it[9] = { 0, 8, 6, 5, 4, 4, 4, 4, 4 };
3399
3400     float r__1;
3401
3402     static int i2, i3, ex, ts, it0, ts2, ts3;
3403     static float rx;
3404
3405 /*  RETURNS Y = 1/X, FOR MP X AND Y.
3406  *  MPROOT (X, -1, Y) HAS THE SAME EFFECT.
3407  *  DIMENSION OF R MUST BE AT LEAST 4*T+10 IN CALLING PROGRAM
3408  *  (BUT Y(1) MAY BE R(3T+9)).
3409  *  NEWTONS METHOD IS USED, SO FINAL ONE OR TWO DIGITS MAY
3410  *  NOT BE CORRECT.
3411  */
3412
3413     --y;             /* Parameter adjustments */
3414     --x;
3415
3416 /* CHECK LEGALITY OF B, T, M AND MXR */
3417
3418     mpchk(&c__4, &c__10);
3419
3420 /* MPADDI REQUIRES 2T+6 WORDS. */
3421
3422     i2 = (MP.t << 1) + 7;
3423     i3 = i2 + MP.t + 2;
3424     if (x[1] != 0) goto L20;
3425
3426     if (v->MPerrors) {
3427         FPRINTF(stderr, "*** ATTEMPTED DIVISION BY ZERO IN CALL TO MPREC ***\n");
3428     }
3429
3430     mperr();
3431     y[1] = 0;
3432     return;
3433
3434 L20:
3435     ex = x[2];
3436
3437 /* TEMPORARILY INCREASE M TO AVOID OVERFLOW */
3438
3439     MP.m += 2;
3440
3441 /* SET EXPONENT TO ZERO SO RX NOT TOO LARGE OR SMALL. */
3442
3443     x[2] = 0;
3444     mpcmr(&x[1], &rx);
3445
3446 /* USE SINGLE-PRECISION RECIPROCAL AS FIRST APPROXIMATION */
3447
3448     r__1 = (float)1. / rx;
3449     mpcrm(&r__1, &MP.r[i2 - 1]);
3450
3451 /* RESTORE EXPONENT */
3452
3453     x[2] = ex;
3454
3455 /* CORRECT EXPONENT OF FIRST APPROXIMATION */
3456
3457     MP.r[i2] -= ex;
3458
3459 /* SAVE T (NUMBER OF DIGITS) */
3460
3461     ts = MP.t;
3462
3463 /* START WITH SMALL T TO SAVE TIME. ENSURE THAT B**(T-1) .GE. 100 */
3464
3465     MP.t = 3;
3466     if (MP.b < 10) MP.t = it[MP.b - 1];
3467     it0 = (MP.t + 4) / 2;
3468     if (MP.t > ts) goto L70;
3469
3470 /* MAIN ITERATION LOOP */
3471
3472 L30:
3473     mpmul(&x[1], &MP.r[i2 - 1], &MP.r[i3 - 1]);
3474     mpaddi(&MP.r[i3 - 1], &c_n1, &MP.r[i3 - 1]);
3475
3476 /* TEMPORARILY REDUCE T */
3477
3478     ts3 = MP.t;
3479     MP.t = (MP.t + it0) / 2;
3480     mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
3481
3482 /* RESTORE T */
3483
3484     MP.t = ts3;
3485     mpsub(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]);
3486     if (MP.t >= ts) goto L50;
3487
3488 /*  FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
3489  *  BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
3490  */
3491
3492     MP.t = ts;
3493
3494 L40:
3495     ts2 = MP.t;
3496     MP.t = (MP.t + it0) / 2;
3497     if (MP.t > ts3) goto L40;
3498
3499     MP.t = min(ts,ts2);
3500     goto L30;
3501
3502 /* RETURN IF NEWTON ITERATION WAS CONVERGING */
3503
3504 L50:
3505     if (MP.r[i3 - 1] == 0 || (MP.r[i2] - MP.r[i3]) << 1 >= MP.t - it0) {
3506         goto L70;
3507     }
3508
3509 /*  THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
3510  *  OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
3511  */
3512
3513     if (v->MPerrors) {
3514         FPRINTF(stderr, "*** ERROR OCCURRED IN MPREC, NEWTON ITERATION NOT CONVERGING PROPERLY ***\n");
3515     }
3516
3517     mperr();
3518
3519 /* MOVE RESULT TO Y AND RETURN AFTER RESTORING T */
3520
3521 L70:
3522     MP.t = ts;
3523     mpstr(&MP.r[i2 - 1], &y[1]);
3524
3525 /* RESTORE M AND CHECK FOR OVERFLOW (UNDERFLOW IMPOSSIBLE) */
3526
3527     MP.m += -2;
3528     if (y[2] <= MP.m) return;
3529
3530     if (v->MPerrors) {
3531         FPRINTF(stderr, "*** OVERFLOW OCCURRED IN MPREC ***\n");
3532     }
3533
3534     mpovfl(&y[1]);
3535 }
3536
3537
3538 static void
3539 mproot(int *x, int *n, int *y)
3540 {
3541
3542 /* Initialized data */
3543
3544     static int it[9] = { 0, 8, 6, 5, 4, 4, 4, 4, 4 };
3545
3546     int i__1;
3547     float r__1;
3548
3549     static int i2, i3, ex, np, ts, it0, ts2, ts3, xes;
3550     static float rx;
3551
3552 /*  RETURNS Y = X**(1/N) FOR INTEGER N, ABS(N) .LE. MAX (B, 64).
3553  *  AND MP NUMBERS X AND Y,
3554  *  USING NEWTONS METHOD WITHOUT DIVISIONS.   SPACE = 4T+10
3555  *  (BUT Y(1) MAY BE R(3T+9))
3556  */
3557
3558     --y;             /* Parameter adjustments */
3559     --x;
3560
3561 /* CHECK LEGALITY OF B, T, M AND MXR */
3562
3563     mpchk(&c__4, &c__10);
3564     if (*n != 1) goto L10;
3565     mpstr(&x[1], &y[1]);
3566     return;
3567
3568 L10:
3569     i2 = (MP.t << 1) + 7;
3570     i3 = i2 + MP.t + 2;
3571     if (*n != 0) goto L30;
3572
3573     if (v->MPerrors) {
3574         FPRINTF(stderr, "*** N = 0 IN CALL TO MPROOT ***\n");
3575     }
3576
3577     goto L50;
3578
3579 L30:
3580     np = C_abs(*n);
3581
3582 /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP .LE. MAX (B, 64) */
3583
3584     if (np <= max(MP.b,64)) goto L60;
3585
3586     if (v->MPerrors) {
3587         FPRINTF(stderr, "*** ABS(N) TOO LARGE IN CALL TO MPROOT ***\n");
3588     }
3589
3590 L50:
3591     mperr();
3592     y[1] = 0;
3593     return;
3594
3595 /* LOOK AT SIGN OF X */
3596
3597 L60:
3598     if (x[1] < 0)  goto L90;
3599     else if (x[1] == 0) goto L70;
3600     else goto L110;
3601
3602 /* X = 0 HERE, RETURN 0 IF N POSITIVE, ERROR IF NEGATIVE */
3603
3604 L70:
3605     y[1] = 0;
3606     if (*n > 0) return;
3607
3608     if (v->MPerrors) {
3609         FPRINTF(stderr, "*** X = 0 AND N NEGATIVE IN CALL TO MPROOT ***\n");
3610     }
3611
3612     goto L50;
3613
3614 /* X NEGATIVE HERE, SO ERROR IF N IS EVEN */
3615
3616 L90:
3617     if (np % 2 != 0) goto L110;
3618
3619     if (v->MPerrors) {
3620         FPRINTF(stderr, "*** X NEGATIVE AND N EVEN IN CALL TO MPROOT ***\n");
3621     }
3622
3623     goto L50;
3624
3625 /* GET EXPONENT AND DIVIDE BY NP */
3626
3627 L110:
3628     xes = x[2];
3629     ex = xes / np;
3630
3631 /* REDUCE EXPONENT SO RX NOT TOO LARGE OR SMALL. */
3632
3633     x[2] = 0;
3634     mpcmr(&x[1], &rx);
3635
3636 /* USE SINGLE-PRECISION ROOT FOR FIRST APPROXIMATION */
3637
3638     r__1 = exp(((float) (np * ex - xes) * log((float) MP.b) -
3639            log((dabs(rx)))) / (float) np);
3640     mpcrm(&r__1, &MP.r[i2 - 1]);
3641
3642 /* SIGN OF APPROXIMATION SAME AS THAT OF X */
3643
3644     MP.r[i2 - 1] = x[1];
3645
3646 /* RESTORE EXPONENT */
3647
3648     x[2] = xes;
3649
3650 /* CORRECT EXPONENT OF FIRST APPROXIMATION */
3651
3652     MP.r[i2] -= ex;
3653
3654 /* SAVE T (NUMBER OF DIGITS) */
3655
3656     ts = MP.t;
3657
3658 /* START WITH SMALL T TO SAVE TIME */
3659
3660     MP.t = 3;
3661
3662 /* ENSURE THAT B**(T-1) .GE. 100 */
3663
3664     if (MP.b < 10) MP.t = it[MP.b - 1];
3665     if (MP.t > ts) goto L160;
3666
3667 /* IT0 IS A NECESSARY SAFETY FACTOR */
3668
3669     it0 = (MP.t + 4) / 2;
3670
3671 /* MAIN ITERATION LOOP */
3672
3673 L120:
3674     mppwr(&MP.r[i2 - 1], &np, &MP.r[i3 - 1]);
3675     mpmul(&x[1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
3676     mpaddi(&MP.r[i3 - 1], &c_n1, &MP.r[i3 - 1]);
3677
3678 /* TEMPORARILY REDUCE T */
3679
3680     ts3 = MP.t;
3681     MP.t = (MP.t + it0) / 2;
3682     mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]);
3683     mpdivi(&MP.r[i3 - 1], &np, &MP.r[i3 - 1]);
3684
3685 /* RESTORE T */
3686
3687     MP.t = ts3;
3688     mpsub(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]);
3689
3690 /*  FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
3691  *  NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
3692  */
3693
3694     if (MP.t >= ts) goto L140;
3695     MP.t = ts;
3696
3697 L130:
3698     ts2 = MP.t;
3699     MP.t = (MP.t + it0) / 2;
3700     if (MP.t > ts3) goto L130;
3701
3702     MP.t = min(ts,ts2);
3703     goto L120;
3704
3705 /*  NOW R(I2) IS X**(-1/NP)
3706  *  CHECK THAT NEWTON ITERATION WAS CONVERGING
3707  */
3708
3709 L140:
3710     if (MP.r[i3 - 1] == 0 || (MP.r[i2] - MP.r[i3]) << 1 >= MP.t - it0) {
3711         goto L160;
3712     }
3713
3714 /*  THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
3715  *  OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
3716  *  IS NOT ACCURATE ENOUGH.
3717  */
3718
3719     if (v->MPerrors) {
3720         FPRINTF(stderr, "*** ERROR OCCURRED IN MPROOT, NEWTON ITERATION NOT CONVERGING PROPERLY ***\n");
3721     }
3722
3723     mperr();
3724
3725 /* RESTORE T */
3726
3727 L160:
3728     MP.t = ts;
3729     if (*n < 0) goto L170;
3730
3731     i__1 = *n - 1;
3732     mppwr(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]);
3733     mpmul(&x[1], &MP.r[i2 - 1], &y[1]);
3734     return;
3735
3736 L170:
3737     mpstr(&MP.r[i2 - 1], &y[1]);
3738 }
3739
3740
3741 void
3742 mpset(int *idecpl, int *itmax2, int *maxdr)
3743 {
3744     int i__1;
3745
3746     static int i, k, w, i2, w2, wn;
3747
3748 /*  SETS BASE (B) AND NUMBER OF DIGITS (T) TO GIVE THE
3749  *  EQUIVALENT OF AT LEAST IDECPL DECIMAL DIGITS.
3750  *  IDECPL SHOULD BE POSITIVE.
3751  *  ITMAX2 IS THE DIMENSION OF ARRAYS USED FOR MP NUMBERS,
3752  *  SO AN ERROR OCCURS IF THE COMPUTED T EXCEEDS ITMAX2 - 2.
3753  *  MPSET ALSO SETS
3754  *        MXR = MAXDR (DIMENSION OF R IN COMMON, .GE. T+4), AND
3755  *          M = (W-1)/4 (MAXIMUM ALLOWABLE EXPONENT),
3756  *  WHERE W IS THE LARGEST INTEGER OF THE FORM 2**K-1 WHICH IS
3757  *  REPRESENTABLE IN THE MACHINE, K .LE. 47
3758  *  THE COMPUTED B AND T SATISFY THE CONDITIONS
3759  *  (T-1)*LN(B)/LN(10) .GE. IDECPL   AND   8*B*B-1 .LE. W .
3760  *  APPROXIMATELY MINIMAL T AND MAXIMAL B SATISFYING
3761  *  THESE CONDITIONS ARE CHOSEN.
3762  *  PARAMETERS IDECPL, ITMAX2 AND MAXDR ARE INTEGERS.
3763  *  BEWARE - MPSET WILL CAUSE AN INTEGER OVERFLOW TO OCCUR
3764  *  ******   IF WORDLENGTH IS LESS THAN 48 BITS.
3765  *           IF THIS IS NOT ALLOWABLE, CHANGE THE DETERMINATION
3766  *           OF W (DO 30 ... TO 30 W = WN) OR SET B, T, M,
3767  *           AND MXR WITHOUT CALLING MPSET.
3768  *  FIRST SET MXR
3769  */
3770
3771     MP.mxr = *maxdr;
3772
3773 /* DETERMINE LARGE REPRESENTABLE INTEGER W OF FORM 2**K - 1 */
3774
3775     w = 0;
3776     k = 0;
3777
3778 /*  ON CYBER 76 HAVE TO FIND K .LE. 47, SO ONLY LOOP
3779  *  47 TIMES AT MOST.  IF GENUINE INTEGER WORDLENGTH
3780  *  EXCEEDS 47 BITS THIS RESTRICTION CAN BE RELAXED.
3781  */
3782
3783     for (i = 1; i <= 47; ++i) {
3784
3785 /*  INTEGER OVERFLOW WILL EVENTUALLY OCCUR HERE
3786  *  IF WORDLENGTH .LT. 48 BITS
3787  */
3788
3789         w2 = w + w;
3790         wn = w2 + 1;
3791
3792 /*  APPARENTLY REDUNDANT TESTS MAY BE NECESSARY ON SOME
3793  *  MACHINES, DEPENDING ON HOW INTEGER OVERFLOW IS HANDLED
3794  */
3795
3796         if (wn <= w || wn <= w2 || wn <= 0) goto L40;
3797         k = i;
3798         w = wn;
3799     }
3800
3801 /* SET MAXIMUM EXPONENT TO (W-1)/4 */
3802
3803 L40:
3804     MP.m = w / 4;
3805     if (*idecpl > 0) goto L60;
3806
3807     if (v->MPerrors) {
3808         FPRINTF(stderr, "*** IDECPL .LE. 0 IN CALL TO MPSET ***\n");
3809     }
3810
3811     mperr();
3812     return;
3813
3814 /* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) .LE. W */
3815
3816 L60:
3817     i__1 = (k - 3) / 2;
3818     MP.b = pow_ii(&c__2, &i__1);
3819
3820 /* 2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT */
3821
3822     MP.t = (int) ((float) (*idecpl) * log((float)10.) / log((float) MP.b) +
3823                   (float) 2.0);
3824
3825 /* SEE IF T TOO LARGE FOR DIMENSION STATEMENTS */
3826
3827     i2 = MP.t + 2;
3828     if (i2 <= *itmax2) goto L80;
3829
3830     if (v->MPerrors) {
3831         FPRINTF(stderr,
3832           "ITMAX2 TOO SMALL IN CALL TO MPSET ***\n*** INCREASE ITMAX2 AND DIMENSIONS OF MP ARRAYS TO AT LEAST %d ***\n",
3833           i2);
3834     }
3835
3836     mperr();
3837
3838 /* REDUCE TO MAXIMUM ALLOWED BY DIMENSION STATEMENTS */
3839
3840     MP.t = *itmax2 - 2;
3841
3842 /* CHECK LEGALITY OF B, T, M AND MXR (AT LEAST T+4) */
3843
3844 L80:
3845     mpchk(&c__1, &c__4);
3846 }
3847
3848
3849 void
3850 mpsin(int *x, int *y)
3851 {
3852     float r__1;
3853
3854     static int i2, i3, ie, xs;
3855     static float rx, ry;
3856
3857 /*  RETURNS Y = SIN(X) FOR MP X AND Y,
3858  *  METHOD IS TO REDUCE X TO (-1, 1) AND USE MPSIN1, SO
3859  *  TIME IS O(M(T)T/LOG(T)).
3860  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 5T+12
3861  *  CHECK LEGALITY OF B, T, M AND MXR
3862  */
3863
3864     --y;               /* Parameter adjustments */
3865     --x;
3866
3867     mpchk(&c__5, &c__12);
3868     i2 = (MP.t << 2) + 11;
3869     if (x[1] != 0) goto L20;
3870
3871 L10:
3872     y[1] = 0;
3873     return;
3874
3875 L20:
3876     xs = x[1];
3877     ie = C_abs(x[2]);
3878     if (ie <= 2) mpcmr(&x[1], &rx);
3879
3880     mpabs(&x[1], &MP.r[i2 - 1]);
3881
3882 /* USE MPSIN1 IF ABS(X) .LE. 1 */
3883
3884     if (mpcmpi(&MP.r[i2 - 1], &c__1) > 0) goto L30;
3885
3886     mpsin1(&MP.r[i2 - 1], &y[1], &c__1);
3887     goto L50;
3888
3889 /*  FIND ABS(X) MODULO 2PI (IT WOULD SAVE TIME IF PI WERE
3890  *  PRECOMPUTED AND SAVED IN COMMON).
3891  *  FOR INCREASED ACCURACY COMPUTE PI/4 USING MPART1
3892  */
3893
3894 L30:
3895     i3 = (MP.t << 1) + 7;
3896     mpart1(&c__5, &MP.r[i3 - 1]);
3897     mpmuli(&MP.r[i3 - 1], &c__4, &MP.r[i3 - 1]);
3898     mpart1(&c__239, &y[1]);
3899     mpsub(&MP.r[i3 - 1], &y[1], &y[1]);
3900     mpdiv(&MP.r[i2 - 1], &y[1], &MP.r[i2 - 1]);
3901     mpdivi(&MP.r[i2 - 1], &c__8, &MP.r[i2 - 1]);
3902     mpcmf(&MP.r[i2 - 1], &MP.r[i2 - 1]);
3903
3904 /* SUBTRACT 1/2, SAVE SIGN AND TAKE ABS */
3905
3906     mpaddq(&MP.r[i2 - 1], &c_n1, &c__2, &MP.r[i2 - 1]);
3907     xs = -xs * MP.r[i2 - 1];
3908     if (xs == 0) goto L10;
3909
3910     MP.r[i2 - 1] = 1;
3911     mpmuli(&MP.r[i2 - 1], &c__4, &MP.r[i2 - 1]);
3912
3913 /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */
3914
3915     if (MP.r[i2] > 0)  {
3916         mpaddi(&MP.r[i2 - 1], &c_n2, &MP.r[i2 - 1]);
3917     }
3918
3919     if (MP.r[i2 - 1] == 0) goto L10;
3920
3921     MP.r[i2 - 1] = 1;
3922     mpmuli(&MP.r[i2 - 1], &c__2, &MP.r[i2 - 1]);
3923
3924 /*  NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE
3925  *  POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT
3926  */
3927
3928     if (MP.r[i2] > 0) goto L40;
3929
3930     mpmul(&MP.r[i2 - 1], &y[1], &MP.r[i2 - 1]);
3931     mpsin1(&MP.r[i2 - 1], &y[1], &c__1);
3932     goto L50;
3933
3934 L40:
3935     mpaddi(&MP.r[i2 - 1], &c_n2, &MP.r[i2 - 1]);
3936     mpmul(&MP.r[i2 - 1], &y[1], &MP.r[i2 - 1]);
3937     mpsin1(&MP.r[i2 - 1], &y[1], &c__0);
3938
3939 L50:
3940     y[1] = xs;
3941     if (ie > 2) return;
3942
3943 /*  CHECK THAT ABSOLUTE ERROR LESS THAN 0.01 IF ABS(X) .LE. 100
3944  *  (IF ABS(X) IS LARGE THEN SINGLE-PRECISION SIN INACCURATE)
3945  */
3946
3947     if (dabs(rx) > (float)100.) return;
3948
3949     mpcmr(&y[1], &ry);
3950     if ((r__1 = ry - sin(rx), dabs(r__1)) < (float) 0.01) return;
3951
3952 /*  THE FOLLOWING MESSAGE MAY INDICATE THAT
3953  *  B**(T-1) IS TOO SMALL.
3954  */
3955
3956     if (v->MPerrors) {
3957         FPRINTF(stderr, "*** ERROR OCCURRED IN MPSIN, RESULT INCORRECT ***\n");
3958     }
3959
3960     mperr();
3961 }
3962
3963
3964 static void
3965 mpsin1(int *x, int *y, int *is)
3966 {
3967     int i__1;
3968
3969     static int i, b2, i2, i3, ts;
3970
3971 /*  COMPUTES Y = SIN(X) IF IS.NE.0, Y = COS(X) IF IS.EQ.0,
3972  *  USING TAYLOR SERIES.   ASSUMES ABS(X) .LE. 1.
3973  *  X AND Y ARE MP NUMBERS, IS AN INTEGER.
3974  *  TIME IS O(M(T)T/LOG(T)).   THIS COULD BE REDUCED TO
3975  *  O(SQRT(T)M(T)) AS IN MPEXP1, BUT NOT WORTHWHILE UNLESS
3976  *  T IS VERY LARGE.  ASYMPTOTICALLY FASTER METHODS ARE
3977  *  DESCRIBED IN THE REFERENCES GIVEN IN COMMENTS
3978  *  TO MPATAN AND MPPIGL.
3979  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 3T+8
3980  *  CHECK LEGALITY OF B, T, M AND MXR
3981  */
3982
3983     --y;               /* Parameter adjustments */
3984     --x;
3985
3986     mpchk(&c__3, &c__8);
3987     if (x[1] != 0) goto L20;
3988
3989 /* SIN(0) = 0, COS(0) = 1 */
3990
3991 L10:
3992     y[1] = 0;
3993     if (*is == 0) mpcim(&c__1, &y[1]);
3994     return;
3995
3996 L20:
3997     i2 = MP.t + 5;
3998     i3 = i2 + MP.t + 2;
3999     b2 = max(MP.b,64) << 1;
4000     mpmul(&x[1], &x[1], &MP.r[i3 - 1]);
4001     if (mpcmpi(&MP.r[i3 - 1], &c__1) <= 0) goto L40;
4002
4003     if (v->MPerrors) {
4004         FPRINTF(stderr, "*** ABS(X) .GT. 1 IN CALL TO MPSIN1 ***\n");
4005     }
4006
4007     mperr();
4008     goto L10;
4009
4010 L40:
4011     if (*is == 0) mpcim(&c__1, &MP.r[i2 - 1]);
4012     if (*is != 0) mpstr(&x[1], &MP.r[i2 - 1]);
4013
4014     y[1] = 0;
4015     i = 1;
4016     ts = MP.t;
4017     if (*is == 0) goto L50;
4018
4019     mpstr(&MP.r[i2 - 1], &y[1]);
4020     i = 2;
4021
4022 /* POWER SERIES LOOP.  REDUCE T IF POSSIBLE */
4023
4024 L50:
4025     MP.t = MP.r[i2] + ts + 2;
4026     if (MP.t <= 2) goto L80;
4027
4028     MP.t = min(MP.t,ts);
4029
4030 /* PUT R(I3) FIRST IN CASE ITS DIGITS ARE MAINLY ZERO */
4031
4032     mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i2 - 1]);
4033
4034 /*  IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
4035  *  DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
4036  */
4037
4038     if (i > b2) goto L60;
4039
4040     i__1 = -i * (i + 1);
4041     mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]);
4042     goto L70;
4043
4044 L60:
4045     i__1 = -i;
4046     mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]);
4047     i__1 = i + 1;
4048     mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]);
4049
4050 L70:
4051     i += 2;
4052     MP.t = ts;
4053     mpadd2(&MP.r[i2 - 1], &y[1], &y[1], &y[1], &c__0);
4054     if (MP.r[i2 - 1] != 0) goto L50;
4055
4056 L80:
4057     MP.t = ts;
4058     if (*is == 0) mpaddi(&y[1], &c__1, &y[1]);
4059 }
4060
4061
4062 void
4063 mpsinh(int *x, int *y)
4064 {
4065     int i__1;
4066
4067     static int i2, i3, xs;
4068
4069 /*  RETURNS Y = SINH(X) FOR MP NUMBERS X AND Y, X NOT TOO LARGE.
4070  *  METHOD IS TO USE MPEXP OR MPEXP1, SPACE = 5T+12
4071  *  SAVE SIGN OF X AND CHECK FOR ZERO, SINH(0) = 0
4072  */
4073
4074     --y;               /* Parameter adjustments */
4075     --x;
4076
4077     xs = x[1];
4078     if (xs != 0) goto L10;
4079
4080     y[1] = 0;
4081     return;
4082
4083 /* CHECK LEGALITY OF B, T, M AND MXR */
4084
4085 L10:
4086     mpchk(&c__5, &c__12);
4087     i3 = (MP.t << 2) + 11;
4088
4089 /* WORK WITH ABS(X) */
4090
4091     mpabs(&x[1], &MP.r[i3 - 1]);
4092     if (MP.r[i3] <= 0) goto L20;
4093
4094 /*  HERE ABS(X) .GE. 1, IF TOO LARGE MPEXP GIVES ERROR MESSAGE
4095  *  INCREASE M TO AVOID OVERFLOW IF SINH(X) REPRESENTABLE
4096  */
4097
4098     MP.m += 2;
4099     mpexp(&MP.r[i3 - 1], &MP.r[i3 - 1]);
4100     mprec(&MP.r[i3 - 1], &y[1]);
4101     mpsub(&MP.r[i3 - 1], &y[1], &y[1]);
4102
4103 /*  RESTORE M.  IF RESULT OVERFLOWS OR UNDERFLOWS, MPDIVI AT
4104  *  STATEMENT 30 WILL ACT ACCORDINGLY.
4105  */
4106
4107     MP.m += -2;
4108     goto L30;
4109
4110 /* HERE ABS(X) .LT. 1 SO USE MPEXP1 TO AVOID CANCELLATION */
4111
4112 L20:
4113     i2 = i3 - (MP.t + 2);
4114     mpexp1(&MP.r[i3 - 1], &MP.r[i2 - 1]);
4115     mpaddi(&MP.r[i2 - 1], &c__2, &MP.r[i3 - 1]);
4116     mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &y[1]);
4117     mpaddi(&MP.r[i2 - 1], &c__1, &MP.r[i3 - 1]);
4118     mpdiv(&y[1], &MP.r[i3 - 1], &y[1]);
4119
4120 /* DIVIDE BY TWO AND RESTORE SIGN */
4121
4122 L30:
4123     i__1 = xs << 1;
4124     mpdivi(&y[1], &i__1, &y[1]);
4125 }
4126
4127
4128 void
4129 mpsqrt(int *x, int *y)
4130 {
4131     static int i, i2, iy3;
4132
4133 /*  RETURNS Y = SQRT(X), USING SUBROUTINE MPROOT IF X .GT. 0.
4134  *  DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 4T+10
4135  *  (BUT Y(1) MAY BE R(3T+9)).  X AND Y ARE MP NUMBERS.
4136  *  CHECK LEGALITY OF B, T, M AND MXR
4137  */
4138
4139     --y;           /* Parameter adjustments */
4140     --x;
4141
4142     mpchk(&c__4, &c__10);
4143
4144 /* MPROOT NEEDS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
4145
4146     i2 = MP.t * 3 + 9;
4147     if (x[1] < 0)  goto L10;
4148     else if (x[1] == 0) goto L30;
4149     else goto L40;
4150
4151 L10:
4152     if (v->MPerrors) {
4153         FPRINTF(stderr, "*** X NEGATIVE IN CALL TO SUBROUTINE MPSQRT ***\n");
4154     }
4155
4156     mperr();
4157
4158 L30:
4159     y[1] = 0;
4160     return;
4161
4162 L40:
4163     mproot(&x[1], &c_n2, &MP.r[i2 - 1]);
4164     i = MP.r[i2 + 1];
4165     mpmul(&x[1], &MP.r[i2 - 1], &y[1]);
4166     iy3 = y[3];
4167     mpext(&i, &iy3, &y[1]);
4168 }
4169
4170
4171 void
4172 mpstr(int *x, int *y)
4173 {
4174     int i__1;
4175
4176     static int i, j, t2;
4177
4178 /*  SETS Y = X FOR MP X AND Y.
4179  *  SEE IF X AND Y HAVE THE SAME ADDRESS (THEY OFTEN DO)
4180  */
4181
4182     --y;                /* Parameter adjustments */
4183     --x;
4184
4185     j = x[1];
4186     y[1] = j + 1;
4187     if (j == x[1]) goto L10;
4188
4189 /* HERE X(1) AND Y(1) MUST HAVE THE SAME ADDRESS */
4190
4191     x[1] = j;
4192     return;
4193
4194 /* HERE X(1) AND Y(1) HAVE DIFFERENT ADDRESSES */
4195
4196 L10:
4197     y[1] = j;
4198
4199 /* NO NEED TO MOVE X(2), ... IF X(1) = 0 */
4200
4201     if (j == 0) return;
4202
4203     t2 = MP.t + 2;
4204     i__1 = t2;
4205     for (i = 2; i <= i__1; ++i) y[i] = x[i];
4206 }
4207
4208
4209 void
4210 mpsub(int *x, int *y, int *z)
4211 {
4212     static int y1[1];
4213
4214 /*  SUBTRACTS Y FROM X, FORMING RESULT IN Z, FOR MP X, Y AND Z.
4215  *  FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING
4216  */
4217
4218     --z;                /* Parameter adjustments */
4219     --y;
4220     --x;
4221
4222     y1[0] = -y[1];
4223     mpadd2(&x[1], &y[1], &z[1], y1, &c__0);
4224 }
4225
4226
4227 void
4228 mptanh(int *x, int *y)
4229 {
4230     float r__1;
4231
4232     static int i2, xs;
4233
4234 /*  RETURNS Y = TANH(X) FOR MP NUMBERS X AND Y,
4235  *  USING MPEXP OR MPEXP1, SPACE = 5T+12
4236  */
4237
4238     --y;             /* Parameter adjustments */
4239     --x;
4240
4241     if (x[1] != 0) goto L10;
4242
4243 /* TANH(0) = 0 */
4244
4245     y[1] = 0;
4246     return;
4247
4248 /* CHECK LEGALITY OF B, T, M AND MXR */
4249
4250 L10:
4251     mpchk(&c__5, &c__12);
4252     i2 = (MP.t << 2) + 11;
4253
4254 /* SAVE SIGN AND WORK WITH ABS(X) */
4255
4256     xs = x[1];
4257     mpabs(&x[1], &MP.r[i2 - 1]);
4258
4259 /* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */
4260
4261     r__1 = (float) MP.t * (float).5 * log((float) MP.b);
4262     mpcrm(&r__1, &y[1]);
4263     if (mpcomp(&MP.r[i2 - 1], &y[1]) <= 0) goto L20;
4264
4265 /* HERE ABS(X) IS VERY LARGE */
4266
4267     mpcim(&xs, &y[1]);
4268     return;
4269
4270 /* HERE ABS(X) NOT SO LARGE */
4271
4272 L20:
4273     mpmuli(&MP.r[i2 - 1], &c__2, &MP.r[i2 - 1]);
4274     if (MP.r[i2] <= 0) goto L30;
4275
4276 /* HERE ABS(X) .GE. 1/2 SO USE MPEXP */
4277
4278     mpexp(&MP.r[i2 - 1], &MP.r[i2 - 1]);
4279     mpaddi(&MP.r[i2 - 1], &c_n1, &y[1]);
4280     mpaddi(&MP.r[i2 - 1], &c__1, &MP.r[i2 - 1]);
4281     mpdiv(&y[1], &MP.r[i2 - 1], &y[1]);
4282     goto L40;
4283
4284 /* HERE ABS(X) .LT. 1/2, SO USE MPEXP1 TO AVOID CANCELLATION */
4285
4286 L30:
4287     mpexp1(&MP.r[i2 - 1], &MP.r[i2 - 1]);
4288     mpaddi(&MP.r[i2 - 1], &c__2, &y[1]);
4289     mpdiv(&MP.r[i2 - 1], &y[1], &y[1]);
4290
4291 /* RESTORE SIGN */
4292
4293 L40:
4294     y[1] = xs * y[1];
4295 }
4296
4297
4298 static void
4299 mpunfl(int *x)
4300 {
4301
4302 /*  CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE
4303  *  EXPONENT OF MP NUMBER X WOULD BE LESS THAN -M.
4304  *  SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
4305  */
4306
4307     --x;             /* Parameter adjustments */
4308
4309     mpchk(&c__1, &c__4);
4310
4311 /*  THE UNDERFLOWING NUMBER IS SET TO ZERO
4312  *  AN ALTERNATIVE WOULD BE TO CALL MPMINR (X) AND RETURN,
4313  *  POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION
4314  *  AFTER A PRESET NUMBER OF UNDERFLOWS.  ACTION COULD EASILY
4315  *  BE DETERMINED BY A FLAG IN LABELLED COMMON.
4316  */
4317
4318     x[1] = 0;
4319 }
4320
4321
4322 static double
4323 mppow_di(double *ap, int *bp)
4324 {
4325     double pow, x;
4326     int n;
4327
4328     pow = 1;
4329     x   = *ap;
4330     n   = *bp;
4331
4332     if (n != 0) {
4333         if (n < 0) {
4334             if (x == 0) return(pow);
4335             n = -n;
4336             x = 1/x;
4337         }
4338         for (;;) {
4339             if (n & 01) pow *= x;
4340             if (n >>= 1) x *= x;
4341             else break;
4342         }
4343     }
4344
4345     return(pow);
4346 }
4347
4348
4349 static int
4350 pow_ii(int *ap, int *bp)
4351 {
4352     int pow, x, n;
4353
4354     pow = 1;
4355     x = *ap;
4356     n = *bp;
4357
4358     if (n > 0) {
4359         for (;;) {
4360             if (n & 01) pow *= x;
4361             if (n >>= 1) x *= x;
4362             else break;
4363         }
4364     }
4365
4366     return(pow);
4367 }
4368
4369
4370 static double
4371 mppow_ri(float *ap, int *bp)
4372 {
4373     double pow, x;
4374     int n;
4375
4376     pow = 1;
4377     x   = *ap;
4378     n   = *bp;
4379
4380     if (n != 0) {
4381         if (n < 0) {
4382             if (x == 0) return(pow);
4383             n = -n;
4384             x = 1/x;
4385         }
4386         for (;;) {
4387             if (n & 01)  pow *= x;
4388             if (n >>= 1) x *= x;
4389             else break;
4390         }
4391     }
4392
4393     return(pow);
4394 }