gcalctool/mp.c

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