C ALGORITHM 719, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 19, NO. 3, SEPTEMBER, 1993, PP. 286-317. README File for the Multiprecision Package and Translator Programs 32-bit computer version Update as of June 30, 1992 The following files are included in this "shar" file: Name Description mpfun.f MPFUN library file testmp.f Test program for MPFUN transmp.f TRANSMP program testran.f Test program for TRANSMP testran.out Reference output of test program These programs should work correctly on any 32-bit system, including those based on IEEE arithmetic. For Crays or other 64-bit systems, request separate versions of these programs from the author. On Unix systems, the following sequence compiles the MPFUN library; compiles, links and executes the test program for MPFUN; compiles the TRANSMP translator program; translates the test program for TRANSMP; compiles, links and executes the translated program; and compares the output of the translated program with the reference output file: Command Notes f77 -c mpfun.f Insert -O2 or -O3 after f77 for an optimized compile. There should be no fatal compiler diagnostics. f77 testmp.f mpfun.o There should be no fatal compiler diagnostics. a.out > testmp.out Check testmp.out to make sure all tests passed. f77 -o transmp transmp.f There should be no fatal compiler diagnostics. transmp < testran.f > tranout.f Check the end of tranout.f to make sure there are no fatal translator errors. f77 tranout.f mpfun.o There should be no fatal compiler diagnostics. a.out > tranout.out This should complete normally. diff tranout.out testran.out There should be no differences here. On IBM workstations, the four f77 lines should be replaced by xlf -O -c mpfun.f xlf testmp.f mpfun.o xlf -qrecur -qcharlen=1600 -o transmp transmp.f xlf tranout.f mpfun.o On HP workstations, comment out the two IMPLICIT AUTOMATIC statements in the transmp.f file before compiling. On DEC workstations, comment out these two lines and compile with the command f77 -automatic -o transmp transmp.f. These codes have been tested quite thoroughly, but a few bugs may remain. If you encounter any, please let me know and I will fix them as soon as possible. David H. Bailey NASA Ames Research Center Mail Stop T045-1 Moffett Field, CA 94035 Tel.: 415-604-4410 Fax: 415-604-3957 E-mail: dbailey@nas.nasa.gov 10 ^ 0 x 3.14159265358979323846264338327950288419716939937510582, 10 ^ 0 x 2.718281828459045235360287471352662497757247093699959574, 10 ^ 8 x 1.5929408, 10 ^ 14 x 1.52779903171104, 10 ^ 7 x 1.798869166210583847918164697465982392760106312722056703, 10 ^ 6 x 4.246554058540655594271278713230977821918624216872862863, 10 ^ 6 x 1.8734, 10 ^ 6 x 3.7468, 10 ^ 0 x 1.78661346435546875, 10 ^ 0 x 2.947257281471990863100813367572530600175591224440627617, 10 ^ 1 x 1.8, 10 ^ 1 x -2.492030753469781077845780729842283661121425243866422996, 10 ^ 0 x 0., C***************************************************************************** C C MPFUN: A MULTIPLE PRECISION FLOATING POINT COMPUTATION PACKAGE C C Standard version C Version Date: June 30, 1992 C C Author: C C David H. Bailey Telephone: 415-604-4410 C NASA Ames Research Center Facsimile: 415-604-3957 C Mail Stop T045-1 Internet: dbailey@nas.nasa.gov C Moffett Field, CA 94035 C USA C C Restrictions: C C This software has now been approved by NASA for unrestricted distribution. C However, usage of this software is subject to the following: C C 1. This software is offered without warranty of any kind, either expressed C or implied. The author would appreciate, however, any reports of bugs C or other difficulties that may be encountered. C 2. If modifications or enhancements to this software are made to this C software by others, NASA Ames reserves the right to obtain this enhanced C software at no cost and with no restrictions on its usage. C 3. The author and NASA Ames are to be acknowledged in any published paper C based on computations using this software. Accounts of practical C applications or other benefits resulting from this software are of C particular interest. Please send a copy of such papers to the author. C C Description: C C The following information is a brief description of this program. For C full details and instructions for usage, see the paper "A Portable High C Performance Multiprecision Package", available from the author. C C This package of Fortran subroutines performs multiprecision floating point C arithmetic. If sufficient main memory is available, the maximum precision C level is at least 16 million digits. The maximum dynamic range is at C least 10^(+-14,000,000). It employs advanced algorithms, including an C FFT-based multiplication routine and some recently discovered C quadratically convergent algorithms for pi, exp and log. The package also C features extensive debug and self-checking facilities, so that it can be C used as a rigorous system integrity test. All of the routines in this C package have been written to facilitate vector and parallel processing. C C For users who do not wish to manually write code that calls these routines, C an automatic translator program is available from the author that converts C ordinary Fortran-77 code into code that calls these routines. Contact the C author for details. C C This package should run correctly on any computer with a Fortran-77 C compiler that meets certain minimal floating point accuracy standards. C Any system based on the IEEE floating point standard, with a 25 bit C mantissa in single precision and a 53 bit mantissa in double precision, C easily meets these requirements. All DEC VAX systems meet these C requirements. All IBM mainframes and workstations meet these requirements. C Cray systems meet all of these requirements with double precision disabled C (i.e. by using only single precision). C C Machine-specific tuning notes may be located by searching for the text C string C> in this program file. It is highly recommended that these notes C be read before running this package on a specific system. If no comment C accompanies a C> string, this indicates that all references to INT in the C next loop may be safely changed to AINT. INT appears to be significantly C faster on many 32 bit systems, but AINT is slightly faster on Crays. Also, C certain vectorizable DO loops that are often not recognized as such by C vectorizing compilers are prefaced with Cray CDIR$ IVDEP directives. On C other vector systems these directives should be replaced by the C appropriate equivalents. C C Instructions for compiling and testing this program are included in the C readme file that accompanies this file. C C***************************************************************************** C BLOCK DATA C C This initializes the parameters in MPCOM1 and the error codes in MPCOM2 C with default values. C> C On IEEE systems and most other 32 bit systems, set BBXC = 4096.D0, C NBTC = 24, NPRC = 32, and MCRC = 7. On Cray systems, set BBXC = 2048.D0, C NBTC = 22, NPRC = 16, and MCRC = 8. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C PARAMETER (BBXC = 4096.D0, NBTC = 24, NPRC = 32, MCRC = 7, $ BDXC = BBXC ** 2, BX2C = BDXC ** 2, RBXC = 1.D0 / BBXC, $ RDXC = RBXC ** 2, RX2C = RDXC ** 2, RXXC = 16.D0 * RX2C) DATA BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR / BBXC, BDXC, $ BX2C, RBXC, RDXC, RX2C, RXXC, NBTC, NPRC/ DATA NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS / $ 16, 0, 6, 0, MCRC, 1, 1, 1, 1024/ DATA KER /72 * 2/ END C SUBROUTINE DPADD (A, NA, B, NB, C, NC) C C This adds the DPE numbers (A, NA) and (B, NB) to yield the sum (C, NC). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION PT(64) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS SAVE PT DATA PT/ 64 * 0.D0/ C IF (IER .NE. 0) THEN C = 0.D0 NC = 0 RETURN ENDIF C C If this is the first call to DPADD, initialize the PT table. C IF (PT(1) .EQ. 0.D0) THEN PT(1) = 0.5D0 C DO 100 I = 2, 64 PT(I) = 0.5D0 * PT(I-1) 100 CONTINUE C ENDIF C C This operation reduces to five cases. C IF (B .EQ. 0.D0) THEN C = A NC = NA ELSE IF (A .EQ. 0.D0) THEN C = B NC = NB ELSE IF (NA .EQ. NB) THEN C = A + B NC = NA ELSE IF (NA .GT. NB) THEN K = NA - NB NC = NA IF (K .GT. 64) THEN C = A ELSE C = A + B * PT(K) ENDIF ELSE K = NB - NA NC = NB IF (K .GT. 64) THEN C = B ELSE C = B + A * PT(K) ENDIF ENDIF IF (C .EQ. 0.D0) THEN NC = 0 GOTO 130 ENDIF C C Normalize the result to a decent range if it is not. C 110 IF (ABS (C) .GE. BDX) THEN C = RDX * C NC = NC + NBT GOTO 110 ENDIF C 120 IF (ABS (C) .LT. 1.D0) THEN C = BDX * C NC = NC - NBT GOTO 120 ENDIF C 130 RETURN END C SUBROUTINE DPDEC (A, NA, B, NB) C C This converts the DPE number (A, NA) to decimal form, i.e. B * 10^NB, C where |B| is between 1 and 10. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (XLT = 0.3010299956639812D0) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C IF (A .NE. 0.D0) THEN T1 = XLT * NA + LOG10 (ABS (A)) NB = T1 IF (T1 .LT. 0.D0) NB = NB - 1 B = SIGN (10.D0 ** (T1 - NB), A) ELSE B = 0.D0 NB = 0 ENDIF C RETURN END C SUBROUTINE DPDIV (A, NA, B, NB, C, NC) C C This divides the DPE number (A, NA) by (B, NB) to yield the quotient C (C, NC). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C IF (IER .NE. 0) THEN C = 0.D0 NC = 0 RETURN ENDIF IF (B .EQ. 0.D0) THEN IF (KER(1) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** DPDIV: Divisor is zero.') IER = 1 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Divide A by B and subtract exponents, unless A is zero. C IF (A .EQ. 0.D0) THEN C = 0.D0 NC = 0 GOTO 120 ELSE C = A / B NC = NA - NB ENDIF C C Normalize the result to a decent range if it is not. C 100 IF (ABS (C) .GE. BDX) THEN C = RDX * C NC = NC + NBT GOTO 100 ENDIF C 110 IF (ABS (C) .LT. 1.D0) THEN C = BDX * C NC = NC - NBT GOTO 110 ENDIF C 120 RETURN END C SUBROUTINE DPMUL (A, NA, B, NB, C, NC) C C This multiplies the DPE number (A, NA) by (B, NB) to yield the product C (C, NC). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS C IF (IER .NE. 0) THEN C = 0.D0 NC = 0 RETURN ENDIF C C Multiply A by B and add exponents, unless either is zero. C IF (A .EQ. 0.D0 .OR. B .EQ. 0.D0) THEN C = 0.D0 NC = 0 GOTO 120 ELSE C = A * B NC = NA + NB ENDIF C C Normalize the result to a decent range if it is not. C 100 IF (ABS (C) .GE. BDX) THEN C = RDX * C NC = NC + NBT GOTO 100 ENDIF C 110 IF (ABS (C) .LT. 1.D0) THEN C = BDX * C NC = NC - NBT GOTO 110 ENDIF C 120 RETURN END C SUBROUTINE DPPWR (A, NA, B, NB, C, NC) C C This raises the DPE number (A, NA) to the (B, NB) power and places the C result in (C, NC). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) PARAMETER (CL2 = 1.4426950408889633D0) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C IF (IER .NE. 0) THEN C = 0.D0 NC = 0 RETURN ENDIF IF (A .LE. 0.D0) THEN IF (KER(2) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** DPPWR: Argument is less than or equal to zero.') IER = 2 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C IF (B .EQ. 0.D0) THEN C = 1.D0 NC = 0 GOTO 120 ENDIF C IF (B .EQ. 1.D0 .AND. NB .EQ. 0) THEN C = A NC = NA GOTO 120 ENDIF C C Compute the base 2 logarithm of A and multiply by B. C AL = CL2 * LOG (A) + NA CALL DPMUL (AL, 0, B, NB, T1, N1) C C Check for possible overflow or underflow. C IF (N1 .GT. 6) THEN IF (T1 .GT. 0.D0) THEN IF (KER(3) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** DPPWR: Overflow') IER = 3 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ELSE C = 0.D0 NC = 0 GOTO 120 ENDIF ENDIF C C Compute 2 raised to the power B * Log_2 (A). C T1 = T1 * 2.D0 ** N1 NC = INT (T1) C = 2.D0 ** (T1 - NC) C C Normalize the result to a decent range if it is not. C 100 IF (ABS (C) .GE. BDX) THEN C = RDX * C NC = NC + NBT GOTO 100 ENDIF C 110 IF (ABS (C) .LT. 1.D0) THEN C = BDX * C NC = NC - NBT GOTO 110 ENDIF C 120 RETURN END C SUBROUTINE DPSQRT (A, NA, B, NB) C C This computes the square root of the DPE number (A, NA) and places the C result in (B, NB). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C IF (IER .NE. 0) THEN B = 0.D0 NB = 0 RETURN ENDIF IF (A .LT. 0.D0) THEN IF (KER(4) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** DPSQRT: Argument is negative.') IER = 4 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C IF (A .EQ. 0.D0) THEN B = 0.D0 NB = 0 GOTO 120 ENDIF C C Divide the exponent of A by two and then take the square root of A. If C NA is not an even number, then we have to multiply A by 10 before taking C the square root. C NB = NA / 2 IF (NA .EQ. 2 * NB) THEN B = SQRT (A) ELSE B = SQRT (2.D0 * A) IF (NA .LT. 0) NB = NB - 1 ENDIF C C Normalize the result to a decent range if it is not. C 100 IF (ABS (B) .GE. BDX) THEN B = RDX * B NB = NB + NBT GOTO 100 ENDIF C 110 IF (ABS (B) .LT. 1.D0) THEN B = BDX * B NB = NB - NBT GOTO 110 ENDIF C 120 RETURN END C SUBROUTINE DPSUB (A, NA, B, NB, C, NC) C C This subtracts the DPE number (B, NB) from (A, NA) to yield the difference C (C, NC). C IMPLICIT DOUBLE PRECISION (A-H, O-Z) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS C IF (IER .NE. 0) THEN C = 0.D0 NC = 0 RETURN ENDIF C BB = -B CALL DPADD (A, NA, BB, NB, C, NC) C RETURN END C SUBROUTINE MPABRT C> C This routine terminates execution. Many users will want to replace the C default STOP with a call to a system routine that provides a traceback. C Examples of code that produce traceback are included here (commented out) C for some systems. C COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS C WRITE (LDB, 1) IER 1 FORMAT ('*** MPABRT: Execution terminated, error code =',I4) C C Use this line on Cray systems. C C CALL ABORT C C Use this line plus the C routine TRACBK (available from author) on C Silicon Graphics IRIS systems. C C CALL TRACBK C C On other systems, merely terminate execution. C STOP END C SUBROUTINE MPADD (A, B, C) C C This routine adds MP numbers A and B to yield the MP sum C. It attempts C to include all significance of A and B in the result, up to the maximum C mantissa length NW. Debug output starts with IDB = 9. C C Max SP space for C: NW + 4 cells. Max DP scratch space: NW + 4 cells. C DOUBLE PRECISION D PARAMETER (NDB = 22) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM4/ D(1024) DIMENSION A(NW+2), B(NW+2), C(NW+4) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. RETURN ENDIF IF (IDB .GE. 9) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPADD I'/(6F12.0)) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 1) (B(I), I = 1, NO) ENDIF C IA = SIGN (1., A(1)) IB = SIGN (1., B(1)) NA = MIN (INT (ABS (A(1))), NW) NB = MIN (INT (ABS (B(1))), NW) C C This first IF block checks for zero inputs. C IF (NA .EQ. 0) THEN C C A is zero -- the result is B. C C(1) = SIGN (NB, IB) C DO 100 I = 2, NB + 2 C(I) = B(I) 100 CONTINUE C GOTO 420 ELSEIF (NB .EQ. 0) THEN C C B is zero -- the result is A. C C(1) = SIGN (NA, IA) C DO 110 I = 1, NA + 2 C(I) = A(I) 110 CONTINUE C GOTO 420 ENDIF MA = A(2) MB = B(2) C C This IF block breaks the problem into different branches depending on C the relative sizes of the exponents of A and B. C IF (MA .EQ. MB) THEN C C A and B have the same exponent. C NM = MIN (NA, NB) NX = MAX (NA, NB) IF (IA .EQ. IB) THEN C C A and B have the same exponent and sign. C D(1) = SIGN (NX, IA) D(2) = MA D(NX+3) = 0.D0 D(NX+4) = 0.D0 C DO 120 I = 3, NM + 2 D(I) = DBLE (A(I)) + DBLE (B(I)) 120 CONTINUE C IF (NA .GT. NB) THEN C C A is longer than B -- include extra words of A in C. C DO 130 I = NM + 3, NA + 2 D(I) = A(I) 130 CONTINUE C ELSEIF (NB .GT. NA) THEN C C B is longer than A -- include extra words of B in C. C DO 140 I = NM + 3, NB + 2 D(I) = B(I) 140 CONTINUE C ENDIF ELSE C C A and B have the same exponent but the opposite sign. It is thus C necessary to scan through each vector until we find an unequal word. C DO 150 I = 3, NM + 2 IF (A(I) .NE. B(I)) GOTO 180 150 CONTINUE C C All words up to the common length are equal. C IF (NA .EQ. NB) THEN C C The length of A is the same as B -- result is zero. C C(1) = 0.D0 C(2) = 0.D0 GOTO 420 ELSEIF (NA .GT. NB) THEN C C A is longer -- thus trailing words of A are shifted to start of C. C NN = NA - NB D(1) = SIGN (NN, IA) D(2) = A(2) - NB D(NN+3) = 0.D0 D(NN+4) = 0.D0 C DO 160 I = 3, NN + 2 D(I) = A(I+NB) 160 CONTINUE C ELSEIF (NB .GT. NA) THEN C C B is longer -- thus trailing words of B are shifted to start of C. C NN = NB - NA D(1) = SIGN (NN, IB) D(2) = B(2) - NA D(NN+3) = 0.D0 D(NN+4) = 0.D0 C DO 170 I = 3, NN + 2 D(I) = B(I+NA) 170 CONTINUE C ENDIF GOTO 410 C C An unequal word was found. C 180 K = I - 3 IF (A(K+3) .GT. B(K+3)) THEN C C A is larger -- subtract B (shifted) from A. C D(1) = SIGN (NX - K, IA) D(2) = A(2) - K D(NX-K+3) = 0.D0 D(NX-K+4) = 0.D0 C DO 190 I = 3, NM - K + 2 D(I) = DBLE (A(I+K)) - DBLE (B(I+K)) 190 CONTINUE C DO 200 I = NB - K + 3, NA - K + 2 D(I) = A(I+K) 200 CONTINUE C DO 210 I = NA - K + 3, NB - K + 2 D(I) = - B(I+K) 210 CONTINUE C ELSE C C B is larger -- subtract A (shifted) from B. C D(1) = SIGN (NX - K, IB) D(2) = B(2) - K D(NX-K+3) = 0.D0 D(NX-K+4) = 0.D0 C DO 220 I = 3, NM - K + 2 D(I) = DBLE (B(I+K)) - DBLE (A(I+K)) 220 CONTINUE C DO 230 I = NB - K + 3, NA - K + 2 D(I) = - A(I+K) 230 CONTINUE C DO 240 I = NA - K + 3, NB - K + 2 D(I) = B(I+K) 240 CONTINUE C ENDIF ENDIF ELSEIF (MA .GT. MB) THEN C C Exponent of A is greater. In other words, A has a larger magnitude. C MC = MA - MB LA = MIN (MC, NA) LB = MIN (MC + NB, NW + 2) LM = MIN (NA, LB) LX = MIN (MAX (NA, LB), NW) D(1) = SIGN (LX, IA) D(2) = A(2) D(LX+3) = 0.D0 D(LX+4) = 0.D0 C DO 250 I = 3, LA + 2 D(I) = A(I) 250 CONTINUE C C If B is shifted NW + 2 or more words to the right of A then C = A. C IF (MC .GE. NW + 2) THEN D(1) = SIGN (NA, IA) GOTO 410 ENDIF IF (MC .GT. NA) THEN C C There is a gap between A and the shifted B. Fill it with zeroes. C DO 260 I = NA + 3, MC + 2 D(I) = 0.D0 260 CONTINUE C LM = MC ENDIF IF (IA .EQ. IB) THEN C C A and B have the same sign -- add common words with B shifted right. C DO 270 I = MC + 3, LM + 2 D(I) = DBLE (A(I)) + DBLE (B(I-MC)) 270 CONTINUE C C Include tail of A or B, whichever is longer after shift. C IF (NA .GT. LB) THEN C DO 280 I = LM + 3, NA + 2 D(I) = A(I) 280 CONTINUE C ELSE C DO 290 I = LM + 3, LB + 2 D(I) = B(I-MC) 290 CONTINUE C ENDIF ELSE C C A and B have different signs -- subtract common words with B shifted right. C DO 300 I = MC + 3, LM + 2 D(I) = DBLE (A(I)) - DBLE (B(I-MC)) 300 CONTINUE C C Include tail of A or B, whichever is longer after shift. C DO 310 I = LM + 3, NA + 2 D(I) = A(I) 310 CONTINUE C DO 320 I = LM + 3, LB + 2 D(I) = - B(I-MC) 320 CONTINUE C ENDIF ELSE C C Exponent of B is greater. In other words, B has a larger magnitude. C MC = MB - MA LB = MIN (MC, NB) LA = MIN (MC + NA, NW + 2) LM = MIN (NB, LA) LX = MIN (MAX (NB, LA), NW) D(1) = SIGN (LX, IB) D(2) = B(2) D(LX+3) = 0.D0 D(LX+4) = 0.D0 C DO 330 I = 3, LB + 2 D(I) = B(I) 330 CONTINUE C C If A is shifted NW + 2 or more words to the right of B then C = B. C IF (MC .GE. NW + 2) THEN D(1) = SIGN (NB, IB) GOTO 410 ENDIF IF (MC .GT. NB) THEN C C There is a gap between B and the shifted A. Fill it with zeroes. C DO 340 I = NB + 3, MC + 2 D(I) = 0.D0 340 CONTINUE C LM = MC ENDIF IF (IB .EQ. IA) THEN C C B and A have the same sign -- add common words with A shifted right. C DO 350 I = MC + 3, LM + 2 D(I) = DBLE (B(I)) + DBLE (A(I-MC)) 350 CONTINUE C C Include tail of B or A, whichever is longer after shift. C DO 360 I = LM + 3, NB + 2 D(I) = B(I) 360 CONTINUE C DO 370 I = LM + 3, LA + 2 D(I) = A(I-MC) 370 CONTINUE C ELSE C C B and A have different signs -- subtract common words with A shifted right. C DO 380 I = MC + 3, LM + 2 D(I) = DBLE (B(I)) - DBLE (A(I-MC)) 380 CONTINUE C C Include tail of B or A, whichever is longer after shift. C DO 390 I = LM + 3, NB + 2 D(I) = B(I) 390 CONTINUE C DO 400 I = LM + 3, LA + 2 D(I) = - A(I-MC) 400 CONTINUE C ENDIF ENDIF C C Fix up result, since some words may be negative or exceed BDX. C 410 CALL MPNORM (C) C 420 IF (IDB .GE. 9) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 2) (C(I), I = 1, NO) 2 FORMAT ('MPADD O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPAGMX (A, B) C C This performs the arithmetic-geometric mean (AGM) iterations. This routine C is called by MPLOGX. It is not intended to be called directly by the user. C C Max SP space for A and B: NW + 4 cells. Max SP scratch space: 5 * NW + 26 C cells. Max DP scratch space: 12 * NW + 6 cells. C COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) DIMENSION A(NW+4), B(NW+4) C IF (IER .NE. 0) THEN A(1) = 0. A(2) = 0. B(1) = 0. B(2) = 0. RETURN ENDIF N4 = NW + 4 NS = 2 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N4 S(K0) = 0. S(K0+1) = 0. L1 = 0 C 100 L1 = L1 + 1 IF (L1 .EQ. 50) THEN IF (KER(5) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPAGMX: Iteration limit exceeded.') IER = 5 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF ENDIF C S1 = S(K0+1) CALL MPADD (A, B, S(K0)) CALL MPMULD (S(K0), 0.5D0, 0, S(K1)) CALL MPMULX (A, B, S(K0)) CALL MPSQRX (S(K0), B) CALL MPEQ (S(K1), A) CALL MPSUB (A, B, S(K0)) C C Check for convergence. C IF (S(K0) .NE. 0. .AND. (S(K0+1) .NE. S1 .OR. S(K0+1) .GE. -2)) $ GOTO 100 C ICS = ISS IF (IDB .GE. 6) WRITE (LDB, 2) L1, S(K0+1) 2 FORMAT ('MPAGMX: Iter., Tol. Achieved =',I5,F8.0) RETURN END C SUBROUTINE MPALER C C This outputs error messages when a single precision scratch space C allocation error is detected. C COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C IF (KER(6) .NE. 0) THEN WRITE (LDB, 1) ICS - 1 1 FORMAT ('*** MPALER: Insufficient single precision scratch ', $ 'space.'/ 'Allocate',I10,' cells in an array in common ', $ 'MPCOM3 of the main '/ 'program and set IMS in common ', $ 'MPCOM1 to this size.') IER = 6 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF C RETURN END C SUBROUTINE MPANG (X, Y, PI, A) C C This computes the MP angle A subtended by the MP pair (X, Y) considered as C a point in the x-y plane. This is more useful than an arctan or arcsin C routine, since it places the result correctly in the full circle, i.e. C -Pi < A <= Pi. PI is the MP value of Pi computed by a previous call to C MPPI. For extra high levels of precision, use MPANGX. The last word of C the result is not reliable. Debug output starts with IDB = 5. C C Max SP space for A: NW + 4 cells. Max SP scratch space: 14 * NW + 81 C cells. Max DP scratch space: NW + 7 cells. C C The Taylor series for Sin converges much more slowly than that of Arcsin. C Thus this routine does not employ Taylor series, but instead computes C Arccos or Arcsin by solving Cos (a) = x or Sin (a) = y using one of the C following Newton iterations, both of which converge to a: C C z_{k+1} = z_k - [x - Cos (z_k)] / Sin (z_k) C z_{k+1} = z_k + [y - Sin (z_k)] / Cos (z_k) C C The first is selected if Abs (x) <= Abs (y); otherwise the second is used. C These iterations are performed with a maximum precision level NW that C is dynamically changed, approximately doubling with each iteration. C See the comment about the parameter NIT in MPDIVX. C DOUBLE PRECISION CL2, CPI, T1, T2, T3 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, CPI = 3.141592653589793D0, $ NIT = 3) DIMENSION A(NW+4), PI(NW+2), X(NW+2), Y(NW+2) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN A(1) = 0. A(2) = 0. RETURN ENDIF IF (IDB .GE. 5) THEN CALL MPDEB ('MPANG I', X) CALL MPDEB ('MPANG I', Y) ENDIF C IX = SIGN (1., X(1)) NX = MIN (INT (ABS (X(1))), NW) IY = SIGN (1., Y(1)) NY = MIN (INT (ABS (Y(1))), NW) C C Check if both X and Y are zero. C IF (NX .EQ. 0 .AND. NY .EQ. 0) THEN IF (KER(7) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPANG: Both arguments are zero.') IER = 7 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if Pi has been precomputed. C CALL MPMDC (PI, T1, N1) IF (N1 .NE. 0 .OR. ABS (T1 - CPI) .GT. RX2) THEN IF (KER(8) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPANG: PI must be precomputed.') IER = 8 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if one of X or Y is zero. C IF (NX .EQ. 0) THEN IF (IY .GT. 0) THEN CALL MPMULD (PI, 0.5D0, 0, A) ELSE CALL MPMULD (PI, -0.5D0, 0, A) ENDIF GOTO 120 ELSEIF (NY .EQ. 0) THEN IF (IX .GT. 0) THEN A(1) = 0. A(2) = 0. ELSE CALL MPEQ (PI, A) ENDIF GOTO 120 ENDIF C N5 = NW + 5 NS = 5 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N5 K2 = K1 + N5 K3 = K2 + N5 K4 = K3 + N5 NWS = NW NW = NW + 1 C C Determine the least integer MQ such that 2 ^ MQ .GE. NW. C T1 = NWS MQ = CL2 * LOG (T1) + 1.D0 - RXX C C Normalize x and y so that x^2 + y^2 = 1. C CALL MPMUL (X, X, S(K0)) CALL MPMUL (Y, Y, S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPSQRT (S(K2), S(K3)) CALL MPDIV (X, S(K3), S(K1)) CALL MPDIV (Y, S(K3), S(K2)) C C Compute initial approximation of the angle. C CALL MPMDC (S(K1), T1, N1) CALL MPMDC (S(K2), T2, N2) N1 = MAX (N1, -66) N2 = MAX (N2, -66) T1 = T1 * 2.D0 ** N1 T2 = T2 * 2.D0 ** N2 T3 = ATAN2 (T2, T1) CALL MPDMC (T3, 0, A) C C The smaller of x or y will be used from now on to measure convergence. C This selects the Newton iteration (of the two listed above) that has the C largest denominator. C IF (ABS (T1) .LE. ABS (T2)) THEN KK = 1 CALL MPEQ (S(K1), S(K0)) ELSE KK = 2 CALL MPEQ (S(K2), S(K0)) ENDIF C NW = 3 IQ = 0 C C Perform the Newton-Raphson iteration described above with a dynamically C changing precision level NW (one greater than powers of two). C DO 110 K = 2, MQ NW = MIN (2 * NW - 2, NWS) + 1 100 CONTINUE CALL MPCSSN (A, PI, S(K1), S(K2)) IF (KK .EQ. 1) THEN CALL MPSUB (S(K0), S(K1), S(K3)) CALL MPDIV (S(K3), S(K2), S(K4)) CALL MPSUB (A, S(K4), S(K1)) ELSE CALL MPSUB (S(K0), S(K2), S(K3)) CALL MPDIV (S(K3), S(K1), S(K4)) CALL MPADD (A, S(K4), S(K1)) ENDIF CALL MPEQ (S(K1), A) IF (K .EQ. MQ - NIT .AND. IQ .EQ. 0) THEN IQ = 1 GOTO 100 ENDIF 110 CONTINUE C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (A) C 120 IF (IDB .GE. 5) CALL MPDEB ('MPANG O', A) C RETURN END C SUBROUTINE MPANGX (X, Y, PI, A) C C This computes the MP angle A subtended by the MP pair (X, Y) considered as C a point in the x-y plane. This is more useful than an arctan or arcsin C routine, since it places the result correctly in the full circle, i.e. C -Pi < A <= Pi. PI is the MP value of Pi computed by a previous call to C MPPI or MPPIX. Before calling MPANGX, the array in MPCOM5 must be C initialized by calling MPINIX. For modest levels of precision, use MPANG. C NW should be a power of two. The last three words of the result are not C reliable. Debug output starts with IDB = 6. C C Max SP space for A: NW + 4 cells. Max SP scratch space: 18 * NW + 78 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine employs a complex arithmetic version of the MPLOGX alogirthm. C DOUBLE PRECISION CPI, T1 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CPI = 3.141592653589793D0) DIMENSION A(NW+4), F0(8), F1(8), F4(8), PI(NW+2), X(NW+2), Y(NW+2) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN A(1) = 0. A(2) = 0. RETURN ENDIF IF (IDB .GE. 6) THEN CALL MPDEB ('MPANGX I', X) CALL MPDEB ('MPANGX I', Y) ENDIF C IX = SIGN (1., X(1)) NX = MIN (INT (ABS (X(1))), NW) IY = SIGN (1., Y(1)) NY = MIN (INT (ABS (Y(1))), NW) NCR = 2 ** MCR C C Check if precision level is too low to justify the advanced routine. C IF (NW .LE. NCR) THEN CALL MPANG (X, Y, PI, A) GOTO 100 ENDIF C C Check if both X and Y are zero. C IF (NX .EQ. 0 .AND. NY .EQ. 0) THEN IF (KER(9) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPANGX: Both arguments are zero.') IER = 9 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if Pi has been precomputed. C CALL MPMDC (PI, T1, N1) IF (N1 .NE. 0 .OR. ABS (T1 - CPI) .GT. RX2) THEN IF (KER(10) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPANGX: PI must be precomputed.') IER = 10 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if one of X or Y is zero. C IF (NX .EQ. 0) THEN IF (IY .GT. 0) THEN CALL MPMULD (PI, 0.5D0, 0, A) ELSE CALL MPMULD (PI, -0.5D0, 0, A) ENDIF GOTO 100 ELSEIF (NY .EQ. 0) THEN IF (IX .GT. 0) THEN A(1) = 0. A(2) = 0. ELSE CALL MPEQ (PI, A) ENDIF GOTO 100 ENDIF C C Define scratch space. C N4 = NW + 4 N42 = 2 * N4 NS = 4 * N42 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N42 K2 = K1 + N42 K3 = K2 + N42 F0(1) = 0. F0(2) = 0. F0(3) = 0. F1(1) = 1. F1(2) = 0. F1(3) = 1. F4(1) = 1. F4(2) = 0. F4(3) = 4. C C Multiply the input by a large power of two. C CALL MPMDC (X, T1, N1) N2 = NBT * (NW / 2 + 2) - N1 TN = N2 CALL MPMULD (X, 1.D0, N2, S(K1)) CALL MPMULD (Y, 1.D0, N2, S(K2)) CALL MPMMPC (S(K1), S(K2), N4, S(K0)) C C Perform AGM iterations. C CALL MPMMPC (F1, F0, N4, S(K1)) CALL MPMMPC (F4, F0, N4, S(K3)) CALL MPCDVX (N4, S(K3), S(K0), S(K2)) CALL MPCAGX (S(K1), S(K2)) C C Compute A = Imag (Pi / (2 * Z)), where Z is the limit of the complex AGM. C CALL MPMULD (S(K1), 2.D0, 0, S(K0)) CALL MPMULD (S(K1+N4), 2.D0, 0, S(K0+N4)) CALL MPMMPC (PI, F0, N4, S(K2)) CALL MPCDVX (N4, S(K2), S(K0), S(K1)) CALL MPEQ (S(K1+N4), A) ICS = ISS C 100 IF (IDB .GE. 6) CALL MPDEB ('MPANGX O', A) C RETURN END C SUBROUTINE MPCADD (L, A, B, C) C C This computes the sum of the MPC numbers A and B and returns the MPC C result in C. L is the offset between real and imaginary parts in A, B C and C. L must be at least NW + 4. Debug output starts with IDB = 9. C C Max SP space for C: 2 * L cells. C DIMENSION A(2*L), B(2*L), C(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. C(L+1) = 0. C(L+2) = 0. RETURN ENDIF IF (IDB .GE. 9) WRITE (LDB, 1) 1 FORMAT ('MPCADD') C IF (L .LT. NW + 4) THEN IF (KER(11) .NE. 0) THEN WRITE (LDB, 2) L, NW + 4 2 FORMAT ('*** MPCADD: Offset parameter is too small',2I8) IER = 11 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C L1 = L + 1 CALL MPADD (A, B, C) CALL MPADD (A(L1), B(L1), C(L1)) C RETURN END C SUBROUTINE MPCAGX (A, B) C C This performs the arithmetic-geometric mean (AGM) iterations. This routine C is called by MPANGX. It is not intended to be called directly by the user. C C Max SP space for A and B: 2*NW + 8 cells. Max SP scratch space: 10*NW + 46 C cells. Max DP scratch space: 12 * NW + 6 cells. C COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) DIMENSION A(2*NW+8), B(2*NW+8) C IF (IER .NE. 0) THEN A(1) = 0. A(2) = 0. B(1) = 0. B(2) = 0. RETURN ENDIF N4 = NW + 4 NS = 4 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + 2 * N4 S(K0) = 0. S(K0+1) = 0. L1 = 0 C 100 L1 = L1 + 1 IF (L1 .EQ. 50) THEN IF (KER(12) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPCAGX: Iteration limit exceeded.') IER = 12 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF ENDIF C S1 = S(K0+1) CALL MPCADD (N4, A, B, S(K0)) CALL MPMULD (S(K0), 0.5D0, 0, S(K1)) CALL MPMULD (S(K0+N4), 0.5D0, 0, S(K1+N4)) CALL MPCMLX (N4, A, B, S(K0)) CALL MPCSQX (N4, S(K0), B) CALL MPCEQ (N4, S(K1), A) CALL MPSUB (A, B, S(K0)) C C Check for convergence. C IF (S(K0) .NE. 0. .AND. (S(K0+1) .NE. S1 .OR. S(K0+1) .GE. -2)) $ GOTO 100 C ICS = ISS IF (IDB .GE. 6) WRITE (LDB, 2) L1, S(K0+1) 2 FORMAT ('MPCAGX: Iter., Tol. Achieved =',I5,F8.0) RETURN END C SUBROUTINE MPCBRT (A, B) C C This computes the cube root of the MP number A and returns the MP result C in B. For extra high levels of precision, use MPCBRX. Debug output C starts with IDB = 7. C C Max SP space for B: NW + 4 cells. Max SP scratch space: 3 * NW + 15 C cells. Max DP scratch space: NW + 5 cells. C C This subroutine employs the following Newton-Raphson iteration, which C converges to A ^ (-2/3): C C X_{n+1} = X_k + (X_k / 3) * (1 - A^2 * X_k^3) C C Multiplying the final approximation to A ^ (-2/3) by A gives the cube C root. These iterations are performed with a maximum precision level NW that C is dynamically changed, approximately doubling with each iteration. C See the comment about the parameter NIT in MPDIVX. C DOUBLE PRECISION CL2, T1, T2 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, NDB = 22, NIT = 3) DIMENSION A(NW+2), B(NW+4), F(8) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. RETURN ENDIF IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPCBRT I'/(6F12.0)) ENDIF C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) C IF (NA .EQ. 0) THEN B(1) = 0. B(2) = 0. GOTO 120 ENDIF IF (IA .LT. 0.D0) THEN IF (KER(13) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPCBRT: Argument is negative.') IER = 13 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N5 = NW + 5 NS = 3 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N5 K2 = K1 + N5 NWS = NW C C Determine the least integer MQ such that 2 ^ MQ .GE. NW. C T1 = NW MQ = CL2 * LOG (T1) + 1.D0 - RXX C C Compute A^2 outside of the iteration loop. C NW = NWS + 1 CALL MPMUL (A, A, S(K0)) C C Compute the initial approximation of A ^ (-2/3). C CALL MPMDC (A, T1, N) N3 = - 2 * N / 3 T2 = (T1 * 2.D0 ** (N + 3.D0 * N3 / 2.D0)) ** (-2.D0 / 3.D0) CALL MPDMC (T2, N3, B) F(1) = 1. F(2) = 0. F(3) = 1. NW = 3 IQ = 0 C C Perform the Newton-Raphson iteration described above with a dynamically C changing precision level NW (one greater than powers of two). C DO 110 K = 2, MQ NW = MIN (2 * NW - 2, NWS) + 1 100 CONTINUE CALL MPMUL (B, B, S(K1)) CALL MPMUL (B, S(K1), S(K2)) CALL MPMUL (S(K0), S(K2), S(K1)) CALL MPSUB (F, S(K1), S(K2)) CALL MPMUL (B, S(K2), S(K1)) CALL MPDIVD (S(K1), 3.D0, 0, S(K2)) CALL MPADD (B, S(K2), S(K1)) CALL MPEQ (S(K1), B) IF (K .EQ. MQ - NIT .AND. IQ .EQ. 0) THEN IQ = 1 GOTO 100 ENDIF 110 CONTINUE C C Multiply by A to give final result. C CALL MPMUL (A, B, S(K1)) CALL MPEQ (S(K1), B) C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (B) C 120 IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 3) (B(I), I = 1, NO) 3 FORMAT ('MPCBRT O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPCBRX (A, B) C C This computes the cube root of the MP number A and returns the MP result C in B. Before calling MPCBRX, the array in MPCOM5 must be initialized by C calling MPINIX. For modest levels of precision, use MPCBRT. NW should be C a power of two. The last three words of the result are not reliable. C Debug output starts with IDB = 6. C C Max SP space for B: NW + 4 cells. Max SP scratch space: 4.5 * NW + 27 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine uses basically the same Newton iteration algorithm as MPCBRT. C In fact, this routine calls MPCBRT to obtain an initial approximation. C See the comment about the parameter NIT in MPDIVX. C DOUBLE PRECISION CL2, T1 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, NDB = 22, NIT = 1) DIMENSION A(NW+2), B(NW+4), F(8) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. RETURN ENDIF IF (IDB .GE. 6) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPCBRX I'/(6F12.0)) ENDIF C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) NCR = 2 ** MCR C IF (NA .EQ. 0) THEN B(1) = 0. B(2) = 0. GOTO 120 ENDIF IF (IA .LT. 0.D0) THEN IF (KER(14) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPCBRX: Argument is negative.') IER = 14 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if precision level is too low to justify the advanced routine. C IF (NW .LE. NCR) THEN CALL MPCBRT (A, B) GOTO 120 ENDIF N4 = NW + 4 NS = 3 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N4 K2 = K1 + N4 NWS = NW C C Determine the least integer MQ such that 2 ^ MQ .GE. NW. C T1 = NW MQ = CL2 * LOG (T1) + 1.D0 - RXX C C Compute A^2 outside of the iteration loop. C CALL MPMULX (A, A, S(K0)) C C Compute the initial approximation of A ^ (-2/3). C NW = NCR CALL MPCBRT (A, S(K1)) CALL MPDIV (S(K1), A, B) F(1) = 1. F(2) = 0. F(3) = 1. IQ = 0 C C Perform the Newton-Raphson iteration described above with a dynamically C changing precision level NW (powers of two). C DO 110 K = MCR + 1, MQ AN = NW NW = MIN (2 * NW, NWS) 100 CONTINUE CALL MPMULX (B, B, S(K1)) CALL MPMULX (B, S(K1), S(K2)) CALL MPMULX (S(K0), S(K2), S(K1)) CALL MPSUB (F, S(K1), S(K2)) S(K2) = MIN (S(K2), AN) CALL MPMULX (B, S(K2), S(K1)) CALL MPDIVD (S(K1), 3.D0, 0, S(K2)) CALL MPADD (B, S(K2), S(K1)) CALL MPEQ (S(K1), B) IF (K .EQ. MQ - NIT .AND. IQ .EQ. 0) THEN IQ = 1 GOTO 100 ENDIF 110 CONTINUE C C Multiply by A to give final result. C CALL MPMULX (A, B, S(K1)) CALL MPEQ (S(K1), B) ICS = ISS C 120 IF (IDB .GE. 6) THEN NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 3) (B(I), I = 1, NO) 3 FORMAT ('MPCBRX O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPCDIV (L, A, B, C) C C This routine divides the MP complex numbers A and B to yield the MPC C quotient C. L is the offset between real and imaginary parts in A, B C and the result C. L must be at least NW + 4. For extra high levels of C precision, use MPCDVX. The last word is not reliable. Debug output C starts with IDB = 7 C C Max SP space for C: 2 * L cells. Max SP scratch space: 5 * NW + 20 C cells. Max DP scratch space: NW + 4 cells. C C This routine employs the formula described in MPCMUL to save multiprecision C multiplications. C PARAMETER (NDB = 22) DIMENSION A(2*L), B(2*L), C(2*L), F(8) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. C(L+1) = 0. C(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 7) THEN WRITE (LDB, 1) L 1 FORMAT ('MPCDIV I',I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCDIV I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 2) (B(I), I = 1, NO) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 2) (B(L+I), I = 1, NO) ENDIF C IF (L .LT. NW + 4) THEN IF (KER(15) .NE. 0) THEN WRITE (LDB, 3) L, NW + 4 3 FORMAT ('*** MPCDIV: Offset parameter is too small',2I8) IER = 15 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C IF (B(1) .EQ. 0. .AND. B(L1) .EQ. 0.) THEN IF (KER(16) .NE. 0) THEN WRITE (LDB, 4) 4 FORMAT ('*** MPCDIV: Divisor is zero.') IER = 16 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N4 = NW + 4 NS = 5 * N4 ISS = ICS ICS = ICS + NS IF (ICS - 1 .GT. IMS) CALL MPALER IHS = MAX (ICS, IHS) K0 = ISS K1 = K0 + N4 K2 = K1 + N4 K3 = K2 + N4 K4 = K3 + N4 F(1) = 1. F(2) = 0. F(3) = 1. C CALL MPMUL (A, B, S(K0)) CALL MPMUL (A(L1), B(L1), S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPSUB (S(K0), S(K1), S(K3)) CALL MPADD (A, A(L1), S(K0)) CALL MPSUB (B, B(L1), S(K1)) CALL MPMUL (S(K0), S(K1), S(K4)) CALL MPSUB (S(K4), S(K3), S(K1)) CALL MPMUL (B, B, S(K0)) CALL MPMUL (B(L1), B(L1), S(K3)) CALL MPADD (S(K0), S(K3), S(K4)) CALL MPDIV (F, S(K4), S(K0)) CALL MPMUL (S(K2), S(K0), C) CALL MPMUL (S(K1), S(K0), C(L1)) ICS = ISS C IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 5) (C(I), I = 1, NO) 5 FORMAT ('MPCDIV O'/(6F12.0)) NO = MIN (INT (ABS (C(L1))), NDB) + 2 WRITE (LDB, 5) (C(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCDVX (L, A, B, C) C C This routine divides the MP complex numbers A and B to yield the MPC C quotient C. L is the offset between real and imaginary parts in A, B C the result C. L must be at least NW + 4. Before calling MPCDVX, the C array in MPCOM5 must be initialized by calling MPINIX. For modest levels C of precision, use MPCDIV. NW should be a power of two. The last two C words are not reliable. Debug output starts with IDB = 7 C C Max SP space for C: 2 * L cells. Max SP scratch space: 7 * NW + 28 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine employs the same scheme as MPCDIV. C PARAMETER (NDB = 22) DIMENSION A(2*L), B(2*L), C(2*L), F(8) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. C(L+1) = 0. C(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 7) THEN WRITE (LDB, 1) L 1 FORMAT ('MPCDVX I',I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCDVX I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 2) (B(I), I = 1, NO) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 2) (B(L+I), I = 1, NO) ENDIF C IF (L .LT. NW + 4) THEN IF (KER(17) .NE. 0) THEN WRITE (LDB, 3) L, NW + 4 3 FORMAT ('*** MPCDVX: Offset parameter is too small',2I8) IER = 17 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C IF (B(1) .EQ. 0. .AND. B(L1) .EQ. 0.) THEN IF (KER(18) .NE. 0) THEN WRITE (LDB, 4) 4 FORMAT ('*** MPCDVX: Divisor is zero.') IER = 18 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N4 = NW + 4 NS = 5 * N4 ISS = ICS ICS = ICS + NS IF (ICS - 1 .GT. IMS) CALL MPALER IHS = MAX (ICS, IHS) K0 = ISS K1 = K0 + N4 K2 = K1 + N4 K3 = K2 + N4 K4 = K3 + N4 F(1) = 1. F(2) = 0. F(3) = 1. C CALL MPMULX (A, B, S(K0)) CALL MPMULX (A(L1), B(L1), S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPSUB (S(K0), S(K1), S(K3)) CALL MPADD (A, A(L1), S(K0)) CALL MPSUB (B, B(L1), S(K1)) CALL MPMULX (S(K0), S(K1), S(K4)) CALL MPSUB (S(K4), S(K3), S(K1)) CALL MPMULX (B, B, S(K0)) CALL MPMULX (B(L1), B(L1), S(K3)) CALL MPADD (S(K0), S(K3), S(K4)) CALL MPDIVX (F, S(K4), S(K0)) CALL MPMUL (S(K2), S(K0), C) CALL MPMUL (S(K1), S(K0), C(L1)) ICS = ISS C IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 5) (C(I), I = 1, NO) 5 FORMAT ('MPCDVX O'/(6F12.0)) NO = MIN (INT (ABS (C(L1))), NDB) + 2 WRITE (LDB, 5) (C(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCEQ (L, A, B) C C This sets the MPC number B equal to the MPC number A. L is the offset C between real and imaginary parts in A and B. Debug output starts with C IDB = 10. C C Max SP space for B: 2 * L cells. C DIMENSION A(2*L), B(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. RETURN ENDIF IF (IDB .GE. 10) WRITE (LDB, 1) 1 FORMAT ('MPCEQ') C I1 = SIGN (1., A(1)) N1 = MIN (INT (ABS (A(1))), NW, L - 2) I2 = SIGN (1., A(L+1)) N2 = MIN (INT (ABS (A(L+1))), NW, L - 2) B(1) = SIGN (N1, I1) B(L+1) = SIGN (N2, I2) C DO 100 I = 2, N1 + 2 B(I) = A(I) 100 CONTINUE C DO 110 I = 2, N2 + 2 B(L+I) = A(L+I) 110 CONTINUE C RETURN END C SUBROUTINE MPCFFT (IS, M, X, Y) C C This routine computes the 2^M -point complex-to-complex FFT of X. See C article by DHB in Intl. J. of Supercomputer Applications, Spring 1988, C p. 82 - 87). X and Y are double precision. X is both the input and the C output array, while Y is a scratch array. Both X and Y must be C dimensioned with 2 * N cells, where N = 2^M. The data in X are assumed C to have real and imaginary parts separated by N cells. A call to MPCFFT C with IS = 1 (or -1) indicates a call to perform a FFT with positive (or C negative) exponentials. M must be at least two. Before calling MPCRFT, C the array in MPCOM5 must be initialized by calling MPINIX. C C In this application, MPCFFT is called by MPRCFT and MPCRFT, which are in C turn called by MPMULX. This routine is not intended to be called directly C by the user. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION X(*), Y(*) C N = 2 ** M C> C For Cray computers, it is most efficient to limit M1 to 6. For most C scalar computers, it is best to limit M1 to 2. Uncomment whichever of the C next two lines is appropriate. C C M1 = MIN (M / 2, 6) M1 = MIN (M / 2, 2) M2 = M - M1 N2 = 2 ** M1 N1 = 2 ** M2 C C Perform one variant of the Stockham FFT. C DO 100 L = 1, M1, 2 CALL MPFFT1 (IS, L, M, X, Y) IF (L .EQ. M1) GOTO 120 CALL MPFFT1 (IS, L + 1, M, Y, X) 100 CONTINUE C C Perform a transposition of X treated as a N2 x N1 x 2 matrix. C CALL MPTRAN (N1, N2, X, Y) C C Perform second variant of the Stockham FFT from Y to X and X to Y. C DO 110 L = M1 + 1, M, 2 CALL MPFFT2 (IS, L, M, Y, X) IF (L .EQ. M) GOTO 160 CALL MPFFT2 (IS, L + 1, M, X, Y) 110 CONTINUE C GOTO 140 C C Perform a transposition of Y treated as a N2 x N1 x 2 matrix. C 120 CALL MPTRAN (N1, N2, Y, X) C C Perform second variant of the Stockham FFT from X to Y and Y to X. C DO 130 L = M1 + 1, M, 2 CALL MPFFT2 (IS, L, M, X, Y) IF (L .EQ. M) GOTO 140 CALL MPFFT2 (IS, L + 1, M, Y, X) 130 CONTINUE C GOTO 160 C C Copy Y to X. C 140 DO 150 I = 1, 2 * N X(I) = Y(I) 150 CONTINUE C 160 RETURN END C SUBROUTINE MPCMLX (L, A, B, C) C C This routine multiplies the MP complex numbers A and B to yield the MPC C product C. L is the offset between real and imaginary parts in A, B and C the result C. L must be at least NW + 4. Before calling MPCMLX, the C array in MPCOM5 must be initialized by calling MPINIX. For modest levels C of precision, use MPCMUL. NW should be a power of two. The last word is C not reliable. Debug output starts with IDB = 7. C C Max SP space for C: 2 * L cells. Max SP scratch space: 4 * NW + 16 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine employs the same scheme as MPCMUL. C PARAMETER (NDB = 22) DIMENSION A(2*L), B(2*L), C(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. C(L+1) = 0. C(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 7) THEN WRITE (LDB, 1) L 1 FORMAT ('MPCMLX I',I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCMLX I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 2) (B(I), I = 1, NO) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 2) (B(L+I), I = 1, NO) ENDIF C IF (L .LT. NW + 4) THEN IF (KER(19) .NE. 0) THEN WRITE (LDB, 3) L, NW + 4 3 FORMAT ('*** MPCMLX: Offset parameter is too small',2I8) IER = 19 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N4 = NW + 4 NS = 4 * N4 ISS = ICS ICS = ICS + NS IF (ICS - 1 .GT. IMS) CALL MPALER IHS = MAX (ICS, IHS) K0 = ISS K1 = K0 + N4 K2 = K1 + N4 K3 = K2 + N4 C CALL MPMULX (A, B, S(K0)) CALL MPMULX (A(L1), B(L1), S(K1)) CALL MPSUB (S(K0), S(K1), C) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPADD (A, A(L1), S(K0)) CALL MPADD (B, B(L1), S(K1)) CALL MPMULX (S(K0), S(K1), S(K3)) CALL MPSUB (S(K3), S(K2), C(L1)) ICS = ISS C IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 4) (C(I), I = 1, NO) 4 FORMAT ('MPCMLX O'/(6F12.0)) NO = MIN (INT (ABS (C(L1))), NDB) + 2 WRITE (LDB, 4) (C(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCMUL (L, A, B, C) C C This routine multiplies the MP complex numbers A and B to yield the MPC C product C. L is the offset between real and imaginary parts in A, B and C the result C. L must be at least NW + 4. For extra high levels of C precision, use MPCMLX. The last word is not reliable. Debug output C starts with IDB = 7. C C Max SP space for C: 2 * L cells. Max SP scratch space: 4 * NW + 16 C cells. Max DP scratch space: NW + 4 cells. C C This routine employs the formula C C (a_1 + a_2 i) (b_1 + b_2 i) = [a_1 b_1 - a_2 b_2] + C [(a_1 + b_1) (a_2 + b_2) - (a_1 b_1 + a_2 b_2)] i C C Note that this formula can be implemented with only three multiplications C whereas the conventional formula requires four. C PARAMETER (NDB = 22) DIMENSION A(2*L), B(2*L), C(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. C(L+1) = 0. C(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE.7) THEN WRITE (LDB, 1) L 1 FORMAT ('MPCMUL I',I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCMUL I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 2) (B(I), I = 1, NO) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 2) (B(L+I), I = 1, NO) ENDIF C IF (L .LT. NW + 4) THEN IF (KER(20) .NE. 0) THEN WRITE (LDB, 3) L, NW + 4 3 FORMAT ('*** MPCMUL: Offset parameter is too small',2I8) IER = 20 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N4 = NW + 4 NS = 4 * N4 ISS = ICS ICS = ICS + NS IF (ICS - 1 .GT. IMS) CALL MPALER IHS = MAX (ICS, IHS) K0 = ISS K1 = K0 + N4 K2 = K1 + N4 K3 = K2 + N4 C CALL MPMUL (A, B, S(K0)) CALL MPMUL (A(L1), B(L1), S(K1)) CALL MPSUB (S(K0), S(K1), C) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPADD (A, A(L1), S(K0)) CALL MPADD (B, B(L1), S(K1)) CALL MPMUL (S(K0), S(K1), S(K3)) CALL MPSUB (S(K3), S(K2), C(L1)) ICS = ISS C IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 4) (C(I), I = 1, NO) 4 FORMAT ('MPCMUL O'/(6F12.0)) NO = MIN (INT (ABS (C(L1))), NDB) + 2 WRITE (LDB, 4) (C(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCPLX (N, LA, A, X1, NX, LX, X) C C This routine finds a complex root of the N-th degree polynomial whose C MPC coefficients are in A by Newton-Raphson iterations, beginning C at the complex DPE value (X1(1), NX(1)) + i (X1(2), NX(2)), and returns C the MPC root in X. The N + 1 coefficients a_0, a_1, ..., a_N are C assumed to start in locations A(1), A(2*LA+1), A(4*LA+1), etc. LA is the C offset between the real and the imaginary parts of each input coefficient. C Typically LA = NW + 4. LX, also an input parameter, is the offset between C the real and the imaginary parts of the result to be stored in X. LX C should be at least NW + 4. Before calling MPCPLX, the array in MPCOM5 C be initialized by calling MPINIX. For modest levels of precision, use C MPCPOL. NW should be a power of two. The last two words of the result C are not reliable. Debug output starts with IDB = 5. C C Max SP space for X: 2 * LX cells. Max SP scratch space: 17.5 * NW + 115 C cells. Max DP scratch space: 12 * NW + 6 cells. C C See the note in MPPOL about repeated roots. C C This routine employs the same scheme as MPCPOL. C CHARACTER*8 CX DOUBLE PRECISION T1, X1 DIMENSION A(2*LA,N+1), NX(2), X(2*LX), X1(2) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN X(1) = 0. X(2) = 0. X(LX+1) = 0. X(LX+2) = 0. ENDIF IF (IDB .GE. 5) THEN WRITE (LDB, 1) N, LX 1 FORMAT ('MPCPLX I',2I6) C DO 100 K = 0, N WRITE (CX, '(I4)') K CALL MPDEB (CX, A(1,K+1)) CALL MPDEB (CX, A(LA+1,K+1)) 100 CONTINUE C WRITE (LDB, 2) X1(2), NX(2) 2 FORMAT ('MPCPLX I',F16.12,' x 10 ^',I6,F20.12,' x 10^',I6) ENDIF C C Check if precision level is too low to justify the advanced routine. C NCR = 2 ** MCR IF (NW .LE. NCR) THEN CALL MPCPOL (N, LA, A, X1, NX, LX, X) L1 = 0 GOTO 150 ENDIF C C Check if the polynomial is proper. C IF (A(1,1) .EQ. 0. .OR. A(1,N+1) .EQ. 0.) THEN IF (KER(21) .NE. 0) THEN WRITE (LDB, 3) 3 FORMAT ('*** MPCPLX: Either the first or last input ', $ 'coefficient is zero.') IER = 21 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N4 = NW + 4 N8 = 2 * N4 NS = 10 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N8 K2 = K1 + N8 K3 = K2 + N8 K4 = K3 + N8 NWS = NW C C Set the initial value. C NW = NCR CALL MPCPOL (N, LA, A, X1, NX, N4, S(K0)) TL = -4. L1 = 0 LS = -10 C C Perform MP Newton-Raphson iterations to solve P(x) = 0. C 110 L1 = L1 + 1 IF (L1 .EQ. 50) THEN IF (KER(22) .NE. 0) THEN WRITE (LDB, 4) 4 FORMAT ('*** MPCPLX: Iteration limit exceeded.') IER = 22 IF (KER(IER) .EQ. 2) CALL MPABRT ICS = ISS NW = NWS RETURN ENDIF ENDIF C C Compute P(x). C CALL MPMMPC (A(1,N+1), A(LA+1,N+1), N4, S(K1)) C DO 120 K = N - 1, 0, -1 CALL MPCMLX (N4, S(K0), S(K1), S(K2)) CALL MPADD (S(K2), A(1,K+1), S(K1)) CALL MPADD (S(K2+N4), A(LA+1,K+1), S(K1+N4)) 120 CONTINUE C C Compute P'(x). C T1 = N CALL MPMULD (A(1,N+1), T1, 0, S(K2)) CALL MPMULD (A(LA+1,N+1), T1, 0, S(K2+N4)) C DO 130 K = N - 1, 1, -1 CALL MPCMLX (N4, S(K0), S(K2), S(K3)) T1 = K CALL MPMULD (A(1,K+1), T1, 0, S(K4)) CALL MPMULD (A(LA+1,K+1), T1, 0, S(K4+N4)) CALL MPCADD (N4, S(K3), S(K4), S(K2)) 130 CONTINUE C C Compute P(x) / P'(x) and update x. C CALL MPCDVX (N4, S(K1), S(K2), S(K3)) CALL MPCSUB (N4, S(K0), S(K3), S(K4)) C IF (IDB .GE. 6) THEN WRITE (LDB, 5) L1 5 FORMAT ('ITERATION',I4) CALL MPDEB ('X', S(K0)) CALL MPDEB (' ', S(K0+N4)) CALL MPDEB ('P(X)', S(K1)) CALL MPDEB (' ', S(K1+N4)) CALL MPDEB ('P''(X)', S(K2)) CALL MPDEB (' ', S(K2+N4)) CALL MPDEB ('CORR', S(K3)) CALL MPDEB (' ', S(K3+N4)) ENDIF CALL MPCEQ (N4, S(K4), S(K0)) C C If this was the second iteration at full precision, there is no need to C continue (the adjusted value of x is correct); otherwise repeat. C IF (L1 .EQ. LS + 1) GOTO 140 IF (S(K3) .NE. 0. .AND. S(K3+1) .GT. TL .OR. S(K3+N4) .NE. 0. $ .AND. S(K3+N4+1) .GT. TL) GOTO 110 C C Newton iterations have converged to current precision. Increase precision C and continue. C IF (NW .EQ. NWS) GOTO 140 NW = MIN (2 * NW, NWS) IF (NW .EQ. NWS) LS = L1 IF (NW .LE. 32) THEN TL = 2 - NW ELSEIF (NW .LE. 256) THEN TL = 3 - NW ELSE TL = 4 - NW ENDIF IF (IDB .GE. 6) THEN WRITE (LDB, 6) NW 6 FORMAT (6X,'New NW =', I8) ENDIF GOTO 110 C 140 CALL MPMMPC (S(K0), S(K0+N4), LX, X) ICS = ISS C 150 IF (IDB .GE. 5) THEN WRITE (LDB, 7) L1 7 FORMAT ('Iteration count:',I5) CALL MPDEB ('MPCPLX O', X) CALL MPDEB (' ', X(LX+1)) ENDIF RETURN END C SUBROUTINE MPCPOL (N, LA, A, X1, NX, LX, X) C C This routine finds a complex root of the N-th degree polynomial whose C MPC coefficients are in A by Newton-Raphson iterations, beginning C at the complex DPE value (X1(1), NX(1)) + i (X1(2), NX(2)), and returns C the MPC root in X. The N + 1 coefficients a_0, a_1, ..., a_N are C assumed to start in locations A(1), A(2*LA+1), A(4*LA+1), etc. LA is the C offset between the real and the imaginary parts of each input coefficient. C Typically LA = NW + 4. LX, also an input parameter, is the offset between C the real and the imaginary parts of the result to be stored in X. LX must C be at least NW + 4. For extra high levels of precision, use MPCPLX. C Debug output starts with IDB = 5. C C Max SP space for X: 2 * LX cells. Max SP scratch space: 15 * NW + 75 C cells. Max DP scratch space: NW + 5 cells. C C See the note about repeated roots in MPPOL. C C This routine employs the complex form of the Newton-Raphson iteration: C C X_{k+1} = X_k - P(X_k) / P'(X_k) C C These iterations are performed with a maximum precision level NW that is C dynamically changed, approximately doubling with each iteration. C CHARACTER*8 CX DOUBLE PRECISION T1, X1 DIMENSION A(2*LA,N+1), NX(2), X(2*LX), X1(2) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN X(1) = 0. X(2) = 0. X(LX+1) = 0. X(LX+2) = 0. ENDIF IF (IDB .GE. 5) THEN WRITE (LDB, 1) N, LX 1 FORMAT ('MPCPOL I',2I6) C DO 100 K = 0, N WRITE (CX, '(I4)') K CALL MPDEB (CX, A(1,K+1)) CALL MPDEB (CX, A(LA+1,K+1)) 100 CONTINUE C WRITE (LDB, 2) X1(1), NX(1), X1(2), NX(2) 2 FORMAT ('MPCPOL I',F16.12,' x 10 ^',I6,F20.12,' x 10^',I6) ENDIF C C Check if the polynomial is proper. C IF (A(1,1) .EQ. 0. .OR. A(1,N+1) .EQ. 0.) THEN IF (KER(23) .NE. 0) THEN WRITE (LDB, 3) 3 FORMAT ('*** MPCPOL: Either the first or last input ', $ 'coefficient is zero.') IER = 23 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N5 = NW + 5 N10 = 2 * N5 NS = 10 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N10 K2 = K1 + N10 K3 = K2 + N10 K4 = K3 + N10 NWS = NW NW = NW + 1 C C Set the initial value. C CALL MPDMC (X1(1), NX(1), S(K0)) CALL MPDMC (X1(2), NX(2), S(K0+N5)) NW = 5 TL = -4. L1 = 0 LS = -10 C C Perform MP Newton-Raphson iterations to solve P(x) = 0. C 110 L1 = L1 + 1 IF (L1 .EQ. 50) THEN IF (KER(24) .NE. 0) THEN WRITE (LDB, 4) 4 FORMAT ('*** MPCPOL: Iteration limit exceeded.') IER = 24 IF (KER(IER) .EQ. 2) CALL MPABRT ICS = ISS NW = NWS RETURN ENDIF ENDIF C C Compute P(x). C CALL MPMMPC (A(1,N+1), A(LA+1,N+1), N5, S(K1)) C DO 120 K = N - 1, 0, -1 CALL MPCMUL (N5, S(K0), S(K1), S(K2)) CALL MPADD (S(K2), A(1,K+1), S(K1)) CALL MPADD (S(K2+N5), A(LA+1,K+1), S(K1+N5)) 120 CONTINUE C C Compute P'(x). C T1 = N CALL MPMULD (A(1,N+1), T1, 0, S(K2)) CALL MPMULD (A(LA+1,N+1), T1, 0, S(K2+N5)) C DO 130 K = N - 1, 1, -1 CALL MPCMUL (N5, S(K0), S(K2), S(K3)) T1 = K CALL MPMULD (A(1,K+1), T1, 0, S(K4)) CALL MPMULD (A(LA+1,K+1), T1, 0, S(K4+N5)) CALL MPCADD (N5, S(K3), S(K4), S(K2)) 130 CONTINUE C C Compute P(x) / P'(x) and update x. C CALL MPCDIV (N5, S(K1), S(K2), S(K3)) CALL MPCSUB (N5, S(K0), S(K3), S(K4)) C IF (IDB .GE. 6) THEN WRITE (LDB, 5) L1 5 FORMAT ('Iteration',I4) CALL MPDEB ('X', S(K0)) CALL MPDEB (' ', S(K0+N5)) CALL MPDEB ('P(X)', S(K1)) CALL MPDEB (' ', S(K1+N5)) CALL MPDEB ('P''(X)', S(K2)) CALL MPDEB (' ', S(K2+N5)) CALL MPDEB ('CORR', S(K3)) CALL MPDEB (' ', S(K3+N5)) ENDIF CALL MPCEQ (N5, S(K4), S(K0)) C C If this was the second iteration at full precision, there is no need to C continue (the adjusted value of x is correct); otherwise repeat. C IF (L1 .EQ. LS + 1) GOTO 140 IF (S(K3) .NE. 0. .AND. S(K3+1) .GT. TL .OR. S(K3+N5) .NE. 0. $ .AND. S(K3+N5+1) .GT. TL) GOTO 110 C C Newton iterations have converged to current precision. Increase precision C and continue. C IF (NW .EQ. NWS + 1) GOTO 140 NW = MIN (2 * NW - 2, NWS) + 1 IF (NW .EQ. NWS + 1) LS = L1 TL = 1 - NW IF (IDB .GE. 6) THEN WRITE (LDB, 6) NW 6 FORMAT (6X,'New NW =', I8) ENDIF GOTO 110 C 140 CALL MPMMPC (S(K0), S(K0+N5), LX, X) C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (X) CALL MPROUN (X(LX+1)) C IF (IDB .GE. 5) THEN WRITE (LDB, 7) L1 7 FORMAT ('Iteration count:',I5) CALL MPDEB ('MPCPOL O', X) CALL MPDEB (' ', X(LX+1)) ENDIF RETURN END C SUBROUTINE MPCPR (A, B, IC) C C This routine compares the MP numbers A and B and returns in IC the value C -1, 0, or 1 depending on whether A < B, A = B, or A > B. It is faster C than merely subtracting A and B and looking at the sign of the result. C Debug output begins with IDB = 9. C DIMENSION A(NW+4), B(NW+4) PARAMETER (NDB = 22) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN IC = 0 RETURN ENDIF IF (IDB .GE. 9) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPCPR I'/(6F12.0)) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 1) (B(I), I = 1, NO) ENDIF IA = SIGN (1., A(1)) IF (A(1) .EQ. 0.) IA = 0 IB = SIGN (1., B(1)) IF (B(1) .EQ. 0.) IB = 0 C C Compare signs. C IF (IA .NE. IB) THEN IC = SIGN (1, IA - IB) GOTO 110 ENDIF C C The signs are the same. Compare exponents. C MA = A(2) MB = B(2) IF (MA .NE. MB) THEN IC = IA * SIGN (1, MA - MB) GOTO 110 ENDIF C C The signs and the exponents are the same. Compare mantissas. C NA = MIN (INT (ABS (A(1))), NW) NB = MIN (INT (ABS (B(1))), NW) C DO 100 I = 3, MIN (NA, NB) + 2 IF (A(I) .NE. B(I)) THEN IC = IA * SIGN (1., A(I) - B(I)) GOTO 110 ENDIF 100 CONTINUE C C The mantissas are the same to the common length. Compare lengths. C IF (NA .NE. NB) THEN IC = IA * SIGN (1, NA - NB) GOTO 110 ENDIF C C The signs, exponents, mantissas and lengths are the same. Thus A = B. C IC = 0 C 110 IF (IDB .GE. 9) WRITE (6, 2) IC 2 FORMAT ('MPCPR O',I4) RETURN END C SUBROUTINE MPCPWR (L, A, N, B) C C This computes the N-th power of the MPC number A and returns the MPC C result C in B. When N is zero, 1 is returned. When N is negative, the C reciprocal of A ^ |N| is returned. L is the offset between real and C imaginary parts in A and B. L should be at least NW + 4. For extra high C levels of precision, use MPCPWX. Debug output starts with IDB = 7. C C Max SP space for B: 2 * L cells. Max SP scratch space: 6 * NW + 30 C cells. Max DP scratch space: NW + 5 cells. C C This routine employs the binary method for exponentiation. C DOUBLE PRECISION CL2, T1 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, NDB = 22) DIMENSION A(2*L), B(2*L), F1(8), F2(8) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 7) THEN WRITE (6, 1) L, N 1 FORMAT ('MPCPWR I',2I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCPWR I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) ENDIF C NA1 = MIN (INT (ABS (A(1))), NW) NA2 = MIN (INT (ABS (A(L1))), NW) IF (NA1 .EQ. 0 .AND. NA2 .EQ. 0) THEN IF (N .GE. 0) THEN B(1) = 0. B(2) = 0. B(L1) = 0. B(L1+1) = 0. GOTO 120 ELSE IF (KER(25) .NE. 0) THEN WRITE (LDB, 3) 3 FORMAT ('*** MPCPWR: Argument is zero and N is negative or', $ ' zero.') IER = 25 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF ENDIF C N5 = NW + 5 NS = 6 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + 2 * N5 K2 = K1 + 2 * N5 NWS = NW NW = NW + 1 NN = ABS (N) F1(1) = 1. F1(2) = 0. F1(3) = 1. F2(1) = 0. F2(2) = 0. CALL MPMMPC (A, A(L1), N5, S(K0)) IF (NN .EQ. 0) THEN CALL MPMMPC (F1, F2, L, B) NW = NWS ICS = ISS GOTO 120 ELSEIF (NN .EQ. 1) THEN CALL MPCEQ (N5, S(K0), S(K2)) GOTO 110 ELSEIF (NN .EQ. 2) THEN CALL MPCMUL (N5, S(K0), S(K0), S(K2)) GOTO 110 ENDIF C C Determine the least integer MN such that 2 ^ MN .GT. NN. C T1 = NN MN = CL2 * LOG (T1) + 1.D0 + RXX CALL MPMMPC (F1, F2, N5, S(K2)) KN = NN C C Compute B ^ N using the binary rule for exponentiation. C DO 100 J = 1, MN KK = KN / 2 IF (KN .NE. 2 * KK) THEN CALL MPCMUL (N5, S(K2), S(K0), S(K1)) CALL MPCEQ (N5, S(K1), S(K2)) ENDIF KN = KK IF (J .LT. MN) THEN CALL MPCMUL (N5, S(K0), S(K0), S(K1)) CALL MPCEQ (N5, S(K1), S(K0)) ENDIF 100 CONTINUE C C Compute reciprocal if N is negative. C 110 IF (N .LT. 0) THEN CALL MPMMPC (F1, F2, N5, S(K1)) CALL MPCDIV (N5, S(K1), S(K2), S(K0)) CALL MPCEQ (N5, S(K0), S(K2)) ENDIF CALL MPMMPC (S(K2), S(N5+K2), L, B) C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (B) CALL MPROUN (B(L1)) C 120 IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 4) (B(I), I = 1, NO) 4 FORMAT ('MPCPWR O'/(6F12.0)) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 4) (B(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCPWX (L, A, N, B) C C This computes the N-th power of the MPC number A and returns the MPC C result C in B. When N is zero, 1 is returned. When N is negative, the C reciprocal of A ^ |N| is returned. L is the offset between real and C imaginary parts in A and B. L should be at least NW + 4. Before calling C MPCPWX, the array in MPCOM5 must be initialized by calling MPINIX. For C modest levels of precision, use MPCPWR. NW should be a power of two. C The last two words of the result are not reliable. Debug output starts C with IDB = 6. C C Max SP space for B: 2 * L cells. Max SP scratch space: 8 * NW + 32 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine employs the binary method for exponentiation. C DOUBLE PRECISION CL2, T1 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, NDB = 22) DIMENSION A(2*L), B(2*L), F1(8), F2(8) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 6) THEN WRITE (6, 1) L, N 1 FORMAT ('MPCPWX I',2I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCPWX I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) ENDIF C NA1 = MIN (INT (ABS (A(1))), NW) NA2 = MIN (INT (ABS (A(L1))), NW) NCR = 2 ** MCR C C Check if precision level of A is too low to justify advanced routine. C IF (NA1 .LE. NCR .AND. NA2 .LE. NCR) THEN CALL MPCPWR (L, A, N, B) GOTO 120 ENDIF IF (NA1 .EQ. 0 .AND. NA2 .EQ. 0) THEN IF (N .GE. 0) THEN B(1) = 0. B(2) = 0. B(L1) = 0. B(L1+1) = 0. GOTO 120 ELSE IF (KER(26) .NE. 0) THEN WRITE (LDB, 3) 3 FORMAT ('*** MPCPWX: Argument is zero and N is negative or', $ ' zero.') IER = 26 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF ENDIF C N4 = NW + 4 NS = 6 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + 2 * N4 K2 = K1 + 2 * N4 NN = ABS (N) F1(1) = 1. F1(2) = 0. F1(3) = 1. F2(1) = 0. F2(2) = 0. CALL MPMMPC (A, A(L1), N4, S(K0)) IF (NN .EQ. 0) THEN CALL MPMMPC (F1, F2, L, B) ICS = ISS GOTO 120 ELSEIF (NN .EQ. 1) THEN CALL MPCEQ (N4, S(K0), S(K2)) GOTO 110 ELSEIF (NN .EQ. 2) THEN CALL MPCMLX (N4, S(K0), S(K0), S(K2)) GOTO 110 ENDIF C C Determine the least integer MN such that 2 ^ MN .GT. NN. C T1 = NN MN = CL2 * LOG (T1) + 1.D0 + RXX CALL MPMMPC (F1, F2, N4, S(K2)) KN = NN C C Compute B ^ N using the binary rule for exponentiation. C DO 100 J = 1, MN KK = KN / 2 IF (KN .NE. 2 * KK) THEN CALL MPCMLX (N4, S(K2), S(K0), S(K1)) CALL MPCEQ (N4, S(K1), S(K2)) ENDIF KN = KK IF (J .LT. MN) THEN CALL MPCMLX (N4, S(K0), S(K0), S(K1)) CALL MPCEQ (N4, S(K1), S(K0)) ENDIF 100 CONTINUE C C Compute reciprocal if N is negative. C 110 IF (N .LT. 0) THEN CALL MPMMPC (F1, F2, N4, S(K1)) CALL MPCDVX (N4, S(K1), S(K2), S(K0)) CALL MPCEQ (N4, S(K0), S(K2)) ENDIF CALL MPMMPC (S(K2), S(N4+K2), L, B) ICS = ISS C 120 IF (IDB .GE. 6) THEN NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 4) (B(I), I = 1, NO) 4 FORMAT ('MPCPWX O'/(6F12.0)) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 4) (B(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCRFT (IS, M, X, Y) C C This performs an N-point complex-to-real FFT, where N = 2^M. X and Y C are double precision arrays. X is both the input and the output data C array, and Y is a scratch array. N/2 + 1 complex pairs are input, with C real and imaginary parts separated by N/2 + 1 locations, and N real C values are output . A call to MPCRFT with IS = 1 (or -1) indicates a call C to perform a complex-to-real FFT with positive (or negative) exponentials. C M must be at least three. The arrays X and Y must be dimensioned with C N + 2 cells. Before calling MPCRFT, the U array in MPCOM5 must be C initialized by calling MPINIX. C C In this application, MPCRFT is called by MPMULX. This routine is not C intended to be called directly by the user. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION X(*), Y(*) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM5/ U(1024) C C Set initial parameters. C K = U(1) MX = MOD (K, 64) NU = K / 64 N = 2 ** M N2 = N / 2 N21 = N2 + 1 N4 = N / 4 KU = N / 2 KN = KU + NU C C Check if input parameters are invalid. C IF ((IS .NE. 1 .AND. IS .NE. -1) .OR. M .LT. 3 .OR. M .GT. MX) $ THEN IF (KER(27) .NE. 0) THEN WRITE (LDB, 1) IS, M, MX 1 FORMAT ('*** MPCRFT: Either U has not been initialized'/ $ 'or else one of the input parameters is invalid', 3I5) IER = 27 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Construct the input to MPCFFT. C Y(1) = 0.5D0 * (X(1) + X(N21)) Y(N2+1) = 0.5D0 * (X(1) - X(N21)) Y(N4+1) = X(N4+1) Y(N4+N2+1) = -IS * X(N4+N2+2) C CDIR$ IVDEP DO 100 K = 2, N4 X11 = X(K) X12 = X(K+N21) X21 = X(N2+2-K) X22 = X(N+3-K) A1 = X11 + X21 A2 = X11 - X21 B1 = X12 + X22 B2 = X12 - X22 U1 = U(K+KU) U2 = IS * U(K+KN) T1 = U1 * B1 + U2 * A2 T2 = U1 * A2 - U2 * B1 Y(K) = 0.5D0 * (A1 - T1) Y(K+N2) = 0.5D0 * (B2 + T2) Y(N2+2-K) = 0.5D0 * (A1 + T1) Y(N+2-K) = 0.5D0 * (-B2 + T2) 100 CONTINUE C C Perform a normal N/2-point FFT on Y. C CALL MPCFFT (IS, M - 1, Y, X) C C Copy Y to X such that Y(k) = X(2k-1) + i X(2k). C CDIR$ IVDEP DO 110 K = 1, N2 X(2*K-1) = Y(K) X(2*K) = Y(K+N2) 110 CONTINUE C RETURN END C SUBROUTINE MPCSHX (A, PI, AL2, X, Y) C C This computes the hyperbolic cosine and sine of the MP number A and C returns the two MP results in X and Y, respectively. PI is the MP value C of Pi computed by a previous call to MPPI or MPPIX. AL2 is the MP value C of Log (10) computed by a previous call to MPLOG or MPLOGX. Before C calling MPCSHX, the array in MPCOM5 must be initialized by calling MPINIX. C For modest levels of precision, use MPCSSH. NW should be a power of two. C The last four words of the result are not reliable. Debug output starts C with IDB = 5. C C Max SP space for X and Y: NW + 4 cells. Max SP scratch space: C 28 * NW + 132 cells. Max DP scratch space: 12 * NX + 6 cells. C DIMENSION A(NW+2), F(8), AL2(NW+2), PI(NW+2), X(NW+4), Y(NW+4) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN X(1) = 0. X(2) = 0. Y(1) = 0. Y(2) = 0. RETURN ENDIF IF (IDB .GE. 5) CALL MPDEB ('MPCSHX I', A) C N4 = NW + 4 NS = 3 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N4 K2 = K1 + N4 F(1) = 1. F(2) = 0. F(3) = 1. C CALL MPEXPX (A, PI, AL2, S(K0)) CALL MPDIVX (F, S(K0), S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPMULD (S(K2), 0.5D0, 0, X) CALL MPSUB (S(K0), S(K1), S(K2)) CALL MPMULD (S(K2), 0.5D0, 0, Y) ICS = ISS C IF (IDB .GE. 5) THEN CALL MPDEB ('MPCSHX O', X) CALL MPDEB ('MPCSHX O', Y) ENDIF RETURN END C SUBROUTINE MPCSQR (L, A, B) C C This routine computes the complex square root of the MPC number C. L is C the offset between real and imaginary parts in A and B. L must be at C least NW + 4. For extra high levels of precision, use MPCSQX. The last C word is not reliable. Debug output starts with IDB = 6. C C Max SP space for B: 2 * L cells. Max SP scratch space: 5 * NW + 22 C cells. Max DP scratch space: NW + 5 cells. C C This routine uses the following formula, where A1 and A2 are the real and C imaginary parts of A, and where R = Sqrt [A1 ^ 2 + A2 ^2]: C C B = Sqrt [(R + A1) / 2] + I Sqrt [(R - A1) / 2] C C If the imaginary part of A is < 0, then the imaginary part of B is also C set to be < 0. C PARAMETER (NDB = 22) DIMENSION A(2*L), B(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 6) THEN WRITE (LDB, 1) L 1 FORMAT ('MPCSQR I',I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCSQR I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) ENDIF C IF (A(1) .EQ. 0. .AND. A(L+1) .EQ. 0.) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. GOTO 100 ENDIF C N4 = NW + 4 NS = 4 * N4 ISS = ICS ICS = ICS + NS IF (ICS - 1 .GT. IMS) CALL MPALER IHS = MAX (ICS, IHS) K0 = ISS K1 = K0 + N4 K2 = K1 + N4 C CALL MPMUL (A, A, S(K0)) CALL MPMUL (A(L1), A(L1), S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPSQRT (S(K2), S(K0)) CALL MPEQ (A, S(K1)) S(K1) = ABS (S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPMULD (S(K2), 0.5D0, 0, S(K1)) CALL MPSQRT (S(K1), S(K0)) CALL MPMULD (S(K0), 2.D0, 0, S(K1)) IF (A(1) .GE. 0.) THEN CALL MPEQ (S(K0), B) CALL MPDIV (A(L1), S(K1), B(L1)) ELSE CALL MPDIV (A(L1), S(K1), B) B(1) = ABS (B(1)) CALL MPEQ (S(K0), B(L1)) B(L1) = SIGN (B(L1), A(L1)) ENDIF ICS = ISS C 100 IF (IDB .GE. 6) THEN NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 3) (B(I), I = 1, NO) 3 FORMAT ('MPCSQR O'/(6F12.0)) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 3) (B(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCSQX (L, A, B) C C This routine computes the complex square root of the MPC number C. L is C the offset between real and imaginary parts in A and B. L must be at C least NW + 4. For modest levels of precision, use MPCSQR. The last two C words are not reliable. Debug output starts with IDB = 5. C C Max SP space for B: 2 * L cells. Max SP scratch space: 6 * NW + 30. C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine uses the same algorithm as MPCSQR. C PARAMETER (NDB = 22) DIMENSION A(2*L), B(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. RETURN ENDIF L1 = L + 1 IF (IDB .GE. 5) THEN WRITE (LDB, 1) L 1 FORMAT ('MPCSQX I',I10) NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 2) (A(I), I = 1, NO) 2 FORMAT ('MPCSQX I'/(6F12.0)) NO = MIN (INT (ABS (A(L1))), NDB) + 2 WRITE (LDB, 2) (A(L+I), I = 1, NO) ENDIF C IF (A(1) .EQ. 0. .AND. A(L+1) .EQ. 0.) THEN B(1) = 0. B(2) = 0. B(L+1) = 0. B(L+2) = 0. GOTO 100 ENDIF C N4 = NW + 4 NS = 4 * N4 ISS = ICS ICS = ICS + NS IF (ICS - 1 .GT. IMS) CALL MPALER IHS = MAX (ICS, IHS) K0 = ISS K1 = K0 + N4 K2 = K1 + N4 C CALL MPMULX (A, A, S(K0)) CALL MPMULX (A(L1), A(L1), S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPSQRX (S(K2), S(K0)) CALL MPEQ (A, S(K1)) S(K1) = ABS (S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPMULD (S(K2), 0.5D0, 0, S(K1)) CALL MPSQRX (S(K1), S(K0)) CALL MPMULD (S(K0), 2.D0, 0, S(K1)) IF (A(1) .GE. 0.) THEN CALL MPEQ (S(K0), B) CALL MPDIVX (A(L1), S(K1), B(L1)) ELSE CALL MPDIVX (A(L1), S(K1), B) B(1) = ABS (B(1)) CALL MPEQ (S(K0), B(L1)) B(L1) = SIGN (B(L1), A(L1)) ENDIF ICS = ISS C 100 IF (IDB .GE. 5) THEN NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 3) (B(I), I = 1, NO) 3 FORMAT ('MPCSQX O'/(6F12.0)) NO = MIN (INT (ABS (B(L1))), NDB) + 2 WRITE (LDB, 3) (B(L+I), I = 1, NO) ENDIF RETURN END C SUBROUTINE MPCSSH (A, AL2, X, Y) C C This computes the hyperbolic cosine and sine of the MP number A and C returns the two MP results in X and Y, respectively. AL2 is the MP value C of Log (10) computed by a previous call to MPLOG. For extra high levels of C precision, use MPCSHX. The last word of the result is not reliable. C Debug output starts with IDB = 5. C C Max SP space for X and Y: NW + 4 cells. Max SP scratch space: 9 * NW + 50 C cells. Max DP scratch space: NW + 6 cells. C DIMENSION A(NW+2), F(8), AL2(NW+2), X(NW+4), Y(NW+4) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN X(1) = 0. X(2) = 0. Y(1) = 0. Y(2) = 0. RETURN ENDIF IF (IDB .GE. 5) CALL MPDEB ('MPCSSH I', A) C N5 = NW + 5 NS = 4 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N5 K2 = K1 + N5 K3 = K2 + N5 NWS = NW NW = NW + 1 F(1) = 1. F(2) = 0. F(3) = 1. C CALL MPEXP (A, AL2, S(K0)) CALL MPDIV (F, S(K0), S(K1)) CALL MPADD (S(K0), S(K1), S(K2)) CALL MPMULD (S(K2), 0.5D0, 0, S(K3)) CALL MPEQ (S(K3), X) CALL MPSUB (S(K0), S(K1), S(K2)) CALL MPMULD (S(K2), 0.5D0, 0, S(K3)) CALL MPEQ (S(K3), Y) C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (X) CALL MPROUN (Y) C IF (IDB .GE. 5) THEN CALL MPDEB ('MPCSSH O', X) CALL MPDEB ('MPCSSH O', Y) ENDIF RETURN END C SUBROUTINE MPCSSN (A, PI, X, Y) C C This computes the cosine and sine of the MP number A and returns the two MP C results in X and Y, respectively. PI is the MP value of Pi computed by a C previous call to MPPI. For extra high levels of precision, use MPCSSX. C The last word of the result is not reliable. Debug output starts with C IDB = 6. C C Max SP space for X and Y: NW + 4 cells. Max SP scratch space: 9 * NW + 47 C cells. Max DP scratch space: NW + 6 cells. C C This routine uses the conventional Taylor's series for Sin (s): C C Sin (s) = s - s^3 / 3! + s^5 / 5! - s^7 / 7! ... C C where s = t - a * pi / 2 - b * pi / 16 and the integers a and b are chosen C to minimize the absolute value of s. We can then compute C C Sin (t) = Sin (s + a * pi / 2 + b * pi / 16) C Cos (t) = Cos (s + a * pi / 2 + b * pi / 16) C C by applying elementary trig identities for sums. The sine and cosine of C b * pi / 16 are of the form 1/2 * Sqrt {2 +- Sqrt [2 +- Sqrt(2)]}. C Reducing t in this manner insures that -Pi / 32 < s <= Pi / 32, which C accelerates convergence in the above series. C DOUBLE PRECISION CPI, T1, T2 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CPI = 3.141592653589793D0) DIMENSION A(NW+2), F(8), PI(NW+2), X(NW+4), Y(NW+4) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN X(1) = 0. X(2) = 0. Y(1) = 0. Y(2) = 0. RETURN ENDIF IF (IDB .GE. 6) CALL MPDEB ('MPCSSN I', A) C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) IF (NA .EQ. 0) THEN X(1) = 1. X(2) = 0. X(3) = 1. Y(1) = 0. Y(2) = 0. L1 = 0 GOTO 120 ENDIF C C Check if Pi has been precomputed. C CALL MPMDC (PI, T1, N1) IF (N1 .NE. 0 .OR. ABS (T1 - CPI) .GT. RX2) THEN IF (KER(28) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPCSSN: PI must be precomputed.') IER = 28 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N5 = NW + 5 NS = 7 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N5 K2 = K1 + N5 K3 = K2 + N5 K4 = K3 + N5 K5 = K4 + N5 K6 = K5 + N5 NWS = NW NW = NW + 1 F(1) = 1. F(2) = 0. F(3) = 1. C C Reduce to between - Pi and Pi. C CALL MPMULD (PI, 2.D0, 0, S(K0)) CALL MPDIV (A, S(K0), S(K1)) CALL MPNINT (S(K1), S(K2)) CALL MPSUB (S(K1), S(K2), S(K3)) C C Determine nearest multiple of Pi / 2, and within a quadrant, the nearest C multiple of Pi / 16. Through most of the rest of this subroutine, KA and C KB are the integers a and b of the algorithm above. C CALL MPMDC (S(K3), T1, N1) IF (N1 .GE. -24) THEN T1 = T1 * 2.D0 ** N1 T2 = 4.D0 * T1 KA = NINT (T2) KB = NINT (8.D0 * (T2 - KA)) ELSE KA = 0 KB = 0 ENDIF T1 = (8 * KA + KB) / 32.D0 CALL MPDMC (T1, 0, S(K1)) CALL MPSUB (S(K3), S(K1), S(K2)) CALL MPMUL (S(K0), S(K2), S(K1)) C C Compute cosine and sine of the reduced argument s using Taylor's series. C IF (S(K1) .EQ. 0.) THEN S(K0) = 0. S(K0+1) = 0. L1 = 0 GOTO 110 ENDIF CALL MPEQ (S(K1), S(K0)) CALL MPMUL (S(K0), S(K0), S(K2)) L1 = 0 C 100 L1 = L1 + 1 IF (L1 .EQ. 10000) THEN IF (KER(29) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPCSSN: Iteration limit exceeded.') IER = 29 IF (KER(IER) .EQ. 2) CALL MPABRT ICS = ISS NW = NWS RETURN ENDIF ENDIF C T2 = - (2.D0 * L1) * (2.D0 * L1 + 1.D0) CALL MPMUL (S(K2), S(K1), S(K3)) CALL MPDIVD (S(K3), T2, 0, S(K1)) CALL MPADD (S(K1), S(K0), S(K3)) CALL MPEQ (S(K3), S(K0)) C C Check for convergence of the series. C IF (S(K1) .NE. 0. .AND. S(K1+1) .GE. S(K0+1) - NW) GOTO 100 C C Compute Cos (s) = Sqrt [1 - Sin^2 (s)]. C 110 CALL MPEQ (S(K0), S(K1)) CALL MPMUL (S(K0), S(K0), S(K2)) CALL MPSUB (F, S(K2), S(K3)) CALL MPSQRT (S(K3), S(K0)) C C Compute cosine and sine of b * Pi / 16. C KC = ABS (KB) F(3) = 2. IF (KC .EQ. 0) THEN S(K2) = 1. S(K2+1) = 0. S(K2+2) = 1. S(K3) = 0. S(K3+1) = 0. ELSE IF (KC .EQ. 1) THEN CALL MPSQRT (F, S(K4)) CALL MPADD (F, S(K4), S(K5)) CALL MPSQRT (S(K5), S(K4)) ELSEIF (KC .EQ. 2) THEN CALL MPSQRT (F, S(K4)) ELSEIF (KC .EQ. 3) THEN CALL MPSQRT (F, S(K4)) CALL MPSUB (F, S(K4), S(K5)) CALL MPSQRT (S(K5), S(K4)) ELSEIF (KC .EQ. 4) THEN S(K4) = 0. S(K4+1) = 0. ENDIF CALL MPADD (F, S(K4), S(K5)) CALL MPSQRT (S(K5), S(K3)) CALL MPMULD (S(K3), 0.5D0, 0, S(K2)) CALL MPSUB (F, S(K4), S(K5)) CALL MPSQRT (S(K5), S(K4)) CALL MPMULD (S(K4), 0.5D0, 0, S(K3)) ENDIF IF (KB .LT. 0) S(K3) = - S(K3) C C Apply the trigonometric summation identities to compute cosine and sine of C s + b * Pi / 16. C CALL MPMUL (S(K0), S(K2), S(K4)) CALL MPMUL (S(K1), S(K3), S(K5)) CALL MPSUB (S(K4), S(K5), S(K6)) CALL MPMUL (S(K1), S(K2), S(K4)) CALL MPMUL (S(K0), S(K3), S(K5)) CALL MPADD (S(K4), S(K5), S(K1)) CALL MPEQ (S(K6), S(K0)) C C This code in effect applies the trigonometric summation identities for C (s + b * Pi / 16) + a * Pi / 2. C IF (KA .EQ. 0) THEN CALL MPEQ (S(K0), X) CALL MPEQ (S(K1), Y) ELSEIF (KA .EQ. 1) THEN CALL MPEQ (S(K1), X) X(1) = - X(1) CALL MPEQ (S(K0), Y) ELSEIF (KA .EQ. -1) THEN CALL MPEQ (S(K1), X) CALL MPEQ (S(K0), Y) Y(1) = - Y(1) ELSEIF (KA .EQ. 2 .OR. KA .EQ. -2) THEN CALL MPEQ (S(K0), X) X(1) = - X(1) CALL MPEQ (S(K1), Y) Y(1) = - Y(1) ENDIF C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (X) CALL MPROUN (Y) C 120 IF (IDB .GE. 6) THEN WRITE (LDB, 3) L1 3 FORMAT ('Iteration count:',I5) CALL MPDEB ('MPCSSN O', X) CALL MPDEB ('MPCSSN O', Y) ENDIF RETURN END C SUBROUTINE MPCSSX (A, PI, X, Y) C C This computes the cosine and sine of the MP number A and returns the two MP C results in X and Y, respectively. PI is the MP value of Pi computed by a C previous call to MPPI or MPPIX. Before calling MPCSSX, the array in C MPCOM5 must be initialized by calling MPINIX. For modest levels of C precision, use MPCSSN. NW should be a power of two. The last four words C of the result are not reliable. Debug output starts with IDB = 5. C C Max SP space for X and Y: NW + 4 cells. Max SP scratch space: 26*NW + 110 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine employs a complex arithmetic version of the scheme used in C MPEXPX. C DOUBLE PRECISION CL2, CPI, T1, T2 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, CPI = 3.141592653589793D0, $ NIT = 1) DIMENSION A(NW+2), F1(8), PI(NW+2), X(NW+4), Y(NW+4) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN X(1) = 0. X(2) = 0. Y(1) = 0. Y(2) = 0. RETURN ENDIF IF (IDB .GE. 5) CALL MPDEB ('MPCSSX I', A) C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) NCR = 2 ** MCR C C Check if precision level is too low to justify advanced routine. C IF (NW .LE. NCR) THEN CALL MPCSSN (A, PI, X, Y) L1 = 0 GOTO 120 ENDIF C C Check if input is zero. C IF (NA .EQ. 0) THEN X(1) = 1. X(2) = 0. X(3) = 1. Y(1) = 0. Y(2) = 0. L1 = 0 GOTO 120 ENDIF C C Check if Pi has been precomputed. C CALL MPMDC (PI, T1, N1) IF (N1 .NE. 0 .OR. ABS (T1 - CPI) .GT. RX2) THEN IF (KER(30) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPCSSX: PI must be precomputed.') IER = 30 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C N4 = NW + 4 N42 = 2 * N4 NS = 4 * N42 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N42 K2 = K1 + N42 K3 = K2 + N42 F1(1) = 1. F1(2) = 0. F1(3) = 1. NWS = NW C C Reduce argument to between - Pi and Pi. C CALL MPMULD (PI, 2.D0, 0, S(K0)) CALL MPDIVX (A, S(K0), S(K1)) CALL MPNINT (S(K1), S(K2)) CALL MPMUL (S(K2), S(K0), S(K1)) CALL MPSUB (A, S(K1), S(K0)) C C Determine the least integer MQ such that 2 ^ MQ .GE. NW. C T2 = NWS MQ = CL2 * LOG (T2) + 1.D0 - RXX CALL MPEQ (F1, S(K2)) C C Compute initial approximation to [Cos (A), Sin (A)]. C NW = NCR CALL MPCSSN (S(K0), PI, S(K3), S(K3+N4)) IQ = 0 C C Perform the Newton-Raphson iteration with a dynamically changing precision C level NW. C DO 110 K = MCR + 1, MQ NW = MIN (2 * NW, NWS) 100 CONTINUE CALL MPANGX (S(K3), S(K3+N4), PI, S(K1)) CALL MPSUB (S(K0), S(K1), S(K2+N4)) CALL MPCMLX (N4, S(K3), S(K2), S(K1)) CALL MPCEQ (N4, S(K1), S(K3)) IF (K .EQ. MQ - NIT .AND. IQ .EQ. 0) THEN IQ = 1 GOTO 100 ENDIF 110 CONTINUE C C The final (cos, sin) result must be normalized to have magnitude 1. C CALL MPMULX (S(K3), S(K3), S(K0)) CALL MPMULX (S(K3+N4), S(K3+N4), S(K0+N4)) CALL MPADD (S(K0), S(K0+N4), S(K1)) CALL MPSQRX (S(K1), S(K2)) CALL MPDIVX (S(K3), S(K2), S(K0)) CALL MPDIVX (S(K3+N4), S(K2), S(K0+N4)) CALL MPMPCM (N4, S(K0), X, Y) ICS = ISS C 120 IF (IDB .GE. 5) THEN CALL MPDEB ('MPCSSX O', X) CALL MPDEB ('MPCSSX O', Y) ENDIF RETURN END C SUBROUTINE MPCSUB (L, A, B, C) C C This subracts the MPC numbers A and B and returns the MPC difference in C C. L is the offset between real and imaginary parts in A, B and C. L C must be at least NW + 4. Debug output starts with IDB = 9. C C Max SP space for C: 2 * L cells. C DIMENSION A(2*L), B(2*L), C(2*L) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. C(L+1) = 0. C(L+2) = 0. RETURN ENDIF IF (IDB .GE. 9) WRITE (LDB, 1) 1 FORMAT ('MPCSUB') C L1 = L + 1 CALL MPSUB (A, B, C) CALL MPSUB (A(L1), B(L1), C(L1)) C RETURN END C SUBROUTINE MPDEB (CS, A) C C This outputs the character string CS, the exponent of the MP number A, and C the first 50 digits of A, all on one line. CS must either be a literal C string not exceeding 12 characters in length or a variable of type C CHARACTER*n, where n does not exceed 12. C CHARACTER*(*) CS CHARACTER*1 B(160) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) DIMENSION A(NW+2) C IF (IER .NE. 0) RETURN IDS = IDB IDB = 0 NWS = NW NW = MIN (NW, 10) CALL MPOUTC (A, B, N) N = MIN (N, 70) WRITE (LDB, 1) CS, ' ', (B(K), K = 1, 4), (B(K), K = 9, N) 1 FORMAT (A12,67A1:/(79A1)) IDB = IDS NW = NWS RETURN END C SUBROUTINE MPDIV (A, B, C) C C This divides the MP number A by the MP number B to yield the MP quotient C. C For extra high levels of precision, use MPDIVX. Debug output starts with C IDB = 8. C C Max SP space for C: NW + 4 cells. Max DP scratch space: NW + 4 cells. C DOUBLE PRECISION D, RB, SS, T1, T2, T3 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (NDB = 22) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM4/ D(1024) DIMENSION A(NW+2), B(NW+2), C(NW+4) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. RETURN ENDIF IF (IDB .GE. 8) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPDIV I'/(6F12.0)) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 1) (B(I), I = 1, NO) ENDIF C IA = SIGN (1., A(1)) IB = SIGN (1., B(1)) NA = MIN (INT (ABS (A(1))), NW) NB = MIN (INT (ABS (B(1))), NW) C C Check if dividend is zero. C IF (NA .EQ. 0) THEN C(1) = 0. C(2) = 0. GOTO 190 ENDIF IF (NB .EQ. 1 .AND. B(3) .EQ. 1.) THEN C C Divisor is 1 or -1 -- result is A or -A. C C(1) = SIGN (NA, IA * IB) C(2) = A(2) - B(2) C DO 100 I = 3, NA + 2 C(I) = A(I) 100 CONTINUE C GOTO 190 ENDIF C C Check if divisor is zero. C IF (NB .EQ. 0) THEN IF (KER(31) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPDIV: Divisor is zero.') IER = 31 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Initialize trial divisor and trial dividend. C T1 = BDX * B(3) IF (NB .GE. 2) T1 = T1 + B(4) IF (NB .GE. 3) T1 = T1 + RDX * B(5) IF (NB .GE. 4) T1 = T1 + RX2 * B(6) RB = 1.D0 / T1 MD = MIN (NA + NB, NW) D(1) = 0.D0 C DO 110 I = 2, NA + 1 D(I) = A(I+1) 110 CONTINUE C DO 120 I = NA + 2, MD + 4 D(I) = 0.D0 120 CONTINUE C C Perform ordinary long division algorithm. First compute only the first C NA words of the quotient. C> DO 140 J = 2, NA + 1 T1 = INT (RB * (BX2 * D(J-1) + BDX * D(J) + D(J+1))) J3 = J - 3 I2 = MIN (NB, NW + 2 - J3) + 2 C> CDIR$ IVDEP DO 130 I = 3, I2 I3 = I + J3 T2 = D(I3) - T1 * B(I) T3 = INT (T2 * RDX) D(I3) = T2 - T3 * BDX D(I3-1) = D(I3-1) + T3 130 CONTINUE C C If the trial divisor was correct, D(J-1) will be zero. If D(J-1) is not C zero, add it (multiplied by the radix) into D(J). C D(J) = D(J) + BDX * D(J-1) D(J-1) = T1 140 CONTINUE C C Compute additional words of the quotient, as long as the remainder C is nonzero. C> DO 160 J = NA + 2, NW + 3 T1 = INT (RB * (BX2 * D(J-1) + BDX * D(J) + D(J+1))) J3 = J - 3 I2 = MIN (NB, NW + 2 - J3) + 2 IJ = I2 + J3 SS = 0.D0 C> CDIR$ IVDEP DO 150 I = 3, I2 I3 = I + J3 T2 = D(I3) - T1 * B(I) T3 = INT (T2 * RDX) D(I3) = T2 - T3 * BDX D(I3-1) = D(I3-1) + T3 SS = SS + ABS (D(I3-1)) 150 CONTINUE C SS = SS + ABS (D(IJ)) D(J) = D(J) + BDX * D(J-1) D(J-1) = T1 IF (SS .EQ. 0.D0) GOTO 170 IF (IJ .LE. NW + 2) D(IJ+2) = 0.D0 160 CONTINUE C C Set sign and exponent, and fix up result. C J = NW + 3 C 170 D(J) = 0.D0 IF (D(1) .EQ. 0.D0) THEN IS = 1 ELSE IS = 2 ENDIF NC = MIN (J - 1, NW) D(NC+3) = 0.D0 D(NC+4) = 0.D0 C DO 180 I = J + 1, 3, -1 D(I) = D(I-IS) 180 CONTINUE C D(1) = SIGN (NC, IA * IB) D(2) = A(2) - B(2) + IS - 2 CALL MPNORM (C) C 190 IF (IDB .GE. 8) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 3) (C(I), I = 1, NO) 3 FORMAT ('MPDIV O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPDIVD (A, B, N, C) C C This routine divides the MP number A by the DPE number (B, N) to yield C the MP quotient C. Debug output starts with IDB = 9. C C Max SP space for C: NW + 4 cells. Max DP space: NW + 4 cells. C DOUBLE PRECISION B, BB, BR, D, DD, T1 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (NDB = 22) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) COMMON /MPCOM4/ D(1024) DIMENSION A(NW+2), C(NW+4), F(8) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. RETURN ENDIF IF (IDB .GE. 9) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPDIVD I'/(6F12.0)) WRITE (LDB, 2) B, N 2 FORMAT ('MPDIVD I',1PD25.15,I10) ENDIF C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) IB = SIGN (1.D0, B) C C Check if dividend is zero. C IF (NA .EQ. 0) THEN C(1) = 0. C(2) = 0. GOTO 150 ENDIF C C Check if divisor is zero. C IF (B .EQ. 0.D0) THEN IF (KER(32) .NE. 0) THEN WRITE (LDB, 3) 3 FORMAT ('*** MPDIVD: Divisor is zero.') IER = 32 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF N1 = N / NBT N2 = N - NBT * N1 BB = ABS (B) * 2.D0 ** N2 C C Reduce BB to within 1 and BDX. C IF (BB .GE. BDX) THEN C DO 100 K = 1, 100 BB = RDX * BB IF (BB .LT. BDX) THEN N1 = N1 + K GOTO 120 ENDIF 100 CONTINUE C ELSEIF (BB .LT. 1.D0) THEN C DO 110 K = 1, 100 BB = BDX * BB IF (BB .GE. 1.D0) THEN N1 = N1 - K GOTO 120 ENDIF 110 CONTINUE C ENDIF C C If B cannot be represented exactly in a single mantissa word, use MPDIV. C 120 IF (BB .NE. AINT (BB)) THEN BB = SIGN (BB, B) CALL MPDMC (BB, N1 * NBT, F) CALL MPDIV (A, F, C) GOTO 150 ENDIF C BR = 1.D0 / BB DD = A(3) C C Perform short division (not vectorizable at present). Continue as long as C the remainder remains nonzero. C> DO 130 J = 2, NW + 3 T1 = INT (BR * DD) D(J+1) = T1 DD = BDX * (DD - T1 * BB) IF (J .LE. NA) THEN DD = DD + A(J+2) ELSE IF (DD .EQ. 0.D0) GOTO 140 ENDIF 130 CONTINUE C C Set sign and exponent of result. C J = NW + 3 C 140 NC = MIN (J - 1, NW) D(1) = SIGN (NC, IA * IB) D(2) = A(2) - N1 IF (J .LE. NW + 2) D(J+2) = 0.D0 IF (J .LE. NW + 1) D(J+3) = 0.D0 CALL MPNORM (C) C 150 IF (IDB .GE. 9) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 4) (C(I), I = 1, NO) 4 FORMAT ('MPDIVD O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPDIVX (A, B, C) C C This divides the MP number A by the MP number B and returns the MP result C in C. Before calling MPDIVX, the array in MPCOM5 must be initialized by C calling MPINIX. For modest levels of precision, use MPDIV. NW should be C a power of two. The last two words of the result are not reliable. Debug C output starts with IDB = 7. C C Max SP space for C: NW + 4 cells. Max SP scratch space: 2 * NW + 8 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This subroutine employs the following Newton-Raphson iteration, which C converges to 1 / B: C C X_{k+1} = X_k + X_k * (1 - B * X_k) C C Multiplying the final approximation to 1 / B by A gives the quotient. C These iterations are performed with a maximum precision level NW that C is dynamically changed, doubling with each iteration. C C One difficulty with this procedure is that errors often accumulate in the C trailing mantissa words. This error can be controlled by repeating one of C the iterations. The iteration that is repeated is controlled by setting C the parameter NIT below: If NIT = 0, the last iteration is repeated (this C is most effective but most expensive). If NIT = 1, then the next-to-last C iteration is repeated, etc. An extra word of precision cannot be used in C this routine (since MPMULX prefers powers of two), so NIT = 0 or 1 is best C unless the user needs maximum speed. C DOUBLE PRECISION CL2, T1 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (CL2 = 1.4426950408889633D0, NDB = 22, NIT = 1) DIMENSION A(NW+2), B(NW+2), C(NW+4), F(8) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN C(1) = 0. C(2) = 0. RETURN ENDIF IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (A(1))), NDB) + 2 WRITE (LDB, 1) (A(I), I = 1, NO) 1 FORMAT ('MPDIVX I'/(6F12.0)) NO = MIN (INT (ABS (B(1))), NDB) + 2 WRITE (LDB, 1) (B(I), I = 1, NO) ENDIF C IA = SIGN (1., A(1)) IB = SIGN (1., B(1)) NA = MIN (INT (ABS (A(1))), NW) NB = MIN (INT (ABS (B(1))), NW) NCR = 2 ** MCR C C Check if dividend is zero. C IF (NA .EQ. 0) THEN C(1) = 0. C(2) = 0. GOTO 120 ENDIF C C Check if divisor is zero. C IF (NB .EQ. 0) THEN IF (KER(33) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPDIVX: Divisor is zero.') IER = 33 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if precision level of divisor is too low to justify the advanced C routine. C IF (NB .LE. NCR) THEN CALL MPDIV (A, B, C) GOTO 120 ENDIF N4 = NW + 4 NS = 2 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N4 NWS = NW C C Determine the least integer MQ such that 2 ^ MQ .GE. NW. C T1 = NW MQ = CL2 * LOG (T1) + 1.D0 - RXX C C Compute the initial approximation of 1 / B to a precision of NCR words. C NW = NCR F(1) = 1. F(2) = 0. F(3) = 1. CALL MPDIV (F, B, C) IQ = 0 C C Perform the Newton-Raphson iterations described above. C DO 110 K = MCR + 1, MQ AN = NW NW = MIN (2 * NW, NWS) 100 CONTINUE CALL MPMULX (B, C, S(K0)) CALL MPSUB (F, S(K0), S(K1)) S(K1) = SIGN (MIN (ABS (S(K1)), AN), S(K1)) CALL MPMULX (C, S(K1), S(K0)) CALL MPADD (C, S(K0), S(K1)) CALL MPEQ (S(K1), C) IF (K .EQ. MQ - NIT .AND. IQ .EQ. 0) THEN IQ = 1 GOTO 100 ENDIF 110 CONTINUE C C Multiply by A to give final result. C CALL MPMULX (A, C, S(K1)) CALL MPEQ (S(K1), C) ICS = ISS C 120 IF (IDB .GE. 7) THEN NO = MIN (INT (ABS (C(1))), NDB) + 2 WRITE (LDB, 3) (C(I), I = 1, NO) 3 FORMAT ('MPDIVX O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPDMC (A, N, B) C C This routine converts the DPE number (A, N) to MP form in B. All bits of C A are recovered in B. However, note for example that if A = 0.1D0 and N C is 0, then B will NOT be the multiprecision equivalent of 1/10. Debug C output starts with IDB = 9. C C Max SP space for B: 8 cells. C DOUBLE PRECISION A, AA DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (NDB = 22) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS DIMENSION B(NW+4) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. RETURN ENDIF IF (IDB .GE. 9) WRITE (LDB, 1) A, N 1 FORMAT ('MPDMC I',1PD25.15,I10) C C Check for zero. C IF (A .EQ. 0.D0) THEN B(1) = 0. B(2) = 0. GOTO 150 ENDIF N1 = N / NBT N2 = N - NBT * N1 AA = ABS (A) * 2.D0 ** N2 C C Reduce AA to within 1 and BDX. C IF (AA .GE. BDX) THEN C DO 100 K = 1, 100 AA = RDX * AA IF (AA .LT. BDX) THEN N1 = N1 + K GOTO 120 ENDIF 100 CONTINUE C ELSEIF (AA .LT. 1.D0) THEN C DO 110 K = 1, 100 AA = BDX * AA IF (AA .GE. 1.D0) THEN N1 = N1 - K GOTO 120 ENDIF 110 CONTINUE C ENDIF C C Store successive sections of AA into B. C 120 B(2) = N1 B(3) = AINT (AA) AA = BDX * (AA - B(3)) B(4) = AINT (AA) AA = BDX * (AA - B(4)) B(5) = AINT (AA) AA = BDX * (AA - B(5)) B(6) = AINT (AA) B(7) = 0. B(8) = 0. C DO 130 I = 6, 3, -1 IF (B(I) .NE. 0.) GOTO 140 130 CONTINUE C 140 AA = I - 2 B(1) = SIGN (AA, A) C 150 IF (IDB .GE. 9) THEN NO = ABS (B(1)) + 2. WRITE (LDB, 2) (B(I), I = 1, NO) 2 FORMAT ('MPDMC O'/(6F12.0)) ENDIF RETURN END C SUBROUTINE MPEQ (A, B) C C This routine sets the MP number B equal to the MP number A. Debug output C starts with IDB = 10. C C Max SP space for B: NW + 2 cells. C C The fact that only NW + 2 cells, and not NW + 4 cells, are copied is C important in some routines that increase the precision level by one. C PARAMETER (NDB = 22) COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS DIMENSION A(NW+2), B(NW+2) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. RETURN ENDIF IF (IDB .GE. 10) WRITE (LDB, 1) 1 FORMAT ('MPEQ') C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) B(1) = SIGN (NA, IA) C DO 100 I = 2, NA + 2 B(I) = A(I) 100 CONTINUE C RETURN END C SUBROUTINE MPEXP (A, AL2, B) C C This computes the exponential function of the MP number A and returns the C MP result in B. AL2 is the MP value of Log(2) produced by a prior call C to MPLOG. For extra high levels of precision, use MPEXPX. The last C word of the result is not reliable. Debug output starts with IDB = 7. C C Max SP space for B: NW + 4 cells. Max SP scratch space: 5 * NW + 25 C cells. Max DP scratch space: NW + 5 cells. C C This routine uses a modification of the Taylor's series for Exp (t): C C Exp (t) = (1 + r + r^2 / 2! + r^3 / 3! + r^4 / 4! ...) ^ q * 2 ^ n C C where q = 256, r = t' / q, t' = t - n Log(2) and where n is chosen so C that -0.5 Log(2) < t' <= 0.5 Log(2). Reducing t mod Log(2) and C dividing by 256 insures that -0.001 < r <= 0.001, which accelerates C convergence in the above series. C DOUBLE PRECISION ALT, T1, T2 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (ALT = 0.693147180559945309D0, NQ = 8) DIMENSION A(NW+2), B(NW+5), AL2(NW+2), F(8) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. RETURN ENDIF IF (IDB .GE. 7) CALL MPDEB ('MPEXP I', A) C IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) CALL MPMDC (A, T1, N1) T1 = T1 * 2.D0 ** N1 C C Unless the argument is near Log (2), Log(2) must be precomputed. This C exception is necessary because MPLOG calls MPEXP to initialize Log (2). C IF (ABS (T1 - ALT) .GT. RDX) THEN CALL MPMDC (AL2, T2, N2) IF (N2 .NE. - NBT .OR. ABS (T2 * 0.5D0 ** NBT - ALT) .GT. RX2) $ THEN IF (KER(34) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPEXP: LOG (2) must be precomputed.') IER = 34 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF ENDIF C C Check for overflows and underflows. C IF (T1 .GE. 1D9) THEN IF (T1 .GT. 0.D0) THEN IF (KER(35) .NE. 0) THEN WRITE (LDB, 2) T1, N1 2 FORMAT ('*** MPEXP: Argument is too large',F12.6,' x 10 ^', $ I8) IER = 35 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ELSE B(1) = 0. B(2) = 0. L1 = 0 GOTO 130 ENDIF ENDIF C N5 = NW + 5 NS = 4 * N5 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N5 K2 = K1 + N5 K3 = K2 + N5 NWS = NW NW = NW + 1 F(1) = 1. F(2) = 0. F(3) = 1. C C Compute the reduced argument A' = A - Log(2) * Nint [A / Log(2)]. Save C NZ = Nint [A / Log(2)] for correcting the exponent of the final result. C IF (ABS (T1 - ALT) .GT. RDX) THEN CALL MPDIV (A, AL2, S(K0)) CALL MPNINT (S(K0), S(K1)) CALL MPMDC (S(K1), T1, N1) NZ = T1 * 2.D0 ** N1 + SIGN (RXX, T1) CALL MPMUL (AL2, S(K1), S(K2)) CALL MPSUB (A, S(K2), S(K0)) ELSE CALL MPEQ (A, S(K0)) NZ = 0 ENDIF TL = S(K0+1) - NW C C Check if the reduced argument is zero. C IF (S(K0) .EQ. 0.D0) THEN S(K0) = 1. S(K0+1) = 0. S(K0+2) = 1. L1 = 0 GOTO 120 ENDIF C C Divide the reduced argument by 2 ^ NQ. C CALL MPDIVD (S(K0), 1.D0, NQ, S(K1)) C C Compute Exp using the usual Taylor series. C CALL MPEQ (F, S(K2)) CALL MPEQ (F, S(K3)) L1 = 0 C 100 L1 = L1 + 1 IF (L1 .EQ. 10000) THEN IF (KER(36) .NE. 0) THEN WRITE (LDB, 3) 3 FORMAT ('*** MPEXP: Iteration limit exceeded.') IER = 36 IF (KER(IER) .EQ. 2) CALL MPABRT NW = NWS ICS = ISS RETURN ENDIF ENDIF C T2 = L1 CALL MPMUL (S(K2), S(K1), S(K0)) CALL MPDIVD (S(K0), T2, 0, S(K2)) CALL MPADD (S(K3), S(K2), S(K0)) CALL MPEQ (S(K0), S(K3)) C C Check for convergence of the series. C IF (S(K2) .NE. 0. .AND. S(K2+1) .GE. TL) GOTO 100 C C Raise to the (2 ^ NQ)-th power. C DO 110 I = 1, NQ CALL MPMUL (S(K0), S(K0), S(K1)) CALL MPEQ (S(K1), S(K0)) 110 CONTINUE C C Multiply by 2 ^ NZ. C 120 CALL MPMULD (S(K0), 1.D0, NZ, S(K1)) CALL MPEQ (S(K1), B) C C Restore original precision level. C NW = NWS ICS = ISS CALL MPROUN (B) C 130 IF (IDB .GE. 7) THEN WRITE (LDB, 4) L1 4 FORMAT ('Iteration count:',I5) CALL MPDEB ('MPEXP O', B) ENDIF RETURN END C SUBROUTINE MPEXPX (A, PI, AL2, B) C C This computes the exponential function of the MP number A and returns the C MP result in B. PI is the MP value of Pi produced by a prior call to MPPI C or MPPIX. AL2 is the MP value of Log(2) produced by a prior call to C MPLOG or MPLOGX. Before calling MPEXPX, the array in MPCOM5 must be C initialized by calling MPINIX. NW should be a power of two. For modest C levels of precision, use MPEXP. The last four words of the result are C not reliable. Debug output starts with IDB = 5. C C Max SP space for B: NW + 4 cells. Max SP scratch space: 12 * NW + 54 C cells. Max DP scratch space: 12 * NW + 6 cells. C C This routine uses the Newton iteration C C b_{k+1} = b_k [a + 1 - log b_k] C C with a dynamically changing level of precision. Logs are performed using C MPLOGX. See the comment about the parameter NIT in MPDIVX. C DOUBLE PRECISION ALT, CL2, CPI, T1, T2 DOUBLE PRECISION BBX, BDX, BX2, RBX, RDX, RX2, RXX PARAMETER (ALT = 0.693147180559945309D0, $ CL2 = 1.4426950408889633D0, CPI = 3.141592653589793238D0, $ NIT = 1) DIMENSION A(NW+2), AL2(NW+2), B(NW+4), F1(8), PI(NW+2) COMMON /MPCOM0/ BBX, BDX, BX2, RBX, RDX, RX2, RXX, NBT, NPR COMMON /MPCOM1/ NW, IDB, LDB, IER, MCR, IRD, ICS, IHS, IMS COMMON /MPCOM2/ KER(72) COMMON /MPCOM3/ S(1024) C IF (IER .NE. 0) THEN B(1) = 0. B(2) = 0. RETURN ENDIF IF (IDB .GE. 5) CALL MPDEB ('MPEXPX I', A) C NCR = 2 ** MCR IA = SIGN (1., A(1)) NA = MIN (INT (ABS (A(1))), NW) CALL MPMDC (A, T1, N1) T1 = T1 * 2.D0 ** N1 C C Check if precision level is too low to justify the advanced routine. C IF (NW .LE. NCR) THEN CALL MPEXP (A, AL2, B) GOTO 120 ENDIF C C Check if Log(2) has been precomputed. C CALL MPMDC (AL2, T2, N2) IF (N2 .NE. - NBT .OR. ABS (T2 * 0.5D0 ** NBT - ALT) .GT. RX2) $ THEN IF (KER(37) .NE. 0) THEN WRITE (LDB, 1) 1 FORMAT ('*** MPEXPX: LOG (2) must be precomputed.') IER = 37 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check if Pi has been precomputed. C CALL MPMDC (PI, T2, N2) IF (N2 .NE. 0 .OR. ABS (T2 - CPI) .GT. RX2) THEN IF (KER(38) .NE. 0) THEN WRITE (LDB, 2) 2 FORMAT ('*** MPEXPX: PI must be precomputed.') IER = 38 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ENDIF C C Check for overflows and underflows. C IF (T1 .GE. 1D9) THEN IF (T1 .GT. 0.D0) THEN IF (KER(39) .NE. 0) THEN WRITE (LDB, 3) T1, N1 3 FORMAT ('*** MPEXPX: Argument is too large',F12.6,' x 10 ^', $ I8) IER = 39 IF (KER(IER) .EQ. 2) CALL MPABRT ENDIF RETURN ELSE B(1) = 0. B(2) = 0. GOTO 120 ENDIF ENDIF C N4 = NW + 4 NS = 3 * N4 ISS = ICS ICS = ICS + NS IHS = MAX (ICS, IHS) IF (ICS - 1 .GT. IMS) CALL MPALER K0 = ISS K1 = K0 + N4 K2 = K1 + N4 NWS = NW F1(1) = 1. F1(2) = 0. F1(3) = 1. C C Determine the least integer MQ such that 2 ^ MQ .GE. NW. C T2 = NWS MQ = CL2 * LOG (T2) + 1.D0 - RXX CALL MPADD (A, F1, S(K0)) C C Compute initial approximation to Exp (A). C NW = NCR CALL MPEXP (A, AL2, B) IQ = 0 C C Perform the Newton-Raphson iteration described above with a dynamically C changing precision level NW. C DO 110 K = MCR + 1, MQ NW = MIN (2 * NW, NWS) 100 CONTINUE CALL MPLOGX (B, PI, AL2, S(K1)) CALL MPSUB (S(K0), S(K1), S(K2)) CALL MPMULX (B, S(K2), S(K1)) CALL MPEQ (S(K1), B) IF (K .EQ. MQ - NIT .AND. IQ .EQ. 0) THEN IQ = 1 GOTO 100 ENDIF 110 CONTINUE C ICS = ISS C 120 IF (IDB .GE. 6) CALL MPDEB ('MPEXPX O', B) RETURN END C SUBROUTINE MPFFT1 (IS, L, M, X, Y) C C Performs the L-th iteration of the first variant of the Stockham FFT. C This routine is called by MPCFFT. It is not intended to be called directly C by the user. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION X(*), Y(*) COMMON /MPCOM5/ U(1024) C C Set initial parameters. C N = 2 ** M K = U(1) NU = K / 64 N1 = N / 2 LK = 2 ** (L - 1) LI = 2 ** (M - L) LJ = 2 * LI KU = LI + 1 KN = KU + NU C DO 100 K = 0, LK - 1 I11 = K * LJ + 1 I12 = I11 + LI I21 = K * LI + 1 I22 = I21 + N1 C CDIR$ IVDEP DO 100 I = 0, LI - 1 U1 = U(KU+I) U2 = IS * U(KN+I) X11 = X(I11+I) X12 = X(I11+I+N) X21 = X(I12+I) X22 = X(I12+I+N) T1 = X11 - X21 T2 = X12 - X22 Y(I21+I) = X11 + X21 Y(I21+I+N) = X12 + X22 Y(I22+I) = U1 * T1 - U2 * T2 Y(I22+I+N) = U1 * T2 + U2 * T1 100 CONTINUE C RETURN END C SUBROUTINE MPFFT2 (IS, L, M, X, Y) C C Performs the L-th iteration of the second variant of the Stockham FFT. C This routine is called by MPCFFT. It is not intended to be called directly C by the user. C IMPLICIT DOUBLE PRECISION (A-H, O-Z) DIMENSION X(*), Y(*) COMMON /MPCOM5/ U(1024) C C Set initial parameters. C N = 2 ** M K = U(1) NU = K / 64 N1 = N / 2 LK = 2 ** (L - 1) LI = 2 ** (M - L) LJ = 2 * LK KU = LI + 1 C DO 100 I = 0, LI - 1 I11 = I * LK + 1 I12 = I11 + N1 I21 = I * LJ + 1 I22 = I21 + LK U1 = U(KU+I) U2 = IS * U(KU+I+NU) C CDIR$ IVDEP DO 100 K = 0, LK - 1 X11 = X(I11+K) X12 = X(I11+K+N) X21 = X(I12+K) X22 = X(I12+K+N) T1 = X11 - X21 T2 = X12 - X22 Y(I21+K) = X11 + X21 Y(I21+K+N) = X12 + X22 Y(I22+K) = U1 * T1 - U2 * T2 Y(I22+K+N) = U1 * T2 + U2 * T1 100 CONTINUE C RETURN END C SUBROUTINE MPINFR (A, B, C) C C Sets B to the integer part of the MP number A and sets C equal to the C fractional part of A. Note that if A = -3.3, then B = -3 and C = -0.3. C Debug output starts with IDB = 9. C C Max SP space for B and C: NW + 4 cells. C PARAMETER (NDB = 22)