C     PROGRAM TO TEST ALOG
C
C     DATA REQUIRED
C
C        NONE
C
C     SUBPROGRAMS REQUIRED FROM THIS PACKAGE
C
C        MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
C                 INFORMATION ON THE FLOATING-POINT ARITHMETIC
C                 SYSTEM.  NOTE THAT THE CALL TO MACHAR CAN
C                 BE DELETED PROVIDED THE FOLLOWING FOUR
C                 PARAMETERS ARE ASSIGNED THE VALUES INDICATED
C
C                 IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM
C                 IT    - THE NUMBER OF BASE-IBETA DIGITS IN THE
C                         SIGNIFICAND OF A FLOATING-POINT NUMBER
C                 XMIN  - THE SMALLEST NON-VANISHING FLOATING-POINT
C                         POWER OF THE RADIX
C                 XMAX  - THE LARGEST FINITE FLOATING-POINT NO.
C
C        REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
C                 NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
C
C
C     STANDARD FORTRAN SUBPROGRAMS REQUIRED
C
C         ABS, ALOG, ALOG10, AMAX1, FLOAT, SIGN, SQRT
C
C
C     LATEST REVISION - DECEMBER 6, 1979
C
C     AUTHOR - W. J. CODY
C              ARGONNE NATIONAL LABORATORY
C
C
      INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,
     1        MAXEXP,MINEXP,N,NEGEP,NGRD
      REAL A,AIT,ALBETA,B,BETA,C,DEL,EIGHT,EPS,EPSNEG,HALF,ONE,
     1     RAN,R6,R7,TENTH,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ
C
      IOUT = 6
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      BETA = FLOAT(IBETA)
      ALBETA = ALOG(BETA)
      AIT = FLOAT(IT)
      J = IT / 3
      ZERO = 0.0E0
      HALF = 0.5E0
      EIGHT = 8.0E0
      TENTH = 0.1E0
      ONE = 1.0E0
      C = ONE
C
      DO 50 I = 1, J
         C = C / BETA
   50 CONTINUE
C
      B = ONE + C
      A = ONE - C
      N = 2000
      XN = FLOAT(N)
      I1 = 0
C-----------------------------------------------------------------
C     RANDOM ARGUMENT ACCURACY TESTS
C-----------------------------------------------------------------
      DO 300 J = 1, 4
         K1 = 0
         K3 = 0
         X1 = ZERO
         R6 = ZERO
         R7 = ZERO
         DEL = (B - A) / XN
         XL = A
C
         DO 200 I = 1, N
            X = DEL * REN(I1) + XL
            IF (J .NE. 1) GO TO 100
            Y = (X - HALF) - HALF
            ZZ = ALOG(X)
            Z = ONE / 3.0E0
            Z = Y * (Z - Y / 4.0E0)
            Z = (Z - HALF) * Y * Y + Y
            GO TO 150
  100       IF (J .NE. 2) GO TO 110
            X = (X + EIGHT) - EIGHT
            Y = X + X / 16.0E0
            Z = ALOG(X)
            ZZ = ALOG(Y) - 7.7746816434842581E-5
            ZZ = ZZ - 31.0E0/512.0E0
            GO TO 150
  110       IF (J .NE. 3) GO TO 120
            X = (X + EIGHT) - EIGHT
            Y = X + X * TENTH
            Z = ALOG10(X)
            ZZ = ALOG10(Y) - 3.7706015822504075E-4
            ZZ = ZZ - 21.0E0/512.0E0
            GO TO 150
  120       Z = ALOG(X*X)
            ZZ = ALOG(X)
            ZZ = ZZ + ZZ
  150       W = ONE
            IF (Z .NE. ZERO) W = (Z - ZZ) / Z
            Z = SIGN(W,Z)
            IF (Z .GT. ZERO) K1 = K1 + 1
            IF (Z .LT. ZERO) K3 = K3 + 1
            W = ABS(W)
            IF (W .LE. R6) GO TO 160
            R6 = W
            X1 = X
  160       R7 = R7 + W*W
            XL = XL + DEL
  200    CONTINUE
C
         K2 = N - K3 - K1
         R7 = SQRT(R7/XN)
         IF (J .EQ. 1) WRITE (IOUT,1000)
         IF (J .EQ. 2) WRITE (IOUT,1001)
         IF (J .EQ. 3) WRITE (IOUT,1005)
         IF (J .EQ. 4) WRITE (IOUT,1002)
         IF (J .EQ. 1) WRITE (IOUT,1009) N,C
         IF (J .NE. 1) WRITE (IOUT,1010) N,A,B
         IF (J .NE. 3) WRITE (IOUT,1011) K1,K2,K3
         IF (J .EQ. 3) WRITE (IOUT,1012) K1,K2,K3
         WRITE (IOUT,1020) IT,IBETA
         W = -999.0E0
         IF (R6 .NE. ZERO) W = ALOG(ABS(R6))/ALBETA
         WRITE (IOUT,1021) R6,IBETA,W,X1
         W = AMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         W = -999.0E0
         IF (R7 .NE. ZERO) W = ALOG(ABS(R7))/ALBETA
         WRITE (IOUT,1023) R7,IBETA,W
         W = AMAX1(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         IF (J .GT. 1) GO TO 230
         A = SQRT(HALF)
         B = 15.0E0 / 16.0E0
         GO TO 300
  230    IF (J .GT. 2) GO TO 240
         A = SQRT(TENTH)
         B = 0.9E0
         GO TO 300
  240    A = 16.0E0
         B = 240.0E0
  300 CONTINUE
C-----------------------------------------------------------------
C     SPECIAL TESTS
C-----------------------------------------------------------------
      WRITE (IOUT,1025)
      WRITE (IOUT,1030)
C
      DO 320 I = 1, 5
         X = REN(I1)
         X = X + X + 15.0E0
         Y = ONE / X
         Z = ALOG(X) + ALOG(Y)
         WRITE (IOUT,1060) X, Z
  320 CONTINUE
C
      WRITE (IOUT,1040)
      X = ONE
      Y = ALOG(X)
      WRITE (IOUT,1041) Y
      X = XMIN
      Y = ALOG(X)
      WRITE (IOUT,1042) X, Y
      X = XMAX
      Y = ALOG(X)
      WRITE (IOUT,1043) X, Y
C-----------------------------------------------------------------
C     TEST OF ERROR RETURNS
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      X = -2.0E0
      WRITE (IOUT,1052) X
      Y = ALOG(X)
      WRITE (IOUT,1055) Y
      X = ZERO
      WRITE (IOUT,1052) X
      Y = ALOG(X)
      WRITE (IOUT,1055) Y
      WRITE (IOUT,1100)
      STOP
 1000 FORMAT(47H1TEST OF ALOG(X) VS T.S. EXPANSION OF ALOG(1+Y)  //)
 1001 FORMAT(44H1TEST OF ALOG(X) VS ALOG(17X/16)-ALOG(17/16)   //)
 1002 FORMAT(32H1TEST OF ALOG(X*X) VS 2 * LOG(X)  //)
 1005 FORMAT(50H1TEST OF ALOG10(X) VS ALOG10(11X/10)-ALOG10(11/10) //)
 1009 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL /
     1 6X,26H(1-EPS,1+EPS), WHERE EPS =, E15.4//)
 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL /
     1 6X,1H(,E15.4,1H,,E15.4,1H)//)
 1011 FORMAT(19H ALOG(X) WAS LARGER,I6,7H TIMES, /
     1     12X,7H AGREED,I6,11H TIMES, AND /
     2     8X,11HWAS SMALLER,I6,7H TIMES.//)
 1012 FORMAT(21H ALOG10(X) WAS LARGER,I6,7H TIMES, /
     1     14X,7H AGREED,I6,11H TIMES, AND /
     2     10X,11HWAS SMALLER,I6,7H TIMES.//)
 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4,
     1    46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER  //)
 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **,
     1  F7.2/4X,16HOCCURRED FOR X =,E17.6)
 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4,
     1  22H SIGNIFICANT DIGITS IS,F7.2//)
 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,E15.4,
     1    3H = ,I4,3H **,F7.2)
 1025 FORMAT(14H1SPECIAL TESTS//)
 1030 FORMAT(52H THE IDENTITY  ALOG(X) = -ALOG(1/X)  WILL BE TESTED.//
     1      8X,1HX,9X,13HF(X) + F(1/X)/)
 1040 FORMAT(//26H TEST OF SPECIAL ARGUMENTS //)
 1041 FORMAT(13H ALOG(1.0) = ,E15.7//)
 1042 FORMAT(19H ALOG(XMIN) = ALOG(,E15.7,4H) = ,E15.7//)
 1043 FORMAT(19H ALOG(XMAX) = ALOG(,E15.7,4H) = ,E15.7//)
 1050 FORMAT(22H1TEST OF ERROR RETURNS//)
 1052 FORMAT(38H ALOG WILL BE CALLED WITH THE ARGUMENT,E15.4/
     1       37H THIS SHOULD TRIGGER AN ERROR MESSAGE//)
 1055 FORMAT(24H ALOG RETURNED THE VALUE,E15.4///)
 1060 FORMAT(2E15.7/)
 1100 FORMAT(25H THIS CONCLUDES THE TESTS )
C     ---------- LAST CARD OF ALOG/ALOG10 TEST PROGRAM ----------
      END
      REAL FUNCTION REN(K)
C
C     RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
C      HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
C      VOL. 8, NO. 10, OCTOBER 1965.
C
C     THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
C      FIXED POINT WORDLENGTH OF AT LEAST 29 BITS.  IT IS
C      BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
C      29 BITS.
C
      INTEGER IY,J,K
      DATA IY/100001/
C
      J = K
      IY = IY * 125
      IY = IY - (IY/2796203) * 2796203
      REN = FLOAT(IY) / 2796203.0E0
      RETURN
C     ---------- LAST CARD OF RAN ----------
      END