C ALGORITHM 567 C C EXTENDED-RANGE ARITHMETIC AND NORMALIZED LEGENDRE POLYNOMIALS C C BY D.W. LOZIER AND J.M. SMITH C C ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE 7,1 (MARCH, 1981) C MAN 10 MAN 20 C *** MAIN PROGRAM *** MAN 30 C THIS IS A DEMONSTRATION PROGRAM FOR THE NORMALIZED LEGENDRE POLY- MAN 40 C NOMIAL SUBROUTINE NORMP, WHICH USES EXTENDED-RANGE ARITHMETIC. SEE THEMAN 50 C COMMENTS IN SETUP, THE SETTING-UP SUBROUTINE FOR EXTENDED-RANGE ARITH-MAN 60 C METIC, FOR FURTHER INFORMATION ON THIS SPECIAL ARITHMETIC AND FOR A MAN 70 C DESCRIPTION OF THE COMPUTER-DEPENDENT PARAMETERS. THESE PARAMETERS AREMAN 80 C THE ONLY INPUT TO THIS PROGRAM. MAN 90 C MAN 100 C SEE THE PROGRAM COMMENTS IN NORMP FOR A DESCRIPTION OF NORMALIZED MAN 110 C LEGENDRE POLYNOMIALS AND FOR A DESCRIPTION OF INPUT TO NORMP. THE MAN 120 C DEGREE (NU), ORDER (MU), AND ARGUMENT (ARG) VALUES ARE SPECIFIED MAN 130 C AS DATA IN THIS PROGRAM. THE USER MAY WISH TO CHANGE THIS DATA. MAN 140 C ARGUMENTS ARE SPECIFIED BY GIVING ANGLES THETA1 AND THETA2 MAN 150 C (WHERE THETA1 + THETA2 = 180) IN DEGREES BETWEEN 0 AND 180. MAN 160 C THE FUNCTION VALUES ARE COMPUTED IN MODE=1 FOR ARG=DCOS(THETA1), MAN 170 C AND IN MODE=2 FOR ARG=THETA1 AND ARG=THETA2. THE REFLECTION MAN 180 C FORMULA IS EMPLOYED WHEN MODE=2 AND ARG=THETA2. THUS EACH MAN 190 C FUNCTION VALUE IS COMPUTED THREE TIMES FROM THREE DIFFERENT MAN 200 C FORMS OF INPUT TO NORMP. MAN 210 C MAN 220 DOUBLE PRECISION ARG, DZERO, PI, PN, RADS1, RADS2, THETA1, THETA2 MAN 230 DIMENSION ARG(3), ISIG(3), PN(51,3), IPN(51,3) MAN 240 DIMENSION NU(1), MU1(2), MU2(2), THETA1(11), THETA2(11) MAN 250 DIMENSION RADS1(11), RADS2(11) MAN 260 C MAN 270 DATA NNU /1/ MAN 280 DATA NU /10000/ MAN 290 C MAN 300 DATA NMU /2/ MAN 310 DATA MU1 /0,9955/ MAN 320 DATA MU2 /45,10000/ MAN 330 C MAN 340 DATA NTHETA /11/ MAN 350 DATA THETA1 /0.1D0,1.0D0,10.0D0,89.0D0,89.9D0,135.0D0,160.0D0, MAN 360 * 176.0D0,178.0D0,179.5D0,179.7D0/ MAN 370 DATA THETA2 /179.9D0,179.0D0,170.0D0,91.0D0,90.1D0,45.0D0,20.0D0, MAN 380 * 4.0D0,2.0D0,0.5D0,0.3D0/ MAN 390 C MAN 400 PI = 4.0D0*DATAN(1.0D0) MAN 410 DO 10 I=1,NTHETA MAN 420 RADS1(I) = THETA1(I)*PI/180.0D0 MAN 430 RADS2(I) = THETA2(I)*PI/180.0D0 MAN 440 10 CONTINUE MAN 450 C MAN 460 C PRINTOUT OF LIBRARY VALUES THAT WILL BE USED IN THIS MAN 470 C DEMONSTRATION PROGRAM OR IN NORMP. MAN 480 C MAN 490 WRITE (6,99999) PI MAN 500 DO 20 I=1,NTHETA MAN 510 WRITE (6,99998) THETA1(I), THETA2(I) MAN 520 ARG(1) = DSIN(RADS1(I)) MAN 530 ARG(2) = DSIN(RADS2(I)) MAN 540 WRITE (6,99997) ARG(1), ARG(2) MAN 550 ARG(1) = DCOS(RADS1(I)) MAN 560 ARG(2) = DCOS(RADS2(I)) MAN 570 WRITE (6,99996) ARG(1), ARG(2) MAN 580 20 CONTINUE MAN 590 C MAN 600 C READ COMPUTER-DEPENDENT PARAMETERS. FOR A DESCRIPTION, SEE MAN 610 C COMMENTS IN SUBROUTINE SETUP. MAN 620 READ (5,99995) IRAD, NRADPL, DZERO, NBITS MAN 630 C MAN 640 C PRINT COMPUTER-DEPENDENT PARAMETERS. MAN 650 WRITE (6,99994) IRAD, NRADPL, DZERO, NBITS MAN 660 C MAN 670 C CALL THE SETTING-UP SUBROUTINE FOR EXTENDED-RANGE ARITHMETIC. MAN 680 CALL SETUP(IRAD, NRADPL, DZERO, NBITS) MAN 690 C MAN 700 C PRINT CONTENTS OF COMMON BLOCKS INITIALIZED BY PRECEDING CALL MAN 710 C TO SETUP. MAN 720 CALL EXRP MAN 730 C MAN 740 DO 110 I=1,NNU MAN 750 DO 100 J=1,NTHETA MAN 760 ARG(1) = DCOS(RADS1(J)) MAN 770 IF (THETA1(J).GT.90.0D0) ARG(1) = -DCOS(RADS2(J)) MAN 780 ARG(2) = RADS1(J) MAN 790 ARG(3) = RADS2(J) MAN 800 DO 90 K=1,NMU MAN 810 C COMPUTE NORMALIZED LEGENDRE POLYNOMIALS FROM THREE DIFFERENT MAN 820 C FORMS OF INPUT. MAN 830 CALL NORMP(NU(I), MU1(K), MU2(K), ARG(1), 1, PN(1,1), MAN 840 * IPN(1,1), ISIG(1)) MAN 850 CALL NORMP(NU(I), MU1(K), MU2(K), ARG(2), 2, PN(1,2), MAN 860 * IPN(1,2), ISIG(2)) MAN 870 CALL NORMP(NU(I), MU1(K), MU2(K), ARG(3), 2, PN(1,3), MAN 880 * IPN(1,3), ISIG(3)) MAN 890 C APPLY REFLECTION FORMULA TO THIRD FORM. MAN 900 MULIM = MU2(K) - MU1(K) + 1 MAN 910 DO 30 L=1,MULIM MAN 920 MU = MU1(K) + L - 1 MAN 930 PN(L,3) = PN(L,3)*DBLE(FLOAT((-1)**(MU+NU(I)))) MAN 940 30 CONTINUE MAN 950 C BEGIN OUTPUT. MAN 960 WRITE (6,99993) NU(I), ARG(1), THETA1(J), THETA2(J) MAN 970 WRITE (6,99992) MAN 980 C FIND COLUMN WITH MINIMUM ESTIMATED LOSS OF SIGNIFICANCE. MAN 990 ISIGM = MIN0(ISIG(1),ISIG(2),ISIG(3)) MAN 1000 IF (ISIG(1).NE.ISIGM) GO TO 40 MAN 1010 M = 1 MAN 1020 M1 = 2 MAN 1030 M2 = 3 MAN 1040 GO TO 60 MAN 1050 40 IF (ISIG(2).NE.ISIGM) GO TO 50 MAN 1060 M = 2 MAN 1070 M1 = 1 MAN 1080 M2 = 3 MAN 1090 GO TO 60 MAN 1100 50 M = 3 MAN 1110 M1 = 1 MAN 1120 M2 = 2 MAN 1130 60 CONTINUE MAN 1140 DO 80 L=1,MULIM MAN 1150 MU = MU1(K) + L - 1 MAN 1160 SIG1 = 1.0 MAN 1170 SIG2 = 1.0 MAN 1180 C CONVERT NORMALIZED LEGENDRE POLYNOMIAL VALUES TO DECIMAL MAN 1190 C FORM FOR PRINTING. MAN 1200 CALL CONVRT(PN(L,1), IPN(L,1)) MAN 1210 CALL CONVRT(PN(L,2), IPN(L,2)) MAN 1220 CALL CONVRT(PN(L,3), IPN(L,3)) MAN 1230 C COMPUTE RELATIVE DIFFERENCES OF THE OTHER TWO COLUMNS COMPARED MAN 1240 C TO THE LEFTMOST COLUMN OF MAXIMUM SIGNIFICANCE, AS MEASURED BY MAN 1250 C THE SIGNIFICANCE LOSS INDICATOR ISIG. SEE OUTPUT PAGES. MAN 1260 IF (PN(L,M).EQ.0.0D0) GO TO 70 MAN 1270 IF (IPN(L,M1).EQ.IPN(L,M)) SIG1 = -1.0D0 + MAN 1280 * PN(L,M1)/PN(L,M) MAN 1290 IF (IPN(L,M2).EQ.IPN(L,M)) SIG2 = -1.0D0 + MAN 1300 * PN(L,M2)/PN(L,M) MAN 1310 C PRINT CURRENT LINE OF OUTPUT TABLE. MAN 1320 70 WRITE (6,99991) MU, (PN(L,N),IPN(L,N),N=1,3), SIG1, SIG2 MAN 1330 80 CONTINUE MAN 1340 C MAN 1350 C END OF OUTPUT TABLE PRINTING. PRINT LOSS OF SIGNIFICANCE INDICATORS.MAN 1360 WRITE (6,99990) ISIG MAN 1370 90 CONTINUE MAN 1380 100 CONTINUE MAN 1390 110 CONTINUE MAN 1400 STOP MAN 1410 99999 FORMAT (47H1LIBRARY VALUES TO BE USED IN THIS COMPUTER RUN/ MAN 1420 * 13H PI =, D26.18) MAN 1430 99998 FORMAT (13H0 THETA1 =, F18.2, 23H DEGREES THETA2 =, MAN 1440 * F18.2, 8H DEGREES) MAN 1450 99997 FORMAT (13H SIN(THETA1)=, D26.18, 15H SIN(THETA2)=, D26.18) MAN 1460 99996 FORMAT (13H COS(THETA1)=, D26.18, 15H COS(THETA2)=, D26.18) MAN 1470 99995 FORMAT (2I4, D10.0, I4) MAN 1480 99994 FORMAT (13H1INPUT VALUES/37H IRAD NRADPL DZERO NB,MAN 1490 * 3HITS/2I8, D16.8, I8/) MAN 1500 99993 FORMAT (54H1NORMALIZED LEGENDRE POLYNOMIALS P(NU,MU,X) FOR FIXED ,MAN 1510 * 59HDEGREE NU AND ARGUMENT X, COMPUTED BY THE SUBROUTINE NORMP./ MAN 1520 * 53H COLUMN 1 .. ARGUMENT X SPECIFIED DIRECTLY. (MODE=1)./ MAN 1530 * 61H COLUMN 2 .. ANGLE THETA1 SPECIFIED, SO THAT X=COS(THETA1). (,MAN 1540 * 8HMODE=2)./50H COLUMN 3 .. ANGLE THETA2 SPECIFIED, SO THAT X=COS,MAN 1550 * 23H(180-THETA2). (MODE=2)./34H REL DIFFS ARE OF OTHER TWO COLS C,MAN 1560 * 61HOMPARED TO LEFTMOST COL OF MAX SIGNIFICANCE, AS MEASURED BY T,MAN 1570 * 36HHE SIGNIFICANCE LOSS INDICATOR ISIG./5H NU =, I8, 6H, X =, MAN 1580 * D26.18, 11H, THETA1 =, F7.2, 19H DEGREES, THETA2 =, F7.2, MAN 1590 * 9H DEGREES./) MAN 1600 99992 FORMAT (7X, 2HMU, 2X, 32H...........COLUMN 1............., 2X, MAN 1610 * 32H...........COLUMN 2............., 2X, 20H...........COLUMN 3.,MAN 1620 * 12H............, 4X, 16HREL. DIFFERENCES) MAN 1630 99991 FORMAT (1X, I8, D26.18, I8, D26.18, I8, D26.18, I8, 2E10.2) MAN 1640 99990 FORMAT (50H0ISIG, THE ESTIMATED NO. OF DECIMAL PLACES LOST .., MAN 1650 * I3, 14H FOR COLUMN 1,, I3, 14H FOR COLUMN 2,, I3, 11H FOR COLUMN,MAN 1660 * 3H 3.) MAN 1670 END MAN 1680 C *** SUBROUTINE NORMP *** NOR 10 SUBROUTINE NORMP(NU, MU1, MU2, ARG, MODE, PN, IPN, ISIG) NOR 20 INTEGER NU, MU1, MU2, MODE, IPN, ISIG DOUBLE PRECISION ARG, PN C C SUBROUTINE TO CALCULATE NORMALIZED LEGENDRE POLYNOMIALS C OF VARYING ORDER AND OF FIXED ARGUMENT AND DEGREE. C THE ORDER MU AND DEGREE NU ARE NONNEGATIVE INTEGERS AND THE C ARGUMENT IS REAL. THE ALGORITHM REQUIRES THE USE OF C NUMBERS OUTSIDE THE NORMAL MACHINE RANGE. THEREFORE, C THIS SUBROUTINE EMPLOYS A SPECIAL ARITHMETIC CALLED C EXTENDED-RANGE ARITHMETIC. (SEE SMITH, J.M., OLVER, C F.W.J., AND LOZIER, D.W., EXTENDED-RANGE ARITHMETIC C AND NORMALIZED LEGENDRE POLYNOMIALS, ACM TRANSACTIONS C ON MATHEMATICAL SOFTWARE, MARCH 1981). C C IN EXTENDED-RANGE ARITHMETIC EACH NUMBER IS REPRESENTED C AS A PAIR (X,IX) WHICH HAS THE VALUE X*RADIX**IX, WHERE C RADIX IS THE BASE IN WHICH CALCULATIONS ARE PERFORMED. FOR C A FULL DESCRIPTION OF EXTENDED-RANGE ARITHMETIC, SEE THE C COMMENTS AT THE BEGINNING OF SUBROUTINE SETUP. C C THE NORMALIZED LEGENDRE POLYNOMIAL IS A MULTIPLE OF THE C ASSOCIATED LEGENDRE POLYNOMIAL OF THE FIRST KIND C WHERE THE NORMALIZING COEFFICIENT IS CHOSEN SO THAT C THE INTEGRAL OF THE SQUARE OF THE FUNCTION FROM -1 C TO +1 IS EQUAL TO 1 (SEE JAHNKE,E. EMDE,F. AND LOSCH,F., C TABLES OF HIGHER FUNCTIONS, MCGRAW HILL, NEW YORK, C 1960, P. 121). C C THE INPUT VALUES TO NORMP ARE NU, MU1, MU2, ARG, AND MODE. C THESE MUST SATISFY: C 1. NU IS A NON-NEGATIVE INTEGER SPECIFYING THE DEGREE. C 2. MU1 AND MU2 ARE NON-NEGATIVE INTEGERS, WITH C MU1.LE.MU2, SPECIFYING THE RANGE OF ORDERS. C 3. MODE IS AN INTEGER, RESTRICTED TO 1 OR 2, C WHICH INDICATES WHICH OF TWO FORMS OF THE C DOUBLE-PRECISION VARIABLE ARG IS INTENDED: C A. IF MODE=1 THEN NORMALIZED LEGENDRE(NU,MU,ARG) IS C COMPUTED FOR ALL MU SUCH THAT MU1 .LE. MU .LE. MU2 . C IN THIS CASE, -1 .LE. ARG .LE. 1 MUST BE SATISFIED. C B. IF MODE=2 THEN NORMALIZED LEGENDRE(NU,MU,COS ARG) IS C COMPUTED FOR ALL MU SUCH THAT MU1 .LE. MU .LE. MU2 . C IN THIS CASE, -PI .LT. ARG .LT. PI MUST BE SATISFIED. C C THE OUTPUT OF SUBROUTINE NORMP CONSISTS OF THE TWO C ARRAYS PN AND IPN AND THE ERROR ESTIMATE ISIG. THE C NORMALIZED LEGENDRE POLYNOMIALS ARE STORED AS EXTENDED- C RANGE NUMBERS WITH DOUBLE-PRECISION PRINCIPAL PARTS C IN ARRAY PN AND AUXILIARY INDICES IN ARRAY IPN. THUS C (PN(1),IPN(1))=NORMALIZED LEGENDRE(NU,MU1,X) C (PN(2),IPN(2))=NORMALIZED LEGENDRE(NU,MU1+1,X) C . C . C (PN(K),IPN(K))=NORMALIZED LEGENDRE(NU,MU2,X) C WHERE K=MU2-MU1+1 AND X=ARG IF MODE=1, OR X=COS ARG IF MODE=2. C FINALLY, ISIG IS AN ESTIMATE OF THE NUMBER OF DECIMAL DIGITS C LOST THROUGH ROUNDING ERRORS IN THE COMPUTATION. IF ARG IS C ACCURATE TO N SIGNIFICANT DECIMALS, THEN THE COMPUTED FUNCTION C VALUES ARE ACCCURATE TO N-ISIG SIGNIFICANT DECIMALS (EXCEPT IN C NEIGHBORHOODS OF ZEROS OF THE FUNCTIONS). C DOUBLE PRECISION C1, C2, P, P1, P2, P3, S, SX, T, TX, X, DK COMMON /EXR1/ LUERR C C ARRAYS PN AND IPN ARE DIMENSIONED 1 IN NORMP. IN THE C CALLING PROGRAM THE DIMENSIONS OF PN AND IPN MUST BE C AT LEAST MU2-MU1+1 TO AVOID IMPROPER STORING OF RESULTS. C DIMENSION PN(1), IPN(1) C C TEST FOR PROPER INPUT VALUES. C IF (NU.LT.0) GO TO 110 IF (MU1.LT.0) GO TO 110 IF (MU1.GT.MU2) GO TO 110 IF (NU.EQ.0) GO TO 90 IF (MODE.LT.1 .OR. MODE.GT.2) GO TO 110 GO TO (10, 20), MODE 10 IF (DABS(ARG).GT.1.0D0) GO TO 120 IF (DABS(ARG).EQ.1.0D0) GO TO 90 X = ARG SX = DSQRT((1.0D0+DABS(X))*((0.5D0-DABS(X))+0.5D0)) TX = X/SX ISIG = DLOG10(2.0D0*DBLE(FLOAT(NU))*(5.0D0+TX**2)) + 0.5D0 GO TO 30 20 IF (DABS(ARG).GT.4.0D0*DATAN(1.0D0)) GO TO 120 IF (ARG.EQ.0.0D0) GO TO 90 X = DCOS(ARG) SX = DABS(DSIN(ARG)) TX = X/SX ISIG = DLOG10(2.0D0*DBLE(FLOAT(NU))*(5.0D0+DABS(ARG*TX))) + 0.5D0 C C BEGIN CALCULATION C 30 MU = MU2 I = MU2 - MU1 + 1 C C IF MU.GT.NU, NORMALIZED LEGENDRE(NU,MU,X)=0. C 40 IF (MU.LE.NU) GO TO 50 PN(I) = 0.0D0 IPN(I) = 0 I = I - 1 MU = MU - 1 GO TO 40 50 MU = NU C C P1 = 0. = NORMALIZED LEGENDRE(NU,NU+1,X) C P1 = 0.0D0 IP1 = 0 C C CALCULATE P2 = NORMALIZED LEGENDRE(NU,NU,X) C P2 = 1.0D0 IP2 = 0 P3 = 0.5D0 DK = 2.0D0 DO 60 J=1,NU P3 = ((DK+1.0D0)/DK)*P3 P2 = P2*SX CALL ADJUST(P2, IP2) DK = DK + 2.0D0 60 CONTINUE P2 = P2*DSQRT(P3) CALL ADJUST(P2, IP2) S = 2.0D0*TX T = 1.0D0/DBLE(FLOAT(NU)) IF (MU2.LT.NU) GO TO 70 PN(I) = P2 IPN(I) = IP2 I = I - 1 C C RECURRENCE PROCESS C 70 P = DBLE(FLOAT(MU))*T C1 = 1.0D0/DSQRT((1.0D0-P+T)*(1.0D0+P)) C2 = S*P*C1*P2 C1 = -DSQRT((1.0D0+P+T)*(1.0D0-P))*C1*P1 CALL ADD(C2, IP2, C1, IP1, P, IP) MU = MU - 1 IF (MU.GT.MU2) GO TO 80 C C STORE IN ARRAY PN FOR RETURN TO CALLING ROUTINE. C PN(I) = P IPN(I) = IP I = I - 1 80 P1 = P2 IP1 = IP2 P2 = P IP2 = IP IF (MU.LE.MU1) GO TO 140 GO TO 70 C C SPECIAL CASE WHEN X=-1 OR +1. C 90 K = MU2 - MU1 + 1 DO 100 I=1,K PN(I) = 0.0D0 IPN(I) = 0 100 CONTINUE ISIG = 0 IF (MU1.GT.0) GO TO 140 PN(1) = DSQRT(DBLE(FLOAT(NU))+0.5D0) IPN(1) = 0 IF (MOD(NU,2).EQ.0) GO TO 140 IF (MODE.EQ.1 .AND. ARG.EQ.1.0D0) GO TO 140 IF (MODE.EQ.2) GO TO 140 PN(1) = -PN(1) GO TO 140 C C ERROR PRINTOUTS AND TERMINATION. C 110 WRITE (LUERR,99999) GO TO 130 120 WRITE (LUERR,99998) 130 WRITE (LUERR,99997) NU, MU1, MU2, ARG, MODE STOP C C RETURN TO CALLING PROGRAM C 140 RETURN 99999 FORMAT (53H0***ERR IN NORMP...NU, MU1, MU2, OR MODE NOT VALID***) 99998 FORMAT (38H0***ERR IN NORMP...ARG OUT OF RANGE***) 99997 FORMAT (9H NU =, I15/9H MU1 =, I15/9H MU2 =, * I15/9H ARG =, D26.18/9H MODE =, I15/) END C *** SUBROUTINE SETUP *** SET 10 C USAGE SET 20 C CALL SETUP(IRAD,NRADPL,DZERO,NBITS) SET 30 C DESCRIPTION SET 40 C SET 50 C SUBROUTINE SETUP MUST BE CALLED PRIOR TO CALLING ANY OTHER SET 60 C EXTENDED-RANGE SUBROUTINE. IT CALCULATES AND STORES SEVERAL SET 70 C MACHINE-DEPENDENT CONSTANTS IN COMMON BLOCKS. THE USER MUST SET 80 C SUPPLY FOUR CONSTANTS THAT PERTAIN TO HIS PARTICULAR COMPUTER. SET 90 C THE CONSTANTS ARE SET 100 C SET 110 C IRAD = THE INTERNAL BASE OF DOUBLE-PRECISION SET 120 C ARITHMETIC IN THE COMPUTER. SET 130 C NRADPL = THE NUMBER OF RADIX PLACES CARRIED IN SET 140 C THE DOUBLE-PRECISION REPRESENTATION. SET 150 C DZERO = THE SMALLEST OF 1/DMIN, DMAX, DMAXLN WHERE SET 160 C DMIN = THE SMALLEST POSITIVE DOUBLE-PRECISION SET 170 C NUMBER OR AN UPPER BOUND TO THIS NUMBER, SET 180 C DMAX = THE LARGEST DOUBLE-PRECISION NUMBER SET 190 C OR A LOWER BOUND TO THIS NUMBER, SET 200 C DMAXLN = THE LARGEST DOUBLE-PRECISION NUMBER SET 210 C SUCH THAT DLOG(DMAXLN) CAN BE COMPUTED BY THE SET 220 C FORTRAN SYSTEM (ON MOST SYSTEMS DMAXLN = DMAX). SET 230 C NBITS = THE NUMBER OF BITS (EXCLUSIVE OF SIGN) IN SET 240 C AN INTEGER COMPUTER WORD. SET 250 C SET 260 SUBROUTINE SETUP(IRAD, NRADPL, DZERO, NBITS) SET 270 INTEGER IRAD, NRADPL, NBITS DOUBLE PRECISION DZERO C C THIS IS THE SETTING-UP SUBROUTINE FOR A PACKAGE OF SUBROUTINES C THAT FACILITATE THE USE OF EXTENDED-RANGE ARITHMETIC. EXTENDED-RANGE C ARITHMETIC ON A PARTICULAR COMPUTER IS DEFINED ON THE SET OF NUMBERS C OF THE FORM C C (X,IX) = X*RADIX**IX C C WHERE X IS A DOUBLE-PRECISION NUMBER CALLED THE PRINCIPAL PART, C IX IS AN INTEGER CALLED THE AUXILIARY INDEX, AND RADIX IS THE C INTERNAL BASE OF THE DOUBLE-PRECISION ARITHMETIC. OBVIOUSLY, C EACH REAL NUMBER IS REPRESENTABLE WITHOUT ERROR BY MORE THAN ONE C EXTENDED-RANGE FORM. CONVERSIONS BETWEEN DIFFERENT FORMS ARE C ESSENTIAL IN CARRYING OUT ARITHMETIC OPERATIONS. WITH THE CHOICE C OF RADIX WE HAVE MADE, AND THE SUBROUTINES WE HAVE WRITTEN, THESE C CONVERSIONS ARE PERFORMED WITHOUT ERROR (AT LEAST ON MOST COMPUTERS). C (SEE SMITH, J.M., OLVER, F.W.J., AND LOZIER, D.W., EXTENDED-RANGE C ARITHMETIC AND NORMALIZED LEGENDRE POLYNOMIALS, ACM TRANSACTIONS ON C MATHEMATICAL SOFTWARE, MARCH 1981). C C AN EXTENDED-RANGE NUMBER (X,IX) IS SAID TO BE IN ADJUSTED FORM IF C X AND IX ARE ZERO OR C C RADIX**(-L) .LE. DABS(X) .LT. RADIX**L C C IS SATISFIED, WHERE L IS A COMPUTER-DEPENDENT INTEGER DEFINED IN THIS C SUBROUTINE. TWO EXTENDED-RANGE NUMBERS IN ADJUSTED FORM CAN BE ADDED, C SUBTRACTED, MULTIPLIED OR DIVIDED (IF THE DIVISOR IS NONZERO) WITHOUT C CAUSING OVERFLOW OR UNDERFLOW IN THE PRINCIPAL PART OF THE RESULT. C WITH PROPER USE OF THE EXTENDED-RANGE SUBROUTINES, THE ONLY OVERFLOW C THAT CAN OCCUR IS INTEGER OVERFLOW IN THE AUXILIARY INDEX. IF THIS C IS DETECTED, THE EXTENDED-RANGE SOFTWARE PRINTS A MESSAGE AND STOPS. C C MULTIPLICATION AND DIVISION IS PERFORMED BY SETTING C C (X,IX)*(Y,IY) = (X*Y,IX+IY) C OR C (X,IX)/(Y,IY) = (X/Y,IX-IY). C C PRE-ADJUSTMENT OF THE OPERANDS IS ESSENTIAL TO AVOID C OVERFLOW OR UNDERFLOW OF THE PRINCIPAL PART. SUBROUTINE C ADJUST (SEE BELOW) MAY BE CALLED TO TRANSFORM ANY EXTENDED- C RANGE NUMBER INTO ADJUSTED FORM. C C ADDITION AND SUBTRACTION REQUIRE THE USE OF SUBROUTINE ADD C (SEE BELOW). THE INPUT OPERANDS NEED NOT BE IN ADJUSTED FORM. C HOWEVER, THE RESULT OF ADDITION OR SUBTRACTION IS RETURNED C IN ADJUSTED FORM. THUS, FOR EXAMPLE, IF (X,IX),(Y,IY), C (U,IU), AND (V,IV) ARE IN ADJUSTED FORM, THEN C C (X,IX)*(Y,IY) + (U,IU)*(V,IV) C C CAN BE COMPUTED AND STORED IN ADJUSTED FORM WITH NO EXPLICIT C CALLS TO ADJUST. C C WHEN AN EXTENDED-RANGE NUMBER IS TO BE PRINTED, IT MUST BE C CONVERTED TO AN EXTENDED-RANGE FORM WITH DECIMAL RADIX. SUBROUTINE C CONVRT IS PROVIDED FOR THIS PURPOSE. C C THE SUBROUTINES CONTAINED IN THIS PACKAGE ARE C C *** SUBROUTINE ADD *** C USAGE C CALL ADD(X,IX,Y,IY,Z,IZ) C DESCRIPTION C FORMS THE EXTENDED-RANGE SUM (Z,IZ) = C (X,IX) + (Y,IY). (Z,IZ) IS ADJUSTED C BEFORE RETURNING. THE INPUT OPERANDS C NEED NOT BE IN ADJUSTED FORM, BUT THEIR C PRINCIPAL PARTS MUST SATISFY C RADIX**(-2L).LE.DABS(X).LE.RADIX**(2L), C RADIX**(-2L).LE.DABS(Y).LE.RADIX**(2L). C C *** SUBROUTINE ADJUST *** C USAGE C CALL ADJUST(X,IX) C DESCRIPTION C TRANSFORMS (X,IX) SO THAT C RADIX**(-L) .LE. DABS(X) .LT. RADIX**L. C ON MOST COMPUTERS THIS TRANSFORMATION DOES C NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS C THE NUMBER BASE OF DOUBLE-PRECISION ARITHMETIC. C C *** SUBROUTINE CNV210 *** C USAGE C CALL CNV210(K,Z,J) C DESCRIPTION C GIVEN K THIS SUBROUTINE COMPUTES J AND Z C SUCH THAT RADIX**K = Z*10**J, WHERE Z IS IN C THE RANGE 1/10 .LE. Z .LT. 1. C THE VALUE OF Z WILL BE ACCURATE TO FULL C DOUBLE PRECISION PROVIDED THE NUMBER C OF DECIMAL PLACES IN THE LARGEST C INTEGER PLUS THE NUMBER OF DECIMAL C PLACES CARRIED IN DOUBLE PRECISION DOES NOT C EXCEED 60. CNV210 IS CALLED BY SUBROUTINE C CONVRT WHEN NECESSARY. THE USER SHOULD C NEVER NEED TO CALL CNV210 DIRECTLY. C C *** SUBROUTINE CONVRT *** C USAGE C CALL CONVRT(X,IX) C DESCRIPTION C CONVERTS (X,IX) = X*RADIX**IX C TO DECIMAL FORM IN PREPARATION FOR C PRINTING, SO THAT (X,IX) = X*10**IX C WHERE 1/10 .LE. DABS(X) .LT. 1 C IS RETURNED, EXCEPT THAT IF C (DABS(X),IX) IS BETWEEN RADIX**(-2L) C AND RADIX**(2L) THEN THE REDUCED C FORM WITH IX = 0 IS RETURNED. C C *** SUBROUTINE REDUCE *** C USAGE C CALL REDUCE(X,IX) C DESCRIPTION C IF C RADIX**(-2L) .LE. (DABS(X),IX) .LE. RADIX**(2L) C THEN REDUCE TRANSFORMS (X,IX) SO THAT IX=0. C IF (X,IX) IS OUTSIDE THE ABOVE RANGE, C THEN REDUCE TAKES NO ACTION. C THIS SUBROUTINE IS USEFUL IF THE C RESULTS OF EXTENDED-RANGE CALCULATIONS C ARE TO BE USED IN SUBSEQUENT ORDINARY C DOUBLE-PRECISION CALCULATIONS. C C LUERR IS THE OUTPUT UNIT NUMBER FOR ERROR MESSAGES, SET C EQUAL TO 6 BELOW. SETUP DETERMINES THE OTHER COMMON VARIABLES C FROM THE INPUT. COMMON /EXR1/ LUERR DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R COMMON /EXR2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX COMMON /EXR3/ NLG102, MLG102, LG102(21) C DIMENSION LOG102(20), LGTEMP(20) C C LOG102 CONTAINS THE FIRST 60 DIGITS OF LOG10(2) FOR USE IN C CONVERSION OF EXTENDED-RANGE NUMBERS TO BASE 10 . DATA LOG102 /301,029,995,663,981,195,213,738,894,724,493,026,768, * 189,881,462,108,541,310,428/ C LUERR = 6 IF (IRAD.EQ.2) GO TO 10 IF (IRAD.EQ.4) GO TO 10 IF (IRAD.EQ.8) GO TO 10 IF (IRAD.EQ.16) GO TO 10 WRITE (LUERR,99999) GO TO 110 10 CONTINUE IF (IRAD.EQ.2) LOG2R = 1 IF (IRAD.EQ.4) LOG2R = 2 IF (IRAD.EQ.8) LOG2R = 3 IF (IRAD.EQ.16) LOG2R = 4 RADIX = IRAD DLG10R = DLOG10(RADIX) L = 0.5D0*DLOG10(0.5D0*DZERO)/DLG10R L2 = 2*L IF (L.GE.4) GO TO 20 WRITE (LUERR,99998) GO TO 110 20 RADIXL = RADIX**L RAD2L = RADIXL**2 IF (15.LE.NBITS .AND. NBITS.LE.60) GO TO 30 WRITE (LUERR,99997) GO TO 110 30 CONTINUE KMAX = 2**(NBITS-1) - L2 NB = (NBITS-1)/2 MLG102 = 2**NB IF (1.LE.NRADPL*LOG2R .AND. NRADPL*LOG2R.LE.120) GO TO 40 WRITE (LUERR,99996) GO TO 110 40 CONTINUE NLG102 = NRADPL*LOG2R/NB + 3 NP1 = NLG102 + 1 C C AFTER COMPLETION OF THE FOLLOWING LOOP, IC CONTAINS C THE INTEGER PART AND LGTEMP CONTAINS THE FRACTIONAL PART C OF LOG10(IRAD) IN RADIX 1000. IC = 0 DO 50 II=1,20 I = 21 - II IT = LOG2R*LOG102(I) + IC IC = IT/1000 LGTEMP(I) = MOD(IT,1000) 50 CONTINUE C C AFTER COMPLETION OF THE FOLLOWING LOOP, LG102 CONTAINS C LOG10(IRAD) IN RADIX MLG102. THE RADIX POINT IS C BETWEEN LG102(1) AND LG102(2). LG102(1) = IC DO 80 I=2,NP1 LG102(I) = 0 DO 70 J=1,NB IC = 0 DO 60 KK=1,20 K = 21 - KK IT = 2*LGTEMP(K) + IC IC = IT/1000 LGTEMP(K) = MOD(IT,1000) 60 CONTINUE LG102(I) = 2*LG102(I) + IC 70 CONTINUE 80 CONTINUE C C CHECK SPECIAL CONDITIONS REQUIRED BY SUBROUTINES... IF (NRADPL.LT.L) GO TO 90 WRITE (LUERR,99995) GO TO 110 90 IF (6*L.LE.KMAX) GO TO 100 WRITE (LUERR,99994) GO TO 110 100 CONTINUE RETURN 110 STOP 99999 FORMAT (44H ***ERR IN SETUP...IMPROPER VALUE OF IRAD***) 99998 FORMAT (45H ***ERR IN SETUP...IMPROPER VALUE OF DZERO***) 99997 FORMAT (45H ***ERR IN SETUP...IMPROPER VALUE OF NBITS***) 99996 FORMAT (46H ***ERR IN SETUP...IMPROPER VALUE OF NRADPL***) 99995 FORMAT (30H ***ERR IN SETUP...NRADPL.GE.L) 99994 FORMAT (30H ***ERR IN SETUP...6*L.GT.KMAX) END C *** SUBROUTINE ADD *** ADD 10 C USAGE ADD 20 C CALL ADD(X,IX,Y,IY,Z,IZ) ADD 30 C DESCRIPTION ADD 40 C FORMS THE EXTENDED-RANGE SUM (Z,IZ) = ADD 50 C (X,IX) + (Y,IY). (Z,IZ) IS ADJUSTED ADD 60 C BEFORE RETURNING. THE INPUT OPERANDS ADD 70 C NEED NOT BE IN ADJUSTED FORM, BUT THEIR ADD 80 C PRINCIPAL PARTS MUST SATISFY ADD 90 C RADIX**(-2L).LE.DABS(X).LE.RADIX**(2L), ADD 100 C RADIX**(-2L).LE.DABS(Y).LE.RADIX**(2L). ADD 110 C ADD 120 SUBROUTINE ADD(X, IX, Y, IY, Z, IZ) ADD 130 DOUBLE PRECISION X, Y, Z INTEGER IX, IY, IZ C DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R COMMON /EXR2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX C DOUBLE PRECISION S, T C C THE CONDITIONS IMPOSED ON L AND KMAX BY THIS SUBROUTINE C ARE C (1) 1 .LT. L .LE. 0.5D0*LOGR(0.5D0*DZERO) C C (2) NRADPL .LT. L .LE. KMAX/6 C C (3) KMAX .LE. (2**NBITS - 4*L - 1)/2 C C THESE CONDITIONS MUST BE MET BY APPROPRIATE CODING C IN SUBROUTINE SETUP. C IF (X.NE.0.0D0) GO TO 10 Z = Y IZ = IY GO TO 220 10 IF (Y.NE.0.0D0) GO TO 20 Z = X IZ = IX GO TO 220 20 CONTINUE IF (IX.GE.0 .AND. IY.GE.0) GO TO 40 IF (IX.LT.0 .AND. IY.LT.0) GO TO 40 IF (IABS(IX).LE.6*L .AND. IABS(IY).LE.6*L) GO TO 40 IF (IX.GE.0) GO TO 30 Z = Y IZ = IY GO TO 220 30 CONTINUE Z = X IZ = IX GO TO 220 40 I = IX - IY IF (I) 80, 50, 90 50 IF (DABS(X).GT.1.0D0 .AND. DABS(Y).GT.1.0D0) GO TO 60 IF (DABS(X).LT.1.0D0 .AND. DABS(Y).LT.1.0D0) GO TO 70 Z = X + Y IZ = IX GO TO 220 60 S = X/RADIXL T = Y/RADIXL Z = S + T IZ = IX + L GO TO 220 70 S = X*RADIXL T = Y*RADIXL Z = S + T IZ = IX - L GO TO 220 80 S = Y IS = IY T = X GO TO 100 90 S = X IS = IX T = Y 100 CONTINUE C C AT THIS POINT, THE ONE OF (X,IX) OR (Y,IY) THAT HAS THE C LARGER AUXILIARY INDEX IS STORED IN (S,IS). THE PRINCIPAL C PART OF THE OTHER INPUT IS STORED IN T. C I1 = IABS(I)/L I2 = MOD(IABS(I),L) IF (DABS(T).GE.RADIXL) GO TO 130 IF (DABS(T).GE.1.0D0) GO TO 120 IF (RADIXL*DABS(T).GE.1.0D0) GO TO 110 J = I1 + 1 T = T*RADIX**(L-I2) GO TO 140 110 J = I1 T = T*RADIX**(-I2) GO TO 140 120 J = I1 - 1 IF (J.LT.0) GO TO 110 T = T*RADIX**(-I2)/RADIXL GO TO 140 130 J = I1 - 2 IF (J.LT.0) GO TO 120 T = T*RADIX**(-I2)/RAD2L 140 CONTINUE C C AT THIS POINT, SOME OR ALL OF THE DIFFERENCE IN THE C AUXILIARY INDICES HAS BEEN USED TO EFFECT A LEFT SHIFT C OF T. THE SHIFTED VALUE OF T SATISFIES C C RADIX**(-2*L) .LE. DABS(T) .LE. 1.0D0 C C AND, IF J=0, NO FURTHER SHIFTING REMAINS TO BE DONE. C IF (J.EQ.0) GO TO 190 IF (DABS(S).GE.RADIXL .OR. J.GT.3) GO TO 150 IF (DABS(S).GE.1.0D0) GO TO (180, 150, 150), J IF (RADIXL*DABS(S).GE.1.0D0) GO TO (180, 170, 150), J GO TO (180, 170, 160), J 150 Z = S IZ = IS GO TO 220 160 S = S*RADIXL 170 S = S*RADIXL 180 S = S*RADIXL 190 CONTINUE C C AT THIS POINT, THE REMAINING DIFFERENCE IN THE C AUXILIARY INDICES HAS BEEN USED TO EFFECT A RIGHT SHIFT C OF S. IF THE SHIFTED VALUE OF S WOULD HAVE EXCEEDED C RADIX**L, THEN (S,IS) IS RETURNED AS THE VALUE OF THE C SUM. C IF (DABS(S).GT.1.0D0 .AND. DABS(T).GT.1.0D0) GO TO 200 IF (DABS(S).LT.1.0D0 .AND. DABS(T).LT.1.0D0) GO TO 210 Z = S + T IZ = IS - J*L GO TO 220 200 S = S/RADIXL T = T/RADIXL Z = S + T IZ = IS - J*L + L GO TO 220 210 S = S*RADIXL T = T*RADIXL Z = S + T IZ = IS - J*L - L 220 CALL ADJUST(Z, IZ) RETURN END C *** SUBROUTINE ADJUST *** ADJ 10 C USAGE ADJ 20 C CALL ADJUST(X,IX) ADJ 30 C DESCRIPTION ADJ 40 C TRANSFORMS (X,IX) SO THAT ADJ 50 C RADIX**(-L) .LE. DABS(X) .LT. RADIX**L. ADJ 60 C ON MOST COMPUTERS THIS TRANSFORMATION DOES ADJ 70 C NOT CHANGE THE MANTISSA OF X PROVIDED RADIX IS ADJ 80 C THE NUMBER BASE OF DOUBLE-PRECISION ARITHMETIC. ADJ 90 C ADJ 100 SUBROUTINE ADJUST(X, IX) ADJ 110 DOUBLE PRECISION X INTEGER IX C COMMON /EXR1/ LUERR DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R COMMON /EXR2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX C C THE CONDITION IMPOSED ON L AND KMAX BY THIS SUBROUTINE C IS C 2*L .LE. KMAX C C THIS CONDITION MUST BE MET BY APPROPRIATE CODING C IN SUBROUTINE SETUP. C IF (X.EQ.0.0D0) GO TO 50 IF (DABS(X).GE.1.0D0) GO TO 20 IF (RADIXL*DABS(X).GE.1.0D0) GO TO 60 X = X*RAD2L IF (IX.LT.0) GO TO 10 IX = IX - L2 GO TO 70 10 IF (IX.LT.-KMAX+L2) GO TO 40 IX = IX - L2 GO TO 70 20 IF (DABS(X).LT.RADIXL) GO TO 60 X = X/RAD2L IF (IX.GT.0) GO TO 30 IX = IX + L2 GO TO 70 30 IF (IX.GT.KMAX-L2) GO TO 40 IX = IX + L2 GO TO 70 40 WRITE (LUERR,99999) STOP 50 IX = 0 60 IF (IABS(IX).GT.KMAX) GO TO 40 70 RETURN 99999 FORMAT (44H0***ERR IN ADJUST...OVERFLOW IN AUX INDEX***) END C ***SUBROUTINE CNV210 *** CNV 10 C USAGE CNV 20 C CALL CNV210(K,Z,J) CNV 30 C DESCRIPTION CNV 40 C GIVEN K THIS SUBROUTINE COMPUTES J AND Z CNV 50 C SUCH THAT RADIX**K = Z*10**J, WHERE Z IS IN CNV 60 C THE RANGE 1/10 .LE. Z .LT. 1. CNV 70 C THE VALUE OF Z WILL BE ACCURATE TO FULL CNV 80 C DOUBLE PRECISION PROVIDED THE NUMBER CNV 90 C OF DECIMAL PLACES IN THE LARGEST CNV 100 C INTEGER PLUS THE NUMBER OF DECIMAL CNV 110 C PLACES CARRIED IN DOUBLE PRECISION DOES NOT CNV 120 C EXCEED 60. CNV210 IS CALLED BY SUBROUTINE CNV 130 C CONVRT WHEN NECESSARY. THE USER SHOULD CNV 140 C NEVER NEED TO CALL CNV210 DIRECTLY. CNV 150 C CNV 160 SUBROUTINE CNV210(K, Z, J) CNV 170 INTEGER K, J DOUBLE PRECISION Z C COMMON /EXR1/ LUERR COMMON /EXR3/ NLG102, MLG102, LG102(21) C C THE CONDITIONS IMPOSED ON NLG102, MLG102, AND LG102 BY C THIS SUBROUTINE ARE C C (1) NLG102 .GE. 2 C C (2) MLG102 .GE. 1 C C (3) 2*MLG102*(MLG102 - 1) .LE. 2**NBITS - 1 C C THESE CONDITIONS MUST BE MET BY APPROPRIATE CODING C IN SUBROUTINE SETUP. C IF (K.EQ.0) GO TO 70 M = MLG102 KA = IABS(K) KA1 = KA/M KA2 = MOD(KA,M) IF (KA1.GE.M) GO TO 60 NM1 = NLG102 - 1 NP1 = NLG102 + 1 IT = KA2*LG102(NP1) IC = IT/M ID = MOD(IT,M) Z = ID IF (KA1.GT.0) GO TO 20 DO 10 II=1,NM1 I = NP1 - II IT = KA2*LG102(I) + IC IC = IT/M ID = MOD(IT,M) Z = Z/DBLE(FLOAT(M)) + DBLE(FLOAT(ID)) 10 CONTINUE JA = KA*LG102(1) + IC GO TO 40 20 CONTINUE DO 30 II=1,NM1 I = NP1 - II IT = KA2*LG102(I) + KA1*LG102(I+1) + IC IC = IT/M ID = MOD(IT,M) Z = Z/DBLE(FLOAT(M)) + DBLE(FLOAT(ID)) 30 CONTINUE JA = KA*LG102(1) + KA1*LG102(2) + IC 40 CONTINUE Z = Z/DBLE(FLOAT(M)) IF (K.GT.0) GO TO 50 J = -JA Z = 10.0D0**(-Z) GO TO 80 50 CONTINUE J = JA + 1 Z = 10.0D0**(Z-1.0D0) GO TO 80 60 CONTINUE C THIS ERROR OCCURS IF K EXCEEDS MLG102**2 - 1 IN MAGNITUDE. C WRITE (LUERR,99999) STOP 70 CONTINUE J = 0 Z = 1.0D0 80 RETURN 99999 FORMAT (34H ***ERR IN CNV210...K TOO LARGE***) END C *** SUBROUTINE CONVRT *** CON 10 C USAGE CON 20 C CALL CONVRT(X,IX) CON 30 C DESCRIPTION CON 40 C CONVERTS (X,IX) = X*RADIX**IX CON 50 C TO DECIMAL FORM IN PREPARATION FOR CON 60 C PRINTING, SO THAT (X,IX) = X*10**IX CON 70 C WHERE 1/10 .LE. DABS(X) .LT. 1 CON 80 C IS RETURNED, EXCEPT THAT IF CON 90 C (DABS(X),IX) IS BETWEEN RADIX**(-2L) CON 100 C AND RADIX**(2L) THEN THE REDUCED CON 110 C FORM WITH IX = 0 IS RETURNED. CON 120 C CON 130 SUBROUTINE CONVRT(X, IX) CON 140 DOUBLE PRECISION X INTEGER IX C C THE CONDITIONS IMPOSED ON L AND KMAX BY THIS SUBROUTINE C ARE C (1) 4 .LE. L .LE. 2**NBITS - 1 - KMAX C C (2) KMAX .LE. ((2**NBITS)-2)/LOG10R - L C C THESE CONDITIONS MUST BE MET BY APPROPRIATE CODING C IN SUBROUTINE SETUP. C DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R COMMON /EXR2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX C DOUBLE PRECISION A, B, Z C DATA ISPACE /1/ C THE PARAMETER ISPACE IS THE INCREMENT USED IN FORM- C ING THE AUXILIARY INDEX OF THE DECIMAL EXTENDED-RANGE C FORM. THE RETURNED VALUE OF IX WILL BE AN INTEGER MULT- C IPLE OF ISPACE. ISPACE MUST SATISFY 1 .LE. ISPACE .LE. C L/2. IF A VALUE GREATER THAN 1 IS TAKEN, THE RETURNED C VALUE OF X WILL SATISFY 10**(-ISPACE) .LE. DABS(X) .LE. 1 C IF (DABS(X),IX) .LT. RADIX**(-2L), OR 1/10 .LE. DABS(X) C .LT. 10**(ISPACE-1) IF (DABS(X),IX) .GT. RADIX**(2L). C CALL REDUCE(X, IX) IF (IX.EQ.0) GO TO 150 CALL ADJUST(X, IX) C C CASE 1 IS WHEN (X,IX) IS LESS THAN RADIX**(-2L) IN MAGNITUDE, C CASE 2 IS WHEN (X,IX) IS GREATER THAN RADIX**(2L) IN MAGNITUDE. ITEMP = 1 ICASE = (3+ISIGN(ITEMP,IX))/2 GO TO (10, 20), ICASE 10 IF (DABS(X).LT.1.0D0) GO TO 30 X = X/RADIXL IX = IX + L GO TO 30 20 IF (DABS(X).GE.1.0D0) GO TO 30 X = X*RADIXL IX = IX - L 30 CONTINUE C C AT THIS POINT, RADIX**(-L) .LE. DABS(X) .LT. 1.0D0 IN CASE 1, C 1.0D0 .LE. DABS(X) .LT. RADIX**L IN CASE 2. I = DLOG10(DABS(X))/DLG10R A = RADIX**I GO TO (40, 60), ICASE 40 IF (A.LE.RADIX*DABS(X)) GO TO 50 I = I - 1 A = A/RADIX GO TO 40 50 IF (DABS(X).LT.A) GO TO 80 I = I + 1 A = A*RADIX GO TO 50 60 IF (A.LE.DABS(X)) GO TO 70 I = I - 1 A = A/RADIX GO TO 60 70 IF (DABS(X).LT.RADIX*A) GO TO 80 I = I + 1 A = A*RADIX GO TO 70 80 CONTINUE C C AT THIS POINT I IS SUCH THAT C RADIX**(I-1) .LE. DABS(X) .LT. RADIX**I IN CASE 1, C RADIX**I .LE. DABS(X) .LT. RADIX**(I+1) IN CASE 2. ITEMP = DBLE(FLOAT(ISPACE))/DLG10R A = RADIX**ITEMP B = 10.0D0**ISPACE 90 IF (A.LE.B) GO TO 100 ITEMP = ITEMP - 1 A = A/RADIX GO TO 90 100 IF (B.LT.A*RADIX) GO TO 110 ITEMP = ITEMP + 1 A = A*RADIX GO TO 100 110 CONTINUE C C AT THIS POINT ITEMP IS SUCH THAT C RADIX**ITEMP .LE. 10**ISPACE .LT. RADIX**(ITEMP+1). IF (ITEMP.GT.0) GO TO 120 C ITEMP = 0 IF, AND ONLY IF, ISPACE = 1 AND RADIX = 16.0D0 X = X*RADIX**(-I) IX = IX + I CALL CNV210(IX, Z, J) X = X*Z IX = J GO TO (130, 140), ICASE 120 CONTINUE I1 = I/ITEMP X = X*RADIX**(-I1*ITEMP) IX = IX + I1*ITEMP C C AT THIS POINT, C RADIX**(-ITEMP) .LE. DABS(X) .LT. 1.0D0 IN CASE 1, C 1.0D0 .LE. DABS(X) .LT. RADIX**ITEMP IN CASE 2. CALL CNV210(IX, Z, J) J1 = J/ISPACE J2 = J - J1*ISPACE X = X*Z*10.0D0**J2 IX = J1*ISPACE C C AT THIS POINT, C 10.0D0**(-2*ISPACE) .LE. DABS(X) .LT. 1.0D0 IN CASE 1, C 10.0D0**-1 .LE. DABS(X) .LT. 10.0D0**(2*ISPACE-1) IN CASE 2. GO TO (130, 140), ICASE 130 IF (B*DABS(X).GE.1.0D0) GO TO 150 X = X*B IX = IX - ISPACE GO TO 130 140 IF (10.0D0*DABS(X).LT.B) GO TO 150 X = X/B IX = IX + ISPACE GO TO 140 150 RETURN END C *** SUBROUTINE REDUCE *** RED 10 C USAGE RED 20 C CALL REDUCE(X,IX) RED 30 C DESCRIPTION RED 40 C IF RED 50 C RADIX**(-2L) .LE. (DABS(X),IX) .LE. RADIX**(2L) RED 60 C THEN REDUCE TRANSFORMS (X,IX) SO THAT IX=0. RED 70 C IF (X,IX) IS OUTSIDE THE ABOVE RANGE, RED 80 C THEN REDUCE TAKES NO ACTION. RED 90 C THIS SUBROUTINE IS USEFUL IF THE RED 100 C RESULTS OF EXTENDED-RANGE CALCULATIONS RED 110 C ARE TO BE USED IN SUBSEQUENT ORDINARY RED 120 C DOUBLE-PRECISION CALCULATIONS. RED 130 C RED 140 SUBROUTINE REDUCE(X, IX) RED 150 DOUBLE PRECISION X INTEGER IX C DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R COMMON /EXR2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX DOUBLE PRECISION XA C IF (X.EQ.0.0D0) GO TO 90 XA = DABS(X) IF (IX.EQ.0) GO TO 70 IXA = IABS(IX) IXA1 = IXA/L2 IXA2 = MOD(IXA,L2) IF (IX.GT.0) GO TO 40 10 CONTINUE IF (XA.GT.1.0D0) GO TO 20 XA = XA*RAD2L IXA1 = IXA1 + 1 GO TO 10 20 XA = XA/RADIX**IXA2 IF (IXA1.EQ.0) GO TO 70 DO 30 I=1,IXA1 IF (XA.LT.1.0D0) GO TO 100 XA = XA/RAD2L 30 CONTINUE GO TO 70 C 40 CONTINUE IF (XA.LT.1.0D0) GO TO 50 XA = XA/RAD2L IXA1 = IXA1 + 1 GO TO 40 50 XA = XA*RADIX**IXA2 IF (IXA1.EQ.0) GO TO 70 DO 60 I=1,IXA1 IF (XA.GT.1.0D0) GO TO 100 XA = XA*RAD2L 60 CONTINUE 70 IF (XA.GT.RAD2L) GO TO 100 IF (XA.GT.1.0D0) GO TO 80 IF (RAD2L*XA.LT.1.0D0) GO TO 100 80 X = DSIGN(XA,X) 90 IX = 0 100 RETURN END C *** SUBROUTINE EXRP *** EXR 10 C SUBROUTINE TO PRINT CONTENTS OF COMMON BLOCKS. EXR 20 SUBROUTINE EXRP EXR 30 COMMON /EXR1/ LUERR DOUBLE PRECISION RADIX, RADIXL, RAD2L, DLG10R COMMON /EXR2/ RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX COMMON /EXR3/ NLG102, MLG102, LG102(21) WRITE (6,99999) LUERR WRITE (6,99998) RADIX, RADIXL, RAD2L, DLG10R, L, L2, KMAX WRITE (6,99997) NLG102, MLG102, (LG102(I),I=1,NLG102), * LG102(NLG102+1) RETURN 99999 FORMAT (/25H CONTENTS OF COMMON/EXR1//9H LUERR =, I26) 99998 FORMAT (/25H CONTENTS OF COMMON/EXR2//9H RADIX =, * D26.18/9H RADIXL =, D26.18/9H RAD2L =, D26.18/9H DLG10R =, * D26.18/9H L =, I26/9H L2 =, I26/9H KMAX =, I26) 99997 FORMAT (/25H CONTENTS OF COMMON/EXR3//9H NLG102 =, I26/8H MLG102 , * 1H=, I26/9H LG102 =, 5I12/(9X, 5I12)) END