cde/programs/dtcalc/mp.c

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