home *** CD-ROM | disk | FTP | other *** search
/ High Voltage Shareware / high1.zip / high1 / DIR9 / ACMALG01.ZIP / ACM665.FOR < prev    next >
Text File  |  1991-02-18  |  13KB  |  291 lines
  1. C      ALGORITHM 665, COLLECTED ALGORITHMS FROM ACM.
  2. C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  3. C      VOL. 14, NO. 4, PP. 303-311.
  4.       SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  5.      1                   MAXEXP,EPS,EPSNEG,XMIN,XMAX)
  6. C-----------------------------------------------------------------------
  7. C  This Fortran 77 subroutine is intended to determine the parameters
  8. C   of the floating-point arithmetic system specified below.  The
  9. C   determination of the first three uses an extension of an algorithm
  10. C   due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
  11. C   but not all, of the improvements suggested by M. Gentleman and S.
  12. C   Marovich, CACM 17 (1974), pp. 276-277.  An earlier version of this
  13. C   program was published in the book Software Manual for the
  14. C   Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
  15. C   Englewood Cliffs, NJ, 1980.
  16. C
  17. C  The program as given here must be modified before compiling.  If
  18. C   a single (double) precision version is desired, change all
  19. C   occurrences of CS (CD) in columns 1 and 2 to blanks.
  20. C
  21. C  Parameter values reported are as follows:
  22. C
  23. C       IBETA   - the radix for the floating-point representation
  24. C       IT      - the number of base IBETA digits in the floating-point
  25. C                 significand
  26. C       IRND    - 0 if floating-point addition chops
  27. C                 1 if floating-point addition rounds, but not in the
  28. C                   IEEE style
  29. C                 2 if floating-point addition rounds in the IEEE style
  30. C                 3 if floating-point addition chops, and there is
  31. C                   partial underflow
  32. C                 4 if floating-point addition rounds, but not in the
  33. C                   IEEE style, and there is partial underflow
  34. C                 5 if floating-point addition rounds in the IEEE style,
  35. C                   and there is partial underflow
  36. C       NGRD    - the number of guard digits for multiplication with
  37. C                 truncating arithmetic.  It is
  38. C                 0 if floating-point arithmetic rounds, or if it
  39. C                   truncates and only  IT  base  IBETA digits
  40. C                   participate in the post-normalization shift of the
  41. C                   floating-point significand in multiplication;
  42. C                 1 if floating-point arithmetic truncates and more
  43. C                   than  IT  base  IBETA  digits participate in the
  44. C                   post-normalization shift of the floating-point
  45. C                   significand in multiplication.
  46. C       MACHEP  - the largest negative integer such that
  47. C                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
  48. C                 MACHEP is bounded below by  -(IT+3)
  49. C       NEGEPS  - the largest negative integer such that
  50. C                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
  51. C                 NEGEPS is bounded below by  -(IT+3)
  52. C       IEXP    - the number of bits (decimal places if IBETA = 10)
  53. C                 reserved for the representation of the exponent
  54. C                 (including the bias or sign) of a floating-point
  55. C                 number
  56. C       MINEXP  - the largest in magnitude negative integer such that
  57. C                 FLOAT(IBETA)**MINEXP is positive and normalized
  58. C       MAXEXP  - the smallest positive power of  BETA  that overflows
  59. C       EPS     - the smallest positive floating-point number such
  60. C                 that  1.0+EPS .NE. 1.0. In particular, if either
  61. C                 IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
  62. C                 Otherwise,  EPS = (FLOAT(IBETA)**MACHEP)/2
  63. C       EPSNEG  - A small positive floating-point number such that
  64. C                 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
  65. C                 or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
  66. C                 Otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
  67. C                 NEGEPS is bounded below by -(IT+3), EPSNEG may not
  68. C                 be the smallest number that can alter 1.0 by
  69. C                 subtraction.
  70. C       XMIN    - the smallest non-vanishing normalized floating-point
  71. C                 power of the radix, i.e.,  XMIN = FLOAT(IBETA)**MINEXP
  72. C       XMAX    - the largest finite floating-point number.  In
  73. C                 particular  XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
  74. C                 Note - on some machines  XMAX  will be only the
  75. C                 second, or perhaps third, largest number, being
  76. C                 too small by 1 or 2 units in the last digit of
  77. C                 the significand.
  78. C
  79. C     Latest revision - April 20, 1987
  80. C
  81. C     Author - W. J. Cody
  82. C              Argonne National Laboratory
  83. C
  84. C-----------------------------------------------------------------------
  85.       INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP,
  86.      1        MINEXP,MX,NEGEP,NGRD,NXRES
  87. CS    REAL A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA,
  88. CS   1     TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
  89. CD    DOUBLE PRECISION A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,
  90. CD   1                 T,TEMP,TEMPA,TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
  91. C-----------------------------------------------------------------------
  92. CS    CONV(I) = REAL(I)
  93. CD    CONV(I) = DBLE(I)
  94.       ONE = CONV(1)
  95.       TWO = ONE + ONE
  96.       ZERO = ONE - ONE
  97. C-----------------------------------------------------------------------
  98. C  Determine IBETA, BETA ala Malcolm.
  99. C-----------------------------------------------------------------------
  100.       A = ONE
  101.    10 A = A + A
  102.          TEMP = A+ONE
  103.          TEMP1 = TEMP-A
  104.          IF (TEMP1-ONE .EQ. ZERO) GO TO 10
  105.       B = ONE
  106.    20 B = B + B
  107.          TEMP = A+B
  108.          ITEMP = INT(TEMP-A)
  109.          IF (ITEMP .EQ. 0) GO TO 20
  110.       IBETA = ITEMP
  111.       BETA = CONV(IBETA)
  112. C-----------------------------------------------------------------------
  113. C  Determine IT, IRND.
  114. C-----------------------------------------------------------------------
  115.       IT = 0
  116.       B = ONE
  117.   100 IT = IT + 1
  118.          B = B * BETA
  119.          TEMP = B+ONE
  120.          TEMP1 = TEMP-B
  121.          IF (TEMP1-ONE .EQ. ZERO) GO TO 100
  122.       IRND = 0
  123.       BETAH = BETA / TWO
  124.       TEMP = A+BETAH
  125.       IF (TEMP-A .NE. ZERO) IRND = 1
  126.       TEMPA = A + BETA
  127.       TEMP = TEMPA+BETAH
  128.       IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2
  129. C-----------------------------------------------------------------------
  130. C  Determine NEGEP, EPSNEG.
  131. C-----------------------------------------------------------------------
  132.       NEGEP = IT + 3
  133.       BETAIN = ONE / BETA
  134.       A = ONE
  135.       DO 200 I = 1, NEGEP
  136.          A = A * BETAIN
  137.   200 CONTINUE
  138.       B = A
  139.   210 TEMP = ONE-A
  140.          IF (TEMP-ONE .NE. ZERO) GO TO 220
  141.          A = A * BETA
  142.          NEGEP = NEGEP - 1
  143.       GO TO 210
  144.   220 NEGEP = -NEGEP
  145.       EPSNEG = A
  146.       IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 300
  147.       A = (A*(ONE+A)) / TWO
  148.       TEMP = ONE-A
  149.       IF (TEMP-ONE .NE. ZERO) EPSNEG = A
  150. C-----------------------------------------------------------------------
  151. C  Determine MACHEP, EPS.
  152. C-----------------------------------------------------------------------
  153.   300 MACHEP = -IT - 3
  154.       A = B
  155.   310 TEMP = ONE+A
  156.          IF (TEMP-ONE .NE. ZERO) GO TO 320
  157.          A = A * BETA
  158.          MACHEP = MACHEP + 1
  159.       GO TO 310
  160.   320 EPS = A
  161.       TEMP = TEMPA+BETA*(ONE+EPS)
  162.       IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 350
  163.       A = (A*(ONE+A)) / TWO
  164.       TEMP = ONE+A
  165.       IF (TEMP-ONE .NE. ZERO) EPS = A
  166. C-----------------------------------------------------------------------
  167. C  Determine NGRD.
  168. C-----------------------------------------------------------------------
  169.   350 NGRD = 0
  170.       TEMP = ONE+EPS
  171.       IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1
  172. C-----------------------------------------------------------------------
  173. C  Determine IEXP, MINEXP, XMIN.
  174. C
  175. C  Loop to determine largest I and K = 2**I such that
  176. C         (1/BETA) ** (2**(I))
  177. C  does not underflow.
  178. C  Exit from loop is signaled by an underflow.
  179. C-----------------------------------------------------------------------
  180.       I = 0
  181.       K = 1
  182.       Z = BETAIN
  183.       T = ONE + EPS
  184.       NXRES = 0
  185.   400 Y = Z
  186.          Z = Y * Y
  187. C-----------------------------------------------------------------------
  188. C  Check for underflow here.
  189. C-----------------------------------------------------------------------
  190.          A = Z * ONE
  191.          TEMP = Z * T
  192.          IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410
  193.          TEMP1 = TEMP * BETAIN
  194.          IF (TEMP1*BETA .EQ. Z) GO TO 410
  195.          I = I + 1
  196.          K = K + K
  197.       GO TO 400
  198.   410 IF (IBETA .EQ. 10) GO TO 420
  199.       IEXP = I + 1
  200.       MX = K + K
  201.       GO TO 450
  202. C-----------------------------------------------------------------------
  203. C  This segment is for decimal machines only.
  204. C-----------------------------------------------------------------------
  205.   420 IEXP = 2
  206.       IZ = IBETA
  207.   430 IF (K .LT. IZ) GO TO 440
  208.          IZ = IZ * IBETA
  209.          IEXP = IEXP + 1
  210.       GO TO 430
  211.   440 MX = IZ + IZ - 1
  212. C-----------------------------------------------------------------------
  213. C  Loop to determine MINEXP, XMIN.
  214. C  Exit from loop is signaled by an underflow.
  215. C-----------------------------------------------------------------------
  216.   450 XMIN = Y
  217.          Y = Y * BETAIN
  218. C-----------------------------------------------------------------------
  219. C  Check for underflow here.
  220. C-----------------------------------------------------------------------
  221.          A = Y * ONE
  222.          TEMP = Y * T
  223.          IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460
  224.          K = K + 1
  225.          TEMP1 = TEMP * BETAIN
  226.          IF (TEMP1*BETA .NE. Y) GO TO 450
  227.       NXRES = 3
  228.       XMIN = Y
  229.   460 MINEXP = -K
  230. C-----------------------------------------------------------------------
  231. C  Determine MAXEXP, XMAX.
  232. C-----------------------------------------------------------------------
  233.       IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500
  234.       MX = MX + MX
  235.       IEXP = IEXP + 1
  236.   500 MAXEXP = MX + MINEXP
  237. C-----------------------------------------------------------------
  238. C  Adjust IRND to reflect partial underflow.
  239. C-----------------------------------------------------------------
  240.       IRND = IRND + NXRES
  241. C-----------------------------------------------------------------
  242. C  Adjust for IEEE-style machines.
  243. C-----------------------------------------------------------------
  244.       IF ((IRND .EQ. 2) .OR. (IRND .EQ. 5)) MAXEXP = MAXEXP - 2
  245. C-----------------------------------------------------------------
  246. C  Adjust for non-IEEE machines with partial underflow.
  247. C-----------------------------------------------------------------
  248.       IF ((IRND .EQ. 3) .OR. (IRND .EQ. 4)) MAXEXP = MAXEXP - IT
  249. C-----------------------------------------------------------------
  250. C  Adjust for machines with implicit leading bit in binary
  251. C  significand, and machines with radix point at extreme
  252. C  right of significand.
  253. C-----------------------------------------------------------------
  254.       I = MAXEXP + MINEXP
  255.       IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1
  256.       IF (I .GT. 20) MAXEXP = MAXEXP - 1
  257.       IF (A .NE. Y) MAXEXP = MAXEXP - 2
  258.       XMAX = ONE - EPSNEG
  259.       IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
  260.       XMAX = XMAX / (BETA * BETA * BETA * XMIN)
  261.       I = MAXEXP + MINEXP + 3
  262.       IF (I .LE. 0) GO TO 520
  263.       DO 510 J = 1, I
  264.           IF (IBETA .EQ. 2) XMAX = XMAX + XMAX
  265.           IF (IBETA .NE. 2) XMAX = XMAX * BETA
  266.   510 CONTINUE
  267.   520 RETURN
  268. C---------- LAST CARD OF MACHAR ----------
  269.       END
  270. C
  271. C   TEST DRIVER FOR MACHAR
  272. C
  273.       INTEGER IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP
  274. CS    REAL EPS,EPSNEG,XMIN,XMAX
  275. CD    DOUBLE PRECISION EPS,EPSNEG,XMIN,XMAX
  276. C
  277.       WRITE(6,1001)
  278.       CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
  279.      1      EPS,EPSNEG,XMIN,XMAX)
  280.       WRITE (6,1000) IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
  281.      1   MAXEXP,EPS,EPSNEG,XMIN,XMAX
  282.       STOP
  283.  1000 FORMAT(19H OUTPUT FROM MACHAR//7H BETA =,I5/7H    T =,I5/
  284.      1   7H  RND =,I5/7H NGRD =,I5/9H MACHEP =,I5/8H NEGEP =,I5/
  285.      2   7H IEXP =,I5/9H MINEXP =,I5/9H MAXEXP =,I5/6H EPS =,
  286.      3   E25.13/9H EPSNEG =,E25.13/7H XMIN =,E25.13/
  287.      4   7H XMAX =,E25.13)
  288.  1001 FORMAT(18H1GOING INTO MACHAR/)
  289. C  ---------- LAST CARD OF DRIVER ----------
  290.       END
  291.