2 * CDE - Common Desktop Environment
4 * Copyright (c) 1993-2012, The Open Group. All rights reserved.
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)
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
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with these librararies and programs; if not, write
20 * to the Free Software Foundation, Inc., 51 Franklin Street, Fifth
21 * Floor, Boston, MA 02110-1301 USA
23 /* $XConsortium: mp.c /main/3 1995/11/01 12:42:08 rswiston $ */
26 * Contains the actual math functions. *
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. *
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.
38 * It has been converted from FORTRAN into C using the freely available
39 * f2c translator, available via netlib on research.att.com.
41 * The subsequently converted C code has then been tidied up, mainly to
42 * remove any dependencies on the libI77 and libF77 support libraries.
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
58 int b, t, m, mxr, r[MP_SIZE] ;
61 /* Table of constant values */
69 static int c__12 = 12 ;
70 static int c_n2 = -2 ;
71 static int c__10 = 10 ;
72 static int c__32 = 32 ;
75 static int c__14 = 14 ;
76 static int c_n1 = -1 ;
77 static int c__239 = 239 ;
79 static int c__16 = 16 ;
81 extern char *mpstrs[] ; /* MP errors (visible with -D option). */
82 extern char *vstrs[] ; /* Various strings. */
84 extern Vars v ; /* Calctool variables and options. */
92 /* SETS Y = ABS(X) FOR MP NUMBERS X AND Y */
94 --y ; /* Parameter adjustments */
108 /* ADDS X AND Y, FORMING RESULT IN Z, WHERE X, Y AND Z ARE MP
109 * NUMBERS. FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING.
112 --z ; /* Parameter adjustments */
116 mpadd2(&x[1], &y[1], &z[1], &y[1], &c__0) ;
122 mpadd2(x, y, z, y1, trunc)
123 int *x, *y, *z, *y1, *trunc ;
128 static int ed, re, rs, med ;
130 /* CALLED BY MPADD, MPSUB ETC.
131 * X, Y AND Z ARE MP NUMBERS, Y1 AND TRUNC ARE INTEGERS.
132 * TO FORCE CALL BY REFERENCE RATHER THAN VALUE/RESULT, Y1 IS
133 * DECLARED AS AN ARRAY, BUT ONLY Y1(1) IS EVER USED.
134 * SETS Z = X + Y1(1)*ABS(Y), WHERE Y1(1) = +- Y(1).
135 * IF TRUNC.EQ.0 R*-ROUNDING IS USED, OTHERWISE TRUNCATION.
136 * R*-ROUNDING IS DEFINED IN KUKI AND CODI, COMM. ACM
137 * 16(1973), 223. (SEE ALSO BRENT, IEEE TC-22(1973), 601.)
138 * CHECK FOR X OR Y ZERO
141 --y1 ; /* Parameter adjustments */
146 if (x[1] != 0) goto L20 ;
148 /* X = 0 OR NEGLIGIBLE, SO RESULT = +-Y */
151 mpstr(&y[1], &z[1]) ;
156 if (y1[1] != 0) goto L40 ;
158 /* Y = 0 OR NEGLIGIBLE, SO RESULT = X */
161 mpstr(&x[1], &z[1]) ;
168 if (C_abs(s) <= 1) goto L60 ;
170 mpchk(&c__1, &c__4) ;
173 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ADD2A]) ;
174 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ADD2B]) ;
181 /* COMPARE EXPONENTS */
186 if (ed < 0) goto L90 ;
187 else if (ed == 0) goto L70 ;
190 /* EXPONENTS EQUAL SO COMPARE SIGNS, THEN FRACTIONS IF NEC. */
193 if (s > 0) goto L100 ;
196 for (j = 1; j <= i__1; ++j)
198 if ((i__2 = x[j + 2] - y[j + 2]) < 0) goto L100 ;
199 else if (i__2 == 0) continue ;
208 /* HERE EXPONENT(Y) .GE. EXPONENT(X) */
211 if (med > MP.t) goto L10 ;
216 mpadd3(&x[1], &y[1], &s, &med, &re) ;
218 /* NORMALIZE, ROUND OR TRUNCATE, AND RETURN */
221 mpnzr(&rs, &re, &z[1], trunc) ;
224 /* ABS(X) .GT. ABS(Y) */
227 if (med > MP.t) goto L30 ;
232 mpadd3(&y[1], &x[1], &s, &med, &re) ;
238 mpadd3(x, y, s, med, re)
239 int *x, *y, *s, *med, *re ;
243 static int c, i, j, i2, i2p, ted ;
245 /* CALLED BY MPADD2, DOES INNER LOOPS OF ADDITION */
247 --y ; /* Parameter adjustments */
255 /* CLEAR GUARD DIGITS TO RIGHT OF X DIGITS */
258 if (i <= ted) goto L20 ;
265 if (*s < 0) goto L130 ;
267 /* HERE DO ADDITION, EXPONENT(Y) .GE. EXPONENT(X) */
269 if (i <= MP.t) goto L40 ;
273 MP.r[i - 1] = x[j + 2] ;
275 if (i > MP.t) goto L30 ;
278 if (i <= *med) goto L60 ;
281 c = y[i + 2] + x[j + 2] + c ;
282 if (c < MP.b) goto L50 ;
284 /* CARRY GENERATED HERE */
286 MP.r[i - 1] = c - MP.b ;
291 /* NO CARRY GENERATED HERE */
300 if (i <= 0) goto L90 ;
303 if (c < MP.b) goto L70 ;
314 /* NO CARRY POSSIBLE HERE */
319 MP.r[i - 1] = y[i + 2] ;
326 /* MUST SHIFT RIGHT HERE AS CARRY OFF END */
330 for (j = 2 ; j <= i__1 ; ++j)
333 MP.r[i] = MP.r[i - 1] ;
339 /* HERE DO SUBTRACTION, ABS(Y) .GT. ABS(X) */
343 MP.r[i - 1] = c - x[j + 2] ;
345 if (MP.r[i - 1] >= 0) goto L120 ;
347 /* BORROW GENERATED HERE */
350 MP.r[i - 1] += MP.b ;
356 if (i > MP.t) goto L110 ;
359 if (i <= *med) goto L160 ;
362 c = y[i + 2] + c - x[j + 2] ;
363 if (c >= 0) goto L150 ;
365 /* BORROW GENERATED HERE */
367 MP.r[i - 1] = c + MP.b ;
372 /* NO BORROW GENERATED HERE */
384 if (c >= 0) goto L70 ;
386 MP.r[i - 1] = c + MP.b ;
398 /* ADDS MULTIPLE-PRECISION X TO INTEGER IY
399 * GIVING MULTIPLE-PRECISION Z.
400 * DIMENSION OF R IN CALLING PROGRAM MUST BE
401 * AT LEAST 2T+6 (BUT Z(1) MAY BE R(T+5)).
402 * DIMENSION R(6) BECAUSE RALPH COMPILER ON UNIVAC 1100 COMPUTERS
403 * OBJECTS TO DIMENSION R(1).
404 * CHECK LEGALITY OF B, T, M AND MXR
407 --z ; /* Parameter adjustments */
410 mpchk(&c__2, &c__6) ;
411 mpcim(iy, &MP.r[MP.t + 4]) ;
412 mpadd(&x[1], &MP.r[MP.t + 4], &z[1]) ;
422 /* ADDS THE RATIONAL NUMBER I/J TO MP NUMBER X, MP RESULT IN Y
423 * DIMENSION OF R MUST BE AT LEAST 2T+6
424 * CHECK LEGALITY OF B, T, M AND MXR
427 --y ; /* Parameter adjustments */
430 mpchk(&c__2, &c__6) ;
431 mpcqm(i, j, &MP.r[MP.t + 4]) ;
432 mpadd(&x[1], &MP.r[MP.t + 4], &y[1]) ;
443 static int i, b2, i2, id, ts ;
445 /* COMPUTES MP Y = ARCTAN(1/N), ASSUMING INTEGER N .GT. 1.
446 * USES SERIES ARCTAN(X) = X - X**3/3 + X**5/5 - ...
447 * DIMENSION OF R IN CALLING PROGRAM MUST BE
449 * CHECK LEGALITY OF B, T, M AND MXR
452 --y ; /* Parameter adjustments */
454 mpchk(&c__2, &c__6) ;
455 if (*n > 1) goto L20 ;
457 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_PART1]) ;
467 /* SET SUM TO X = 1/N */
469 mpcqm(&c__1, n, &y[1]) ;
471 /* SET ADDITIVE TERM TO X */
473 mpstr(&y[1], &MP.r[i2 - 1]) ;
477 /* ASSUME AT LEAST 16-BIT WORD. */
480 if (*n < b2) id = b2 * 7 * b2 / (*n * *n) ;
482 /* MAIN LOOP. FIRST REDUCE T IF POSSIBLE */
485 MP.t = ts + 2 + MP.r[i2] - y[2] ;
486 if (MP.t < 2) goto L60 ;
488 MP.t = min(MP.t,ts) ;
490 /* IF (I+2)*N**2 IS NOT REPRESENTABLE AS AN INTEGER THE DIVISION
491 * FOLLOWING HAS TO BE PERFORMED IN SEVERAL STEPS.
494 if (i >= id) goto L40 ;
497 i__2 = (i + 2) * *n * *n ;
498 mpmulq(&MP.r[i2 - 1], &i__1, &i__2, &MP.r[i2 - 1]) ;
504 mpmulq(&MP.r[i2 - 1], &i__1, &i__2, &MP.r[i2 - 1]) ;
505 mpdivi(&MP.r[i2 - 1], n, &MP.r[i2 - 1]) ;
506 mpdivi(&MP.r[i2 - 1], n, &MP.r[i2 - 1]) ;
515 /* ADD TO SUM, USING MPADD2 (FASTER THAN MPADD) */
517 mpadd2(&MP.r[i2 - 1], &y[1], &y[1], &y[1], &c__0) ;
518 if (MP.r[i2 - 1] != 0) goto L30 ;
534 /* RETURNS Y = ARCSIN(X), ASSUMING ABS(X) .LE. 1,
535 * FOR MP NUMBERS X AND Y.
536 * Y IS IN THE RANGE -PI/2 TO +PI/2.
537 * METHOD IS TO USE MPATAN, SO TIME IS O(M(T)T).
538 * DIMENSION OF R MUST BE AT LEAST 5T+12
539 * CHECK LEGALITY OF B, T, M AND MXR
542 --y ; /* Parameter adjustments */
545 mpchk(&c__5, &c__12) ;
546 i3 = (MP.t << 2) + 11 ;
547 if (x[1] == 0) goto L30 ;
549 if (x[2] <= 0) goto L40 ;
551 /* HERE ABS(X) .GE. 1. SEE IF X = +-1 */
553 mpcim(&x[1], &MP.r[i3 - 1]) ;
554 if (mpcomp(&x[1], &MP.r[i3 - 1]) != 0) goto L10 ;
556 /* X = +-1 SO RETURN +-PI/2 */
559 i__1 = MP.r[i3 - 1] << 1 ;
560 mpdivi(&y[1], &i__1, &y[1]) ;
564 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ASIN]) ;
572 /* HERE ABS(X) .LT. 1 SO USE ARCTAN(X/SQRT(1 - X**2)) */
575 i2 = i3 - (MP.t + 2) ;
576 mpcim(&c__1, &MP.r[i2 - 1]) ;
577 mpstr(&MP.r[i2 - 1], &MP.r[i3 - 1]) ;
578 mpsub(&MP.r[i2 - 1], &x[1], &MP.r[i2 - 1]) ;
579 mpadd(&MP.r[i3 - 1], &x[1], &MP.r[i3 - 1]) ;
580 mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]) ;
581 mproot(&MP.r[i3 - 1], &c_n2, &MP.r[i3 - 1]) ;
582 mpmul(&x[1], &MP.r[i3 - 1], &y[1]) ;
583 mpatan(&y[1], &y[1]) ;
595 static int i, q, i2, i3, ie, ts ;
596 static float rx, ry ;
598 /* RETURNS Y = ARCTAN(X) FOR MP X AND Y, USING AN O(T.M(T)) METHOD
599 * WHICH COULD EASILY BE MODIFIED TO AN O(SQRT(T)M(T))
600 * METHOD (AS IN MPEXP1). Y IS IN THE RANGE -PI/2 TO +PI/2.
601 * FOR AN ASYMPTOTICALLY FASTER METHOD, SEE - FAST MULTIPLE-
602 * PRECISION EVALUATION OF ELEMENTARY FUNCTIONS
603 * (BY R. P. BRENT), J. ACM 23 (1976), 242-251,
604 * AND THE COMMENTS IN MPPIGL.
605 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 5T+12
606 * CHECK LEGALITY OF B, T, M AND MXR
609 --y ; /* Parameter adjustments */
612 mpchk(&c__5, &c__12) ;
622 mpstr(&x[1], &MP.r[i3 - 1]) ;
624 if (ie <= 2) mpcmr(&x[1], &rx) ;
628 /* REDUCE ARGUMENT IF NECESSARY BEFORE USING SERIES */
631 if (MP.r[i3] < 0) goto L30 ;
633 if (MP.r[i3] == 0 && MP.r[i3 + 1] + 1 << 1 <= MP.b)
637 mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], &y[1]) ;
638 mpaddi(&y[1], &c__1, &y[1]) ;
639 mpsqrt(&y[1], &y[1]) ;
640 mpaddi(&y[1], &c__1, &y[1]) ;
641 mpdiv(&MP.r[i3 - 1], &y[1], &MP.r[i3 - 1]) ;
644 /* USE POWER SERIES NOW ARGUMENT IN (-0.5, 0.5) */
647 mpstr(&MP.r[i3 - 1], &y[1]) ;
648 mpmul(&MP.r[i3 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]) ;
652 /* SERIES LOOP. REDUCE T IF POSSIBLE. */
655 MP.t = ts + 2 + MP.r[i3] ;
656 if (MP.t <= 2) goto L50 ;
658 MP.t = min(MP.t,ts) ;
659 mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i3 - 1]) ;
662 mpmulq(&MP.r[i3 - 1], &i__1, &i__2, &MP.r[i3 - 1]) ;
665 mpadd(&y[1], &MP.r[i3 - 1], &y[1]) ;
666 if (MP.r[i3 - 1] != 0) goto L40 ;
668 /* RESTORE T, CORRECT FOR ARGUMENT REDUCTION, AND EXIT */
672 mpmuli(&y[1], &q, &y[1]) ;
674 /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS EXPONENT
675 * OF X IS LARGE (WHEN ATAN MIGHT NOT WORK)
681 if ((r__1 = ry - atan(rx), dabs(r__1)) < dabs(ry) * (float).01)
684 /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL. */
686 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ATAN]) ;
700 static int i, k, i2, ib, ie, re, tp, rs ;
701 static double db, dj ;
703 /* CONVERTS DOUBLE-PRECISION NUMBER DX TO MULTIPLE-PRECISION Z.
704 * SOME NUMBERS WILL NOT CONVERT EXACTLY ON MACHINES
705 * WITH BASE OTHER THAN TWO, FOUR OR SIXTEEN.
706 * THIS ROUTINE IS NOT CALLED BY ANY OTHER ROUTINE IN MP,
707 * SO MAY BE OMITTED IF DOUBLE-PRECISION IS NOT AVAILABLE.
708 * CHECK LEGALITY OF B, T, M AND MXR
711 --z ; /* Parameter adjustments */
713 mpchk(&c__1, &c__4) ;
718 if (*dx < 0.) goto L20 ;
719 else if (*dx == 0) goto L10 ;
722 /* IF DX = 0D0 RETURN 0 */
745 if (dj < 1.) goto L60 ;
747 /* INCREASE IE AND DIVIDE DJ BY 16. */
754 if (dj >= .0625) goto L70 ;
760 /* NOW DJ IS DY DIVIDED BY SUITABLE POWER OF 16
767 /* DB = DFLOAT(B) IS NOT ANSI STANDARD SO USE FLOAT AND DBLE */
769 db = (double) ((float) MP.b) ;
771 /* CONVERSION LOOP (ASSUME DOUBLE-PRECISION OPS. EXACT) */
774 for (i = 1; i <= i__1; ++i)
777 MP.r[i - 1] = (int) dj ;
778 dj -= (double) ((float) MP.r[i - 1]) ;
781 /* NORMALIZE RESULT */
783 mpnzr(&rs, &re, &z[1], &c__0) ;
787 i__1 = MP.b * 7 * MP.b ;
788 ib = max(i__1, 32767) / 16 ;
791 /* NOW MULTIPLY BY 16**IE */
793 if (ie < 0) goto L90 ;
794 else if (ie == 0) goto L130 ;
800 for (i = 1; i <= i__1; ++i)
803 if (tp <= ib && tp != MP.b && i < k) continue ;
804 mpdivi(&z[1], &tp, &z[1]) ;
811 for (i = 1; i <= i__1; ++i)
814 if (tp <= ib && tp != MP.b && i < ie) continue ;
815 mpmuli(&z[1], &tp, &z[1]) ;
830 /* CHECKS LEGALITY OF B, T, M AND MXR.
831 * THE CONDITION ON MXR (THE DIMENSION OF R IN COMMON) IS THAT
835 /* CHECK LEGALITY OF B, T AND M */
837 if (MP.b > 1) goto L40 ;
841 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKC]) ;
842 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKD]) ;
847 if (MP.t > 1) goto L60 ;
851 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKE]) ;
852 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKF]) ;
857 if (MP.m > MP.t) goto L80 ;
861 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKG]) ;
862 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKH]) ;
866 /* 8*B*B-1 SHOULD BE REPRESENTABLE, IF NOT WILL OVERFLOW
867 * AND MAY BECOME NEGATIVE, SO CHECK FOR THIS
871 ib = (MP.b << 2) * MP.b - 1 ;
872 if (ib > 0 && (ib << 1) + 1 > 0) goto L100 ;
875 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKI]) ;
879 /* CHECK THAT SPACE IN COMMON IS SUFFICIENT */
882 mx = *i * MP.t + *j ;
883 if (MP.mxr >= mx) return ;
885 /* HERE COMMON IS TOO SMALL, SO GIVE ERROR MESSAGE. */
891 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKJ]) ;
892 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CHKL]) ;
893 msg = (char *) XtMalloc(strlen( mpstrs[(int) MP_CHKM]) + 200);
894 sprintf(msg, mpstrs[(int) MP_CHKM], *i, *j, mx) ;
895 _DtSimpleError (v->appname, DtWarning, NULL, msg);
896 sprintf(msg, mpstrs[(int) MP_CHKN], MP.mxr, MP.t) ;
897 _DtSimpleError (v->appname, DtWarning, NULL, msg);
914 /* CONVERTS INTEGER IX TO MULTIPLE-PRECISION Z.
915 * CHECK LEGALITY OF B, T, M AND MXR
918 --z ; /* Parameter adjustments */
920 mpchk(&c__1, &c__4) ;
922 if (n < 0) goto L20 ;
923 else if (n == 0) goto L10 ;
938 /* SET EXPONENT TO T */
946 for (i = 2; i <= i__1; ++i) z[i + 1] = 0 ;
952 /* NORMALIZE BY CALLING MPMUL2 */
954 mpmul2(&z[1], &c__1, &z[1], &c__1) ;
968 static double db, dz2 ;
970 /* CONVERTS MULTIPLE-PRECISION X TO DOUBLE-PRECISION DZ.
971 * ASSUMES X IS IN ALLOWABLE RANGE FOR DOUBLE-PRECISION
972 * NUMBERS. THERE IS SOME LOSS OF ACCURACY IF THE
974 * CHECK LEGALITY OF B, T, M AND MXR
977 --x ; /* Parameter adjustments */
979 mpchk(&c__1, &c__4) ;
981 if (x[1] == 0) return ;
983 /* DB = DFLOAT(B) IS NOT ANSI STANDARD, SO USE FLOAT AND DBLE */
985 db = (double) ((float) MP.b) ;
987 for (i = 1; i <= i__1; ++i)
989 *dz = db * *dz + (double) ((float) x[i + 2]) ;
992 /* CHECK IF FULL DOUBLE-PRECISION ACCURACY ATTAINED */
996 /* TEST BELOW NOT ALWAYS EQUIVALENT TO - IF (DZ2.LE.DZ) GO TO 20,
997 * FOR EXAMPLE ON CYBER 76.
999 if (dz2 - *dz <= 0.) goto L20 ;
1002 /* NOW ALLOW FOR EXPONENT */
1006 *dz *= mppow_di(&db, &i__1) ;
1008 /* CHECK REASONABLENESS OF RESULT. */
1010 if (*dz <= 0.) goto L30 ;
1012 /* LHS SHOULD BE .LE. 0.5 BUT ALLOW FOR SOME ERROR IN DLOG */
1014 if ((d__1 = (double) ((float) x[2]) - (log(*dz) / log((double)
1015 ((float) MP.b)) + .5), C_abs(d__1)) > .6)
1018 if (x[1] < 0) *dz = -(*dz) ;
1021 /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
1022 * TRY USING MPCMDE INSTEAD.
1026 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CMD]) ;
1040 static int x2, il, ip, xs ;
1042 /* FOR MP X AND Y, RETURNS FRACTIONAL PART OF X IN Y,
1043 * I.E., Y = X - INTEGER PART OF X (TRUNCATED TOWARDS 0).
1046 --y ; /* Parameter adjustments */
1049 if (x[1] != 0) goto L20 ;
1051 /* RETURN 0 IF X = 0 */
1060 /* RETURN 0 IF EXPONENT SO LARGE THAT NO FRACTIONAL PART */
1062 if (x2 >= MP.t) goto L10 ;
1064 /* IF EXPONENT NOT POSITIVE CAN RETURN X */
1066 if (x2 > 0) goto L30 ;
1068 mpstr(&x[1], &y[1]) ;
1071 /* CLEAR ACCUMULATOR */
1075 for (i = 1; i <= i__1; ++i) MP.r[i - 1] = 0 ;
1079 /* MOVE FRACTIONAL PART OF X TO ACCUMULATOR */
1082 for (i = il; i <= i__1; ++i) MP.r[i - 1] = x[i + 2] ;
1084 for (i = 1; i <= 4; ++i)
1091 /* NORMALIZE RESULT AND RETURN */
1093 mpnzr(&xs, &x2, &y[1], &c__1) ;
1104 static int i, j, k, j1, x2, kx, xs, izs ;
1106 /* CONVERTS MULTIPLE-PRECISION X TO INTEGER IZ,
1107 * ASSUMING THAT X NOT TOO LARGE (ELSE USE MPCMIM).
1108 * X IS TRUNCATED TOWARDS ZERO.
1109 * IF INT(X)IS TOO LARGE TO BE REPRESENTED AS A SINGLE-
1110 * PRECISION INTEGER, IZ IS RETURNED AS ZERO. THE USER
1111 * MAY CHECK FOR THIS POSSIBILITY BY TESTING IF
1112 * ((X(1).NE.0).AND.(X(2).GT.0).AND.(IZ.EQ.0)) IS TRUE ON
1113 * RETURN FROM MPCMI.
1116 --x ; /* Parameter adjustments */
1120 if (xs == 0) return ;
1122 if (x[2] <= 0) return ;
1126 for (i = 1; i <= i__1; ++i)
1130 if (i <= MP.t) *iz += x[i + 2] ;
1132 /* CHECK FOR SIGNS OF INTEGER OVERFLOW */
1134 if (*iz <= 0 || *iz <= izs) goto L30 ;
1137 /* CHECK THAT RESULT IS CORRECT (AN UNDETECTED OVERFLOW MAY
1143 for (i = 1; i <= i__1; ++i)
1148 if (k <= MP.t) kx = x[k + 2] ;
1149 if (kx != j - MP.b * j1) goto L30 ;
1152 if (j != 0) goto L30 ;
1154 /* RESULT CORRECT SO RESTORE SIGN AND RETURN */
1159 /* HERE OVERFLOW OCCURRED (OR X WAS UNNORMALIZED), SO
1177 /* RETURNS Y = INTEGER PART OF X (TRUNCATED TOWARDS 0), FOR MP X AND Y.
1178 * USE IF Y TOO LARGE TO BE REPRESENTABLE AS A SINGLE-PRECISION INTEGER.
1179 * (ELSE COULD USE MPCMI).
1180 * CHECK LEGALITY OF B, T, M AND MXR
1183 --y ; /* Parameter adjustments */
1186 mpchk(&c__1, &c__4) ;
1187 mpstr(&x[1], &y[1]) ;
1188 if (y[1] == 0) return ;
1192 /* IF EXPONENT LARGE ENOUGH RETURN Y = X */
1194 if (il > MP.t) return ;
1196 /* IF EXPONENT SMALL ENOUGH RETURN ZERO */
1198 if (il > 1) goto L10 ;
1203 /* SET FRACTION TO ZERO */
1207 for (i = il; i <= i__1; ++i) y[i + 2] = 0 ;
1218 /* COMPARES MP NUMBER X WITH INTEGER I, RETURNING
1222 * DIMENSION OF R IN COMMON AT LEAST 2T+6
1223 * CHECK LEGALITY OF B, T, M AND MXR
1226 --x ; /* Parameter adjustments */
1228 mpchk(&c__2, &c__6) ;
1230 /* CONVERT I TO MULTIPLE-PRECISION AND COMPARE */
1232 mpcim(i, &MP.r[MP.t + 4]) ;
1233 ret_val = mpcomp(&x[1], &MP.r[MP.t + 4]) ;
1245 /* COMPARES MP NUMBER X WITH REAL NUMBER RI, RETURNING
1249 * DIMENSION OF R IN COMMON AT LEAST 2T+6
1250 * CHECK LEGALITY OF B, T, M AND MXR
1253 --x ; /* Parameter adjustments */
1255 mpchk(&c__2, &c__6) ;
1257 /* CONVERT RI TO MULTIPLE-PRECISION AND COMPARE */
1259 mpcrm(ri, &MP.r[MP.t + 4]) ;
1260 ret_val = mpcomp(&x[1], &MP.r[MP.t + 4]) ;
1274 static float rb, rz2 ;
1276 /* CONVERTS MULTIPLE-PRECISION X TO SINGLE-PRECISION RZ.
1277 * ASSUMES X IN ALLOWABLE RANGE. THERE IS SOME LOSS OF
1278 * ACCURACY IF THE EXPONENT IS LARGE.
1279 * CHECK LEGALITY OF B, T, M AND MXR
1282 --x ; /* Parameter adjustments */
1284 mpchk(&c__1, &c__4) ;
1286 if (x[1] == 0) return ;
1290 for (i = 1; i <= i__1; ++i)
1292 *rz = rb * *rz + (float) x[i + 2] ;
1295 /* CHECK IF FULL SINGLE-PRECISION ACCURACY ATTAINED */
1297 rz2 = *rz + (float) 1.0 ;
1298 if (rz2 <= *rz) goto L20 ;
1301 /* NOW ALLOW FOR EXPONENT */
1305 *rz *= mppow_ri(&rb, &i__1) ;
1307 /* CHECK REASONABLENESS OF RESULT */
1309 if (*rz <= (float)0.) goto L30 ;
1311 /* LHS SHOULD BE .LE. 0.5, BUT ALLOW FOR SOME ERROR IN ALOG */
1313 if ((r__1 = (float) x[2] - (log(*rz) / log((float) MP.b) + (float).5),
1314 dabs(r__1)) > (float).6)
1317 if (x[1] < 0) *rz = -(double)(*rz) ;
1320 /* FOLLOWING MESSAGE INDICATES THAT X IS TOO LARGE OR SMALL -
1321 * TRY USING MPCMRE INSTEAD.
1325 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CMR]) ;
1336 int ret_val, i__1, i__2 ;
1340 /* COMPARES THE MULTIPLE-PRECISION NUMBERS X AND Y,
1341 * RETURNING +1 IF X .GT. Y,
1343 * AND 0 IF X .EQ. Y.
1346 --y ; /* Parameter adjustments */
1349 if ((i__1 = x[1] - y[1]) < 0) goto L10 ;
1350 else if (i__1 == 0) goto L30 ;
1365 /* SIGN(X) = SIGN(Y), SEE IF ZERO */
1368 if (x[1] != 0) goto L40 ;
1375 /* HAVE TO COMPARE EXPONENTS AND FRACTIONS */
1380 for (i = 2; i <= i__1; ++i)
1382 if ((i__2 = x[i] - y[i]) < 0) goto L60 ;
1383 else if (i__2 == 0) continue ;
1392 /* ABS(X) .LT. ABS(Y) */
1398 /* ABS(X) .GT. ABS(Y) */
1412 /* RETURNS Y = COS(X) FOR MP X AND Y, USING MPSIN AND MPSIN1.
1413 * DIMENSION OF R IN COMMON AT LEAST 5T+12.
1416 --y ; /* Parameter adjustments */
1419 if (x[1] != 0) goto L10 ;
1423 mpcim(&c__1, &y[1]) ;
1426 /* CHECK LEGALITY OF B, T, M AND MXR */
1429 mpchk(&c__5, &c__12) ;
1430 i2 = MP.t * 3 + 12 ;
1432 /* SEE IF ABS(X) .LE. 1 */
1434 mpabs(&x[1], &y[1]) ;
1435 if (mpcmpi(&y[1], &c__1) <= 0) goto L20 ;
1437 /* HERE ABS(X) .GT. 1 SO USE COS(X) = SIN(PI/2 - ABS(X)),
1438 * COMPUTING PI/2 WITH ONE GUARD DIGIT.
1442 mppi(&MP.r[i2 - 1]) ;
1443 mpdivi(&MP.r[i2 - 1], &c__2, &MP.r[i2 - 1]) ;
1445 mpsub(&MP.r[i2 - 1], &y[1], &y[1]) ;
1446 mpsin(&y[1], &y[1]) ;
1449 /* HERE ABS(X) .LE. 1 SO USE POWER SERIES */
1452 mpsin1(&y[1], &y[1], &c__0) ;
1463 /* RETURNS Y = COSH(X) FOR MP NUMBERS X AND Y, X NOT TOO LARGE.
1464 * USES MPEXP, DIMENSION OF R IN COMMON AT LEAST 5T+12
1467 --y ; /* Parameter adjustments */
1470 if (x[1] != 0) goto L10 ;
1474 mpcim(&c__1, &y[1]) ;
1477 /* CHECK LEGALITY OF B, T, M AND MXR */
1480 mpchk(&c__5, &c__12) ;
1481 i2 = (MP.t << 2) + 11 ;
1482 mpabs(&x[1], &MP.r[i2 - 1]) ;
1484 /* IF ABS(X) TOO LARGE MPEXP WILL PRINT ERROR MESSAGE
1485 * INCREASE M TO AVOID OVERFLOW WHEN COSH(X) REPRESENTABLE
1489 mpexp(&MP.r[i2 - 1], &MP.r[i2 - 1]) ;
1490 mprec(&MP.r[i2 - 1], &y[1]) ;
1491 mpadd(&MP.r[i2 - 1], &y[1], &y[1]) ;
1493 /* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, MPDIVI WILL
1498 mpdivi(&y[1], &c__2, &y[1]) ;
1509 /* CONVERTS THE RATIONAL NUMBER I/J TO MULTIPLE PRECISION Q. */
1511 --q ; /* Parameter adjustments */
1516 if (j1 < 0) goto L30 ;
1517 else if (j1 == 0) goto L10 ;
1521 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_CQM]) ;
1533 if (j1 != 1) mpdivi(&q[1], &j1, &q[1]) ;
1545 static int i, k, i2, ib, ie, re, tp, rs ;
1546 static float rb, rj ;
1548 /* CONVERTS SINGLE-PRECISION NUMBER RX TO MULTIPLE-PRECISION Z.
1549 * SOME NUMBERS WILL NOT CONVERT EXACTLY ON MACHINES
1550 * WITH BASE OTHER THAN TWO, FOUR OR SIXTEEN.
1551 * CHECK LEGALITY OF B, T, M AND MXR
1554 --z ; /* Parameter adjustments */
1556 mpchk(&c__1, &c__4) ;
1561 if (*rx < (float) 0.0) goto L20 ;
1562 else if (*rx == 0) goto L10 ;
1565 /* IF RX = 0E0 RETURN 0 */
1575 rj = -(double)(*rx) ;
1588 if (rj < (float)1.0) goto L60 ;
1590 /* INCREASE IE AND DIVIDE RJ BY 16. */
1593 rj *= (float) 0.0625 ;
1597 if (rj >= (float).0625) goto L70 ;
1603 /* NOW RJ IS DY DIVIDED BY SUITABLE POWER OF 16.
1611 /* CONVERSION LOOP (ASSUME SINGLE-PRECISION OPS. EXACT) */
1614 for (i = 1; i <= i__1; ++i)
1617 MP.r[i - 1] = (int) rj ;
1618 rj -= (float) MP.r[i - 1] ;
1621 /* NORMALIZE RESULT */
1623 mpnzr(&rs, &re, &z[1], &c__0) ;
1627 i__1 = MP.b * 7 * MP.b ;
1628 ib = max(i__1, 32767) / 16 ;
1631 /* NOW MULTIPLY BY 16**IE */
1633 if (ie < 0) goto L90 ;
1634 else if (ie == 0) goto L130 ;
1640 for (i = 1; i <= i__1; ++i)
1643 if (tp <= ib && tp != MP.b && i < k) continue ;
1644 mpdivi(&z[1], &tp, &z[1]) ;
1651 for (i = 1; i <= i__1; ++i)
1654 if (tp <= ib && tp != MP.b && i < ie) continue ;
1655 mpmuli(&z[1], &tp, &z[1]) ;
1668 static int i, i2, ie, iz3 ;
1670 /* SETS Z = X/Y, FOR MP X, Y AND Z.
1671 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 4T+10
1672 * (BUT Z(1) MAY BE R(3T+9)).
1673 * CHECK LEGALITY OF B, T, M AND MXR
1676 --z ; /* Parameter adjustments */
1680 mpchk(&c__4, &c__10) ;
1682 /* CHECK FOR DIVISION BY ZERO */
1684 if (y[1] != 0) goto L20 ;
1686 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_DIVA]) ;
1692 /* SPACE USED BY MPREC IS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
1697 /* CHECK FOR X = 0 */
1699 if (x[1] != 0) goto L30 ;
1704 /* INCREASE M TO AVOID OVERFLOW IN MPREC */
1709 /* FORM RECIPROCAL OF Y */
1711 mprec(&y[1], &MP.r[i2 - 1]) ;
1713 /* SET EXPONENT OF R(I2) TO ZERO TO AVOID OVERFLOW IN MPMUL */
1721 mpmul(&x[1], &MP.r[i2 - 1], &z[1]) ;
1723 mpext(&i, &iz3, &z[1]) ;
1725 /* RESTORE M, CORRECT EXPONENT AND RETURN */
1729 if (z[2] >= -MP.m) goto L40 ;
1731 /* UNDERFLOW HERE */
1737 if (z[2] <= MP.m) return ;
1741 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_DIVB]) ;
1754 static int c, i, j, k, b2, c2, i2, j1, j2, r1 ;
1755 static int j11, kh, re, iq, ir, rs, iqj ;
1757 /* DIVIDES MP X BY THE SINGLE-PRECISION INTEGER IY GIVING MP Z.
1758 * THIS IS MUCH FASTER THAN DIVISION BY AN MP NUMBER.
1761 --z ; /* Parameter adjustments */
1766 if (j < 0) goto L30 ;
1767 else if (j == 0) goto L10 ;
1771 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_DIVIA]) ;
1782 /* CHECK FOR ZERO DIVIDEND */
1784 if (rs == 0) goto L120 ;
1786 /* CHECK FOR DIVISION BY B */
1788 if (j != MP.b) goto L50 ;
1789 mpstr(&x[1], &z[1]) ;
1790 if (re <= -MP.m) goto L240 ;
1795 /* CHECK FOR DIVISION BY 1 OR -1 */
1798 if (j != 1) goto L60 ;
1799 mpstr(&x[1], &z[1]) ;
1808 /* IF J*B NOT REPRESENTABLE AS AN INTEGER HAVE TO SIMULATE
1809 * LONG DIVISION. ASSUME AT LEAST 16-BIT WORD.
1814 i__1 = MP.b << 3, i__2 = 32767 / MP.b ;
1815 b2 = max(i__1,i__2) ;
1816 if (j >= b2) goto L130 ;
1818 /* LOOK FOR FIRST NONZERO DIGIT IN QUOTIENT */
1823 if (i <= MP.t) c += x[i + 2] ;
1825 if (r1 < 0) goto L210 ;
1826 else if (r1 == 0) goto L70 ;
1829 /* ADJUST EXPONENT AND GET T+4 DIGITS IN QUOTIENT */
1834 c = MP.b * (c - j * r1) ;
1836 if (i >= MP.t) goto L100 ;
1839 for (k = 2; k <= i__1; ++k)
1843 MP.r[k - 1] = c / j ;
1844 c = MP.b * (c - j * MP.r[k - 1]) ;
1846 if (c < 0) goto L210 ;
1851 for (k = kh; k <= i__1; ++k)
1853 MP.r[k - 1] = c / j ;
1854 c = MP.b * (c - j * MP.r[k - 1]) ;
1856 if (c < 0) goto L210 ;
1858 /* NORMALIZE AND ROUND RESULT */
1861 mpnzr(&rs, &re, &z[1], &c__0) ;
1864 /* HERE NEED SIMULATED DOUBLE-PRECISION DIVISION */
1869 j2 = j - j1 * MP.b ;
1872 /* LOOK FOR FIRST NONZERO DIGIT */
1878 if (i <= MP.t) c2 = x[i + 2] ;
1879 if ((i__1 = c - j1) < 0) goto L140 ;
1880 else if (i__1 == 0) goto L150 ;
1884 if (c2 < j2) goto L140 ;
1886 /* COMPUTE T+4 QUOTIENT DIGITS */
1893 /* MAIN LOOP FOR LARGE ABS(IY) CASE */
1897 if (k > i2) goto L120 ;
1900 /* GET APPROXIMATE QUOTIENT FIRST */
1905 /* NOW REDUCE SO OVERFLOW DOES NOT OCCUR */
1908 if (iq < b2) goto L190 ;
1910 /* HERE IQ*B WOULD POSSIBLY OVERFLOW SO INCREASE IR */
1916 iq = iq * MP.b - ir * j2 ;
1917 if (iq >= 0) goto L200 ;
1919 /* HERE IQ NEGATIVE SO IR WAS TOO LARGE */
1925 if (i <= MP.t) iq += x[i + 2] ;
1928 /* R(K) = QUOTIENT, C = REMAINDER */
1930 MP.r[k - 1] = iqj + ir ;
1932 if (c >= 0) goto L170 ;
1934 /* CARRY NEGATIVE SO OVERFLOW MUST HAVE OCCURRED */
1937 mpchk(&c__1, &c__4) ;
1939 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_DIVIB]) ;
1946 /* UNDERFLOW HERE */
1960 /* RETURNS LOGICAL VALUE OF (X .EQ. Y) FOR MP X AND Y. */
1962 --y ; /* Parameter adjustments */
1965 ret_val = mpcomp(&x[1], &y[1]) == 0 ;
1974 /* THIS ROUTINE IS CALLED WHEN AN ERROR CONDITION IS ENCOUNTERED, AND
1975 * AFTER A MESSAGE HAS BEEN WRITTEN TO STDERR.
1978 doerr(vstrs[(int) V_ERROR]) ;
1989 static int i, i2, i3, ie, ix, ts, xs, tss ;
1990 static float rx, ry, rlb ;
1992 /* RETURNS Y = EXP(X) FOR MP X AND Y.
1993 * EXP OF INTEGER AND FRACTIONAL PARTS OF X ARE COMPUTED
1994 * SEPARATELY. SEE ALSO COMMENTS IN MPEXP1.
1995 * TIME IS O(SQRT(T)M(T)).
1996 * DIMENSION OF R MUST BE AT LEAST 4T+10 IN CALLING PROGRAM
1997 * CHECK LEGALITY OF B, T, M AND MXR
2000 --y ; /* Parameter adjustments */
2003 mpchk(&c__4, &c__10) ;
2004 i2 = (MP.t << 1) + 7 ;
2005 i3 = i2 + MP.t + 2 ;
2007 /* CHECK FOR X = 0 */
2009 if (x[1] != 0) goto L10 ;
2010 mpcim(&c__1, &y[1]) ;
2013 /* CHECK IF ABS(X) .LT. 1 */
2016 if (x[2] > 0) goto L20 ;
2018 /* USE MPEXP1 HERE */
2020 mpexp1(&x[1], &y[1]) ;
2021 mpaddi(&y[1], &c__1, &y[1]) ;
2024 /* SEE IF ABS(X) SO LARGE THAT EXP(X) WILL CERTAINLY OVERFLOW
2025 * OR UNDERFLOW. 1.01 IS TO ALLOW FOR ERRORS IN ALOG.
2029 rlb = log((float) MP.b) * (float)1.01 ;
2030 r__1 = -(double)((float) (MP.m + 1)) * rlb ;
2031 if (mpcmpr(&x[1], &r__1) >= 0) goto L40 ;
2033 /* UNDERFLOW SO CALL MPUNFL AND RETURN */
2040 r__1 = (float) MP.m * rlb ;
2041 if (mpcmpr(&x[1], &r__1) <= 0) goto L70 ;
2046 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_EXPA]) ;
2051 /* NOW SAFE TO CONVERT X TO REAL */
2056 /* SAVE SIGN AND WORK WITH ABS(X) */
2059 mpabs(&x[1], &MP.r[i3 - 1]) ;
2061 /* IF ABS(X) .GT. M POSSIBLE THAT INT(X) OVERFLOWS,
2065 if (dabs(rx) > (float) MP.m)
2066 mpdivi(&MP.r[i3 - 1], &c__32, &MP.r[i3 - 1]) ;
2068 /* GET FRACTIONAL AND INTEGER PARTS OF ABS(X) */
2070 mpcmi(&MP.r[i3 - 1], &ix) ;
2071 mpcmf(&MP.r[i3 - 1], &MP.r[i3 - 1]) ;
2073 /* ATTACH SIGN TO FRACTIONAL PART AND COMPUTE EXP OF IT */
2075 MP.r[i3 - 1] = xs * MP.r[i3 - 1] ;
2076 mpexp1(&MP.r[i3 - 1], &y[1]) ;
2077 mpaddi(&y[1], &c__1, &y[1]) ;
2079 /* COMPUTE E-2 OR 1/E USING TWO EXTRA DIGITS IN CASE ABS(X) LARGE
2080 * (BUT ONLY ONE EXTRA DIGIT IF T .LT. 4)
2085 if (MP.t < 4) ts = MP.t + 1 ;
2088 i3 = i2 + MP.t + 2 ;
2090 mpcim(&xs, &MP.r[i2 - 1]) ;
2093 /* LOOP FOR E COMPUTATION. DECREASE T IF POSSIBLE. */
2099 i__1 = ts, i__2 = ts + 2 + MP.r[i2] ;
2100 MP.t = min(i__1,i__2) ;
2101 if (MP.t <= 2) goto L90 ;
2104 mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]) ;
2106 mpadd2(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i3 - 1], &
2107 MP.r[i2 - 1], &c__0) ;
2108 if (MP.r[i2 - 1] != 0) goto L80 ;
2110 /* RAISE E OR 1/E TO POWER IX */
2115 mpaddi(&MP.r[i3 - 1], &c__2, &MP.r[i3 - 1]) ;
2116 mppwr(&MP.r[i3 - 1], &ix, &MP.r[i3 - 1]) ;
2122 /* MULTIPLY EXPS OF INTEGER AND FRACTIONAL PARTS */
2124 mpmul(&y[1], &MP.r[i3 - 1], &y[1]) ;
2126 /* MUST CORRECT RESULT IF DIVIDED BY 32 ABOVE. */
2128 if (dabs(rx) <= (float) MP.m || y[1] == 0) goto L110 ;
2130 for (i = 1; i <= 5; ++i)
2133 /* SAVE EXPONENT TO AVOID OVERFLOW IN MPMUL */
2137 mpmul(&y[1], &y[1], &y[1]) ;
2140 /* CHECK FOR UNDERFLOW AND OVERFLOW */
2142 if (y[2] < -MP.m) goto L30 ;
2143 if (y[2] > MP.m) goto L50 ;
2146 /* CHECK THAT RELATIVE ERROR LESS THAN 0.01 UNLESS ABS(X) LARGE
2147 * (WHEN EXP MIGHT OVERFLOW OR UNDERFLOW)
2151 if (dabs(rx) > (float)10.0) return ;
2154 if ((r__1 = ry - exp(rx), dabs(r__1)) < ry * (float)0.01) return ;
2156 /* THE FOLLOWING MESSAGE MAY INDICATE THAT
2157 * B**(T-1) IS TOO SMALL, OR THAT M IS TOO SMALL SO THE
2158 * RESULT UNDERFLOWED.
2161 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_EXPB]) ;
2174 static int i, q, i2, i3, ib, ic, ts ;
2177 /* ASSUMES THAT X AND Y ARE MP NUMBERS, -1 .LT. X .LT. 1.
2178 * RETURNS Y = EXP(X) - 1 USING AN O(SQRT(T).M(T)) ALGORITHM
2179 * DESCRIBED IN - R. P. BRENT, THE COMPLEXITY OF MULTIPLE-
2180 * PRECISION ARITHMETIC (IN COMPLEXITY OF COMPUTATIONAL PROBLEM
2181 * SOLVING, UNIV. OF QUEENSLAND PRESS, BRISBANE, 1976, 126-165).
2182 * ASYMPTOTICALLY FASTER METHODS EXIST, BUT ARE NOT USEFUL
2183 * UNLESS T IS VERY LARGE. SEE COMMENTS TO MPATAN AND MPPIGL.
2184 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 3T+8
2185 * CHECK LEGALITY OF B, T, M AND MXR
2188 --y ; /* Parameter adjustments */
2191 mpchk(&c__3, &c__8) ;
2193 i3 = i2 + MP.t + 2 ;
2195 /* CHECK FOR X = 0 */
2197 if (x[1] != 0) goto L20 ;
2203 /* CHECK THAT ABS(X) .LT. 1 */
2206 if (x[2] <= 0) goto L40 ;
2208 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_EXP1]) ;
2214 mpstr(&x[1], &MP.r[i2 - 1]) ;
2215 rlb = log((float) MP.b) ;
2217 /* COMPUTE APPROXIMATELY OPTIMAL Q (AND DIVIDE X BY 2**Q) */
2219 q = (int) (sqrt((float) MP.t * (float).48 * rlb) + (float) x[2] *
2220 (float)1.44 * rlb) ;
2224 if (q <= 0) goto L60 ;
2228 for (i = 1; i <= i__1; ++i)
2231 if (ic < ib && ic != MP.b && i < q) continue ;
2232 mpdivi(&MP.r[i2 - 1], &ic, &MP.r[i2 - 1]) ;
2237 if (MP.r[i2 - 1] == 0) goto L10 ;
2238 mpstr(&MP.r[i2 - 1], &y[1]) ;
2239 mpstr(&MP.r[i2 - 1], &MP.r[i3 - 1]) ;
2243 /* SUM SERIES, REDUCING T WHERE POSSIBLE */
2246 MP.t = ts + 2 + MP.r[i3] - y[2] ;
2247 if (MP.t <= 2) goto L80 ;
2249 MP.t = min(MP.t,ts) ;
2250 mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]) ;
2252 mpdivi(&MP.r[i3 - 1], &i, &MP.r[i3 - 1]) ;
2254 mpadd2(&MP.r[i3 - 1], &y[1], &y[1], &y[1], &c__0) ;
2255 if (MP.r[i3 - 1] != 0) goto L70 ;
2259 if (q <= 0) return ;
2261 /* APPLY (X+1)**2 - 1 = X(2 + X) FOR Q ITERATIONS */
2264 for (i = 1; i <= i__1; ++i)
2266 mpaddi(&y[1], &c__2, &MP.r[i2 - 1]) ;
2267 mpmul(&MP.r[i2 - 1], &y[1], &y[1]) ;
2279 /* ROUTINE CALLED BY MPDIV AND MPSQRT TO ENSURE THAT
2280 * RESULTS ARE REPRESENTED EXACTLY IN T-2 DIGITS IF THEY
2281 * CAN BE. X IS AN MP NUMBER, I AND J ARE INTEGERS.
2284 --x ; /* Parameter adjustments */
2286 if (x[1] == 0 || MP.t <= 2 || *i == 0) return ;
2288 /* COMPUTE MAXIMUM POSSIBLE ERROR IN THE LAST PLACE */
2290 q = (*j + 1) / *i + 1 ;
2291 s = MP.b * x[MP.t + 1] + x[MP.t + 2] ;
2292 if (s > q) goto L10 ;
2294 /* SET LAST TWO DIGITS TO ZERO */
2301 if (s + q < MP.b * MP.b) return ;
2305 x[MP.t + 1] = MP.b - 1 ;
2306 x[MP.t + 2] = MP.b ;
2308 /* NORMALIZE X (LAST DIGIT B IS OK IN MPMULI) */
2310 mpmuli(&x[1], &c__1, &x[1]) ;
2319 static int i, j, is, js ;
2321 /* RETURNS K = K/GCD AND L = L/GCD, WHERE GCD IS THE
2322 * GREATEST COMMON DIVISOR OF K AND L.
2323 * SAVE INPUT PARAMETERS IN LOCAL VARIABLES
2330 if (js == 0) goto L30 ;
2332 /* EUCLIDEAN ALGORITHM LOOP */
2336 if (is == 0) goto L20 ;
2338 if (js != 0) goto L10 ;
2341 /* HERE JS IS THE GCD OF I AND J */
2348 /* IF J = 0 RETURN (1, 0) UNLESS I = 0, THEN (0, 0) */
2352 if (is == 0) *k = 0 ;
2364 /* RETURNS LOGICAL VALUE OF (X .GE. Y) FOR MP X AND Y. */
2366 --y ; /* Parameter adjustments */
2369 ret_val = mpcomp(&x[1], &y[1]) >= 0 ;
2380 /* RETURNS LOGICAL VALUE OF (X .GT. Y) FOR MP X AND Y. */
2382 --y ; /* Parameter adjustments */
2385 ret_val = mpcomp(&x[1], &y[1]) > 0 ;
2396 /* RETURNS LOGICAL VALUE OF (X .LE. Y) FOR MP X AND Y. */
2398 --y ; /* Parameter adjustments */
2401 ret_val = mpcomp(&x[1], &y[1]) <= 0 ;
2412 static int e, k, i2, i3 ;
2413 static float rx, rlx ;
2415 /* RETURNS Y = LN(X), FOR MP X AND Y, USING MPLNS.
2416 * RESTRICTION - INTEGER PART OF LN(X) MUST BE REPRESENTABLE
2417 * AS A SINGLE-PRECISION INTEGER. TIME IS O(SQRT(T).M(T)).
2418 * FOR SMALL INTEGER X, MPLNI IS FASTER.
2419 * ASYMPTOTICALLY FASTER METHODS EXIST (EG THE GAUSS-SALAMIN
2420 * METHOD, SEE MPLNGS), BUT ARE NOT USEFUL UNLESS T IS LARGE.
2421 * SEE COMMENTS TO MPATAN, MPEXP1 AND MPPIGL.
2422 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 6T+14.
2423 * CHECK LEGALITY OF B, T, M AND MXR
2426 --y ; /* Parameter adjustments */
2429 mpchk(&c__6, &c__14) ;
2430 i2 = (MP.t << 2) + 11 ;
2431 i3 = i2 + MP.t + 2 ;
2433 /* CHECK THAT X IS POSITIVE */
2434 if (x[1] > 0) goto L20 ;
2436 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_LNA]) ;
2442 /* MOVE X TO LOCAL STORAGE */
2445 mpstr(&x[1], &MP.r[i2 - 1]) ;
2449 /* LOOP TO GET APPROXIMATE LN(X) USING SINGLE-PRECISION */
2452 mpaddi(&MP.r[i2 - 1], &c_n1, &MP.r[i3 - 1]) ;
2454 /* IF POSSIBLE GO TO CALL MPLNS */
2456 if (MP.r[i3 - 1] == 0 || MP.r[i3] + 1 <= 0) goto L50 ;
2458 /* REMOVE EXPONENT TO AVOID FLOATING-POINT OVERFLOW */
2462 mpcmr(&MP.r[i2 - 1], &rx) ;
2464 /* RESTORE EXPONENT AND COMPUTE SINGLE-PRECISION LOG */
2467 rlx = log(rx) + (float) e * log((float) MP.b) ;
2468 r__1 = -(double)rlx ;
2469 mpcrm(&r__1, &MP.r[i3 - 1]) ;
2471 /* UPDATE Y AND COMPUTE ACCURATE EXP OF APPROXIMATE LOG */
2473 mpsub(&y[1], &MP.r[i3 - 1], &y[1]) ;
2474 mpexp(&MP.r[i3 - 1], &MP.r[i3 - 1]) ;
2476 /* COMPUTE RESIDUAL WHOSE LOG IS STILL TO BE FOUND */
2478 mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]) ;
2480 /* MAKE SURE NOT LOOPING INDEFINITELY */
2482 if (k < 10) goto L30 ;
2484 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_LNB]) ;
2489 /* COMPUTE FINAL CORRECTION ACCURATELY USING MPLNS */
2492 mplns(&MP.r[i3 - 1], &MP.r[i3 - 1]) ;
2493 mpadd(&y[1], &MP.r[i3 - 1], &y[1]) ;
2504 static int i2, i3, i4, ts, it0, ts2, ts3 ;
2506 /* RETURNS MP Y = LN(1+X) IF X IS AN MP NUMBER SATISFYING THE
2507 * CONDITION ABS(X) .LT. 1/B, ERROR OTHERWISE.
2508 * USES NEWTONS METHOD TO SOLVE THE EQUATION
2509 * EXP1(-Y) = X, THEN REVERSES SIGN OF Y.
2510 * (HERE EXP1(Y) = EXP(Y) - 1 IS COMPUTED USING MPEXP1).
2511 * TIME IS O(SQRT(T).M(T)) AS FOR MPEXP1, SPACE = 5T+12.
2512 * CHECK LEGALITY OF B, T, M AND MXR
2515 --y ; /* Parameter adjustments */
2518 mpchk(&c__5, &c__12) ;
2519 i2 = (MP.t << 1) + 7 ;
2520 i3 = i2 + MP.t + 2 ;
2521 i4 = i3 + MP.t + 2 ;
2523 /* CHECK FOR X = 0 EXACTLY */
2525 if (x[1] != 0) goto L10 ;
2529 /* CHECK THAT ABS(X) .LT. 1/B */
2532 if (x[2] + 1 <= 0) goto L30 ;
2534 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_LNSA]) ;
2540 /* SAVE T AND GET STARTING APPROXIMATION TO -LN(1+X) */
2544 mpstr(&x[1], &MP.r[i3 - 1]) ;
2545 mpdivi(&x[1], &c__4, &MP.r[i2 - 1]) ;
2546 mpaddq(&MP.r[i2 - 1], &c_n1, &c__3, &MP.r[i2 - 1]) ;
2547 mpmul(&x[1], &MP.r[i2 - 1], &MP.r[i2 - 1]) ;
2548 mpaddq(&MP.r[i2 - 1], &c__1, &c__2, &MP.r[i2 - 1]) ;
2549 mpmul(&x[1], &MP.r[i2 - 1], &MP.r[i2 - 1]) ;
2550 mpaddi(&MP.r[i2 - 1], &c_n1, &MP.r[i2 - 1]) ;
2551 mpmul(&x[1], &MP.r[i2 - 1], &y[1]) ;
2553 /* START NEWTON ITERATION USING SMALL T, LATER INCREASE */
2557 i__1 = 5, i__2 = 13 - (MP.b << 1) ;
2558 MP.t = max(i__1,i__2) ;
2559 if (MP.t > ts) goto L80 ;
2561 it0 = (MP.t + 5) / 2 ;
2564 mpexp1(&y[1], &MP.r[i4 - 1]) ;
2565 mpmul(&MP.r[i3 - 1], &MP.r[i4 - 1], &MP.r[i2 - 1]) ;
2566 mpadd(&MP.r[i4 - 1], &MP.r[i2 - 1], &MP.r[i4 - 1]) ;
2567 mpadd(&MP.r[i3 - 1], &MP.r[i4 - 1], &MP.r[i4 - 1]) ;
2568 mpsub(&y[1], &MP.r[i4 - 1], &y[1]) ;
2569 if (MP.t >= ts) goto L60 ;
2571 /* FOLLOWING LOOP COMPUTES NEXT VALUE OF T TO USE.
2572 * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE,
2573 * WE CAN ALMOST DOUBLE T EACH TIME.
2581 MP.t = (MP.t + it0) / 2 ;
2582 if (MP.t > ts3) goto L50 ;
2587 /* CHECK THAT NEWTON ITERATION WAS CONVERGING AS EXPECTED */
2590 if (MP.r[i4 - 1] == 0 || MP.r[i4] << 1 <= it0 - MP.t)
2595 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_LNSB]) ;
2596 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_LNSC]) ;
2601 /* REVERSE SIGN OF Y AND RETURN */
2616 /* RETURNS LOGICAL VALUE OF (X .LT. Y) FOR MP X AND Y. */
2618 --y ; /* Parameter adjustments */
2621 ret_val = mpcomp(&x[1], &y[1]) < 0 ;
2634 /* SETS X TO THE LARGEST POSSIBLE POSITIVE MP NUMBER
2635 * CHECK LEGALITY OF B, T, M AND MXR
2638 --x ; /* Parameter adjustments */
2640 mpchk(&c__1, &c__4) ;
2643 /* SET FRACTION DIGITS TO B-1 */
2646 for (i = 1; i <= i__1; ++i) x[i + 2] = it ;
2648 /* SET SIGN AND EXPONENT */
2658 int *u, *v, *w, *j ;
2664 /* PERFORMS INNER MULTIPLICATION LOOP FOR MPMUL
2665 * NOTE THAT CARRIES ARE NOT PROPAGATED IN INNER LOOP,
2666 * WHICH SAVES TIME AT THE EXPENSE OF SPACE.
2669 --v ; /* Parameter adjustments */
2673 for (i = 1; i <= i__1; ++i) u[i] += *w * v[i] ;
2682 int i__1, i__2, i__3, i__4 ;
2684 static int c, i, j, i2, j1, re, ri, xi, rs, i2p ;
2686 /* MULTIPLIES X AND Y, RETURNING RESULT IN Z, FOR MP X, Y AND Z.
2687 * THE SIMPLE O(T**2) ALGORITHM IS USED, WITH
2688 * FOUR GUARD DIGITS AND R*-ROUNDING.
2689 * ADVANTAGE IS TAKEN OF ZERO DIGITS IN X, BUT NOT IN Y.
2690 * ASYMPTOTICALLY FASTER ALGORITHMS ARE KNOWN (SEE KNUTH,
2691 * VOL. 2), BUT ARE DIFFICULT TO IMPLEMENT IN FORTRAN IN AN
2692 * EFFICIENT AND MACHINE-INDEPENDENT MANNER.
2693 * IN COMMENTS TO OTHER MP ROUTINES, M(T) IS THE TIME
2694 * TO PERFORM T-DIGIT MP MULTIPLICATION. THUS
2695 * M(T) = O(T**2) WITH THE PRESENT VERSION OF MPMUL,
2696 * BUT M(T) = O(T.LOG(T).LOG(LOG(T))) IS THEORETICALLY POSSIBLE.
2697 * CHECK LEGALITY OF B, T, M AND MXR
2700 --z ; /* Parameter adjustments */
2704 mpchk(&c__1, &c__4) ;
2708 /* FORM SIGN OF PRODUCT */
2711 if (rs != 0) goto L10 ;
2713 /* SET RESULT TO ZERO */
2718 /* FORM EXPONENT OF PRODUCT */
2723 /* CLEAR ACCUMULATOR */
2726 for (i = 1; i <= i__1; ++i) MP.r[i - 1] = 0 ;
2728 /* PERFORM MULTIPLICATION */
2732 for (i = 1; i <= i__1; ++i)
2736 /* FOR SPEED, PUT THE NUMBER WITH MANY ZEROS FIRST */
2738 if (xi == 0) continue ;
2742 i__3 = MP.t, i__4 = i2 - i ;
2743 i__2 = min(i__3,i__4) ;
2744 mpmlp(&MP.r[i], &y[3], &xi, &i__2) ;
2746 if (c > 0) continue ;
2748 /* CHECK FOR LEGAL BASE B DIGIT */
2750 if (xi < 0 || xi >= MP.b) goto L90 ;
2752 /* PROPAGATE CARRIES AT END AND EVERY EIGHTH TIME,
2753 * FASTER THAN DOING IT EVERY TIME.
2757 for (j = 1; j <= i__2; ++j)
2760 ri = MP.r[j1 - 1] + c ;
2761 if (ri < 0) goto L70 ;
2763 MP.r[j1 - 1] = ri - MP.b * c ;
2765 if (c != 0) goto L90 ;
2769 if (c == 8) goto L60 ;
2770 if (xi < 0 || xi >= MP.b) goto L90 ;
2773 for (j = 1; j <= i__1; ++j)
2776 ri = MP.r[j1 - 1] + c ;
2777 if (ri < 0) goto L70 ;
2779 MP.r[j1 - 1] = ri - MP.b * c ;
2781 if (c != 0) goto L90 ;
2783 /* NORMALIZE AND ROUND RESULT */
2786 mpnzr(&rs, &re, &z[1], &c__0) ;
2790 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_MULA]) ;
2797 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_MULB]) ;
2798 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_MULC]) ;
2809 mpmul2(x, iy, z, trunc)
2810 int *x, *iy, *z, *trunc ;
2814 static int c, i, j, c1, c2, j1, j2 ;
2815 static int t1, t3, t4, ij, re, ri, is, ix, rs ;
2817 /* MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z.
2818 * MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER
2819 * EVEN IF SOME DIGITS ARE GREATER THAN B-1.
2820 * RESULT IS ROUNDED IF TRUNC.EQ.0, OTHERWISE TRUNCATED.
2823 --z ; /* Parameter adjustments */
2827 if (rs == 0) goto L10 ;
2829 if (j < 0) goto L20 ;
2830 else if (j == 0) goto L10 ;
2843 /* CHECK FOR MULTIPLICATION BY B */
2845 if (j != MP.b) goto L50 ;
2846 if (x[2] < MP.m) goto L40 ;
2848 mpchk(&c__1, &c__4) ;
2850 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_MUL2A]) ;
2856 mpstr(&x[1], &z[1]) ;
2861 /* SET EXPONENT TO EXPONENT(X) + 4 */
2866 /* FORM PRODUCT IN ACCUMULATOR */
2873 /* IF J*B NOT REPRESENTABLE AS AN INTEGER WE HAVE TO SIMULATE
2874 * DOUBLE-PRECISION MULTIPLICATION.
2879 i__1 = MP.b << 3, i__2 = 32767 / MP.b ;
2880 if (j >= max(i__1,i__2)) goto L110 ;
2883 for (ij = 1; ij <= i__1; ++ij)
2886 ri = j * x[i + 2] + c ;
2888 MP.r[i + 3] = ri - MP.b * c ;
2891 /* CHECK FOR INTEGER OVERFLOW */
2893 if (ri < 0) goto L130 ;
2895 /* HAVE TO TREAT FIRST FOUR WORDS OF R SEPARATELY */
2897 for (ij = 1; ij <= 4; ++ij)
2902 MP.r[i - 1] = ri - MP.b * c ;
2904 if (c == 0) goto L100 ;
2906 /* HAVE TO SHIFT RIGHT HERE AS CARRY OFF END */
2910 for (ij = 1; ij <= i__1; ++ij)
2913 MP.r[i] = MP.r[i - 1] ;
2917 MP.r[0] = ri - MP.b * c ;
2919 if (c < 0) goto L130 ;
2920 else if (c == 0) goto L100 ;
2923 /* NORMALIZE AND ROUND OR TRUNCATE RESULT */
2926 mpnzr(&rs, &re, &z[1], trunc) ;
2929 /* HERE J IS TOO LARGE FOR SINGLE-PRECISION MULTIPLICATION */
2933 j2 = j - j1 * MP.b ;
2938 for (ij = 1; ij <= i__1; ++ij)
2941 c2 = c - MP.b * c1 ;
2944 if (i > 0) ix = x[i + 2] ;
2947 c = j1 * ix + c1 + is ;
2948 MP.r[i + 3] = ri - MP.b * is ;
2951 if (c < 0) goto L130 ;
2952 else if (c == 0) goto L100 ;
2955 /* CAN ONLY GET HERE IF INTEGER OVERFLOW OCCURRED */
2958 mpchk(&c__1, &c__4) ;
2960 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_MUL2B]) ;
2972 /* MULTIPLIES MP X BY SINGLE-PRECISION INTEGER IY GIVING MP Z.
2973 * THIS IS FASTER THAN USING MPMUL. RESULT IS ROUNDED.
2974 * MULTIPLICATION BY 1 MAY BE USED TO NORMALIZE A NUMBER
2975 * EVEN IF THE LAST DIGIT IS B.
2978 --z ; /* Parameter adjustments */
2981 mpmul2(&x[1], iy, &z[1], &c__0) ;
2988 int *x, *i, *j, *y ;
2994 /* MULTIPLIES MP X BY I/J, GIVING MP Y */
2996 --y ; /* Parameter adjustments */
2999 if (*j != 0) goto L20 ;
3000 mpchk(&c__1, &c__4) ;
3002 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_MULQ]) ;
3008 if (*i != 0) goto L40 ;
3014 /* REDUCE TO LOWEST TERMS */
3020 if (C_abs(is) == 1) goto L50 ;
3022 mpdivi(&x[1], &js, &y[1]) ;
3023 mpmul2(&y[1], &is, &y[1], &c__0) ;
3030 mpdivi(&x[1], &i__1, &y[1]) ;
3040 /* SETS Y = -X FOR MP NUMBERS X AND Y */
3042 --y ; /* Parameter adjustments */
3045 mpstr(&x[1], &y[1]) ;
3052 mpnzr(rs, re, z, trunc)
3053 int *rs, *re, *z, *trunc ;
3057 static int i, j, k, b2, i2, is, it, i2m, i2p ;
3059 /* ASSUMES LONG (I.E. (T+4)-DIGIT) FRACTION IN
3060 * R, SIGN = RS, EXPONENT = RE. NORMALIZES,
3061 * AND RETURNS MP RESULT IN Z. INTEGER ARGUMENTS RS AND RE
3062 * ARE NOT PRESERVED. R*-ROUNDING IS USED IF TRUNC.EQ.0
3065 --z ; /* Parameter adjustments */
3068 if (*rs != 0) goto L20 ;
3070 /* STORE ZERO IN Z */
3076 /* CHECK THAT SIGN = +-1 */
3079 if (C_abs(*rs) <= 1) goto L40 ;
3083 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_NZRA]) ;
3084 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_NZRB]) ;
3090 /* LOOK FOR FIRST NONZERO DIGIT */
3094 for (i = 1; i <= i__1; ++i)
3097 if (MP.r[i - 1] > 0) goto L60 ;
3105 if (is == 0) goto L90 ;
3112 for (j = 1; j <= i__1; ++j)
3115 MP.r[j - 1] = MP.r[k - 1] ;
3119 for (j = i2p; j <= i__1; ++j) MP.r[j - 1] = 0 ;
3121 /* CHECK TO SEE IF TRUNCATION IS DESIRED */
3124 if (*trunc != 0) goto L150 ;
3126 /* SEE IF ROUNDING NECESSARY
3127 * TREAT EVEN AND ODD BASES DIFFERENTLY
3131 if (b2 << 1 != MP.b) goto L130 ;
3133 /* B EVEN. ROUND IF R(T+1).GE.B2 UNLESS R(T) ODD AND ALL ZEROS
3137 if ((i__1 = MP.r[MP.t] - b2) < 0) goto L150 ;
3138 else if (i__1 == 0) goto L100 ;
3142 if (MP.r[MP.t - 1] % 2 == 0) goto L110 ;
3144 if (MP.r[MP.t + 1] + MP.r[MP.t + 2] +
3145 MP.r[MP.t + 3] == 0)
3152 for (j = 1; j <= i__1; ++j)
3156 if (MP.r[i - 1] < MP.b) goto L150 ;
3160 /* EXCEPTIONAL CASE, ROUNDED UP TO .10000... */
3166 /* ODD BASE, ROUND IF R(T+1)... .GT. 1/2 */
3169 for (i = 1; i <= 4; ++i)
3172 if ((i__1 = MP.r[it - 1] - b2) < 0) goto L150 ;
3173 else if (i__1 == 0) continue ;
3177 /* CHECK FOR OVERFLOW */
3180 if (*re <= MP.m) goto L170 ;
3182 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_NZRC]) ;
3187 /* CHECK FOR UNDERFLOW */
3190 if (*re < -MP.m) goto L190 ;
3192 /* STORE RESULT IN Z */
3197 for (i = 1; i <= i__1; ++i) z[i + 2] = MP.r[i - 1] ;
3200 /* UNDERFLOW HERE */
3213 /* CALLED ON MULTIPLE-PRECISION OVERFLOW, IE WHEN THE
3214 * EXPONENT OF MP NUMBER X WOULD EXCEED M.
3215 * AT PRESENT EXECUTION IS TERMINATED WITH AN ERROR MESSAGE
3216 * AFTER CALLING MPMAXR(X), BUT IT WOULD BE POSSIBLE TO RETURN,
3217 * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION AFTER
3218 * A PRESET NUMBER OF OVERFLOWS. ACTION COULD EASILY BE DETERMINED
3219 * BY A FLAG IN LABELLED COMMON.
3220 * M MAY HAVE BEEN OVERWRITTEN, SO CHECK B, T, M ETC.
3223 --x ; /* Parameter adjustments */
3225 mpchk(&c__1, &c__4) ;
3227 /* SET X TO LARGEST POSSIBLE POSITIVE NUMBER */
3231 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_OVFL]) ;
3233 /* TERMINATE EXECUTION BY CALLING MPERR */
3249 /* SETS MP X = PI TO THE AVAILABLE PRECISION.
3250 * USES PI/4 = 4.ARCTAN(1/5) - ARCTAN(1/239).
3252 * DIMENSION OF R MUST BE AT LEAST 3T+8
3253 * CHECK LEGALITY OF B, T, M AND MXR
3256 --x ; /* Parameter adjustments */
3258 mpchk(&c__3, &c__8) ;
3260 /* ALLOW SPACE FOR MPART1 */
3262 i2 = (MP.t << 1) + 7 ;
3263 mpart1(&c__5, &MP.r[i2 - 1]) ;
3264 mpmuli(&MP.r[i2 - 1], &c__4, &MP.r[i2 - 1]) ;
3265 mpart1(&c__239, &x[1]) ;
3266 mpsub(&MP.r[i2 - 1], &x[1], &x[1]) ;
3267 mpmuli(&x[1], &c__4, &x[1]) ;
3269 /* RETURN IF ERROR IS LESS THAN 0.01 */
3272 if ((r__1 = rx - (float)3.1416, dabs(r__1)) < (float) 0.01) return ;
3274 /* FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL */
3276 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_PI]) ;
3287 static int i2, n2, ns ;
3289 /* RETURNS Y = X**N, FOR MP X AND Y, INTEGER N, WITH 0**0 = 1.
3290 * R MUST BE DIMENSIONED AT LEAST 4T+10 IN CALLING PROGRAM
3291 * (2T+6 IS ENOUGH IF N NONNEGATIVE).
3294 --y ; /* Parameter adjustments */
3299 if (n2 < 0) goto L20 ;
3300 else if (n2 == 0) goto L10 ;
3303 /* N = 0, RETURN Y = 1. */
3306 mpcim(&c__1, &y[1]) ;
3312 mpchk(&c__4, &c__10) ;
3314 if (x[1] != 0) goto L60 ;
3318 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_PWRA]) ;
3319 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_PWRB]) ;
3328 mpchk(&c__2, &c__6) ;
3329 if (x[1] != 0) goto L60 ;
3331 /* X = 0, N .GT. 0, SO Y = 0 */
3340 mpstr(&x[1], &y[1]) ;
3342 /* IF N .LT. 0 FORM RECIPROCAL */
3344 if (*n < 0) mprec(&y[1], &y[1]) ;
3345 mpstr(&y[1], &MP.r[i2 - 1]) ;
3347 /* SET PRODUCT TERM TO ONE */
3349 mpcim(&c__1, &y[1]) ;
3351 /* MAIN LOOP, LOOK AT BITS OF N2 FROM RIGHT */
3356 if (n2 << 1 != ns) mpmul(&y[1], &MP.r[i2 - 1], &y[1]) ;
3357 if (n2 <= 0) return ;
3359 mpmul(&MP.r[i2 - 1], &MP.r[i2 - 1], &MP.r[i2 - 1]) ;
3370 /* RETURNS Z = X**Y FOR MP NUMBERS X, Y AND Z, WHERE X IS
3371 * POSITIVE (X .EQ. 0 ALLOWED IF Y .GT. 0). SLOWER THAN
3372 * MPPWR AND MPQPWR, SO USE THEM IF POSSIBLE.
3373 * DIMENSION OF R IN COMMON AT LEAST 7T+16
3374 * CHECK LEGALITY OF B, T, M AND MXR
3377 --z ; /* Parameter adjustments */
3381 mpchk(&c__7, &c__16) ;
3382 if (x[1] < 0) goto L10 ;
3383 else if (x[1] == 0) goto L30 ;
3387 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_PWR2A]) ;
3391 /* HERE X IS ZERO, RETURN ZERO IF Y POSITIVE, OTHERWISE ERROR */
3394 if (y[1] > 0) goto L60 ;
3396 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_PWR2B]) ;
3401 /* RETURN ZERO HERE */
3407 /* USUAL CASE HERE, X POSITIVE
3408 * USE MPLN AND MPEXP TO COMPUTE POWER
3412 i2 = MP.t * 6 + 15 ;
3413 mpln(&x[1], &MP.r[i2 - 1]) ;
3414 mpmul(&y[1], &MP.r[i2 - 1], &z[1]) ;
3416 /* IF X**Y IS TOO LARGE, MPEXP WILL PRINT ERROR MESSAGE */
3418 mpexp(&z[1], &z[1]) ;
3428 /* Initialized data */
3430 static int it[9] = { 0, 8, 6, 5, 4, 4, 4, 4, 4 } ;
3434 static int i2, i3, ex, ts, it0, ts2, ts3 ;
3437 /* RETURNS Y = 1/X, FOR MP X AND Y.
3438 * MPROOT (X, -1, Y) HAS THE SAME EFFECT.
3439 * DIMENSION OF R MUST BE AT LEAST 4*T+10 IN CALLING PROGRAM
3440 * (BUT Y(1) MAY BE R(3T+9)).
3441 * NEWTONS METHOD IS USED, SO FINAL ONE OR TWO DIGITS MAY
3445 --y ; /* Parameter adjustments */
3448 /* CHECK LEGALITY OF B, T, M AND MXR */
3450 mpchk(&c__4, &c__10) ;
3452 /* MPADDI REQUIRES 2T+6 WORDS. */
3454 i2 = (MP.t << 1) + 7 ;
3455 i3 = i2 + MP.t + 2 ;
3456 if (x[1] != 0) goto L20 ;
3458 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_RECA]) ;
3467 /* TEMPORARILY INCREASE M TO AVOID OVERFLOW */
3471 /* SET EXPONENT TO ZERO SO RX NOT TOO LARGE OR SMALL. */
3476 /* USE SINGLE-PRECISION RECIPROCAL AS FIRST APPROXIMATION */
3478 r__1 = (float)1. / rx ;
3479 mpcrm(&r__1, &MP.r[i2 - 1]) ;
3481 /* RESTORE EXPONENT */
3485 /* CORRECT EXPONENT OF FIRST APPROXIMATION */
3489 /* SAVE T (NUMBER OF DIGITS) */
3493 /* START WITH SMALL T TO SAVE TIME. ENSURE THAT B**(T-1) .GE. 100 */
3496 if (MP.b < 10) MP.t = it[MP.b - 1] ;
3497 it0 = (MP.t + 4) / 2 ;
3498 if (MP.t > ts) goto L70 ;
3500 /* MAIN ITERATION LOOP */
3503 mpmul(&x[1], &MP.r[i2 - 1], &MP.r[i3 - 1]) ;
3504 mpaddi(&MP.r[i3 - 1], &c_n1, &MP.r[i3 - 1]) ;
3506 /* TEMPORARILY REDUCE T */
3509 MP.t = (MP.t + it0) / 2 ;
3510 mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]) ;
3515 mpsub(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]) ;
3516 if (MP.t >= ts) goto L50 ;
3518 /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE
3519 * BECAUSE NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
3526 MP.t = (MP.t + it0) / 2 ;
3527 if (MP.t > ts3) goto L40 ;
3529 MP.t = min(ts,ts2) ;
3532 /* RETURN IF NEWTON ITERATION WAS CONVERGING */
3535 if (MP.r[i3 - 1] == 0 || MP.r[i2] - MP.r[i3] << 1 >=
3539 /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
3540 * OR THAT THE STARTING APPROXIMATION IS NOT ACCURATE ENOUGH.
3545 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_RECB]) ;
3546 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_RECC]) ;
3551 /* MOVE RESULT TO Y AND RETURN AFTER RESTORING T */
3555 mpstr(&MP.r[i2 - 1], &y[1]) ;
3557 /* RESTORE M AND CHECK FOR OVERFLOW (UNDERFLOW IMPOSSIBLE) */
3560 if (y[2] <= MP.m) return ;
3562 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_RECD]) ;
3574 /* Initialized data */
3576 static int it[9] = { 0, 8, 6, 5, 4, 4, 4, 4, 4 } ;
3581 static int i2, i3, ex, np, ts, it0, ts2, ts3, xes ;
3584 /* RETURNS Y = X**(1/N) FOR INTEGER N, ABS(N) .LE. MAX (B, 64).
3585 * AND MP NUMBERS X AND Y,
3586 * USING NEWTONS METHOD WITHOUT DIVISIONS. SPACE = 4T+10
3587 * (BUT Y(1) MAY BE R(3T+9))
3590 --y ; /* Parameter adjustments */
3593 /* CHECK LEGALITY OF B, T, M AND MXR */
3595 mpchk(&c__4, &c__10) ;
3596 if (*n != 1) goto L10 ;
3597 mpstr(&x[1], &y[1]) ;
3601 i2 = (MP.t << 1) + 7 ;
3602 i3 = i2 + MP.t + 2 ;
3603 if (*n != 0) goto L30 ;
3605 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ROOTA]) ;
3612 /* LOSS OF ACCURACY IF NP LARGE, SO ONLY ALLOW NP .LE. MAX (B, 64) */
3614 if (np <= max(MP.b,64)) goto L60 ;
3616 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ROOTB]) ;
3623 /* LOOK AT SIGN OF X */
3626 if (x[1] < 0) goto L90 ;
3627 else if (x[1] == 0) goto L70 ;
3630 /* X = 0 HERE, RETURN 0 IF N POSITIVE, ERROR IF NEGATIVE */
3634 if (*n > 0) return ;
3636 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ROOTC]) ;
3640 /* X NEGATIVE HERE, SO ERROR IF N IS EVEN */
3643 if (np % 2 != 0) goto L110 ;
3645 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ROOTD]) ;
3649 /* GET EXPONENT AND DIVIDE BY NP */
3655 /* REDUCE EXPONENT SO RX NOT TOO LARGE OR SMALL. */
3660 /* USE SINGLE-PRECISION ROOT FOR FIRST APPROXIMATION */
3662 r__1 = exp(((float) (np * ex - xes) * log((float) MP.b) - log((dabs(
3663 rx)))) / (float) np) ;
3664 mpcrm(&r__1, &MP.r[i2 - 1]) ;
3666 /* SIGN OF APPROXIMATION SAME AS THAT OF X */
3668 MP.r[i2 - 1] = x[1] ;
3670 /* RESTORE EXPONENT */
3674 /* CORRECT EXPONENT OF FIRST APPROXIMATION */
3678 /* SAVE T (NUMBER OF DIGITS) */
3682 /* START WITH SMALL T TO SAVE TIME */
3686 /* ENSURE THAT B**(T-1) .GE. 100 */
3688 if (MP.b < 10) MP.t = it[MP.b - 1] ;
3689 if (MP.t > ts) goto L160 ;
3691 /* IT0 IS A NECESSARY SAFETY FACTOR */
3693 it0 = (MP.t + 4) / 2 ;
3695 /* MAIN ITERATION LOOP */
3698 mppwr(&MP.r[i2 - 1], &np, &MP.r[i3 - 1]) ;
3699 mpmul(&x[1], &MP.r[i3 - 1], &MP.r[i3 - 1]) ;
3700 mpaddi(&MP.r[i3 - 1], &c_n1, &MP.r[i3 - 1]) ;
3702 /* TEMPORARILY REDUCE T */
3705 MP.t = (MP.t + it0) / 2 ;
3706 mpmul(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i3 - 1]) ;
3707 mpdivi(&MP.r[i3 - 1], &np, &MP.r[i3 - 1]) ;
3712 mpsub(&MP.r[i2 - 1], &MP.r[i3 - 1], &MP.r[i2 - 1]) ;
3714 /* FOLLOWING LOOP ALMOST DOUBLES T (POSSIBLE BECAUSE
3715 * NEWTONS METHOD HAS 2ND ORDER CONVERGENCE).
3718 if (MP.t >= ts) goto L140 ;
3723 MP.t = (MP.t + it0) / 2 ;
3724 if (MP.t > ts3) goto L130 ;
3726 MP.t = min(ts,ts2) ;
3729 /* NOW R(I2) IS X**(-1/NP)
3730 * CHECK THAT NEWTON ITERATION WAS CONVERGING
3734 if (MP.r[i3 - 1] == 0 || MP.r[i2] - MP.r[i3] << 1 >=
3738 /* THE FOLLOWING MESSAGE MAY INDICATE THAT B**(T-1) IS TOO SMALL,
3739 * OR THAT THE INITIAL APPROXIMATION OBTAINED USING ALOG AND EXP
3740 * IS NOT ACCURATE ENOUGH.
3745 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ROOTE]) ;
3746 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_ROOTF]) ;
3755 if (*n < 0) goto L170 ;
3758 mppwr(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]) ;
3759 mpmul(&x[1], &MP.r[i2 - 1], &y[1]) ;
3763 mpstr(&MP.r[i2 - 1], &y[1]) ;
3769 mpset(idecpl, itmax2, maxdr)
3770 int *idecpl, *itmax2, *maxdr ;
3774 static int i, k, w, i2, w2, wn ;
3776 /* SETS BASE (B) AND NUMBER OF DIGITS (T) TO GIVE THE
3777 * EQUIVALENT OF AT LEAST IDECPL DECIMAL DIGITS.
3778 * IDECPL SHOULD BE POSITIVE.
3779 * ITMAX2 IS THE DIMENSION OF ARRAYS USED FOR MP NUMBERS,
3780 * SO AN ERROR OCCURS IF THE COMPUTED T EXCEEDS ITMAX2 - 2.
3782 * MXR = MAXDR (DIMENSION OF R IN COMMON, .GE. T+4), AND
3783 * M = (W-1)/4 (MAXIMUM ALLOWABLE EXPONENT),
3784 * WHERE W IS THE LARGEST INTEGER OF THE FORM 2**K-1 WHICH IS
3785 * REPRESENTABLE IN THE MACHINE, K .LE. 47
3786 * THE COMPUTED B AND T SATISFY THE CONDITIONS
3787 * (T-1)*LN(B)/LN(10) .GE. IDECPL AND 8*B*B-1 .LE. W .
3788 * APPROXIMATELY MINIMAL T AND MAXIMAL B SATISFYING
3789 * THESE CONDITIONS ARE CHOSEN.
3790 * PARAMETERS IDECPL, ITMAX2 AND MAXDR ARE INTEGERS.
3791 * BEWARE - MPSET WILL CAUSE AN INTEGER OVERFLOW TO OCCUR
3792 * ****** IF WORDLENGTH IS LESS THAN 48 BITS.
3793 * IF THIS IS NOT ALLOWABLE, CHANGE THE DETERMINATION
3794 * OF W (DO 30 ... TO 30 W = WN) OR SET B, T, M,
3795 * AND MXR WITHOUT CALLING MPSET.
3801 /* DETERMINE LARGE REPRESENTABLE INTEGER W OF FORM 2**K - 1 */
3806 /* ON CYBER 76 HAVE TO FIND K .LE. 47, SO ONLY LOOP
3807 * 47 TIMES AT MOST. IF GENUINE INTEGER WORDLENGTH
3808 * EXCEEDS 47 BITS THIS RESTRICTION CAN BE RELAXED.
3811 for (i = 1; i <= 47; ++i)
3814 /* INTEGER OVERFLOW WILL EVENTUALLY OCCUR HERE
3815 * IF WORDLENGTH .LT. 48 BITS
3821 /* APPARENTLY REDUNDANT TESTS MAY BE NECESSARY ON SOME
3822 * MACHINES, DEPENDING ON HOW INTEGER OVERFLOW IS HANDLED
3825 if (wn <= w || wn <= w2 || wn <= 0) goto L40 ;
3830 /* SET MAXIMUM EXPONENT TO (W-1)/4 */
3834 if (*idecpl > 0) goto L60 ;
3836 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_SETB]) ;
3841 /* B IS THE LARGEST POWER OF 2 SUCH THAT (8*B*B-1) .LE. W */
3844 i__1 = (k - 3) / 2 ;
3845 MP.b = pow_ii(&c__2, &i__1) ;
3847 /* 2E0 BELOW ENSURES AT LEAST ONE GUARD DIGIT */
3849 MP.t = (int) ((float) (*idecpl) * log((float)10.) / log((float)
3850 MP.b) + (float)2.0) ;
3852 /* SEE IF T TOO LARGE FOR DIMENSION STATEMENTS */
3855 if (i2 <= *itmax2) goto L80 ;
3861 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_SETC]) ;
3862 _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_SETD]) ;
3863 msg = (char *) XtMalloc(strlen( mpstrs[(int) MP_SETE]) + 200);
3864 sprintf(msg, mpstrs[(int) MP_SETE], i2) ;
3865 _DtSimpleError (v->appname, DtWarning, NULL, msg);
3871 /* REDUCE TO MAXIMUM ALLOWED BY DIMENSION STATEMENTS */
3873 MP.t = *itmax2 - 2 ;
3875 /* CHECK LEGALITY OF B, T, M AND MXR (AT LEAST T+4) */
3878 mpchk(&c__1, &c__4) ;
3889 static int i2, i3, ie, xs ;
3890 static float rx, ry ;
3892 /* RETURNS Y = SIN(X) FOR MP X AND Y,
3893 * METHOD IS TO REDUCE X TO (-1, 1) AND USE MPSIN1, SO
3894 * TIME IS O(M(T)T/LOG(T)).
3895 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 5T+12
3896 * CHECK LEGALITY OF B, T, M AND MXR
3899 --y ; /* Parameter adjustments */
3902 mpchk(&c__5, &c__12) ;
3903 i2 = (MP.t << 2) + 11 ;
3904 if (x[1] != 0) goto L20 ;
3913 if (ie <= 2) mpcmr(&x[1], &rx) ;
3915 mpabs(&x[1], &MP.r[i2 - 1]) ;
3917 /* USE MPSIN1 IF ABS(X) .LE. 1 */
3919 if (mpcmpi(&MP.r[i2 - 1], &c__1) > 0) goto L30 ;
3921 mpsin1(&MP.r[i2 - 1], &y[1], &c__1) ;
3924 /* FIND ABS(X) MODULO 2PI (IT WOULD SAVE TIME IF PI WERE
3925 * PRECOMPUTED AND SAVED IN COMMON).
3926 * FOR INCREASED ACCURACY COMPUTE PI/4 USING MPART1
3930 i3 = (MP.t << 1) + 7 ;
3931 mpart1(&c__5, &MP.r[i3 - 1]) ;
3932 mpmuli(&MP.r[i3 - 1], &c__4, &MP.r[i3 - 1]) ;
3933 mpart1(&c__239, &y[1]) ;
3934 mpsub(&MP.r[i3 - 1], &y[1], &y[1]) ;
3935 mpdiv(&MP.r[i2 - 1], &y[1], &MP.r[i2 - 1]) ;
3936 mpdivi(&MP.r[i2 - 1], &c__8, &MP.r[i2 - 1]) ;
3937 mpcmf(&MP.r[i2 - 1], &MP.r[i2 - 1]) ;
3939 /* SUBTRACT 1/2, SAVE SIGN AND TAKE ABS */
3941 mpaddq(&MP.r[i2 - 1], &c_n1, &c__2, &MP.r[i2 - 1]) ;
3942 xs = -xs * MP.r[i2 - 1] ;
3943 if (xs == 0) goto L10 ;
3946 mpmuli(&MP.r[i2 - 1], &c__4, &MP.r[i2 - 1]) ;
3948 /* IF NOT LESS THAN 1, SUBTRACT FROM 2 */
3951 mpaddi(&MP.r[i2 - 1], &c_n2, &MP.r[i2 - 1]) ;
3953 if (MP.r[i2 - 1] == 0) goto L10 ;
3956 mpmuli(&MP.r[i2 - 1], &c__2, &MP.r[i2 - 1]) ;
3958 /* NOW REDUCED TO FIRST QUADRANT, IF LESS THAN PI/4 USE
3959 * POWER SERIES, ELSE COMPUTE COS OF COMPLEMENT
3962 if (MP.r[i2] > 0) goto L40 ;
3964 mpmul(&MP.r[i2 - 1], &y[1], &MP.r[i2 - 1]) ;
3965 mpsin1(&MP.r[i2 - 1], &y[1], &c__1) ;
3969 mpaddi(&MP.r[i2 - 1], &c_n2, &MP.r[i2 - 1]) ;
3970 mpmul(&MP.r[i2 - 1], &y[1], &MP.r[i2 - 1]) ;
3971 mpsin1(&MP.r[i2 - 1], &y[1], &c__0) ;
3975 if (ie > 2) return ;
3977 /* CHECK THAT ABSOLUTE ERROR LESS THAN 0.01 IF ABS(X) .LE. 100
3978 * (IF ABS(X) IS LARGE THEN SINGLE-PRECISION SIN INACCURATE)
3981 if (dabs(rx) > (float)100.) return ;
3984 if ((r__1 = ry - sin(rx), dabs(r__1)) < (float) 0.01) return ;
3986 /* THE FOLLOWING MESSAGE MAY INDICATE THAT
3987 * B**(T-1) IS TOO SMALL.
3990 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_SIN]) ;
4003 static int i, b2, i2, i3, ts ;
4005 /* COMPUTES Y = SIN(X) IF IS.NE.0, Y = COS(X) IF IS.EQ.0,
4006 * USING TAYLOR SERIES. ASSUMES ABS(X) .LE. 1.
4007 * X AND Y ARE MP NUMBERS, IS AN INTEGER.
4008 * TIME IS O(M(T)T/LOG(T)). THIS COULD BE REDUCED TO
4009 * O(SQRT(T)M(T)) AS IN MPEXP1, BUT NOT WORTHWHILE UNLESS
4010 * T IS VERY LARGE. ASYMPTOTICALLY FASTER METHODS ARE
4011 * DESCRIBED IN THE REFERENCES GIVEN IN COMMENTS
4012 * TO MPATAN AND MPPIGL.
4013 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 3T+8
4014 * CHECK LEGALITY OF B, T, M AND MXR
4017 --y ; /* Parameter adjustments */
4020 mpchk(&c__3, &c__8) ;
4021 if (x[1] != 0) goto L20 ;
4023 /* SIN(0) = 0, COS(0) = 1 */
4027 if (*is == 0) mpcim(&c__1, &y[1]) ;
4032 i3 = i2 + MP.t + 2 ;
4033 b2 = max(MP.b,64) << 1 ;
4034 mpmul(&x[1], &x[1], &MP.r[i3 - 1]) ;
4035 if (mpcmpi(&MP.r[i3 - 1], &c__1) <= 0) goto L40 ;
4037 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_SIN1]) ;
4043 if (*is == 0) mpcim(&c__1, &MP.r[i2 - 1]) ;
4044 if (*is != 0) mpstr(&x[1], &MP.r[i2 - 1]) ;
4049 if (*is == 0) goto L50 ;
4051 mpstr(&MP.r[i2 - 1], &y[1]) ;
4054 /* POWER SERIES LOOP. REDUCE T IF POSSIBLE */
4057 MP.t = MP.r[i2] + ts + 2 ;
4058 if (MP.t <= 2) goto L80 ;
4060 MP.t = min(MP.t,ts) ;
4062 /* PUT R(I3) FIRST IN CASE ITS DIGITS ARE MAINLY ZERO */
4064 mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &MP.r[i2 - 1]) ;
4066 /* IF I*(I+1) IS NOT REPRESENTABLE AS AN INTEGER, THE FOLLOWING
4067 * DIVISION BY I*(I+1) HAS TO BE SPLIT UP.
4070 if (i > b2) goto L60 ;
4072 i__1 = -i * (i + 1) ;
4073 mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]) ;
4078 mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]) ;
4080 mpdivi(&MP.r[i2 - 1], &i__1, &MP.r[i2 - 1]) ;
4085 mpadd2(&MP.r[i2 - 1], &y[1], &y[1], &y[1], &c__0) ;
4086 if (MP.r[i2 - 1] != 0) goto L50 ;
4090 if (*is == 0) mpaddi(&y[1], &c__1, &y[1]) ;
4101 static int i2, i3, xs ;
4103 /* RETURNS Y = SINH(X) FOR MP NUMBERS X AND Y, X NOT TOO LARGE.
4104 * METHOD IS TO USE MPEXP OR MPEXP1, SPACE = 5T+12
4105 * SAVE SIGN OF X AND CHECK FOR ZERO, SINH(0) = 0
4108 --y ; /* Parameter adjustments */
4112 if (xs != 0) goto L10 ;
4117 /* CHECK LEGALITY OF B, T, M AND MXR */
4120 mpchk(&c__5, &c__12) ;
4121 i3 = (MP.t << 2) + 11 ;
4123 /* WORK WITH ABS(X) */
4125 mpabs(&x[1], &MP.r[i3 - 1]) ;
4126 if (MP.r[i3] <= 0) goto L20 ;
4128 /* HERE ABS(X) .GE. 1, IF TOO LARGE MPEXP GIVES ERROR MESSAGE
4129 * INCREASE M TO AVOID OVERFLOW IF SINH(X) REPRESENTABLE
4133 mpexp(&MP.r[i3 - 1], &MP.r[i3 - 1]) ;
4134 mprec(&MP.r[i3 - 1], &y[1]) ;
4135 mpsub(&MP.r[i3 - 1], &y[1], &y[1]) ;
4137 /* RESTORE M. IF RESULT OVERFLOWS OR UNDERFLOWS, MPDIVI AT
4138 * STATEMENT 30 WILL ACT ACCORDINGLY.
4144 /* HERE ABS(X) .LT. 1 SO USE MPEXP1 TO AVOID CANCELLATION */
4147 i2 = i3 - (MP.t + 2) ;
4148 mpexp1(&MP.r[i3 - 1], &MP.r[i2 - 1]) ;
4149 mpaddi(&MP.r[i2 - 1], &c__2, &MP.r[i3 - 1]) ;
4150 mpmul(&MP.r[i3 - 1], &MP.r[i2 - 1], &y[1]) ;
4151 mpaddi(&MP.r[i2 - 1], &c__1, &MP.r[i3 - 1]) ;
4152 mpdiv(&y[1], &MP.r[i3 - 1], &y[1]) ;
4154 /* DIVIDE BY TWO AND RESTORE SIGN */
4158 mpdivi(&y[1], &i__1, &y[1]) ;
4167 static int i, i2, iy3 ;
4169 /* RETURNS Y = SQRT(X), USING SUBROUTINE MPROOT IF X .GT. 0.
4170 * DIMENSION OF R IN CALLING PROGRAM MUST BE AT LEAST 4T+10
4171 * (BUT Y(1) MAY BE R(3T+9)). X AND Y ARE MP NUMBERS.
4172 * CHECK LEGALITY OF B, T, M AND MXR
4175 --y ; /* Parameter adjustments */
4178 mpchk(&c__4, &c__10) ;
4180 /* MPROOT NEEDS 4T+10 WORDS, BUT CAN OVERLAP SLIGHTLY. */
4183 if (x[1] < 0) goto L10 ;
4184 else if (x[1] == 0) goto L30 ;
4188 if (v->MPerrors) _DtSimpleError (v->appname, DtWarning, NULL, mpstrs[(int) MP_SQRT]) ;
4197 mproot(&x[1], &c_n2, &MP.r[i2 - 1]) ;
4199 mpmul(&x[1], &MP.r[i2 - 1], &y[1]) ;
4201 mpext(&i, &iy3, &y[1]) ;
4212 static int i, j, t2 ;
4214 /* SETS Y = X FOR MP X AND Y.
4215 * SEE IF X AND Y HAVE THE SAME ADDRESS (THEY OFTEN DO)
4218 --y ; /* Parameter adjustments */
4223 if (j == x[1]) goto L10 ;
4225 /* HERE X(1) AND Y(1) MUST HAVE THE SAME ADDRESS */
4230 /* HERE X(1) AND Y(1) HAVE DIFFERENT ADDRESSES */
4235 /* NO NEED TO MOVE X(2), ... IF X(1) = 0 */
4237 if (j == 0) return ;
4241 for (i = 2; i <= i__1; ++i) y[i] = x[i] ;
4252 /* SUBTRACTS Y FROM X, FORMING RESULT IN Z, FOR MP X, Y AND Z.
4253 * FOUR GUARD DIGITS ARE USED, AND THEN R*-ROUNDING
4256 --z ; /* Parameter adjustments */
4261 mpadd2(&x[1], &y[1], &z[1], y1, &c__0) ;
4274 /* RETURNS Y = TANH(X) FOR MP NUMBERS X AND Y,
4275 * USING MPEXP OR MPEXP1, SPACE = 5T+12
4278 --y ; /* Parameter adjustments */
4281 if (x[1] != 0) goto L10 ;
4288 /* CHECK LEGALITY OF B, T, M AND MXR */
4291 mpchk(&c__5, &c__12) ;
4292 i2 = (MP.t << 2) + 11 ;
4294 /* SAVE SIGN AND WORK WITH ABS(X) */
4297 mpabs(&x[1], &MP.r[i2 - 1]) ;
4299 /* SEE IF ABS(X) SO LARGE THAT RESULT IS +-1 */
4301 r__1 = (float) MP.t * (float).5 * log((float) MP.b) ;
4302 mpcrm(&r__1, &y[1]) ;
4303 if (mpcomp(&MP.r[i2 - 1], &y[1]) <= 0) goto L20 ;
4305 /* HERE ABS(X) IS VERY LARGE */
4310 /* HERE ABS(X) NOT SO LARGE */
4313 mpmuli(&MP.r[i2 - 1], &c__2, &MP.r[i2 - 1]) ;
4314 if (MP.r[i2] <= 0) goto L30 ;
4316 /* HERE ABS(X) .GE. 1/2 SO USE MPEXP */
4318 mpexp(&MP.r[i2 - 1], &MP.r[i2 - 1]) ;
4319 mpaddi(&MP.r[i2 - 1], &c_n1, &y[1]) ;
4320 mpaddi(&MP.r[i2 - 1], &c__1, &MP.r[i2 - 1]) ;
4321 mpdiv(&y[1], &MP.r[i2 - 1], &y[1]) ;
4324 /* HERE ABS(X) .LT. 1/2, SO USE MPEXP1 TO AVOID CANCELLATION */
4327 mpexp1(&MP.r[i2 - 1], &MP.r[i2 - 1]) ;
4328 mpaddi(&MP.r[i2 - 1], &c__2, &y[1]) ;
4329 mpdiv(&MP.r[i2 - 1], &y[1], &y[1]) ;
4344 /* CALLED ON MULTIPLE-PRECISION UNDERFLOW, IE WHEN THE
4345 * EXPONENT OF MP NUMBER X WOULD BE LESS THAN -M.
4346 * SINCE M MAY HAVE BEEN OVERWRITTEN, CHECK B, T, M ETC.
4349 --x ; /* Parameter adjustments */
4351 mpchk(&c__1, &c__4) ;
4353 /* THE UNDERFLOWING NUMBER IS SET TO ZERO
4354 * AN ALTERNATIVE WOULD BE TO CALL MPMINR (X) AND RETURN,
4355 * POSSIBLY UPDATING A COUNTER AND TERMINATING EXECUTION
4356 * AFTER A PRESET NUMBER OF UNDERFLOWS. ACTION COULD EASILY
4357 * BE DETERMINED BY A FLAG IN LABELLED COMMON.
4381 if (x == 0) return(pow) ;
4387 if (n & 01) pow *= x;
4388 if (n >>= 1) x *= x ;
4409 if (n & 01) pow *= x ;
4410 if (n >>= 1) x *= x ;
4433 if (x == 0) return(pow) ;
4439 if (n & 01) pow *= x ;
4440 if (n >>= 1) x *= x ;