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