C  Program to test CLOG
C
C  Accuracy tests are based on the identity
C
C            log(z) = log(z*z)/2
C     and
C            log(z) = log(17z/16) - log(17/16)
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     RESET, TABLAT, REPORT - programs to report results.
C
C
C  Generic Fortran subprograms required:
C
C         CMPLX (or DCMPLX), LOG, REAL (or DBLE), SQRT
C
C  The program as given here must be modified before compiling.
C   If a single (double) precision version is desired, change all
C   occurrences of CS (CD) in columns 1 and 2 to blanks.  For a
C   calibration version, change all occurrences of CS and of CC
C   in columns 1 and 2 to blanks.
C
C  Latest revision - November 20, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
C  Declare named common
C----------------------------------------------------------------------
      COMMON /ELEFUN/ CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,
     1       CZ,CZZ,AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN,
     2       IBETA,IT,KX1,KX2,KX3,KY1,KY2,KY3,K1,K2,K3,N,
     3       IFLAG,RFLAG,PFLAG
      LOGICAL IFLAG,RFLAG,PFLAG
      INTEGER IBETA,IT,KX1,KX2,KX3,KY1,KY2,KY3,K1,K2,K3,N
CS    REAL
      DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
      COMPLEX*16
     1       CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,CZ,CZZ
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
      INTEGER I,IEXP,IOUT,IRND,I1,J,MACHEP,MAXEXP,MINEXP,NDUM,
     1        NEGEP,NGRD
CS    REAL
      DOUBLE PRECISION
     1   BETA,CONV,DELX,DELY,EPS,EPSNEG,HALF,ONE,REN,SCALE,TEN,
     2   THOUS,W,X,XL,XMAX,XMIN,Y,Z,ZERO
CS    COMPLEX
      COMPLEX*16
     1   CONVC,CY
CC   2   ,FLOG
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,HALF,ONE,TEN,THOUS/0.0E0,0.5E0,1.0E0,10.0E0,1.0E3/
      DATA ZERO,HALF,ONE,TEN,THOUS/0.0D0,0.5D0,1.0D0,10.0D0,1.0D3/
C----------------------------------------------------------------------
C  Statement functions for conversion between integer or complex
C  and float
C----------------------------------------------------------------------
CS    CONV(NDUM) = REAL(NDUM)
CS    CONVC(X,Y) = CMPLX(X,Y)
      CONV(NDUM) = DBLE(NDUM)
      CONVC(X,Y) = DCMPLX(X,Y)
C----------------------------------------------------------------------
      CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1            MAXEXP,EPS,EPSNEG,XMIN,XMAX)
      BETA = CONV(IBETA)
      ALBETA = LOG(BETA)
      AIT = CONV(IT)
      I1 = IT/2 + 4
      SCALE = BETA ** I1
      AX = ONE + ONE
      BX = TEN
      AY = ZERO
      BY = BX
      N = 2000
      XN = CONV(N)
      I1 = 0
      RFLAG = .TRUE.
      IFLAG = .TRUE.
      PFLAG = .FALSE.
C---------------------------------------------------------------------
C  Random argument accuracy tests
C---------------------------------------------------------------------
      DO 300 J = 1, 3
         CALL RESET
         DELX = (BX - AX) / XN
         DELY = (BY - AY)
         XL = AX
         DO 200 I = 1, N
            X = DELX * REN(I1) + XL
            Y = DELY * REN(I1) + AY
C-----------------------------------------------------------------
C  Purify argument and evaluate identities
C-----------------------------------------------------------------
            Z = X * SCALE
            W = Z + X
            X = W - Z
            Z = Y * SCALE
            W = Z + Y
            Y = W - Z
            CX = CONVC(X,Y)
            Z = X*X - Y*Y
            W = X*Y
            CY = CONVC(Z,W+W)
            CZ = LOG(CX)
            CZZ = LOG(CY)
CC          CZ = FLOG(CX)
CC          CZZ = FLOG(CY)
            CZZ = CZZ * HALF
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
            CALL TABLAT
            XL = XL + DELX
  200    CONTINUE
C-----------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------
         WRITE (IOUT,1000)
         CALL REPORT(IOUT)
         IF (J .EQ. 1) THEN
               AX = THOUS
               BX = AX + AX
               AY = -AX
               BY = -BX - BX
            ELSE
               AX = EPS
               BX = HALF*HALF
               AY = -AX
               BY = -BX
         END IF
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1025)
      WRITE (IOUT,1030)
      DO 320 I = 1, 5
         X = REN(I1) * TEN
         Y = REN(I1) * TEN
         CX = CONVC(X,Y)
         CY = ONE / CX
         CZ = LOG(CX) + LOG(CY)
         WRITE (IOUT,1031) CX,CZ
  320 CONTINUE
      WRITE (IOUT,1040)
      CX = CONVC(ONE,ZERO)
      CZ = LOG(CX)
      WRITE (IOUT,1041) CZ
      WRITE (IOUT,1050)
      CX = CONVC(XMIN,XMIN)
      WRITE (IOUT,1061) CX
      CZ = LOG(CX)
      WRITE (IOUT,1051) CX,CZ
      CX = CONVC(XMAX,XMAX)
      WRITE (IOUT,1061) CX
      CZ = LOG(CX)
      WRITE (IOUT,1052) CX,CZ
C-----------------------------------------------------------------
C  Test of error returns
C-----------------------------------------------------------------
      WRITE (IOUT,1060)
      CX = CONVC(ZERO,ZERO)
      WRITE (IOUT,1062) CX
      CZ = LOG(CX)
      WRITE (IOUT,1065) CZ
      WRITE (IOUT,1100)
      STOP
C---------------------------------------------------------------------
 1000 FORMAT(' Test of LOG(Z) vs LOG(Z*Z)/2'//)
 1025 FORMAT(' Special tests'//)
 1030 FORMAT(' The identity  LOG(Z) = -LOG(1/Z)  will be tested.'//
     1      15X,'Z',22X,'LOG(Z) + LOG(1/Z)'/)
 1031 FORMAT(2(E15.7,E15.7)/)
 1040 FORMAT(//' Test of special arguments'//)
 1041 FORMAT(' LOG(1.0,0.0) = (',(E15.7,E15.7),')'//)
 1050 FORMAT(//' Tests near extreme arguments.  These should not,',
     1         ' but may'/' trigger error messages'//)
 1051 FORMAT(' LOG(XMIN,XMIN) = LOG(',(E15.7,E15.7),') ='/40X,'(',
     1         (E15.7,E15.7),')'/)
 1052 FORMAT(' LOG(XMAX,XMAX) = LOG(',(E15.7,E15.7),') ='/40X,'(',
     1         (E15.7,E15.7),')'/)
 1060 FORMAT(' Test of error returns'/)
 1061 FORMAT(/' LOG will be called with the argument (',
     1          (E15.4,E15.4),')'//)
 1062 FORMAT(/' LOG will be called with the argument (',
     1          (E15.4,E15.4),')'/
     2        ' This should trigger an error message'//)
 1065 FORMAT(' LOG returned the value (',(E15.4,E15.4),')'//)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CLOG Test Program ----------
      END
CC    COMPLEX FUNCTION FLOG(CX)
C----------------------------------------------------------------------
C  Calibration function for clog tests.
C
C  Latest revision - October 31, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
CC    COMPLEX CX
CC    COMPLEX*16 DCX
C----------------------------------------------------------------------
CC    DCX = CX
CC    FLOG = LOG(DCX)
CC    RETURN
C---------- Last line of FLOG ----------
CC    END
