API tools faq
paste
Advertisement
Guest User

NASA/Cray PITEST PI4.f Double - Quartically Convergent

a guest
May 20th, 2018
280
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 63.24 KB | None | 0 0
raw download clone embed print report
  1. C*****************************************************************************
  2. C
  3. C   MULTIPRECISION FLOATING POINT COMPUTATION PACKAGE
  4. C   DOUBLE PRECISION RADIX 2 COMPLEX FFT VERSION
  5. C
  6. C   This version has been modified for higher performance when the actual level
  7. C   of precision of input numbers (i.e., the index of the last nonzero words)
  8. C   varies widely.  The advanced algorithms in MPMULX, MPDIVX, and MPSQRX
  9. C   are only performed when the actual precision levels of the input numbers
  10. C   are above a certain level.  For this reason, it is recommended that users
  11. C   always call MPMULX, MPDIVX, or MPSQRX, and not directly call MPMUL, MPDIV,
  12. C   or MPSQRT.
  13. C
  14. C   Version Date:  June 23, 1988
  15. C
  16. C   Author:
  17. C      David H. Bailey                 Telephone:   415-694-4410
  18. C      NASA Ames Research Center       E-mail:  dbailey@ames-nas.arpa
  19. C      Mail Stop 258-5
  20. C      Moffett Field, CA  94035
  21. C
  22. C   Introduction:
  23. C
  24. C       This package of subroutines performs multiprecision floating point
  25. C   arithmetic.  It employs the most advanced algorithms known as of the
  26. C   above date, and all programs have been written to effectively utilize
  27. C   the vector processing capabilities of a supercomputer.  This package
  28. C   assumes that all double precision variables represent eight bytes of
  29. C   storage, with at least 46 mantissa bits.  No other machine dependent
  30. C   assumptions have been made.
  31. C
  32. C   IMPORTANT NOTE: ALL floating-point variables in the following package are
  33. C   assumed to be double-precision (64 bit) floating-point variables.  In
  34. C   addition, double precision data is to be inferred in ALL references to
  35. C   "cells", "locations", etc.  In particular, ALL scratch space requirements
  36. C   are to be interpreted in units of 64-bit storage locations.
  37. C
  38. C       This package includes the following subroutines  (SP denotes
  39. C   ordinary single precision, MP denotes multiprecision (NW mantissa words),
  40. C   and DMP denotes double multiprecision (2 * NW mantissa words):
  41. C
  42. C   MPSMC        Performs SP to MP conversion.
  43. C   MPMSC        Performs MP to SP conversion.
  44. C   MPEQ         Sets one MP number equal to another.
  45. C   MPADD        Adds two MP numbers to yield MP sum.
  46. C   MPSUB        Subtracts two MP numbers to yield MP difference.
  47. C   MPMUL1       Multiplies a MP number by a SP number to yield MP product.
  48. C   MPDIV1       Divides a MP number by a MP divisor to yield SP quotient.
  49. C   MPDIV2       Divides a MP number by a SP divisor to yield MP quotient.
  50. C   MPNORM       Normalizes a MP number to standard MP form.
  51. C   MPRAND       Returns a pseudo-random MP number between 0 and 1.
  52. C   MPINFR       Returns the integer and fractional part of a MP number.
  53. C   MPFMT        Formats a MP number for output.
  54. C   MPCFFT       Initializes the array U in common MPCOM2 and performs FFT.
  55. C   MPMULX       Multiplies two MP numbers to yield DMP product.
  56. C   MPDIVX       Divides a MP number by a MP divisor to yield MP quotient.
  57. C   MPSQRX       Computes the square root of a MP number.
  58. C
  59. C   A number of other subroutines are part of the package but are not listed
  60. C   here because they are not intended to be directly called by users.
  61. C
  62. C        The format of the multiprecision MP numbers is as follows:  the
  63. C   first word contains the sign (1., 0. or -1.), and the second word
  64. C   contains the exponent (powers of 10^6).  Words numbered 3 up to NW + 2
  65. C   contain the mantissa.  If the sign is zero, then all other words must
  66. C   be zero.  If the sign is nonzero, then word 3 must be nonzero.  The
  67. C   decimal point is assumed after the first mantissa word for numbers with
  68. C   with zero exponent.  The radix of this representation is 10^6, so that
  69. C   each mantissa word contains a floating whole number between 0 and
  70. C   999999.  The decimal radix 10^6 has been used instead of a binary radix
  71. C   so that expensive binary to decimal conversion is not necessary for output.
  72. C
  73. C        The format of the double multiprecision DMP numbers is the same
  74. C   as the MP numbers except that DMP numbers have twice as many mantissa
  75. C   words, for a total of 2 * NW + 2 words.  Note that only MPMULX explicitly
  76. C   computes with DMP numbers.  Other operations may be easily performed on
  77. C   DMP numbers by merely doubling NW and, if necessary, increasing MW by one.
  78. C   If this is done be sure that the scratch arrays in calls to MPMULX,
  79. C   MPDIVX, and MPSQRX have sufficient scratch space, and remember to restore
  80. C   NW and MW to their previous values when the DMP computations are complete.
  81. C
  82. C        To initialize the package, the variables ND, NW, MW, IDB, LDB, and NDB
  83. C   must be set in the common block MPCOM1.  ND is the requested number of
  84. C   digits of precision.  NW is the number of words of precision, which should
  85. C   be at least ND / 6 + 1.  All MP variables must be dimensioned at least
  86. C   NW + 2 and DMP variables must be dimensioned at least 2 * NW + 2 (twice
  87. C   the amount for MP variables will suffice).  MW is used in subroutines
  88. C   MPMULX, MPDIVX, and MPSQRX.  When ND is less than about 200 (depending on
  89. C   the system), MW may be set to zero.  For higher levels of precision, the
  90. C   only permissible values of of NW are 2^(MW - 2), and both MW and NW must
  91. C   be set in common MPCOM1 by the user.  In addition, prior to calling
  92. C   MPMULX, MPDIVX, or MPSQRX for these higher levels of precision, the user
  93. C   must first allocate at least 2 * 2^MX cells in the array U of common block
  94. C   MPCOM2 in the user's main program, and the user must initialize the FFT
  95. C   routine by calling MPCFFT with 0 as the first argument.  Here MX is the
  96. C   largest value of MW that will be used.
  97. C
  98. C        The variable IDB is a debug flag and ordinarily should be set to 0.
  99. C   Setting IDB to an integer between 5 and 9 (or greater) produces debug
  100. C   printouts in varying degrees from the MP subroutines.  Values of IDB
  101. C   between 1 and 4 are available for use as debug flags in the user's
  102. C   calling program if desired.  LDB is the logical unit number for output of
  103. C   these debug messages, and NDB is the number of words output in the debug
  104. C   printout of a MP number.  Typically LDB is set to 6 and NDB is set to 16.
  105. C
  106. C        Common data has been divided into two separate common blocks to
  107. C   facilitate efficient multitasking.  The parameters in MPCOM1 may be
  108. C   frequently changed by a user program and thus should be declared TASK
  109. C   COMMON or the equivalent.  The array U in common MPCOM2, however,
  110. C   is not intended to be altered after it is first initialized.  For this
  111. C   reason, and since U may be a very large array, it makes more sense for
  112. C   MPCOM2 to be an ordinary global common.  Except for the data in these
  113. C   two common blocks, no variables are global or need to be "saved".
  114. C
  115. C*****************************************************************************
  116. C
  117.       SUBROUTINE MPSMC (A, N, B)
  118. C
  119. C   Converts the SP number  A * 10^N  to MP form in  B.
  120. C
  121. C   Note:  The result should be exact provided N is 0 or divisible by IB
  122. C   and B is an exact binary fraction in the range RX^2 .LT. |B| .LT. BX^2.
  123. C
  124.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  125.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, RX2 = RX * RX,
  126.      $  RXX = 0.5D0 * RX)
  127.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  128.       DIMENSION B(1)
  129. C
  130.       NX = NW + 2
  131.       IF (IDB .GE. 8)  WRITE (LDB, 1)  A, N, NX
  132. 1     FORMAT (/'MPSMC I',1PD25.15,2I10)
  133. C
  134. C   Check for zero.
  135. C
  136.       IF (A .EQ. 0.D0)  THEN
  137.         DO 100 I = 1, NX
  138.           B(I) = 0.D0
  139. 100     CONTINUE
  140.         GOTO 200
  141.       ENDIF
  142. C
  143. C   Adjust exponent value.
  144. C
  145.       AA = ABS (A)
  146.       AL = N + LOG10 (AA) * (1.D0 + RX2)
  147.       M = AL / IB
  148.       IF (AL .LT. 0.D0 .AND. AL .NE. M * IB)  M = M - 1
  149.       AM = AA * 10.D0 ** (IB - IB * M + N)
  150. C
  151. C   Convert constant to 2-word precision.
  152. C
  153.       B(1) = SIGN (1.D0, A)
  154.       B(2) = M
  155.       B(3) = AINT (RX * AM + RXX)
  156.       B(4) = AINT (AM - BX * B(3))
  157. C
  158.       DO 110 I = 5, NX
  159.         B(I) = 0.D0
  160. 110   CONTINUE
  161. C
  162. 200   IF (IDB .GE. 8)  WRITE (LDB, 2)  (B(I),I=1,NDB)
  163. 2     FORMAT ('MPSMC O'/(8F9.0))
  164.       RETURN
  165.       END
  166. C
  167.       SUBROUTINE MPMSC (A, B, N)
  168. C
  169. C   Converts the MP number  A  to the SP form  B * 10^N , accurate to about
  170. C   12 decimal places.
  171. C
  172.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  173.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, RX2 = RX * RX)
  174.       DIMENSION A(1)
  175.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  176. C
  177.       NX = NW + 2
  178.       IF (IDB .GE. 8)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  179. 1     FORMAT (/'MPMSC I',I10/(8F9.0))
  180. C
  181.       IF (A(1) .EQ. 0.D0)  THEN
  182.         B = 0.D0
  183.         N = 0
  184.         GOTO 200
  185.       ENDIF
  186. C
  187.       AA = A(3) + A(4) * RX + A(5) * RX2
  188. C
  189. C   Adjust results so that  B  is in the range  1 .LE. |B| .LT. 10.
  190. C
  191.       N2 = LOG10 (AA) + RX2
  192.       B = SIGN (AA * 10.D0 ** (-N2), A(1))
  193.       N = N2 + IB * A(2)
  194. C
  195. 200   IF (IDB .GE. 8)  WRITE (LDB, 2)  B, N
  196. 2     FORMAT ('MPMSC O',F10.0,I10)
  197.       RETURN
  198.       END
  199. C
  200.       SUBROUTINE MPEQ (A, B)
  201. C
  202. C   Sets the MP number  B  equal to the MP number  A.
  203. C
  204.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  205.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  206.       DIMENSION A(1), B(1)
  207. C
  208.       IF (IDB .GE. 8)  WRITE (LDB, 1)
  209. 1     FORMAT (/'MPEQ')
  210. C
  211.       DO 100 I = 1, NW + 2
  212.         B(I) = A(I)
  213. 100   CONTINUE
  214.       RETURN
  215.       END
  216. C
  217.       SUBROUTINE MPADD (A, B, C)
  218. C
  219. C   Adds MP numbers  A  and  B  to yield the MP sum  C.
  220. C
  221.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  222.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  223.       DIMENSION A(1), B(1), C(1)
  224. C
  225.       NX = NW + 2
  226.       IF (IDB .GE. 7)  THEN
  227.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  228. 1       FORMAT (/'MPADD I',I10/(8F9.0))
  229.         WRITE (LDB, 1)  NX, (B(I),I=1,NDB)
  230.       ENDIF
  231. C
  232.       IF (A(1) * B(1) .EQ. 0.D0)  THEN
  233. C
  234. C   A or B is zero.
  235. C
  236.         IF (A(1) .EQ. 0.D0)  THEN
  237.           DO 80 I = 1, NX
  238.             C(I) = B(I)
  239. 80        CONTINUE
  240.           GOTO 300
  241.         ELSE
  242.           DO 90 I = 1, NX
  243.             C(I) = A(I)
  244. 90        CONTINUE
  245.           GOTO 300
  246.         ENDIF
  247. C
  248.       ELSEIF (A(1) * B(1) .GT. 0.D0)  THEN
  249. C
  250. C   A  and  B  have the same sign or one is zero -- perform addition
  251. C   after shifting to match exponents.
  252. C
  253.         C(1) = A(1)
  254.         C(2) = MAX (A(2), B(2))
  255. C
  256.         IF (A(2) .GE. B(2))  THEN
  257.           N = MIN (INT (A(2) - B(2)), NW)
  258.           DO 100 I = 3, N+2
  259.             C(I) = A(I)
  260. 100       CONTINUE
  261.           DO 110 I = N+3, NX
  262.             C(I) = A(I) + B(I-N)
  263. 110       CONTINUE
  264. C
  265.         ELSE
  266.           N = MIN (INT (B(2) - A(2)), NW)
  267.           DO 120 I = 3, N+2
  268.             C(I) = B(I)
  269. 120       CONTINUE
  270.           DO 130 I = N+3, NX
  271.             C(I) = A(I-N) + B(I)
  272. 130       CONTINUE
  273.         ENDIF
  274. C
  275.       ELSE
  276. C
  277. C   A  and  B  have different nonzero signs -- perform subtraction.
  278. C
  279.         IF (A(2) .EQ. B(2))  THEN
  280. C
  281. C   A  and  B  have the same exponent -- search for unequal cells.
  282. C
  283.           C(2) = 0.D0
  284.           DO 140 I = 3, NX
  285.             IF (A(I) .NE. B(I))  GOTO 160
  286. 140       CONTINUE
  287. C
  288. C   ABS (A)  =  ABS (B).  The result is zero.
  289. C
  290.           DO 150 I = 1, NX
  291.             C(I) = 0.D0
  292. 150       CONTINUE
  293.           GOTO 300
  294. C
  295. 160       K = I
  296. C
  297. C   ABS (A)  .NE.  ABS (B)  but exponents are equal -- subtract and normalize.
  298. C
  299.           IF (A(K) .GT. B(K))  THEN
  300.             K = K - 3
  301.             C(1) = A(1)
  302.             C(2) = A(2) - K
  303.             DO 170 I = 3, NX - K
  304.               C(I) = A(I+K) - B(I+K)
  305. 170         CONTINUE
  306. C
  307.           ELSE
  308.             K = K - 3
  309.             C(1) = B(1)
  310.             C(2) = B(2) - K
  311.             DO 180 I = 3, NX - K
  312.               C(I) = B(I+K) - A(I+K)
  313. 180         CONTINUE
  314.           ENDIF
  315. C
  316.           DO 190 I = NX - K + 1, NX
  317.             C(I) = 0.D0
  318. 190       CONTINUE
  319. C
  320.         ELSE
  321. C
  322. C   Exponents are different -- subtract with smaller shifted right.
  323. C
  324.           IF (A(2) .GT. B(2))  THEN
  325.             C(1) = A(1)
  326.             C(2) = A(2)
  327.             N = MIN (INT (A(2) - B(2)), NW)
  328.             DO 200 I = 3, N+2
  329.               C(I) = A(I)
  330. 200         CONTINUE
  331.             DO 210 I = N+3, NX
  332.               C(I) = A(I) - B(I-N)
  333. 210         CONTINUE
  334. C
  335.           ELSE
  336.             C(1) = B(1)
  337.             C(2) = B(2)
  338.             N = MIN (INT (B(2) - A(2)), NW)
  339.             DO 220 I = 3, N+2
  340.               C(I) = B(I)
  341. 220         CONTINUE
  342.             DO 230 I = N+3, NX
  343.               C(I) = B(I) - A(I-N)
  344. 230         CONTINUE
  345. C
  346.           ENDIF
  347.         ENDIF
  348.       ENDIF
  349. C
  350. C   Fix up result (possible negative numbers or overflow).
  351. C
  352.       CALL MPNORM (C)
  353. C
  354. 300   IF (IDB .GE. 7)  WRITE (LDB, 2)  (C(I),I=1,NDB)
  355. 2     FORMAT ('MPADD O'/(8F9.0))
  356.       RETURN
  357.       END
  358. C
  359.       SUBROUTINE MPSUB (A, B, C)
  360. C
  361. C   Subtracts MP numbers  A  and  B  to yield the MP difference  C.
  362. C
  363. C   Note:  This does not work if  A  and  B  are the same variable, i.e.
  364. C          CALL MPSUB (X, X, Y)  .
  365. C
  366.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  367.       DIMENSION  A(1), B(1), C(1)
  368.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  369. C
  370.       IF (IDB .GE. 7)  WRITE (LDB, 1)
  371. 1     FORMAT (/'MPSUB')
  372. C
  373. C   Negate  B  and add.
  374. C
  375.       B(1) = - B(1)
  376.       CALL MPADD (A, B, C)
  377. C
  378. C   Restore the sign of  B.
  379. C
  380.       B(1) = - B(1)
  381.       RETURN
  382.       END
  383. C
  384.       SUBROUTINE MPMUL (A, B, C, NA, NB)
  385. C
  386. C   Multiplies two MP numbers  A  and  B  to yield the exact DMP product  C.
  387. C   Note that the array  C  must contain at least  2 * NW  cells.  NA and NB
  388. C   are the indices of the last nonzero words of A and B.  It is recommended
  389. C   that MPMUL is only called by MPMULX, not directly by the user.
  390. C
  391.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  392.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6)
  393.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  394.       DIMENSION A(1), B(1), C(1)
  395. C
  396.       NX = NW + 2
  397.       IF (IDB .GE. 6)  THEN
  398.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  399. 1       FORMAT (/'MPMUL I',I10/(8F9.0))
  400.         WRITE (LDB, 1)  NX, (B(I),I=1,NDB)
  401.       ENDIF
  402. C
  403.       NWS = NW
  404.       NW2 = 2 * NW
  405.       NX2 = NW2 + 2
  406. C
  407.       DO 100 I = 1, NX2
  408.         C(I) = 0.D0
  409. 100   CONTINUE
  410.       IF (A(1) * B(1) .EQ. 0.D0)  GOTO 200
  411. C
  412. C   Ordinary long multiplication algorithm.
  413. C
  414.       DO 170 J = 3, NA
  415.         AJ = A(J)
  416.         J3 = J - 3
  417. C
  418.         DO 150 I = 3, NB
  419.           C(I+J3) = C(I+J3) + AJ * B(I)
  420. 150     CONTINUE
  421. C
  422. C   Release carries periodically.
  423. C
  424.         IF (J .EQ. NA .OR. J .EQ. 64 * (J / 64)) THEN
  425.           K1 = MAX (3, J - 64)
  426.           K2 = NB + J3
  427. C
  428.           DO 160 I = K1, K2
  429.             T = C(I)
  430.             R = AINT (T * RX)
  431.             C(I) = T - R * BX
  432.             C(I-1) = C(I-1) + R
  433. 160       CONTINUE
  434.         ENDIF
  435. 170   CONTINUE
  436. C
  437.       R = C(2)
  438.       C(1) = A(1) * B(1)
  439.       C(2) = A(2) + B(2)
  440.       C(3) = C(3) + R * BX
  441. C
  442. C   Fix up result (some words may exceed BX).
  443. C
  444.       NW = NW2
  445.       CALL MPNORM (C)
  446.       NW = NWS
  447. C
  448. 200   IF (IDB .GE. 6)  WRITE (LDB, 2)  (C(I),I=1,NDB)
  449. 2     FORMAT ('MPMUL O'/(8F9.0))
  450.       RETURN
  451.       END
  452. C
  453.       SUBROUTINE MPMUL1 (A, B, N, C)
  454. C
  455. C   Multiplies the MP number A by the SP number B * 10^N to yield the
  456. C   MP product C.
  457. C
  458. C   Note:  The result should be exact provided N is 0 or divisible by IB
  459. C   and B is an exact binary fraction in the range RX^2 .LT. |B| .LT. BX^2.
  460. C
  461.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  462.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, RX2 = RX * RX)
  463.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  464.       DIMENSION A(1), C(1), BB(2)
  465. C
  466.       NX = NW + 2
  467.       IF (IDB .GE. 7)  THEN
  468.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  469. 1       FORMAT (/'MPMUL1 I',I10/(8F9.0))
  470.         WRITE (LDB, 2)  B, N
  471. 2       FORMAT ('MPMUL1 I',1PD25.15,I10)
  472.       ENDIF
  473. C
  474. C   Check for zero inputs.
  475. C
  476.       DO 100 I = 1, NX
  477.         C(I) = 0.D0
  478. 100   CONTINUE
  479.       IF (A(1) * B .EQ. 0.D0)  GOTO 200
  480. C
  481. C   Adjust exponent value.
  482. C
  483.       BA = ABS (B)
  484.       BL = N + LOG10 (BA) * (1.D0 + RX2)
  485.       M = BL / IB
  486.       IF (BL .LT. 0.D0 .AND. BL .NE. M * IB)  M = M - 1
  487.       BM = BA * 10.D0 ** (IB - IB * M + N)
  488.       BB(1) = AINT (RX * BM)
  489.       BB(2) = AINT (BM - BX * BB(1))
  490. C
  491. C   Multiply with 2-word precision.
  492. C
  493.       DO 120 J = 1, 2
  494.         BJ = BB(J)
  495.         J1 = J - 1
  496. C
  497.         DO 110 I = 3, NX - J1
  498.           C(I+J1) = C(I+J1) + A(I) * BJ
  499. 110     CONTINUE
  500. 120   CONTINUE
  501. C
  502. C   Release the large carries.
  503. C
  504.       DO 130 I = 3, NX
  505.         T = C(I)
  506.         R = AINT (T * RX)
  507.         C(I) = T - R * BX
  508.         C(I-1) = C(I-1) + R
  509. 130   CONTINUE
  510. C
  511.       R = C(2)
  512.       C(1) = A(1) * SIGN (1.D0, B)
  513.       C(2) = A(2) + M
  514.       C(3) = C(3) + R * BX
  515.       CALL MPNORM (C)
  516. C
  517. 200   IF (IDB .GE. 7)  WRITE (LDB, 3)  (C(I),I=1,NDB)
  518. 3     FORMAT ('MPMUL1 O'/(8F9.0))
  519.       RETURN
  520.       END
  521. C
  522.       SUBROUTINE MPDIV (A, B, C, D)
  523. C
  524. C   Divides the MP number  A  by the MP number  B  to yield the MP
  525. C   quotient  C.   D  is a scratch array with at least  NW  locations.
  526. C
  527.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  528.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, BX2 = BX * BX)
  529.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  530.       DIMENSION A(1), B(1), C(1), D(1)
  531. C
  532.       NX = NW + 2
  533.       IF (IDB .GE. 6)  THEN
  534.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  535. 1       FORMAT (/'MPDIV I',I10/(8F9.0))
  536.         WRITE (LDB, 1)  NX, (B(I),I=1,NDB)
  537.       ENDIF
  538. C
  539. C   Check if dividend is zero.
  540. C
  541.       IF (A(1) .EQ. 0.D0)  THEN
  542.         DO 100 I = 1, NX
  543.           C(I) = 0.D0
  544. 100     CONTINUE
  545.         GOTO 200
  546.       ENDIF
  547. C
  548. C   Check if divisor is zero.
  549. C
  550.       IF (B(1) .EQ. 0.D0)  THEN
  551.         WRITE (LDB, 2)
  552. 2       FORMAT ('*** MPDIV: ZERO DIVISOR ***')
  553.         CALL XABORT
  554.       ENDIF
  555. C
  556. C   Initialize trial dividend.
  557. C
  558.       D(1) = 0.D0
  559.       DO 110 I = 2, NW + 1
  560.         D(I) = A(I+1)
  561. 110   CONTINUE
  562.       D(NX) = 0.D0
  563. C
  564. C   Compute trial divisor.
  565. C
  566.       RB = 1.D0 / (BX * B(3) + B(4) + RX * B(5))
  567. C
  568.       DO 130 J = 2, NX - 1
  569. C
  570. C   Compute current quotient.
  571. C
  572.         S = AINT (RB * (BX2 * D(J-1) + BX * D(J) + D(J+1)))
  573.         C(J) = S
  574.         J3 = J - 3
  575.         IN = NX - J3
  576.         IF (J .EQ. 2)  IN = IN - 1
  577. C
  578.         DO 120 I = 3, IN
  579.           IJ = I + J3
  580.           T = D(IJ) - S * B(I)
  581.           R = AINT (T * RX - 1.D0)
  582.           D(IJ) = T - R * BX
  583.           D(IJ-1) = D(IJ-1) + R
  584. 120     CONTINUE
  585.         D(J) = D(J) + BX * D(J-1)
  586. 130   CONTINUE
  587. C
  588. C   Set sign and exponent of result.
  589. C
  590.       C(NX) = AINT (RB * (BX2 * D(NX-1) + BX * D(NX)))
  591.       C(3) = BX * C(2) + C(3)
  592.       C(1) = A(1) * B(1)
  593.       C(2) = A(2) - B(2) - 1
  594.       CALL MPNORM (C)
  595. C
  596. 200   IF (IDB .GE. 6)  WRITE (LDB, 3)  (C(I),I=1,NDB)
  597. 3     FORMAT ('MPDIV O'/(8F9.0))
  598.       RETURN
  599.       END
  600. C
  601.       SUBROUTINE MPDIV1 (A, B, C, N)
  602. C
  603. C   Divides MP  A  by MP  B  to yield the SP quotient  C * 10^N,  accurate
  604. C   to about 12 decimal places.
  605. C
  606.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  607.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, RX2 = RX * RX)
  608.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  609.       DIMENSION A(1), B(1)
  610. C
  611.       NX = NW + 2
  612.       IF (IDB .GE. 7)  THEN
  613.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  614. 1       FORMAT (/'MPDIV1 I',I10/(8F9.0))
  615.         WRITE (LDB, 1)  NX, (B(I),I=1,NDB)
  616.       ENDIF
  617. C
  618. C   Check for zero dividend.
  619. C
  620.       IF (A(1) .EQ. 0.D0)  THEN
  621.         C = 0.D0
  622.         N = 0
  623.         GOTO 200
  624.       ENDIF
  625. C
  626. C   Check for zero divisor.
  627. C
  628.       IF (B(1) .EQ. 0.D0)  THEN
  629.         WRITE (LDB, 2)
  630. 2       FORMAT ('*** MPDIV1: ZERO DIVISOR ***')
  631.         CALL XABORT
  632.       ENDIF
  633. C
  634.       AA = SIGN (A(3) + A(4) * RX + A(5) * RX2, A(1))
  635.       BB = SIGN (B(3) + B(4) * RX + B(5) * RX2, B(1))
  636.       C = AA / BB
  637.       N = IB * (A(2) - B(2))
  638. C
  639. C   Adjust results so  C  is in the range  1 .LE. |C| .LT. 10.
  640. C
  641.       CA = ABS (C)
  642.       IF (CA .GE. 1.D0)  THEN
  643.         N2 = LOG10 (CA) + RX2
  644.       ELSE
  645.         N2 = LOG10 (CA) - 1.D0 - RX2
  646.       ENDIF
  647.       C = C * 10.D0 ** (-N2)
  648.       N = N + N2
  649. C
  650. 200   IF (IDB .GE. 7)  WRITE (LDB, 3)  C, N
  651. 3     FORMAT ('MPDIV1 O',F10.0,I10)
  652.       RETURN
  653.       END
  654. C
  655.       SUBROUTINE MPDIV2 (A, B, N, C)
  656. C
  657. C   Divides the MP number  A  by the SP number  B * 10^N  to obtain the MP
  658. C   quotient  C.
  659. C
  660. C   Note:  The result should be exact provided N is 0 or divisible by IB
  661. C   and B is an exact binary fraction in the range RX^2 .LT. |B| .LT. BX^2.
  662. C
  663.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  664.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, BX2 = BX * BX,
  665.      $  RX2 = RX * RX)
  666.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  667.       DIMENSION A(1), C(1)
  668. C
  669.       NX = NW + 2
  670.       IF (IDB .GE. 7)  THEN
  671.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  672. 1       FORMAT (/'MPDIV2 I',I10/(8F9.0))
  673.         WRITE (LDB, 2)  B, N
  674. 2       FORMAT ('MPDIV2 I',1PD25.15,I10)
  675.       ENDIF
  676. C
  677. C   Check for zero dividend.
  678. C
  679.       IF (A(1) .EQ. 0.D0)  THEN
  680.         DO 100 I = 1, NX
  681.           C(I) = 0.D0
  682. 100     CONTINUE
  683.         D = 0.D0
  684.         GOTO 200
  685.       ENDIF
  686. C
  687. C   Check for zero divisor.
  688. C
  689.       IF (B .EQ. 0.D0)  THEN
  690.         WRITE (LDB, 3)  B
  691. 3       FORMAT ('*** MPDIV2:  ZERO DIVISOR ***', 1PD15.6)
  692.         CALL XABORT
  693.       ENDIF
  694. C
  695. C   Adjust exponent value.
  696. C
  697.       BA = ABS (B)
  698.       BL = N + LOG10 (BA) * (1.D0 + RX2)
  699.       M = BL / IB
  700.       IF (BL .LT. 0.D0 .AND. BL .NE. M * IB)  M = M - 1
  701.       BM = BA * 10.D0 ** (IB - IB * M + N)
  702.       B1 = AINT (RX * BM)
  703.       B2 = AINT (BM - BX * B1)
  704.       RB = 1.D0 / (BX * B1 + B2)
  705.       D1 = 0.D0
  706.       D2 = A(3)
  707.       D3 = A(4)
  708. C
  709. C   Perform short division (not vectorizable at present).
  710. C
  711.       DO 110 J = 2, NX - 1
  712.         S = AINT (RB * (BX2 * D1 + BX * D2 + D3))
  713.         C(J) = S
  714.         T = D3 - S * B2
  715.         R = AINT ((T + 1) * RX - 1.D0)
  716.         D3 = T - R * BX
  717.         T = R + D2 - S * B1
  718.         R = AINT ((T + 1) * RX - 1.D0)
  719.         D2 = T - R * BX
  720.         D1 = (R + D1) * BX + D2
  721.         D2 = D3
  722.         D3 = A(J+3)
  723.         IF (J .GE. NW)  D3 = 0.D0
  724. 110   CONTINUE
  725. C
  726. C   Set sign and exponent of result.
  727. C
  728.       C(NX) = AINT (RB * (BX2 * D1 + BX * D2))
  729.       C(3) = BX * C(2) + C(3)
  730.       C(1) = A(1) * SIGN (1.D0, B)
  731.       C(2) = A(2) - M - 1
  732.       CALL MPNORM (C)
  733. C
  734. 200   IF (IDB .GE. 7)  WRITE (LDB, 4)  (C(I),I=1,NDB)
  735. 4     FORMAT ('MPDIV2 O'/(8F9.0))
  736.       RETURN
  737.       END
  738. C
  739.       SUBROUTINE MPNORM (A)
  740. C
  741. C   Converts the MP number  A  to the standard normalized form.  Many
  742. C   MP routines call MPNORM to clean up results, because occasionally these
  743. C   operations leave negative numbers or values exceeding the radix  BX  in
  744. C   result arrays.  The user need not call MPNORM directly.
  745. C
  746.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  747.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, RXZ = 2.D0 + 0.5D0 * RX)
  748.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  749.       DIMENSION A(1)
  750. C
  751.       NX = NW + 2
  752.       IF (IDB .GE. 9)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  753. 1     FORMAT (/'MPNORM I',I10/(4F18.0))
  754. C
  755.       IF (A(1) .EQ. 0.D0)  GOTO 200
  756.       A2 = A(2)
  757.       A(2) = 0.D0
  758. C
  759. C   Try vectorized fixup loop three times (this should fix 99%).
  760. C   This loop should be commented out for execution on scalar computers.
  761. C
  762. C      DO 110 K = 1, 3
  763. C        SR = 0.D0
  764. C        DO 100 I = 3, NX
  765. C          R = AINT (A(I) * RX + RXZ) - 2.D0
  766. C          A(I) = A(I) - R * BX
  767. C          A(I-1) = A(I-1) + R
  768. C          SR = SR + ABS (R)
  769. C100     CONTINUE
  770. C        IF (SR .EQ. 0.D0)  GOTO 130
  771. C110   CONTINUE
  772. C
  773. C   Still not fixed - use recursive (unvectorizable) loop.
  774. C
  775.       R = 0.D0
  776.       DO 120 I = NX, 3, -1
  777.         T = R + A(I)
  778.         R = AINT (T * RX + RXZ) - 2.D0
  779.         A(I) = T - R * BX
  780. 120   CONTINUE
  781.       A(2) = A(2) + R
  782. C
  783. C   Check for overflow.
  784. C
  785. 130   IF (A(2) .NE. 0.D0)  THEN
  786.         DO 140 I = NX, 3, -1
  787.           A(I) = A(I-1)
  788. 140     CONTINUE
  789.         A2 = A2 + 1.D0
  790.       ENDIF
  791. C
  792. C   Check if first few mantissa words are zero.
  793. C
  794.       IF (A(3) .EQ. 0.D0)  THEN
  795.         DO 150 I = 4, NX
  796.           IF (A(I) .NE. 0.D0)  GOTO 160
  797. 150     CONTINUE
  798. 160     K = I - 3
  799. CDIR$ IVDEP
  800. CDEC$ IVDEP
  801.         DO 170 I = 3, NX - K
  802.           A(I) = A(I+K)
  803. 170     CONTINUE
  804.         DO 180 I = NX - K + 1, NX
  805.           A(I) = 0.D0
  806. 180     CONTINUE
  807.         A2 = A2 - K
  808.       ENDIF
  809. C
  810.       A(2) = A2
  811. C
  812. 200   IF (IDB .GE. 9)  WRITE (LDB, 2)  (A(I),I=1,NDB)
  813. 2     FORMAT ('MPNORM O'/(8F9.0))
  814.       RETURN
  815.       END
  816. C
  817.       SUBROUTINE MPSQRT (A, B, C)
  818. C
  819. C   Computes the square root of the MP number  A  and returns the MP result
  820. C   in  B.  C  is a scratch array with at least 2 * NW + 4 cells (this is
  821. C   double the amount required for MP numbers).
  822. C
  823. C   This subroutine employs a Newton-Raphson iteration that converges to the
  824. C   square root of A:
  825. C
  826. C          x(n+1) = .5 * (x(n) + a / x(n))
  827. C
  828. C   Note that this subroutine dynamically changes the precision level, doubling
  829. C   it with each iteration.
  830. C
  831.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  832.       DIMENSION A(1), B(1), C(1)
  833.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  834.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, CL2 = 1.442695041D0)
  835. C
  836.       NX = NW + 2
  837.       IF (IDB .GE. 5)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  838. 1     FORMAT (/'MPSQRT I',I10/(8F9.0))
  839.       IDS = IDB
  840.       IF (IDB .LE. 7)  IDB = 0
  841. C
  842. C   Check for zero input.
  843. C
  844.       IF (A(1) .EQ. 0.D0)  THEN
  845.         DO 100 I = 1, NX
  846.           B(I) = 0.D0
  847. 100     CONTINUE
  848.         GOTO 200
  849.       ENDIF
  850. C
  851.       IF (A(1) .LT. 0.D0)  THEN
  852.         WRITE (LDB, 2)
  853. 2       FORMAT ('*** MPSQRT: NEGATIVE INPUT ***')
  854.         CALL XABORT
  855.       ENDIF
  856. C
  857.       NWS = NW
  858.       NX1 = NX + 1
  859.       L = CL2 * LOG (DBLE (NW)) + 1.D0 + RX
  860. C
  861. C   Compute initial approximation.
  862. C
  863.       SA = A(3) + A(4) * RX + A(5) * RX ** 2
  864.       NA = IB * A(2)
  865.       N2 = NA / 2
  866.       IF (NA .NE. 2 * N2)  THEN
  867.         IF (N2 .LT. 0)  N2 = N2 - 1
  868.         SA = SA * 10.D0
  869.       ENDIF
  870.       S2 = SQRT (SA)
  871.       CALL MPSMC (S2, N2, B)
  872.       NW = 2
  873. C
  874. C   Repeat Newton-Raphson iterations.
  875. C
  876.       DO 110 I = 1, L
  877. C
  878. C   Repeat the next-to-last step to insure maximum accuracy.
  879. C
  880.         IF (I .NE. L - 1)  NW = MIN (2 * NW, NWS)
  881.         CALL MPDIV (A, B, C, C(NX1))
  882.         CALL MPADD (B, C, C(NX1))
  883.         CALL MPMUL1 (C(NX1), 0.5D0, 0, B)
  884. 110   CONTINUE
  885. C
  886. 200   IDB = IDS
  887.       IF (IDB .GE. 5)  WRITE (LDB, 3)  (B(I),I=1,NDB)
  888. 3     FORMAT ('MPSQRT O'/(8F9.0))
  889. C
  890.       RETURN
  891.       END
  892. C
  893.       SUBROUTINE MPRAND (A)
  894. C
  895. C   Returns a pseudo-random MP number  A  between  0  and  1.
  896. C
  897.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  898.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6)
  899.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  900.       DIMENSION A(1)
  901. C
  902.       NX = NW + 2
  903.       A(1) = 1.D0
  904.       A(2) = -1.D0
  905. C
  906.       DO 100 I = 3, NX
  907.         A(I) = AINT (BX * XRAND (0))
  908. 100   CONTINUE
  909.  
  910.       IF (IDB .GE. 7)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  911. 1     FORMAT (/'MPRAND',I10/(8F9.0))
  912.       RETURN
  913.       END
  914. C
  915.       SUBROUTINE MPINFR (A, B, C)
  916. C
  917. C   Sets  B  equal to the integer part of the MP number  A
  918. C   and sets  C  equal to the fractional part of  A.
  919. C
  920.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  921.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  922.       DIMENSION  A(1), B(1), C(1)
  923. C
  924.       NX = NW + 2
  925.       IF (IDB .GE. 7)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  926. 1     FORMAT (/'MPINFR I',I10/(8F9.0))
  927. C
  928. C   Check if  A  is zero.
  929. C
  930.       IF (A(1) .EQ. 0.D0)  THEN
  931.         DO 100 I = 1, NX
  932.           B(I) = 0.D0
  933.           C(I) = 0.D0
  934. 100     CONTINUE
  935.         GOTO 200
  936.       ENDIF
  937. C
  938. C   Find word that precedes decimal point.
  939. C
  940.       N = MIN (INT (A(2)) + 3, NX)
  941.       IF (N .LT. 3)  N = 0
  942. C
  943. C   Place integer part in  B.
  944. C
  945.       DO 110 I = 1, N
  946.         B(I) = A(I)
  947. 110   CONTINUE
  948.       DO 120 I = N + 1, NX
  949.         B(I) = 0.D0
  950. 120   CONTINUE
  951. C
  952. C   Integer part is zero,  fractional part is  A.
  953. C
  954.       IF (N .EQ. 0)  THEN
  955.         C(1) = A(1)
  956.         C(2) = A(2)
  957.         GOTO 150
  958.       ENDIF
  959. C
  960. C   Check if initial cells of fractional part are zero.
  961. C
  962.       DO 130 I = N + 1, NX
  963.         IF (A(I) .NE. 0.D0)  GOTO 140
  964. 130   CONTINUE
  965. C
  966. C   The fractional part is zero.
  967. C
  968.       C(1) = 0.D0
  969.       C(2) = 0.D0
  970.       N = NX
  971.       GOTO 150
  972. C
  973. C   Store the fractional part in  C.
  974. C
  975. 140   C(1) = A(1)
  976.       C(2) = N - I
  977.       N = I - 3
  978. 150   DO 160 I = 3, NX - N
  979.         C(I) = A(I+N)
  980. 160   CONTINUE
  981.       DO 170 I = NX - N + 1, NX
  982.         C(I) = 0.D0
  983. 170   CONTINUE
  984. C
  985. 200   IF (IDB .GE. 7)  THEN
  986.         WRITE (LDB, 2)  (B(I),I=1,NDB)
  987. 2       FORMAT ('MPINFR O'/(8F9.0))
  988.         WRITE (LDB, 2)  (C(I),I=1,NDB)
  989.       ENDIF
  990.       RETURN
  991.       END
  992. C
  993.       SUBROUTINE MPFMT (A, B)
  994. C
  995. C   Formats the MP number  A  into character form of length  ND + 20  in the
  996. C   character array  B.  It is similar to the Fortran exponential format
  997. C   (E format), except that the exponent is placed first.  B may be output
  998. C   using A1 format.
  999. C
  1000.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1001.       CHARACTER*1 B
  1002.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, RX2 = RX * RX)
  1003.       DIMENSION  A(1), B(1)
  1004.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1005. C
  1006.       NX = NW + 2
  1007.       IF (IDB .GE. 8)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  1008. 1     FORMAT (/'MPFMT I',I10/(8F9.0))
  1009. C
  1010. C   Determine exact power of ten for exponent.
  1011. C
  1012.       K = 0
  1013.       IF (A(3) .NE. 0)  K = LOG10 (A(3)) + RX2
  1014.       AN = IB * A(2) + K
  1015. C
  1016. C   Place exponent first instead of at the very end as in Fortran.
  1017. C
  1018.       B(1) = '1'
  1019.       B(2) = '0'
  1020.       B(3) = ' '
  1021.       B(4) = '^'
  1022.       CALL MPENC (1, AN, B(5), 10, 1)
  1023.       B(15) = ' '
  1024.       B(16) = 'x'
  1025.       B(17) = ' '
  1026. C
  1027. C
  1028. C   Insert sign.
  1029. C
  1030.       B(18) = ' '
  1031.       IF (A(1) .LT. 0.)  B(18) = '-'
  1032. C
  1033. C   Insert first mantissa digit of number, followed by a period.
  1034. C
  1035.       AN = AINT (A(3) * 10.D0 ** (-K) + RX2)
  1036.       CALL MPENC (1, AN, B(19), 1, 0)
  1037.       B(20) = '.'
  1038. C
  1039. C   Insert the remaining digits of the first word.
  1040. C
  1041.       AN = AINT (A(3) - AN * 10.D0 ** K + RX2)
  1042.       CALL MPENC (1, AN, B(21), K, 0)
  1043. C
  1044. C   Insert the digits of the remaining words.
  1045. C
  1046.       CALL MPENC (ND / IB, A(4), B(K+21), IB, 0)
  1047. C
  1048. 200   RETURN
  1049.       END
  1050. C
  1051.       SUBROUTINE MPENC (N, A, B, LEN, IBL)
  1052. C
  1053. C   Encodes  N  floating whole numbers in the array  A  into the character
  1054. C   array  B.  the length of each number is  LEN  characters.  If  IBL  is
  1055. C   zero, then the numbers are left-filled with zeroes.  Otherwise the numbers
  1056. C   are left-filled with blanks.
  1057. C
  1058.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1059.       CHARACTER*1 B, DIG
  1060.       DIMENSION A(1), B(1), DIG(0:9)
  1061.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1062.       DATA DIG/'0', '1', '2', '3', '4', '5', '6', '7', '8', '9'/
  1063. C
  1064.       IF (LEN .EQ. 0)  GOTO 200
  1065.       IX = 0
  1066.       DO 130 J = 1, N
  1067.         AJ = A(J)
  1068.         XJ = ABS (AJ)
  1069.         IF (AJ .GE. 0.)  THEN
  1070.           L1 = 1
  1071.         ELSE
  1072.           L1 = 2
  1073.           B(IX+1) = ' '
  1074.         ENDIF
  1075. C
  1076. C   Extract digits one at a time, beginning with units digit.
  1077. C
  1078.         DO 100 I = LEN, L1, -1
  1079.           XI = AINT (0.1D0 * XJ + 0.01D0)
  1080.           K = XJ - 10. * XI
  1081.           XJ = XI
  1082.           B(IX+I) = DIG(K)
  1083. 100     CONTINUE
  1084. C
  1085. C   Check if input is too large for field.
  1086. C
  1087.         IF (XI .NE. 0.)  THEN
  1088.           WRITE (LDB, 1)  J, AJ, LEN
  1089. 1         FORMAT ('*** MPENC: INPUT TOO LARGE ***', I5, F15.1, I5)
  1090.           CALL XABORT
  1091.         ENDIF
  1092. C
  1093. C   Left-fill with blanks if requested.
  1094. C
  1095.         IF (IBL .EQ. 1)  THEN
  1096.           DO 110 I = L1, LEN - 1
  1097.             IF (B(IX+I) .NE. '0')  GOTO 120
  1098.             B(IX+I) = ' '
  1099. 110       CONTINUE
  1100.           I = LEN
  1101. C
  1102. C   Insert sign if required to the left of first nonzero digit.
  1103. C
  1104. 120       IF (AJ .LT. 0)  B(IX+I-1) = '-'
  1105.         ENDIF
  1106.         IX = IX + LEN
  1107. 130   CONTINUE
  1108. C
  1109. 200   RETURN
  1110.       END
  1111. C
  1112.       SUBROUTINE MPMULX (A, B, C, D)
  1113. C
  1114. C   Multiplies MP numbers A and B to yield the DMP product C.  The array D is
  1115. C   a scratch array with at least 12 * NW + 6 cells (twelve times the amount
  1116. C   allocated for MP numbers will suffice).  The value of NW must be equal to
  1117. C   2^(MW - 2).  Subroutine MPCFFT must have previously been called to
  1118. C   initialize the array U in common block MPCOM2.
  1119. C
  1120. C   This version checks the input numbers A and B to determine the locations of
  1121. C   the last nonzero words.  If this location is low, so that the actual
  1122. C   precision of the number is low, then MPMUL is called instead, and these
  1123. C   locations and passed to it.
  1124. C
  1125. C   This subroutine uses a special technique involving a discrete convolution
  1126. C   and a fast Fourier transform.  The work factor for this algorithm is
  1127. C   proportional to NW * log (NW), which is much faster than the NW^2 algorithm
  1128. C   used in MPMUL for large values of NW.  The parameter MI is the crossover
  1129. C   point for MW to apply the advanced algorithm.  On Cray systems MI should
  1130. C   be set to 8.  On other systems another value may be slightly better.  On
  1131. C   scalar computers, MI probably should be set to 5.  Change MI here and also
  1132. C   in MPDIVX and MPSQRX.
  1133. C
  1134.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1135.       DIMENSION  A(1), B(1), C(1), D(1)
  1136.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1137.       PARAMETER (BX = 1D6, RX = 1D-6, IB = 6, MI = 5, BB = 1D3,
  1138.      $   RB = 1D-3)
  1139. C
  1140.       NX = NW + 2
  1141.       IF (IDB .GE. 6)  THEN
  1142.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  1143. 1       FORMAT (/'MPMULX I',I10/(8F9.0))
  1144.         WRITE (LDB, 1)  NX, (B(I),I=1,NDB)
  1145.       ENDIF
  1146. C
  1147. C   Check if precision level is too low to justify advanced algorithm.  The
  1148. C   following code searches backwards from the end of the input arrays to find
  1149. C   the indices of the last nonzero words.
  1150. C
  1151.       DO 100 I = NX, 3, -1
  1152.         IF (A(I) .NE. 0.D0)  GOTO 110
  1153. 100   CONTINUE
  1154. 110   NA = I
  1155. C
  1156.       DO 120 I = NX, 3, -1
  1157.         IF (B(I) .NE. 0.D0)  GOTO 130
  1158. 120   CONTINUE
  1159. 130   NB = I
  1160. C
  1161.       IF (MW .LE. MI .OR. NA .LT. 64 .OR. NB .LT. 64)  THEN
  1162.         CALL MPMUL (A, B, C, NA, NB)
  1163.         GOTO 200
  1164.       ENDIF
  1165. C
  1166. C   Check for mismatch between MW and NW.
  1167. C
  1168.       N1 = 2 ** (MW - 2)
  1169.       IF (N1 .NE. NW)  THEN
  1170.         WRITE (LDB, 2)  MW, NW
  1171. 2       FORMAT ('*** MPMULX:  MISMATCHED VALUES OF MW AND NW--', 2I10)
  1172.         CALL XABORT
  1173.       ENDIF
  1174. C
  1175.       N2 = 2 * N1
  1176.       N21 = N2 + 1
  1177.       N4 = 4 * N1
  1178.       N42 = N4 + 2
  1179.       N6 = 6 * N1
  1180.       N63 = N6 + 3
  1181.       N8 = 8 * N1
  1182.       N84 = N8 + 4
  1183. C
  1184. C   Check for zero inputs.
  1185. C
  1186.       IF (A(1) * B(1) .EQ. 0.D0)  THEN
  1187.         DO 140 I = 1, N2 + 2
  1188.           C(I) = 0.D0
  1189. 140     CONTINUE
  1190.         GOTO 200
  1191.       ENDIF
  1192. C
  1193. C   Split input data in A and B into words with half as many bits and place in
  1194. C   (D(I), I=1,N4) and (D(I), I=N4+3,N8+2).  The second half of the two vectors
  1195. C   must be filled with zeroes for the FFT multiply technique to work properly.
  1196. C
  1197. CDIR$ IVDEP
  1198. CDEC$ IVDEP
  1199.       DO 150 I = 1, N1
  1200.         T = A(I+2)
  1201.         A1 = AINT (RB * T)
  1202.         A2 = T - BB * A1
  1203.         D(2*I-1) = A1
  1204.         D(2*I) = A2
  1205.         T = B(I+2)
  1206.         B1 = AINT (RB * T)
  1207.         B2 = T - BB * B1
  1208.         D(2*I-1+N42) = B1
  1209.         D(2*I+N42) = B2
  1210. 150   CONTINUE
  1211. C
  1212. CDIR$ IVDEP
  1213. CDEC$ IVDEP
  1214.       DO 160 I = N2 + 1, N42
  1215.         D(I) = 0.D0
  1216.         D(I+N42) = 0.D0
  1217. 160   CONTINUE
  1218. C
  1219. C   Perform forward real-to-complex FFTs on the two vectors in D.
  1220. C
  1221.       CALL MPRCFT (1, D, D(N84+1))
  1222.       CALL MPRCFT (1, D(N42+1), D(N84+1))
  1223. C
  1224. C   Multiply resulting complex arrays.
  1225. C
  1226. CDIR$ IVDEP
  1227. CDEC$ IVDEP
  1228.       DO 170 I = 1, N21
  1229.         D11 = D(I)
  1230.         D12 = D(I+N21)
  1231.         D21 = D(I+N42)
  1232.         D22 = D(I+N63)
  1233.         D(I) = D11 * D21 - D12 * D22
  1234.         D(I+N21) = D11 * D22 + D12 * D21
  1235. 170   CONTINUE
  1236. C        
  1237. C   Perform reverse FFT.
  1238. C
  1239.       CALL MPCRFT (-1, D, D(N84+1))
  1240. C
  1241. C   Round, recombine into regular words, divide by 2^MW, and release carrys.
  1242. C
  1243.       AN = 1.D0 / N8
  1244.       T1 = AN * D(1)
  1245.       A1 = AINT (T1 + 0.5D0)
  1246.       D(1) = ABS (A1 - T1)
  1247.       C(2) = A1
  1248. C
  1249. CDIR$ IVDEP
  1250. CDEC$ IVDEP
  1251.       DO 180 I = 1, N2 - 1
  1252.         T1 = AN * D(2*I)
  1253.         T2 = AN * D(2*I+1)
  1254.         A1 = AINT (T1 + 0.5D0)
  1255.         A2 = AINT (T2 + 0.5D0)
  1256.         D(2*I) = ABS (A1 - T1)
  1257.         D(2*I+1) = ABS (A2 - T2)
  1258.         C1 = AINT (RX * A1)
  1259.         C2 = A1 - BX * C1
  1260.         C3 = AINT (RX * A2)
  1261.         C4 = A2 - BX * C3
  1262.         C(I+2) = BB * C2 + C4
  1263.         C(I+1) = C(I+1) + BB * C1 + C3
  1264. 180   CONTINUE
  1265. C
  1266.       C(N2+2) = 0.D0
  1267. C
  1268. C   Find the largest rounding error in the loop above.  This code may be
  1269. C   commented out unless very high precision (millions of digits) is used.
  1270. C
  1271. C      T = 0.D0
  1272. C
  1273. C      DO 190 I = 1, N4 - 1
  1274. C        T = MAX (T, D(I))
  1275. C190   CONTINUE
  1276. C
  1277. C      IF (T .GT. 0.25D0)  THEN
  1278. C        WRITE (LDB, 3)  T
  1279. C3       FORMAT ('*** MPMULX:  EXCESSIVE ROUNDING ERROR',F8.4)
  1280. C        CALL XABORT
  1281. C      ENDIF
  1282. C
  1283. C   Adjust sign and exponent.
  1284. C
  1285.       C(1) = A(1) * B(1)
  1286.       C(3) = C(3) + BX * C(2)
  1287.       C(2) = A(2) + B(2)
  1288.       NW = N2
  1289.       CALL MPNORM (C)
  1290.       NW = N1
  1291. C
  1292. 200   IF (IDB .GE. 6)  WRITE (LDB, 4)  (C(I),I=1,NDB)
  1293. 4     FORMAT ('MPMULX O'/(8F9.0))
  1294.       RETURN
  1295.       END
  1296. C
  1297.       SUBROUTINE MPRCFT (IS, X, Y)
  1298. C
  1299. C   Performs an N-point real-to-complex FFT, where N = 2^MW.  X is both the
  1300. C   input and the output data array, and Y is a scratch array.  N real values
  1301. C   are input in X, and N/2 + 1 complex values are output in X, with real and
  1302. C   imaginary parts separated by N/2 + 1 locations.  Before calling MPRCFT,
  1303. C   the U array must be initialized by calling MPCFFT with IS = 0.  A call to
  1304. C   MPRCFT with IS = 1 (or -1) indicates a call to perform a real-to-complex
  1305. C   FFT with positive (or negative) exponentials.  The arrays X and Y must be
  1306. C   dimensioned with N/2 + 1 complex or N + 2 real cells.  U must be
  1307. C   dimensioned the same as in MPCFFT.  The output values from MPRCFT are
  1308. C   twice as large as the results of a complex-to-complex transform on real
  1309. C   data.  MW must be at least three.
  1310. C
  1311. C   This subroutine employs a technique that converts a real-to-complex FFT
  1312. C   to a complex-to-complex FFT.  See Brigham, "The Fast Fourier Transform",
  1313. C   p. 169, and Hockney & Jesshope, "Parallel Computers", p. 303, although
  1314. C   neither reference is complete and error-free.
  1315. C
  1316. C   David H. Bailey    May 9, 1988
  1317. C
  1318.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1319.       DIMENSION X(1), Y(1)
  1320.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1321.       COMMON /MPCOM2/ U(1)
  1322. C
  1323. C   Set initial parameters.
  1324. C
  1325.       K = U(1)
  1326.       MX = MOD (K, 64)
  1327.       NU = K / 64
  1328.       N = 2 ** MW
  1329.       N2 = N / 2
  1330.       N21 = N2 + 1
  1331.       N4 = N / 4
  1332.       KU = N / 2
  1333.       KN = KU + NU
  1334. C
  1335. C   Check if input parameters are invalid.
  1336. C
  1337.       IF ((IS .NE. 1 .AND. IS .NE. -1) .OR. MW .LT. 3 .OR. MW .GT. MX)
  1338.      $  THEN
  1339.         WRITE (LDB, 1)  IS, MW, MX
  1340. 1       FORMAT ('*** MPRCFT ERROR: EITHER U HAS NOT BEEN INITIALIZED'/
  1341.      $    'OR ELSE ONE OF THE INPUT PARAMETERS IS INVALID --', 3I5)
  1342.         CALL XABORT
  1343.       ENDIF
  1344. C
  1345. C   Copy X to Y such that Y(k) = X(2k-1) + i X(2k).
  1346. C
  1347. CDIR$ IVDEP
  1348. CDEC$ IVDEP
  1349.       DO 100 K = 1, N2
  1350.         Y(K) = X(2*K-1)
  1351.         Y(K+N2) = X(2*K)
  1352. 100   CONTINUE
  1353. C
  1354. C   Perform a normal N/2-point FFT on Y.
  1355. C
  1356.       CALL MPCFFT (IS, Y, X)
  1357. C
  1358. C   Reconstruct the FFT of X.
  1359. C
  1360.       X(1) = 2.D0 * (Y(1) + Y(N21))
  1361.       X(N21+1) = 0.D0
  1362.       X(N4+1) = 2.D0 * Y(N4+1)
  1363.       X(N4+1+N21) = 2.D0 * IS * Y(N4+N2+1)
  1364.       X(N21) = 2.D0 * (Y(1) - Y(N21))
  1365.       X(N+2) = 0.D0
  1366. C
  1367. CDIR$ IVDEP
  1368. CDEC$ IVDEP
  1369.       DO 110 K = 2, N4
  1370.         Y11 = Y(K)
  1371.         Y12 = Y(K+N2)
  1372.         Y21 = Y(N2+2-K)
  1373.         Y22 = Y(N+2-K)
  1374.         A1 = Y11 + Y21
  1375.         A2 = Y11 - Y21
  1376.         B1 = Y12 + Y22
  1377.         B2 = Y12 - Y22
  1378.         U1 = U(K+KU)
  1379.         U2 = IS * U(K+KN)
  1380.         T1 = U1 * B1 + U2 * A2
  1381.         T2 = - U1 * A2 + U2 * B1
  1382.         X(K) = A1 + T1
  1383.         X(K+N21) = B2 + T2
  1384.         X(N2+2-K) = A1 - T1
  1385.         X(N+3-K) = -B2 + T2
  1386. 110   CONTINUE
  1387. C
  1388.       RETURN
  1389.       END
  1390. C
  1391.       SUBROUTINE MPCRFT (IS, X, Y)
  1392. C
  1393. C   Performs an N-point complex-to-real FFT, where N = 2^MW.  X is both the
  1394. C   input and the output data array, and Y is a scratch array.  N/2 + 1
  1395. C   complex values are input in X, with real and imaginary parts separated by
  1396. C   N/2 + 1 locations, and N real values are output.  Before calling MPCRFT,
  1397. C   the U array must be initialized by calling MPCFFT with IS = 0.  A call to
  1398. C   MPCRFT with IS = 1 (or -1) indicates a call to perform a complex-to-real
  1399. C   FFT with positive (or negative) exponentials.  The arrays X and Y must be
  1400. C   dimensioned with N/2 + 1 complex or N + 2 real cells.  U must be
  1401. C   dimensioned the same as in MPCFFT.  MW must be at least three.
  1402. C
  1403. C   This subroutine employs a technique that converts a complex-to-real FFT
  1404. C   to a complex-to-complex FFT.  See Brigham, "The Fast Fourier Transform",
  1405. C   p. 169, and Hockney & Jesshope, "Parallel Computers", p. 303, although
  1406. C   neither reference is complete and error-free.
  1407. C
  1408. C   David H. Bailey    May 9, 1988
  1409. C
  1410.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1411.       DIMENSION X(1), Y(1)
  1412.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1413.       COMMON /MPCOM2/ U(1)
  1414. C
  1415. C   Set initial parameters.
  1416. C
  1417.       K = U(1)
  1418.       MX = MOD (K, 64)
  1419.       NU = K / 64
  1420.       N = 2 ** MW
  1421.       N2 = N / 2
  1422.       N21 = N2 + 1
  1423.       N4 = N / 4
  1424.       KU = N / 2
  1425.       KN = KU + NU
  1426. C
  1427. C   Check if input parameters are invalid.
  1428. C
  1429.       IF ((IS .NE. 1 .AND. IS .NE. -1) .OR. MW .LT. 3 .OR. MW .GT. MX)
  1430.      $  THEN
  1431.         WRITE (LDB, 1)  IS, MW, MX
  1432. 1       FORMAT ('*** MPCRFT ERROR: EITHER U HAS NOT BEEN INITIALIZED'/
  1433.      $    'OR ELSE ONE OF THE INPUT PARAMETERS IS INVALID --', 3I5)
  1434.         CALL XABORT
  1435.       ENDIF
  1436. C
  1437. C   Construct the input to MPCFFT.
  1438. C
  1439.       Y(1) = 0.5D0 * (X(1) + X(N21))
  1440.       Y(N2+1) = 0.5D0 * (X(1) - X(N21))
  1441.       Y(N4+1) = X(N4+1)
  1442.       Y(N4+N2+1) = -IS * X(N4+N2+2)
  1443. C      
  1444. CDIR$ IVDEP
  1445. CDEC$ IVDEP
  1446.       DO 100 K = 2, N4
  1447.         X11 = X(K)
  1448.         X12 = X(K+N21)
  1449.         X21 = X(N2+2-K)
  1450.         X22 = X(N+3-K)
  1451.         A1 = X11 + X21
  1452.         A2 = X11 - X21
  1453.         B1 = X12 + X22
  1454.         B2 = X12 - X22
  1455.         U1 = U(K+KU)
  1456.         U2 = IS * U(K+KN)
  1457.         T1 = U1 * B1 + U2 * A2
  1458.         T2 = U1 * A2 - U2 * B1
  1459.         Y(K) = 0.5D0 * (A1 - T1)
  1460.         Y(K+N2) = 0.5D0 * (B2 + T2)
  1461.         Y(N2+2-K) = 0.5D0 * (A1 + T1)
  1462.         Y(N+2-K) = 0.5D0 * (-B2 + T2)
  1463. 100   CONTINUE
  1464. C
  1465. C   Perform a normal N/2-point FFT on Y.
  1466. C
  1467.       CALL MPCFFT (IS, Y, X)
  1468. C
  1469. C   Copy Y to X such that Y(k) = X(2k-1) + i X(2k).
  1470. C
  1471. CDIR$ IVDEP
  1472. CDEC$ IVDEP
  1473.       DO 110 K = 1, N2
  1474.         X(2*K-1) = Y(K)
  1475.         X(2*K) = Y(K+N2)
  1476. 110   CONTINUE
  1477. C
  1478.       RETURN
  1479.       END
  1480. C
  1481.       SUBROUTINE MPCFFT (IS, X, Y)
  1482. C
  1483. C   Computes the 2^M-point complex-to-complex FFT of X using an algorithm due
  1484. C   to Swarztrauber, coupled with some fast methods for performing power-of-
  1485. C   two matrix transpositions (see article by DHB in Intl. J. of Supercomputer
  1486. C   Applications, Spring 1988, p. 82 - 87). This is the radix 2 version.
  1487. C   X is both the input and the output array, while Y is a scratch array.
  1488. C   Both X and Y must be dimensioned with 2 * N real cells, where N = 2^M.
  1489. C   The data in X are assumed to have real and imaginary parts separated
  1490. C   by N cells.  M in the above is actually MW - 1, since in this application
  1491. C   MPCFFT is called by MPRCFT and MPCRFT on half-size complex arrays.
  1492. C
  1493. C   Before calling MPRCFT or MPCRFT to perform an FFT, the array U must be
  1494. C   initialized by calling MPCFFT with IS set to 0 and MW set to MX, where MX
  1495. C   is the maximum value of MW for any subsequent call.  U must be dimensioned
  1496. C   with at least 2 * NX real cells, where NX = 2^MX.  MW must be at least two.
  1497. C
  1498. C   This routine and the next three have been adapted for use in the
  1499. C   multiprecision package -- parameters MW and U are passed in common.
  1500. C
  1501. C   David H. Bailey     May 9, 1988
  1502. C
  1503.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1504.       PARAMETER (PI = 3.141592653589793238D0)
  1505.       DIMENSION X(1), Y(1)
  1506.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1507.       COMMON /MPCOM2/ U(1)
  1508. C
  1509.       IF (IS .EQ. 0)  THEN
  1510. C
  1511. C   Initialize the U array with sines and cosines in a manner that permits
  1512. C   stride one access at each FFT iteration.
  1513. C
  1514.         M = MW
  1515.         N = 2 ** M
  1516.         NU = N
  1517.         U(1) = 64 * N + M
  1518.         KU = 2
  1519.         KN = KU + NU
  1520.         LN = 1
  1521. C
  1522.         DO 110 J = 1, M
  1523.           T = PI / LN
  1524. CDIR$ IVDEP
  1525. CDEC$ IVDEP
  1526.           DO 100 I = 0, LN - 1
  1527.             TI = I * T
  1528.             U(I+KU) = COS (TI)
  1529.             U(I+KN) = SIN (TI)
  1530. 100       CONTINUE
  1531. C
  1532.           KU = KU + LN
  1533.           KN = KU + NU
  1534.           LN = 2 * LN
  1535. 110     CONTINUE
  1536. C
  1537.         RETURN
  1538.       ENDIF
  1539. C
  1540. C   A normal call to MPCFFT starts here.  M1 is the number of the first variant
  1541. C   radix-2 Stockham iterations to be performed.  The second variant is faster
  1542. C   on most machines after the first few iterations, since in the second
  1543. C   variant it is not necessary to access roots of unity in the inner DO loop.
  1544. C   Thus it is most efficient to limit M1 to some value.  For the Cray-2,
  1545. C   it is most efficient to limit M1 to 6.  On other systems this limit will
  1546. C   be different.  For scalar systems, M1 probably should be limited to 2.
  1547. C
  1548.       M = MW - 1
  1549.       N = 2 ** M
  1550.       M1 = MIN0 (M / 2, 2)
  1551.       M2 = M - M1
  1552.       N2 = 2 ** M1
  1553.       N1 = 2 ** M2
  1554. C
  1555. C   Perform one variant of the Stockham FFT.
  1556. C
  1557.       DO 120 L = 1, M1, 2
  1558.         CALL MPFFT1 (IS, L, X, Y)
  1559.         IF (L .EQ. M1) GOTO 140
  1560.         CALL MPFFT1 (IS, L + 1, Y, X)
  1561. 120   CONTINUE
  1562. C
  1563. C   Perform a transposition of X treated as a N2 x N1 x 2 matrix.
  1564. C
  1565.       CALL MPTRAN (N1, N2, X, Y)
  1566. C
  1567. C   Perform second variant of the Stockham FFT from Y to X and X to Y.
  1568. C
  1569.       DO 130 L = M1 + 1, M, 2
  1570.         CALL MPFFT2 (IS, L, Y, X)
  1571.         IF (L .EQ. M) GOTO 180
  1572.         CALL MPFFT2 (IS, L + 1, X, Y)
  1573. 130   CONTINUE
  1574. C
  1575.       GOTO 160
  1576. C
  1577. C   Perform a transposition of Y treated as a N2 x N1 x 2 matrix.
  1578. C
  1579. 140   CALL MPTRAN (N1, N2, Y, X)
  1580. C
  1581. C   Perform second variant of the Stockham FFT from X to Y and Y to X.
  1582. C
  1583.       DO 150 L = M1 + 1, M, 2
  1584.         CALL MPFFT2 (IS, L, X, Y)
  1585.         IF (L .EQ. M) GOTO 160
  1586.         CALL MPFFT2 (IS, L + 1, Y, X)
  1587. 150   CONTINUE
  1588. C
  1589.       GOTO 180
  1590. C
  1591. C   Copy Y to X.
  1592. C
  1593. 160   DO 170 I = 1, 2 * N
  1594.         X(I) = Y(I)
  1595. 170   CONTINUE
  1596. C
  1597. 180   CONTINUE
  1598.       RETURN
  1599.       END
  1600. C
  1601.       SUBROUTINE MPFFT1 (IS, L, X, Y)
  1602. C
  1603. C   Performs the L-th iteration of the first variant of the Stockham FFT.
  1604. C
  1605.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1606.       DIMENSION X(1), Y(1)
  1607.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1608.       COMMON /MPCOM2/ U(1)
  1609. C
  1610. C   Set initial parameters.
  1611. C
  1612.       M = MW - 1
  1613.       N = 2 ** M
  1614.       K = U(1)
  1615.       MX = MOD (K, 64)
  1616.       NU = K / 64
  1617.       N1 = N / 2
  1618.       LK = 2 ** (L - 1)
  1619.       LI = 2 ** (M - L)
  1620.       LJ = 2 * LI
  1621.       KU = LI + 1
  1622.       KN = KU + NU
  1623. C
  1624.       DO 100 K = 0, LK - 1
  1625.         I11 = K * LJ + 1
  1626.         I12 = I11 + LI
  1627.         I21 = K * LI + 1
  1628.         I22 = I21 + N1
  1629. C
  1630. CDIR$ IVDEP
  1631. CDEC$ IVDEP
  1632.         DO 100 I = 0, LI - 1
  1633.           U1 = U(KU+I)
  1634.           U2 = IS * U(KN+I)
  1635.           X11 = X(I11+I)
  1636.           X12 = X(I11+I+N)
  1637.           X21 = X(I12+I)
  1638.           X22 = X(I12+I+N)
  1639.           T1 = X11 - X21
  1640.           T2 = X12 - X22
  1641.           Y(I21+I) = X11 + X21
  1642.           Y(I21+I+N) = X12 + X22
  1643.           Y(I22+I) = U1 * T1 - U2 * T2
  1644.           Y(I22+I+N) = U1 * T2 + U2 * T1
  1645. 100   CONTINUE
  1646. C
  1647.       RETURN
  1648.       END
  1649. C
  1650.       SUBROUTINE MPFFT2 (IS, L, X, Y)
  1651. C
  1652. C   Performs the L-th iteration of the second variant of the Stockham FFT.
  1653. C
  1654.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1655.       DIMENSION X(1), Y(1)
  1656.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1657.       COMMON /MPCOM2/ U(1)
  1658. C
  1659. C   Set initial parameters.
  1660. C
  1661.       M = MW - 1
  1662.       N = 2 ** M
  1663.       K = U(1)
  1664.       MX = MOD (K, 64)
  1665.       NU = K / 64
  1666.       N1 = N / 2
  1667.       LK = 2 ** (L - 1)
  1668.       LI = 2 ** (M - L)
  1669.       LJ = 2 * LK
  1670.       KU = LI + 1
  1671. C
  1672.       DO 100 I = 0, LI - 1
  1673.         I11 = I * LK + 1
  1674.         I12 = I11 + N1
  1675.         I21 = I * LJ + 1
  1676.         I22 = I21 + LK
  1677.         U1 = U(KU+I)
  1678.         U2 = IS * U(KU+I+NU)
  1679. C
  1680. CDIR$ IVDEP
  1681. CDEC$ IVDEP
  1682.         DO 100 K = 0, LK - 1
  1683.           X11 = X(I11+K)
  1684.           X12 = X(I11+K+N)
  1685.           X21 = X(I12+K)
  1686.           X22 = X(I12+K+N)
  1687.           T1 = X11 - X21
  1688.           T2 = X12 - X22
  1689.           Y(I21+K) = X11 + X21
  1690.           Y(I21+K+N) = X12 + X22
  1691.           Y(I22+K) = U1 * T1 - U2 * T2
  1692.           Y(I22+K+N) = U1 * T2 + U2 * T1
  1693. 100   CONTINUE
  1694. C
  1695.       RETURN
  1696.       END
  1697. C
  1698.       SUBROUTINE MPTRAN (N1, N2, X, Y)
  1699. C
  1700. C   Performs a transpose of the vector X, returning the result in Y.  X is
  1701. C   treated as a N1 x N2 complex matrix, and Y is treated as a N2 x N1 complex
  1702. C   matrix.  The complex data is assumed stored with real and imaginary parts
  1703. C   separated by N1 x N2 locations.  If this routine is to be used for an
  1704. C   application involving only real data, then the second line of all inner DO
  1705. C   loops may be deleted.
  1706. C
  1707.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1708.       DIMENSION X(1), Y(1)
  1709. C
  1710.       N = N1 * N2
  1711. C
  1712. C   Perform one of three techniques, depending on N.  The best strategy varies
  1713. C   with the computer system.  The following strategy is best for the Cray-2.
  1714. C   For scalar computers, the outer IF block should be commented out.
  1715. C
  1716. C      IF (N1 .LT. 32 .OR. N2 .LT. 32) THEN
  1717.         IF (N1 .GE. N2) THEN
  1718.           GOTO 100
  1719.         ELSE
  1720.           GOTO 120
  1721.         ENDIF
  1722. C      ELSE
  1723. C        GOTO 140
  1724. C      ENDIF
  1725. C
  1726. C   Scheme 1:  Perform a simple transpose in the usual way.  This is usually
  1727. C   the best on vector computers if N2 is odd, or if both N1 and N2 are small,
  1728. C   and N1 is larger than N2.
  1729. C
  1730. 100   DO 110 J = 0, N2 - 1
  1731. CDIR$ IVDEP
  1732. CDEC$ IVDEP
  1733.         DO 110 I = 0, N1 - 1
  1734.             Y(I*N2+J+1) = X(I+J*N1+1)
  1735.             Y(I*N2+J+1+N) = X(I+J*N1+1+N)
  1736. 110   CONTINUE
  1737. C
  1738.       GOTO 200
  1739. C
  1740. C   Scheme 2:  Perform a simple transpose with the loops reversed.  This is
  1741. C   usually the best on vector computers if N1 is odd, or if both N1 and N2 are
  1742. C   small, and N2 is larger than N1.
  1743. C
  1744. 120   DO 130 I = 0, N1 - 1
  1745. CDIR$ IVDEP
  1746. CDEC$ IVDEP
  1747.         DO 130 J = 0, N2 - 1
  1748.             Y(J+I*N2+1) = X(J*N1+I+1)
  1749.             Y(J+I*N2+1+N) = X(J*N1+I+1+N)
  1750. 130   CONTINUE
  1751. C
  1752.       GOTO 200
  1753. C
  1754. C   Scheme 3:  Perform the transpose along diagonals to insure odd strides.
  1755. C   This works well on moderate vector, variable stride computers, when both
  1756. C   N1 and N2 are divisible by reasonably large powers of two (32 or larger on
  1757. C   Cray computers).
  1758. C
  1759. 140   N11 = N1 + 1
  1760.       N21 = N2 + 1
  1761.       IF (N1 .GE. N2) THEN
  1762.         K1 = N1
  1763.         K2 = N2
  1764.         I11 = N1
  1765.         I12 = 1
  1766.         I21 = 1
  1767.         I22 = N2
  1768.       ELSE
  1769.         K1 = N2
  1770.         K2 = N1
  1771.         I11 = 1
  1772.         I12 = N2
  1773.         I21 = N1
  1774.         I22 = 1
  1775.       ENDIF
  1776. C
  1777.       DO 150 J = 0, K2 - 1
  1778.         J1 = J * I11 + 1
  1779.         J2 = J * I12 + 1
  1780. CDIR$ IVDEP
  1781. CDEC$ IVDEP
  1782.         DO 150 I = 0, K2 - 1 - J
  1783.           Y(N21*I+J2) = X(N11*I+J1)
  1784.           Y(N21*I+J2+N) = X(N11*I+J1+N)
  1785. 150   CONTINUE
  1786. C
  1787.       DO 160 J = 1, K1 - K2 - 1
  1788.         J1 = J * I21 + 1
  1789.         J2 = J * I22 + 1
  1790. CDIR$ IVDEP
  1791. CDEC$ IVDEP
  1792.         DO 160 I = 0, K2 - 1
  1793.           Y(N21*I+J2) = X(N11*I+J1)
  1794.           Y(N21*I+J2+N) = X(N11*I+J1+N)
  1795. 160   CONTINUE
  1796. C
  1797.       DO 170 J = K1 - K2, K1 - 1
  1798.         J1 = J * I21 + 1
  1799.         J2 = J * I22 + 1
  1800. CDIR$ IVDEP
  1801. CDEC$ IVDEP
  1802.         DO 170 I = 0, K1 - 1 - J
  1803.           Y(N21*I+J2) = X(N11*I+J1)
  1804.           Y(N21*I+J2+N) = X(N11*I+J1+N)
  1805. 170   CONTINUE
  1806. C
  1807.       GOTO 200
  1808. C
  1809. C   Scheme 4:  Explicitly compute the transposition permutation (with an odd
  1810. C   stride offset) and utilize gather-scatter hardware.  This scheme works
  1811. C   well on long vector machines when N1 and N2 are both powers of two.
  1812. C
  1813. 180   N11 = N1 + 1
  1814. C
  1815.       DO 190 I = 0, N - 1
  1816.         J = N11 * I
  1817.         K = J - N * (J / N)
  1818.         J = K / N1
  1819.         L = J + N2 * (K - N1 * J)
  1820.         Y(L+1) = X(K+1)
  1821.         Y(L+1+N) = X(K+1+N)
  1822. 190   CONTINUE
  1823. C
  1824. 200   RETURN
  1825.       END
  1826. C
  1827.       SUBROUTINE MPDIVX (A, B, C, D)
  1828. C
  1829. C   Divides the MP dividend A by the MP divisor B to yield the MP quotient C.
  1830. C   The array D is a scratch array with at least 15 * NW + 12 cells (fifteen
  1831. C   times the amount allocated for MP numbers will suffice).  Since this
  1832. C   subroutine calls MPMULX, NW must be equal to 2^(MW - 2), and MPCFFT
  1833. C   must have been previously called to initialize the array U in common
  1834. C   MPCOM2.
  1835. C
  1836. C   This subroutine calls MPDIV is the value of MW is less than 8.
  1837. C
  1838. C   This subroutine employs the following Newton-Raphson iteration to compute
  1839. C   the reciprocal of B:
  1840. C
  1841. C          x(n+1) = x(n) * (2 - b * x(n))
  1842. C
  1843. C   Multiplying the result by A yields the quotient.  Note that this subroutine
  1844. C   dynamically changes the precision level, doubling with each iteration.
  1845. C
  1846.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1847.       DIMENSION  A(1), B(1), C(1), D(1)
  1848.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1849.       PARAMETER (BX = 1D6, RX = 1D-6, MI = 5)
  1850. C
  1851.       NX = NW + 2
  1852.       IF (IDB .GE. 6)  THEN
  1853.         WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  1854. 1       FORMAT (/'MPDIVX I',I10/(8F9.0))
  1855.         WRITE (LDB, 1)  NX, (B(I),I=1,NDB)
  1856.       ENDIF
  1857.       IDS = IDB
  1858.       IF (IDB .LE. 7)  IDB = 0
  1859. C
  1860. C   Check if precision level is too low to justify the advanced algorithm.
  1861. C
  1862.       IF (MW .LE. MI)  THEN
  1863.         CALL MPDIV (A, B, C, D)
  1864.         GOTO 200
  1865.       ENDIF
  1866. C
  1867. C   Check for zero dividend.
  1868. C
  1869.       IF (A(1) .EQ. 0.D0)  THEN
  1870.         DO 100 I = 1, NX
  1871.           D(I) = 0.D0
  1872. 100     CONTINUE
  1873.         GOTO 200
  1874.       ENDIF
  1875. C
  1876.       NWS = NW
  1877.       NX1 = NX + 1
  1878.       NX2 = 2 * NX + 1
  1879.       NX3 = 3 * NX + 1
  1880.       MWS = MW
  1881. C
  1882. C   Compute initial approximation to two extra words precision.
  1883. C
  1884.       CALL MPSMC (0.D0, 0, C)
  1885.       MW = MI
  1886.       NW = 2 ** (MI - 2) + 2
  1887.       CALL MPSMC (1.D0, 0, D)
  1888.       CALL MPDIV (D, B, C, D(NX1))
  1889.       NW = NW - 2
  1890. C
  1891. C   Repeat Newton-Raphson iterations.
  1892. C
  1893.       DO 110 I = MI, MWS
  1894. C
  1895. C   Repeat the next-to-last step to insure maximum accuracy.
  1896. C
  1897.         IF (I .NE. MWS - 1)  THEN
  1898.           MW = MW + 1
  1899.           NW = 2 * NW
  1900.         ENDIF
  1901. C
  1902.         CALL MPMULX (B, C, D(NX1), D(NX3))
  1903.         CALL MPSMC (2.D0, 0, D(NX2))
  1904.         CALL MPSUB (D(NX2), D(NX1), D)
  1905.         CALL MPMULX (C, D, D(NX1), D(NX3))
  1906.         CALL MPEQ (D(NX1), C)
  1907. 110   CONTINUE
  1908. C
  1909. C   Multiply by A to complete the division calculation.
  1910. C
  1911.       CALL MPMULX (A, C, D, D(NX2))
  1912.       CALL MPEQ (D, C)
  1913. C
  1914. 200   IDB = IDS
  1915.       IF (IDB .GE. 6)  WRITE (LDB, 2)  (C(I),I=1,NDB)
  1916. 2     FORMAT ('MPDIVX O'/(8F9.0))
  1917.       RETURN
  1918.       END
  1919. C
  1920.       SUBROUTINE MPSQRX (A, B, C)
  1921. C
  1922. C   Computes the square root of the MP number  A  and returns the MP result
  1923. C   in  B.  C  is a scratch array with at least 15 * NW + 12 cells (fifteen
  1924. C   times the amount allocated for MP numbers will suffice).  The value
  1925. C   of NW must be equal to 2^(MW - 2).  Subroutine MPCFFT must have been
  1926. C   previously called to initialize the array U in common block MPCOM2.
  1927. C   As an aside, this subroutine also returns the reciprocal of the square
  1928. C   root of  A  beginning in location  C(2*NW+3).
  1929. C
  1930. C   This subroutine calls MPSQRT if MW is less than 8.
  1931. C
  1932. C   This subroutine employs a Newton-Raphson iteration that converges to the
  1933. C   reciprocal of the square root of A:
  1934. C
  1935. C          x(n+1) = .5 * x(n) * (3 - a * x(n)^2)
  1936. C
  1937. C   Multiplying the result by A yields the square root.  This method was
  1938. C   selected instead of the usual one because this iteration avoids divisions,
  1939. C   which are more expensive than multiplications in extra-high precision mode.
  1940. C   Note that this subroutine dynamically changes the precision level, doubling
  1941. C   it with each iteration.
  1942. C
  1943.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1944.       DIMENSION A(1), B(1), C(1)
  1945.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  1946.       PARAMETER (BX = 1D6, RX = 1D-6, MI = 5)
  1947. C
  1948.       NX = NW + 2
  1949.       IF (IDB .GE. 5)  WRITE (LDB, 1)  NX, (A(I),I=1,NDB)
  1950. 1     FORMAT (/'MPSQRX I',I10/(8F9.0))
  1951.       IDS = IDB
  1952.       IF (IDB .LE. 7)  IDB = 0
  1953. C
  1954. C   Check for zero input.
  1955. C
  1956.       IF (A(1) .EQ. 0.D0)  THEN
  1957.         DO 100 I = 1, NX
  1958.           B(I) = 0.D0
  1959. 100     CONTINUE
  1960.         GOTO 200
  1961.       ENDIF
  1962. C
  1963.       IF (A(1) .LT. 0.D0)  THEN
  1964.         WRITE (LDB, 2)
  1965. 2       FORMAT ('*** MPSQRX: NEGATIVE INPUT ***')
  1966.         CALL XABORT
  1967.       ENDIF
  1968. C
  1969.       NWS = NW
  1970.       NX1 = NX + 1
  1971.       NX2 = 2 * NX + 1
  1972.       NX3 = 3 * NX + 1
  1973.       MWS = MW
  1974. C
  1975. C   Check if precision level is too low to justify advanced algorithm.
  1976. C
  1977.       IF (MW .LE. MI)  THEN
  1978.         CALL MPSQRT (A, B, C)
  1979.         CALL MPSMC (1.D0, 0, C)
  1980.         CALL MPDIV (C, B, C(NX1), C(NX2))
  1981.         GOTO 200
  1982.       ENDIF
  1983. C
  1984. C   Compute initial approximation to two extra words precision.
  1985. C
  1986.       CALL MPSMC (0.D0, 0, B)
  1987.       MW = MI
  1988.       NW = 2 ** (MI - 2) + 2
  1989.       CALL MPSQRT (A, C, C(NX1))
  1990.       CALL MPSMC (1.D0, 0, C(NX1))
  1991.       CALL MPDIV (C(NX1), C, B, C(NX2))
  1992.       NW = NW - 2
  1993. C
  1994. C   Repeat Newton-Raphson iterations.
  1995. C
  1996.       DO 110 I = MI, MWS
  1997. C
  1998. C   Repeat the next-to-last step to insure maximum accuracy.
  1999. C
  2000.         IF (I .NE. MWS - 1)  THEN
  2001.           MW = MW + 1
  2002.           NW = 2 * NW
  2003.         ENDIF
  2004. C
  2005.         CALL MPMULX (B, B, C, C(NX2))
  2006.         CALL MPMULX (A, C, C(NX1), C(NX3))
  2007.         CALL MPSMC (3.D0, 0, C(NX2))
  2008.         CALL MPSUB (C(NX2), C(NX1), C)
  2009.         CALL MPMULX (B, C, C(NX1), C(NX3))
  2010.         CALL MPMUL1 (C(NX1), 0.5D0, 0, B)
  2011. 110   CONTINUE
  2012. C
  2013. C   Multiply result by A to complete square root calculation.
  2014. C
  2015.       CALL MPMULX (A, B, C, C(NX2))
  2016.       CALL MPEQ (B, C(NX1))
  2017.       CALL MPEQ (C, B)
  2018. C
  2019. 200   IDB = IDS
  2020.       IF (IDB .GE. 5)  WRITE (LDB, 3)  (B(I),I=1,NDB)
  2021. 3     FORMAT ('MPSQRX O'/(8F9.0))
  2022.       RETURN
  2023.       END
  2024. C
  2025.       FUNCTION XRAND (IS)
  2026. C
  2027. C   Returns a pseudorandom DP floating number between 0 and 1 if IS is zero.
  2028. C   If IS is not zero, then the seed is set with IS.  2^26 pseudorandom numbers
  2029. C   with 28 bits each are returned before repeating.
  2030. C
  2031. C  This subroutine is included for completeness only.  If a library pseudo-
  2032. C  random number generator is available, it should be used instead.
  2033. C
  2034.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  2035.       PARAMETER (F7 = 5.D0 ** 7,  R28 = 2.D0 ** (-28),
  2036.      $   S0 = 31415927.D0 * R28)
  2037.       SAVE S
  2038.       DATA S/S0/
  2039. C
  2040.       IF (IS .NE. 0)  S = IS * R28
  2041.       XRAND = S
  2042.       T = F7 * S
  2043.       S = T - INT (T)
  2044. C
  2045.       RETURN
  2046.       END
  2047. C
  2048.       FUNCTION XTIME (TX)
  2049. C
  2050. C   This subroutine is included for completeness only.  Most Fortran systems
  2051. C   have a system call, often called SECOND, that returns the cumulative job
  2052. C   CPU time in seconds.  Calls to the appropriate system routine should be
  2053. C   immediately substituted for the calls to XTIME in the above program.
  2054. C
  2055.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  2056.       XTIME = SECNDS(0.0)
  2057.       RETURN
  2058.       END
  2059. C
  2060.       SUBROUTINE XABORT
  2061. C
  2062. C   This subroutine is included for completeness only.  Most Fortran systems
  2063. C   have a system call, often called ABORT, that terminates execution with
  2064. C   a traceback message.  Calls to the appropriate system routine should be
  2065. C   immediately substituted for the calls to XABORT in the above program.
  2066. C
  2067.       STOP
  2068.       END
  2069.  
  2070.       PROGRAM PI4
  2071. C
  2072. C   Computes PI to any requested number of digits, up to a maximum of about
  2073. C   ten million, depending on the floating-point precision available on the
  2074. C   computer being used.  This program is intended for use as a supercomputer
  2075. C   system check.  This is the quartically convergent, double-precision,
  2076. C   radix-2, complex FFT version.
  2077. C
  2078. C   This program is designed for high performance on vector computer systems.
  2079. C   Vectorizable loops that are frequently not recognized as such are
  2080. C   preceded by the Cray compiler directive CDIR$ IVDEP.  These directives
  2081. C   should be replaced by their equivalent on other vector computer systems.
  2082. C   For scalar computers, search for the word 'scalar' in comments and make
  2083. C   the recommended modifications for best performance.
  2084. C
  2085. C   David H. Bailey    May 9, 1988
  2086. C
  2087. C   The algorithm that is used for computing PI is as follows:
  2088. C
  2089. C   Set   A(0) = 6 - 4 * SQRT(2),  Y(0) = SQRT (2) - 1
  2090. C
  2091. C   and then iterate the following operations:
  2092. C
  2093. C      Y(N+1)  =  [1 - (1 - Y(N)^4) ^ .25] / [1 + (1 - Y(N)^4) ^ .25]
  2094. C      A(N+1)  =  A(N) * (1 + Y(N+1)) ^ 4  -  2 * 4^(N+1)
  2095. C                 * Y(N+1) * (1 + Y(N+1) + Y(N+1)^2)
  2096. C
  2097. C   Then  A(N)  converges quartically to  1 / PI.  To be specific,
  2098. C
  2099. C      | A(N) - 1/PI |  <  32 * 4^N * EXP (-2 * PI * 4^N)
  2100. C
  2101. C   Reference:  Borwein & Borwein, "Elliptic Integrals and Approximations
  2102. C               to PI", to appear
  2103. C
  2104. C   Set the parameter MX in the following PARAMETER statement to set the
  2105. C   level of precision.  NDX is the number of digits of precision.
  2106. C
  2107. C   JDB 20090806
  2108. C   The Intel CDEC$ IVDEP directive was added after each equivilant
  2109. C   Cray directive.
  2110. C
  2111. C   JDB 20090806
  2112. C   Modified the time reporting to breakdown to hours, mins, secs, and
  2113. C   subseconds.  Added newlines to the initial banner, and to the end
  2114. C   of the time display to make standard output more readable.
  2115. C
  2116.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  2117.       CHARACTER*1 PID
  2118.       PARAMETER  (MX = 23, NX = 2 ** (MX-2) + 2, NDX = 6 * (NX - 3),
  2119.      $  NU = 8 * NX)
  2120.       DIMENSION  A(NX), Y(NX), C1(NX), X1(NX), X2(NX), X3(NX), X4(2*NX),
  2121.      $  S(15*NX), PID(NDX+100)
  2122.       COMMON /MPCOM1/ ND, NW, MW, IDB, LDB, NDB
  2123.       COMMON /MPCOM2/ U(NU)
  2124. C
  2125.       T_HOURS = 0.0
  2126.       T_MINS = 0.0
  2127.       T_SSECS = 0.0
  2128.  
  2129.       ND = NDX
  2130.       MW = MX
  2131.       NW = NX - 2
  2132.       IDB = 0
  2133.       LDB = 0
  2134.       NDB = 16
  2135.       NI = (MX + 1) / 2
  2136.       CALL MPCFFT (0, A, Y)
  2137.       WRITE (LDB, 1)  MW, NW, ND
  2138. 1     FORMAT (/'PI4 COMPUTATION TEST -- DP COMPLEX FFT MP VERSION'//
  2139.      $  'MW =', I3, 3X, 'NW =', I10, 3X, 'ND =', I10/)
  2140. C
  2141.       TA = XTIME (TX)
  2142.       CALL MPSMC (1.D0, 0, C1)
  2143.       CALL MPSMC (2.D0, 0, X1)
  2144.       CALL MPSQRX (X1, X2, S)
  2145.       CALL MPMUL1 (X2, 4.D0, 0, X1)
  2146.       CALL MPSMC (6.D0, 0, X3)
  2147.       CALL MPSUB (X3, X1, A)
  2148.       CALL MPSUB (X2, C1, Y)
  2149. C
  2150.       DO 100 K = 1, NI
  2151.         WRITE (LDB, 2)  K, NI
  2152. 2       FORMAT ('ITERATION', 2I4)
  2153. C
  2154. C   Update Y.
  2155. C
  2156.         CALL MPMULX (Y, Y, X4, S)
  2157.         CALL MPEQ (X4, X1)
  2158.         CALL MPMULX (X1, X1, X4, S)
  2159.         CALL MPSUB (C1, X4, X1)
  2160.         CALL MPSQRX (X1, X2, S)
  2161.         CALL MPSQRX (X2, X1, S)
  2162.         CALL MPSUB (C1, X1, X2)
  2163.         CALL MPADD (C1, X1, X3)
  2164.         CALL MPDIVX (X2, X3, Y, S)
  2165. C
  2166. C   Update A.
  2167. C
  2168.         CALL MPADD (C1, Y, X1)
  2169.         CALL MPMULX (X1, X1, X4, S)
  2170.         CALL MPEQ (X4, X2)
  2171.         CALL MPMULX (X2, X2, X4, S)
  2172.         CALL MPEQ (X4, X2)
  2173.         CALL MPMULX (A, X2, X4, S)
  2174.         CALL MPEQ (X4, X2)
  2175.         CALL MPMULX (Y, Y, X4, S)
  2176.         CALL MPADD (X1, X4, X3)
  2177.         CALL MPMULX (Y, X3, X4, S)
  2178.         CX = 2.D0 * 4.D0 ** K
  2179.         CALL MPMUL1 (X4, CX, 0, X3)
  2180.         CALL MPSUB (X2, X3, A)
  2181. 100   CONTINUE
  2182. C
  2183. C   Output results -- the reciprocal of A is the final output value.
  2184. C
  2185. 110   CALL MPDIVX (C1, A, X1, S)
  2186.  
  2187.       TB = XTIME (TX)
  2188.       T_SECS = TB - TA
  2189.  
  2190.       IF (T_SECS .GE. 3600) THEN
  2191.          T_HOURS = INT(T_SECS) / 3600
  2192.          T_SECS = T_SECS - (T_HOURS * 3600)
  2193.       ENDIF
  2194.  
  2195.       IF (T_SECS .GE. 60) THEN
  2196.          T_MINS = INT(T_SECS) / 60
  2197.          T_SECS = T_SECS - (T_MINS * 60)
  2198.       ENDIF
  2199.  
  2200.       IF (T_SECS .GE. 1) THEN
  2201.          T_SSECS = T_SECS - INT(T_SECS)
  2202.          T_SECS = INT(T_SECS)
  2203.       ELSE
  2204.          T_SSECS = T_SECS
  2205.       ENDIF
  2206.  
  2207.       WRITE (LDB, 3) INT(T_HOURS), INT(T_MINS), INT(T_SECS), T_SSECS
  2208.       WRITE (MX, 3) INT(T_HOURS), INT(T_MINS), INT(T_SECS), T_SSECS
  2209.  
  2210. 3     FORMAT (/'CPU TIME = ',I2.2,':',I2.2,':',I2.2,F9.8/)
  2211.  
  2212.       CALL MPFMT (X1, PID)
  2213.       WRITE (MX, 4) (PID(I), I = 1, 20)
  2214. 4     FORMAT (10(80A1/))
  2215.       WRITE (MX, 4) (PID(I), I = 21, ND+20)
  2216. C
  2217.       STOP
  2218.       END
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!