C      ALGORITHM 714, COLLECTED ALGORITHMS FROM ACM.
C      THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
C      VOL. 19, NO. 1, MARCH, 1993, PP. 1-21.
C----------------------------------------------------------------------
C  Program to test CABS
C
C  Accuracy tests are based on the identities
C
C            cabs((3x,4x)) = 5x
C     and
C            cabs((5x,12x)) = 13x
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 six
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  Generic Fortran subprograms required:
C
C     ABS, LOG, MAX, 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 7, 1990
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,NDUM,NEGEP,NGRD
CS    REAL
CD    DOUBLE PRECISION
     1   A,AIT,ALL9,ALBETA,B,BETA,CONV,DEL,EIGHT,EPS,EPSNEG,FIVE,FOUR,
     2   ONE,REN,R6,R7,SIXTEN,THREE,TWELVE,TWENTY,W,X,XL,XMAX,XMIN,XN,
     3   X1,Y,Z,ZERO
CS    COMPLEX
CD    COMPLEX*16
     1   CONVC,CX
CC   2   ,FABS
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ONE,ZERO,TWENTY,ALL9/1.0E0,0.0E0,2.0E1,-999.0E0/
CS    DATA THREE,FOUR,FIVE,TWELVE/3.0E0,4.0E0,5.0E0,12.0E0/
CS    DATA EIGHT,SIXTEN/8.0E0,16.0E0/
CD    DATA ONE,ZERO,TWENTY,ALL9/1.0D0,0.0D0,2.0D1,-999.0D0/
CD    DATA THREE,FOUR,FIVE,TWELVE/3.0D0,4.0D0,5.0D0,12.0D0/
CD    DATA EIGHT,SIXTEN/8.0D0,16.0D0/
C----------------------------------------------------------------------
C  Statement functions for conversion between integer and float
C----------------------------------------------------------------------
CS    CONV(NDUM) = REAL(NDUM)
CS    CONVC(X,Y) = CMPLX(X,Y)
CD    CONV(NDUM) = DBLE(NDUM)
CD    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)
      A = ONE
      B = TWENTY
      N = 2000
      XN = CONV(N)
      DEL = (B - A) / XN
      I1 = 0
C-----------------------------------------------------------------
C  Random argument accuracy tests
C-----------------------------------------------------------------
      DO 300 J = 1, 2
         K1 = 0
         K3 = 0
         X1 = ZERO
         R6 = ZERO
         R7 = ZERO
         XL = A
         DO 200 I = 1, N
            X = DEL * REN(I1) + XL
            IF (J .EQ. 1) THEN
                  Y = X * EIGHT
                  Z = X + Y
                  X = Z - Y
                  Y = X * FIVE
                  CX = CONVC(THREE*X,FOUR*X)
               ELSE
                  Y = X * (SIXTEN + SIXTEN)
                  Z = X + Y
                  X = Z - Y
                  Y = X * TWELVE + X
                  CX = CONVC(FIVE*X,TWELVE*X)
            END IF
            Z = ABS(CX)
CC          Z = FABS(CX)
            W = (Z - Y) / Z
C-----------------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------------
            IF (W .GT. ZERO) THEN
                  K1 = K1 + 1
               ELSE IF (W .LT. ZERO) THEN
                  K3 = K3 + 1
            END IF
            W = ABS(W)
            IF (W .GT. R6) THEN
               R6 = W
               X1 = X
            END IF
            R7 = R7 + W * W
            XL = XL + DEL
  200    CONTINUE
C-----------------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------------
         K2 = N - K1 - K3
         R7 = SQRT(R7/XN)
         IF (J .EQ. 1) THEN
               WRITE (IOUT,1000)
            ELSE
               WRITE (IOUT,1005)
         END IF
         WRITE (IOUT,1010) N,A,B
         WRITE (IOUT,1011) K1,K2,K3
         WRITE (IOUT,1020) IT,IBETA
         IF (R6 .NE. ZERO) THEN
               W = LOG(ABS(R6))/ALBETA
            ELSE
               W = ALL9
         END IF
         WRITE (IOUT,1021) R6,IBETA,W,X1
         W = MAX(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
         IF (R7 .NE. ZERO) THEN
               W = LOG(ABS(R7))/ALBETA
            ELSE
               W = ALL9
         END IF
         WRITE (IOUT,1023) R7,IBETA,W
         W = MAX(AIT+W,ZERO)
         WRITE (IOUT,1022) IBETA,W
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1040)
      X = XMIN
      CX = CONVC(THREE*X,FOUR*X)
      Y = FIVE*X
      WRITE (IOUT,1041)
      Z = ABS(CX)
      WRITE (IOUT,1042) Z,Y
CC    Z = FABS(CX)
C-----------------------------------------------------------------
C  Test of error returns
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      X = XMAX/SIXTEN
      CX = CONVC(FIVE*X,TWELVE*X)
      WRITE (IOUT,1051) CX
      Z = ABS(CX)
CC    Z = FABS(CX)
      WRITE (IOUT,1055) Z
      X = XMAX/FOUR
      CX = CONVC(THREE*X,XMAX)
      WRITE (IOUT,1052) CX
      Z = ABS(CX)
CC    Z = FABS(CX)
      WRITE (IOUT,1055) Z
      WRITE (IOUT,1100)
      STOP
C-----------------------------------------------------------------
 1000 FORMAT('1Test of CABS((3X,4X)) vs 5X'//)
 1005 FORMAT(//' Test of CABS((5X,12X)) vs 13X'//)
 1010 FORMAT(I7,' Random values of X were drawn from the interval '/
     1 6X,'(',E15.4,',',E15.4,')'//)
 1011 FORMAT(3X,' CABS was larger',I6,' times,'/
     1     12X,' agreed',I6,' times, and'/
     2     8X,'was smaller',I6,' times.'//)
 1020 FORMAT(' There are',I4,' base',I4,
     1    ' significant digits in a floating-point number'//)
 1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **',
     1  F7.2/4X,'occurred for X =',E17.6)
 1022 FORMAT(' The estimated loss of base',I4,
     1  ' significant digits is',F7.2//)
 1023 FORMAT(' The root mean square relative error was',E15.4,
     1    ' = ',I4,' **',F7.2)
 1040 FORMAT(//' Test of special arguments'//)
 1041 FORMAT(' Test of |(3*XMIN,4*XMIN)| vs 5*XMIN.'/
     1       ' This should not trigger an error message'//)
 1042 FORMAT(' CABS((3*XMIN,4*XMIN)) = ',E15.7/
     2       ' 5*XMIN                = ',E15.7//)
 1050 FORMAT(' Test of error returns'//)
 1051 FORMAT(' CABS will be called with the argument',(E15.4,E15.4)/
     1       ' This should not trigger an error message'//)
 1052 FORMAT(' CABS will be called with the argument',(E15.4,E15.4)/
     1       ' This should trigger an error message'//)
 1055 FORMAT(' CABS returned the value ',E15.4///)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CABS Test Program ----------
      END
 
CC    COMPLEX FUNCTION FABS(CX)
C----------------------------------------------------------------------
C  Calibration function for cabs 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    FABS = ABS(DCX)
CC    RETURN
C---------- Last line of FABS ----------
CC    END
 
C----------------------------------------------------------------------
C  Program to test CEXP
C
C  Accuracy tests are based on the identity
C
C            exp(z) = exp(z-w)*exp(w)
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         AIMAG (or DIMAG), AINT, CMPLX (or DCMPLX), EXP,
C
C         LOG, REAL (or DBLE), SQRT
C
C
C  Latest revision - October 31, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    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
CD    DOUBLE PRECISION
     1   BETA,CONV,CONVI,CONVR,DELX,DELY,EPS,EPSNEG,
     2   HALF,ONE,ONEP6,ONE16,REN,SIXTEN,TEN,THREE,
     3   W,X,XL,XMAX,XMIN,Y,ZERO
CS    COMPLEX
CD    COMPLEX*16
     1   CDEL,CONVC,CY,EXPDEL,ZDUM
CC   2   ,FEXP
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,ONE16,HALF,ONE/0.0E0,0.0625E0,0.5E0,1.0E0/
CS    DATA ONEP6,THREE,TEN,SIXTEN/1.625E0,3.0E0,10.0E0,16.0E0/
CS    DATA EXPDEL/(6.2416044877018563681E-2,6.6487597751003112768E-2)/
CS    DATA CDEL/(0.0625E0,0.0625E0)/
CD    DATA ZERO,ONE16,HALF,ONE/0.0D0,0.0625D0,0.5D0,1.0D0/
CD    DATA ONEP6,THREE,TEN,SIXTEN/1.625D0,3.0D0,10.0D0,16.0D0/
CD    DATA EXPDEL/(6.2416044877018563681D-2,6.6487597751003112768D-2)/
CD    DATA CDEL/(0.0625D0,0.0625D0)/
C----------------------------------------------------------------------
C  Statement functions for conversion between integer or complex
C  and float
C----------------------------------------------------------------------
CS    CONV(NDUM) = REAL(NDUM)
CS    CONVR(ZDUM) = REAL(ZDUM)
CS    CONVI(ZDUM) = AIMAG(ZDUM)
CS    CONVC(X,Y) = CMPLX(X,Y)
CD    CONV(NDUM) = DBLE(NDUM)
CD    CONVR(ZDUM) = DBLE(ZDUM)
CD    CONVI(ZDUM) = DIMAG(ZDUM)
CD    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)
      AX = ZERO
      BX = ONE
      AY = ONE16
      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  Carefully purify arguments and evaluate identities
C-----------------------------------------------------------------
            CX = CONVC(X,Y)
            CY = CX - CDEL
            CZ = EXP(CX)
            CZZ = EXP(CY)
CC          CZ = FEXP(CX)
CC          CZZ = FEXP(CY)
            CZZ = CZZ + CZZ*EXPDEL
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 = ONEP6
               BX = THREE
            ELSE
               AX = SIXTEN
               BX = AX + ONE
         END IF
         AY = AX
         BY = BX
  300 CONTINUE
C---------------------------------------------------------------------
C  Special tests
C---------------------------------------------------------------------
      WRITE (IOUT,1025)
      WRITE (IOUT,1030)
      DO 320 I = 1, 5
         X = REN(I1) * BETA
         Y = REN(I1) * BETA
         CX = CONVC(X,Y)
         CY = -CX
         CZ = EXP(CX) * EXP(CY) - ONE
         WRITE (IOUT,1070) CX, CZ
  320 CONTINUE
      WRITE (IOUT,1040)
      CX = CONVC(ZERO,ZERO)
      CZ = EXP(CX) - ONE
      WRITE (IOUT,1041) CZ
      WRITE (IOUT,1045)
      CX = CONVC(ZERO,TEN)
      CZ = EXP(CX)
      X = COS(TEN)
      Y = SIN(TEN)
      W = (CONVR(CZ)-X)/X
      WRITE (IOUT,1046) W
      W = (CONVI(CZ)-Y)/Y
      WRITE (IOUT,1047) W
      CX = CONVC(TEN,ZERO)
      CZ = EXP(CX)
      X = EXP(TEN)
      W = (CONVR(CZ)-X)/X
      WRITE (IOUT,1048) W
      WRITE (IOUT,1050)
      WRITE (IOUT,1052)
      X = AINT(LOG(XMIN))
      CX = CONVC(X,ONE)
      CZ = EXP(CX)
      WRITE (IOUT,1056) CX,CZ
      Y = SQRT(XMAX)
      X = AINT(LOG(XMAX))
      CX = CONVC(X,ONE)
      CZ = EXP(CX)
      WRITE (IOUT,1056) CX,CZ
      CX = CONVC (X*HALF,TEN)
      CY = CX * HALF
      CZ = EXP(CX)
      CZZ = EXP (CY)
      CZZ = CZZ * CZZ
      WRITE (IOUT,1057) CX, CZ, CY, CZZ
C---------------------------------------------------------------------
C  Test of error returns
C---------------------------------------------------------------------
      WRITE (IOUT,1060)
      X = -ONE / SQRT(XMIN)
      CX = CONVC(X,ONE)
      WRITE (IOUT,1061) CX
      WRITE (IOUT,1062)
      CZ = EXP(CX)
      WRITE (IOUT,1065) CZ
      Y = SQRT(XMAX)
      CX = CONVC(TEN,Y)
      WRITE (IOUT,1061) CX
      WRITE (IOUT,1063)
      CZ = EXP(CX)
      WRITE (IOUT,1065) CZ
      CX = CONVC(-X,-ONE)
      WRITE (IOUT,1061) CX
      WRITE (IOUT,1064)
      CZ = EXP(CX)
      WRITE (IOUT,1065) CZ
      WRITE (IOUT,1100)
      STOP
C---------------------------------------------------------------------
 1000 FORMAT(' Test of EXP(Z) vs EXP(Z-DEL)EXP(DEL), DEL = '
     1       '(1/16,1/16)'//)
 1025 FORMAT(' Special tests'//)
 1030 FORMAT(' The identity  EXP(Z)*EXP(-Z) = 1.0  will be tested.'//
     1      8X,'Z',24X,'F(Z)*F(-Z) - 1'/)
 1040 FORMAT(//' Test of special arguments'//)
 1041 FORMAT(' EXP(0.0,0.0) - 1.0E0 = (',(E15.7,E15.7),')'/)
 1045 FORMAT(//' Test cos(y), sin(y), and exp(x)'//)
 1046 FORMAT(' [REAL(EXP(0,10)) - COS(10)]/COS(10) = ',E15.7/)
 1047 FORMAT(' [IMAG(EXP(0,10)) - SIN(10)]/SIN(10) = ',E15.7/)
 1048 FORMAT(' [REAL(EXP(10,0)) - EXP(10)]/EXP(10) = ',E15.7/)
 1050 FORMAT(//' Tests near extreme arguments'//)
 1052 FORMAT(' The real part of the first result below may underflow.'/)
 1056 FORMAT(' EXP(',(E13.6,E13.6),') = (',(E13.6,E13.6),')'/)
 1057 FORMAT(/' If EXP(',(E13.6,E13.6),') = (',(E13.6,E13.6),
     1 ') is not about'/' EXP(',(E13.6,E13.6),')**2 = (',
     2 (E13.6,E13.6),'), there is an'/' argument reduction error'//)
 1060 FORMAT(' Test of error returns'/)
 1061 FORMAT(/' EXP will be called with the argument (',
     1          (E15.4,E15.4),')'/)
 1062 FORMAT(' The real component of the argument is so negative'/
     1        ' that the real exponential should underflow to zero.'/
     2        ' Thus, the result should be (0,0).'/)
 1063 FORMAT(' The imaginary component of the argument is too large'/
     1        ' for the sin and cos computations to make sense.'/)
 1064 FORMAT(' The real component of the argument is so large'/
     1        ' that the real exponential should overflow.  The'/
     2        ' result should be (inf,inf) with an error return.'/)
 1065 FORMAT(' EXP returned the value (',(E15.4,E15.4),')'//)
 1070 FORMAT(2(E15.7,E15.7)/)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CEXP Test Program ----------
      END
 
CC    COMPLEX FUNCTION FEXP(CX)
C----------------------------------------------------------------------
C  Calibration function for cexp tests.
C
C  Latest revision - July 12, 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    FEXP = EXP(DCX)
CC    RETURN
C---------- Last line of FEXP ----------
CC    END
C----------------------------------------------------------------------
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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    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
CD    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
CD    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/
CD    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)
CD    CONV(NDUM) = DBLE(NDUM)
CD    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
 
 
C----------------------------------------------------------------------
C  Program to test complex power (exponentiation)
C
C  Accuracy tests are based on the identity
C
C            z = z ** (1,0)
C
C            z * z = Z ** (2,0)
C     and
C            z ** w = (z*z) ** (w/2)
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              MAXEXP - the smallest integer power of FLOAT(IBETA)
C                       that causes overflow;
C              XMIN   - the smallest non-vanishing floating-point
C                       power of the radix.
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         ABS, CMPLX (or DCMPLX), LOG, REAL (or DBLE)
C
C  Latest revision - November 29, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    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
CD    DOUBLE PRECISION
     1   BETA,CONV,CONVR,C1,C2,DELX,DELY,EPS,EPSNEG,FOUR,HALF,
     2   ONE,REN,SCALE,TEN,W,X,XL,XMAX,XMIN,Y,ZERO
CS    COMPLEX
CD    COMPLEX*16
     1   CONE,CONVC,CTWO,CY
CC   2   ,FPOW
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,HALF,ONE,FOUR,TEN/0.0E0,0.5E0,1.0E0,4.0E0,10.0E0/
CS    DATA CONE,CTWO/(1.0E0,0.0E0),(2.0E0,0.0E0)/
CS    DATA C1,C2/0.95E0,1.05E0/
CD    DATA ZERO,HALF,ONE,FOUR,TEN/0.0D0,0.5D0,1.0D0,4.0D0,10.0D0/
CD    DATA CONE,CTWO/(1.0D0,0.0D0),(2.0D0,0.0D0)/
CD    DATA C1,C2/0.95D0,1.05D0/
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)
CS    CONVR(CY) = REAL(CY)
CD    CONV(NDUM) = DBLE(NDUM)
CD    CONVC(X,Y) = DCMPLX(X,Y)
CD    CONVR(CY) = DBLE(CY)
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
      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  Carefully purify arguments and evaluate identities
C-----------------------------------------------------------------
            CX = CONVC(X,Y)
            IF (J .EQ. 1) THEN
                  CX = CONVC(X,Y)
                  CZ = CX ** CONE
CC                CZ = FPOW(CX,CONE)
                  CZZ = CX
               ELSE IF (J .EQ. 2) THEN
  100             W = ABS(Y/X)
                  IF ((W .GT. C1) .AND. (W .LT. C2)) THEN
                     Y = DELY * REN(I1) + AY
                     GO TO 100
                  END IF
                  CX = CONVC(X,Y)
                  CZ = CX * SCALE
                  CY = CX + CZ
                  CX = CY - CZ
                  CZ = CX ** CTWO
CC                CZ = FPOW(CX,CTWO)
                  CZZ = CX * CX
               ELSE IF (J .EQ. 3) THEN
                  CX = CONVC(X,Y)
                  X = DELY * REN(I1) + AX
                  Y = DELY * REN(I1) + AY
                  CZ = CX * SCALE
                  CY = CX + CZ
                  CX = CY - CZ
                  CY = CONVC(X,Y)
                  CV = CY
                  CZ = CX ** CY
                  CZZ = (CX*CX) ** (CY*HALF)
CC                CZ = FPOW(CX,CY)
CC                CZZ = FPOW(CX*CX,CY*HALF)
            END IF
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
            CALL TABLAT
            XL = XL + DELX
  200    CONTINUE
C-----------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------
         IF (J .EQ. 1) THEN
               WRITE (IOUT,1000)
            ELSE IF (J .EQ. 2) THEN
               WRITE (IOUT,1001)
            ELSE
               WRITE (IOUT,1002)
               WRITE (IOUT,1003)
         END IF
         CALL REPORT(IOUT)
         IF (J .EQ. 2) THEN
            PFLAG = .TRUE.
            AX = FOUR
            AY = AX
         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)
         X = REN(I1) * TEN
         Y = REN(I1) * TEN
         CZ = CONVC(X,Y)
         CY = ONE / CX
         CZZ = (CX**CZ - CY**(-CZ))/ABS(CX**CZ)
CC       CZZ = (FPOW(CX,CZ) - FPOW(CY,-CZ))/ABS(FPOW(CX,CZ))
         WRITE (IOUT,1031) CX,CZ,CZZ
  320 CONTINUE
      WRITE (IOUT,1040)
      Y = BETA * XMIN
      CX = CONVC(BETA,Y)
      W = CONV(MAXEXP-1)
      CY = CONVC(W,ZERO)
      CZ = CX ** CY
CC    CZ = FPOW(CX,CY)
      X = BETA
      DO 330 I = 3, MAXEXP
         X = X * BETA
  330 CONTINUE
      CZZ = X * CONVC(ONE,XMIN*W)
      WRITE (IOUT,1041) CZ,CZZ
      WRITE (IOUT,1042) CZ-CZZ
C-----------------------------------------------------------------
C  Test of error returns
C-----------------------------------------------------------------
      WRITE (IOUT,1060)
      WRITE (IOUT,1061)
      CX = CONVC(ZERO,ZERO)
      CZ = CX ** CX
      WRITE (IOUT,1065) CZ
      WRITE (IOUT,1100)
      STOP
C---------------------------------------------------------------------
 1000 FORMAT(' Test of Z ** (1,0) vs Z'//)
 1001 FORMAT(' Test of Z ** (2,0) vs Z*Z'//)
 1002 FORMAT(' Test of Z**W vs (Z*Z) ** (W/2)'//)
 1003 FORMAT(' *****************************************************'/
     1       ' * Very large MRE values are the norm for this test; *'/
     2       ' * the frequency of exact agreement is meaningful.   *'/
     3       ' *****************************************************'//)
 1025 FORMAT(' Special tests'//)
 1030 FORMAT(' The identity  Z**W = (1/Z)**(-W)  will be tested.'//
     1      11X,'Z',21X,'W',10X,'(Z**W - (1/Z)**(-W))/CABS(Z**W)'/)
 1031 FORMAT(2(E11.4,E11.4),(E15.7,E15.7)/)
 1040 FORMAT(//' Test near an extreme argument.  This should not,'
     1         ' but may'/' trigger an error message'//)
 1041 FORMAT(' (BETA,BETA*XMIN)**(MAXEXP-1,0)       = ',
     1       '(',(E15.7,E15.7),')'/
     2       ' BETA**(MAXEXP-1)*(1,(MAXEXP-1)*XMIN) = ',
     3       '(',(E15.7,E15.7),')')
 1042 FORMAT(' Difference ',26X,'= ','(',(E15.7,E15.7),')'//)
 1060 FORMAT(' Test of error returns'/)
 1061 FORMAT(/' Test of (0,0) ** (0,0).'/
     1        ' This should trigger an error message'//)
 1065 FORMAT(' Test returned the value (',(E15.4,E15.4),')'//)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CPOW Test Program ----------
      END
 
CC    COMPLEX FUNCTION FPOW(CX,CY)
C----------------------------------------------------------------------
C  Calibration function for complex power tests.
C
C  Latest revision - October 31, 1990
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C----------------------------------------------------------------------
CC    COMPLEX CX,CY
CC    COMPLEX*16 DCX,DCY
C----------------------------------------------------------------------
CC    DCX = CX
CC    DCY = CY
CC    FPOW = DCX ** DCY
CC    RETURN
C---------- Last line of FPOW ----------
CC    END
 
C---------------------------------------------------------------------
C  Program to test CSIN/CCOS
C
C  Accuracy tests are based on the identities
C
C            sin(z) = sin(z-w)*cos(w) + cos(z-w)*sin(w)
C     and
C            cos(z) = cos(z-w)*cos(w) - sin(z-w)*sin(w)
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 five
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              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         ABS, AIMAG (or DIMAG), CMPLX (or DCMPLX), COS,
C
C         LOG, MAX, REAL (or DBLE), SIN
C
C
C  Latest revision - October 31, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    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
CD    DOUBLE PRECISION
     1   ALXMAX,BETA,BETAP,C,CONV,CONVR,CONVI,DELX,DELY,ELEVEN,EPS,
     2   EPSNEG,HALF,ONE,REN,SIXTEN,TEN,W,X,XL,XMAX,XMIN,Y,ZERO
CS    COMPLEX
CD    COMPLEX*16
     1   CC,CONVC,CDEL,SINDEL,COSDEL,CS,CY,SIN11A,SIN11B,SIN22A,
     2   SIN22B,ZDUM
CC   3   ,FSIN,FCOS
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,HALF,ONE/0.0E0,0.5E0,1.0E0/
CS    DATA TEN,ELEVEN,SIXTEN/10.0E0,11.0E0,16.0E0/
CS    DATA SINDEL/(6.2581348413276935585E-2,6.2418588008436587236E-2)/
CS    DATA COSDEL/(-2.5431314180235545803E-6,-3.9062493377261771826E-3)/
CS    DATA SIN11A/(-134167.0E0,+593.0E0)/
CS    DATA SIN22A/(-1187.0E0,-134163.0E0)/
CS    DATA SIN11B/(-3.2928849557976510994E-1,7.8989452897413630729E-1)/
CS    DATA SIN22B/(-5.6815858848429427210E-1,-3.8738909081835398666E-1)/
CS    DATA CDEL/(0.0625,0.0625)/
CD    DATA ZERO,HALF,ONE/0.0D0,0.5D0,1.0D0/
CD    DATA TEN,ELEVEN,SIXTEN/10.0D0,11.0D0,16.0D0/
CD    DATA SINDEL/(6.2581348413276935585D-2,6.2418588008436587236D-2)/
CD    DATA COSDEL/(-2.5431314180235545803D-6,-3.9062493377261771826D-3)/
CD    DATA SIN11A/(-134167.0D0,+593.0D0)/
CD    DATA SIN22A/(-1187.0D0,-134163.0D0)/
CD    DATA SIN11B/(-3.2928849557976510994D-1,7.8989452897413630729D-1)/
CD    DATA SIN22B/(-5.6815858848429427210D-1,-3.8738909081835398666D-1)/
CD    DATA CDEL/(0.0625,0.0625)/
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)
CS    CONVR(ZDUM) = REAL(ZDUM)
CS    CONVI(ZDUM) = AIMAG(ZDUM)
CD    CONV(NDUM) = DBLE(NDUM)
CD    CONVC(X,Y) = DCMPLX(X,Y)
CD    CONVR(ZDUM) = DBLE(ZDUM)
CD    CONVI(ZDUM) = DIMAG(ZDUM)
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)
      AX = ONE/SIXTEN
      BX = ONE
      AY = ZERO
      BY = ONE
      N = 2000
      XN = CONV(N)
      I1 = 0
      IFLAG = .FALSE.
      RFLAG = .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
            CX = CONVC(X,Y)
C----------------------------------------------------------------
C  Carefully purify arguments and evaluate identities
C----------------------------------------------------------------
            CY = CX - CDEL
            CS = SIN(CY)
            CC = COS(CY)
CC          CS = FSIN(CY)
CC          CC = FCOS(CY)
            IF (J .NE. 3) THEN
                  CZ = SIN(CX)
CC                CZ = FSIN(CX)
                  CZZ = CS + (CS*COSDEL+CC*SINDEL)
               ELSE
                  CZ = COS(CX)
CC                CZ = FCOS(CX)
                  CZZ = CC + (CC*COSDEL-CS*SINDEL)
            END IF
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
            CALL TABLAT
            XL = XL + DELX
  200    CONTINUE
C---------------------------------------------------------------------
C  Process and output statistics
C---------------------------------------------------------------------
         IF (J .NE. 3) THEN
               WRITE (IOUT,1000)
            ELSE
               WRITE (IOUT,1005)
         END IF
         CALL REPORT(IOUT)
         IF (J .EQ. 1) THEN
               AX = SIXTEN
               BX = AX+ONE
         END IF
         AY = AX
         BY = BX
         IFLAG = .TRUE.
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1025)
      C = ONE / BETA ** (IT/2)
      X = ELEVEN
      Y = X + ONE + HALF
      X = X + X
      CX = CONVC(X+C,Y)
      CY = CONVC(X-C,Y)
      CZ = (SIN(CX) - SIN(CY)) / (C + C)
      WRITE (IOUT,1026) CZ
      WRITE (IOUT,1030)
      DO 320 I = 1, 5
         X = REN(I1) * TEN
         Y = REN(I1) * TEN
         CX = CONVC(X,Y)
         CZ = SIN(CX) + SIN(-CX)
         WRITE (IOUT,1060) CX, CZ
  320 CONTINUE
      WRITE (IOUT,1040)
      X = ELEVEN
      Y = X + ONE + HALF
      CX = CONVC(X,Y)
      CZ = SIN (CX)
      CS = (CZ-SIN11A)-SIN11B
      C = CONVI(CS)/CONVI(CZ)
      IF (C .NE. ZERO) THEN
            W = LOG(ABS(C))/ALBETA+AIT
         ELSE
            W = ZERO
      END IF
      W = MAX(W,ZERO)
      WRITE (IOUT,1042) CX,CZ,IBETA,W
      IF (W .LT. ONE) THEN
            WRITE (IOUT,1046)
         ELSE
            WRITE (IOUT,1048)
      END IF
      X = X + X
      CX = CONVC(X,Y)
      CZ = SIN (CX)
      CS = (CZ-SIN22A)-SIN22B
      C = CONVR(CS)/CONVR(CZ)
      IF (C .NE. ZERO) THEN
            W = LOG(ABS(C))/ALBETA+AIT
         ELSE
            W = ZERO
      END IF
      W = MAX(W,ZERO)
      WRITE (IOUT,1044) CX,CZ,IBETA,W
      IF (W .LT. ONE) THEN
            WRITE (IOUT,1046)
         ELSE
            WRITE (IOUT,1048)
      END IF
C-----------------------------------------------------------------
C  Test of error returns
C-----------------------------------------------------------------
      WRITE (IOUT,1050)
      ALXMAX = LOG(XMAX)-HALF
      CX = CONVC(ONE,ALXMAX)
      WRITE (IOUT,1053) CX
      CY = SIN(CX)
      WRITE (IOUT,1055) CY
      BETAP = BETA ** IT
      CX = CONVC(BETAP,ONE)
      WRITE (IOUT,1052) CX
      CY = SIN(CX)
      WRITE (IOUT,1055) CY
      CX = CONVC(ONE,BETAP)
      WRITE (IOUT,1052) CX
      CY = SIN(CX)
      WRITE (IOUT,1055) CY
      WRITE (IOUT,1100)
      STOP
C----------------------------------------------------------------
 1000 FORMAT(' Test of SIN(X) vs SIN(W+D), W = X-D, D = (1/16,1/16)'//)
 1005 FORMAT(' Test of COS(X) vs COS(W+D), W = X-D, D = (1/16,1/16)'//)
 1025 FORMAT(' Special Tests'//)
 1026 FORMAT(' If (',(F9.1,F7.1),') is not almost (-134163,1188),'/
     1   14X,25HSIN has the wrong period.  //)
 1030 FORMAT(' The identity   SIN(-Z) = -SIN(Z)   will be tested.'//
     1      16X,'X',23X,'F(Z) + F(-Z)'/)
 1040 FORMAT(/' Test the significance in the real or imaginary'
     1  ' components near'/' the zeros of the real sin'
     2  ' or cos.'/)
 1042 FORMAT(' SIN(',(E15.7,E15.7),') = (',(E15.7,E15.7),')'/
     1  ' The estimated loss of base ',I4,' digits in the'
     2  ' imaginary component is ',F7.2)
 1044 FORMAT(' SIN(',(E15.7,E15.7),') = (',(E15.7,E15.7),')'/
     1  ' The estimated loss of base ',I4,' digits in the'
     2  ' real component is ',F7.2)
 1046 FORMAT(' This is acceptable'/)
 1048 FORMAT(' This is too large'/)
 1050 FORMAT(' Test of error returns'//)
 1052 FORMAT(' SIN will be called with the argument (',(E15.4,E15.4),
     1       ')'/' This should trigger an error message'//)
 1053 FORMAT(' SIN will be called with the argument (',(E15.4,E15.4),
     1       ')'/' This should not trigger an error message'//)
 1055 FORMAT(' SIN returned the value (',(E15.4,E15.4),')'///)
 1060 FORMAT(2(E15.7,E15.7)/)
 1100 FORMAT(' This concludes the tests')
C---------- Last card of CSIN/CCOS Test Program ----------
      END
 
CC    COMPLEX FUNCTION FSIN(CX)
C----------------------------------------------------------------------
C  Calibration function for csin 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    FSIN = SIN(DCX)
CC    RETURN
C---------- Last line of FSIN ----------
CC    END
 
CC    COMPLEX FUNCTION FCOS(CX)
C----------------------------------------------------------------------
C  Calibration function for ccos 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    FCOS = COS(DCX)
CC    RETURN
C---------- Last line of FCOS ----------
CC    END
 
C----------------------------------------------------------------------
C  Program to test CSQRT
C
C  Accuracy tests are based on the identity
C
C            z = sqrt(z*z)
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         ABS, 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 15, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    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
CD    DOUBLE PRECISION
     1   BETA,CONV,DELX,DELY,EPS,EPSNEG,ONE,REN,SCALE,TEN,
     2   W,X,XL,XMAX,XMIN,Y,Z,ZERO
CS    COMPLEX
CD    COMPLEX*16
     1   CONVC,CY
CC   2   ,FSQRT
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
      DATA IOUT/6/
CS    DATA ZERO,ONE,TEN/0.0E0,1.0E0,10.0E0/
CD    DATA ZERO,ONE,TEN/0.0D0,1.0D0,10.0D0/
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)
CD    CONV(NDUM) = DBLE(NDUM)
CD    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 = ZERO
      BX = TEN
      AY = AX
      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, 2
         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)
C-----------------------------------------------------------------
C  Carefully purify arguments 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)
            CZZ = CX
            IF (J .EQ. 1) THEN
                  CZ = SQRT(CY)
CC                CZ = FSQRT(CY)
               ELSE
                  CZ = -SQRT(CY)
CC                CZ = -FSQRT(CY)
            END IF
C-----------------------------------------------------------------
C  Accumulate results
C-----------------------------------------------------------------
            CALL TABLAT
            XL = XL + DELX
  200    CONTINUE
C-----------------------------------------------------------------
C  Process and output statistics
C-----------------------------------------------------------------
         IF (J .EQ. 1) THEN
               WRITE (IOUT,1000)
            ELSE
               WRITE (IOUT,1001)
         END IF
         CALL REPORT(IOUT)
         BX = -BX * BX
         BY = -BX
  300 CONTINUE
C-----------------------------------------------------------------
C  Special tests
C-----------------------------------------------------------------
      WRITE (IOUT,1040)
      CX = CONVC(ONE,ZERO)
      CZ = SQRT(CX)
      WRITE (IOUT,1041) CZ
      WRITE (IOUT,1050)
      CX = CONVC(XMIN,XMIN)
      WRITE (IOUT,1061) CX
      CZ = SQRT(CX)
      WRITE (IOUT,1051) CX,CZ
      CX = CONVC(XMAX,XMAX)
      WRITE (IOUT,1061) CX
      CZ = SQRT(CX)
      WRITE (IOUT,1052) CX,CZ
      WRITE (IOUT,1100)
      STOP
C---------------------------------------------------------------------
 1000 FORMAT(' Test of SQRT(Z*Z) vs Z'//)
 1001 FORMAT(' Test of -SQRT(Z*Z) vs Z'//)
 1040 FORMAT(//' Test of special arguments'//)
 1041 FORMAT(' SQRT(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(' SQRT(XMIN,XMIN) = SQRT(',(E15.7,E15.7),') ='/40X,'(',
     1         (E15.7,E15.7),')'/)
 1052 FORMAT(' SQRT(XMAX,XMAX) = SQRT(',(E15.7,E15.7),') ='/40X,'(',
     1         (E15.7,E15.7),')'/)
 1061 FORMAT(/' SQRT will be called with the argument (',
     1          (E15.4,E15.4),')'//)
 1100 FORMAT(' This concludes the tests')
C--------- Last card of CSQRT Test Program ----------
      END
 
CC    COMPLEX FUNCTION FSQRT(CX)
C----------------------------------------------------------------------
C  Calibration function for csqrt tests.
C
C  Latest revision - October 10, 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    FSQRT = SQRT(DCX)
CC    RETURN
C---------- Last line of FSQRT ----------
CC    END
 
      SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,
     1                   MAXEXP,EPS,EPSNEG,XMIN,XMAX)
C-----------------------------------------------------------------------
C  This Fortran 77 subroutine is intended to determine the parameters
C   of the floating-point arithmetic system specified below.  The
C   determination of the first three uses an extension of an algorithm
C   due to M. Malcolm, CACM 15 (1972), pp. 949-951, incorporating some,
C   but not all, of the improvements suggested by M. Gentleman and S.
C   Marovich, CACM 17 (1974), pp. 276-277.  An earlier version of this
C   program was published in the book Software Manual for the
C   Elementary Functions by W. J. Cody and W. Waite, Prentice-Hall,
C   Englewood Cliffs, NJ, 1980.
C
C  The program as given here must be modified before compiling.  If
C   a single (double) precision version is desired, change all
C   occurrences of CS (CD) in columns 1 and 2 to blanks.
C
C  Parameter values reported are as follows:
C
C       IBETA   - the radix for the floating-point representation
C       IT      - the number of base IBETA digits in the floating-point
C                 significand
C       IRND    - 0 if floating-point addition chops
C                 1 if floating-point addition rounds, but not in the
C                   IEEE style
C                 2 if floating-point addition rounds in the IEEE style
C                 3 if floating-point addition chops, and there is
C                   partial underflow
C                 4 if floating-point addition rounds, but not in the
C                   IEEE style, and there is partial underflow
C                 5 if floating-point addition rounds in the IEEE style,
C                   and there is partial underflow
C       NGRD    - the number of guard digits for multiplication with
C                 truncating arithmetic.  It is
C                 0 if floating-point arithmetic rounds, or if it
C                   truncates and only  IT  base  IBETA digits
C                   participate in the post-normalization shift of the
C                   floating-point significand in multiplication;
C                 1 if floating-point arithmetic truncates and more
C                   than  IT  base  IBETA  digits participate in the
C                   post-normalization shift of the floating-point
C                   significand in multiplication.
C       MACHEP  - the largest negative integer such that
C                 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, except that
C                 MACHEP is bounded below by  -(IT+3)
C       NEGEPS  - the largest negative integer such that
C                 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, except that
C                 NEGEPS is bounded below by  -(IT+3)
C       IEXP    - the number of bits (decimal places if IBETA = 10)
C                 reserved for the representation of the exponent
C                 (including the bias or sign) of a floating-point
C                 number
C       MINEXP  - the largest in magnitude negative integer such that
C                 FLOAT(IBETA)**MINEXP is positive and normalized
C       MAXEXP  - the smallest positive power of  BETA  that overflows
C       EPS     - the smallest positive floating-point number such
C                 that  1.0+EPS .NE. 1.0. In particular, if either
C                 IBETA = 2  or  IRND = 0, EPS = FLOAT(IBETA)**MACHEP.
C                 Otherwise,  EPS = (FLOAT(IBETA)**MACHEP)/2
C       EPSNEG  - A small positive floating-point number such that
C                 1.0-EPSNEG .NE. 1.0. In particular, if IBETA = 2
C                 or  IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS.
C                 Otherwise,  EPSNEG = (IBETA**NEGEPS)/2.  Because
C                 NEGEPS is bounded below by -(IT+3), EPSNEG may not
C                 be the smallest number that can alter 1.0 by
C                 subtraction.
C       XMIN    - the smallest non-vanishing normalized floating-point
C                 power of the radix, i.e.,  XMIN = FLOAT(IBETA)**MINEXP
C       XMAX    - the largest finite floating-point number.  In
C                 particular  XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP
C                 Note - on some machines  XMAX  will be only the
C                 second, or perhaps third, largest number, being
C                 too small by 1 or 2 units in the last digit of
C                 the significand.
C
C     Latest revision - December 4, 1987
C
C     Author - W. J. Cody
C              Argonne National Laboratory
C
C-----------------------------------------------------------------------
      INTEGER I,IBETA,IEXP,IRND,IT,ITEMP,IZ,J,K,MACHEP,MAXEXP,
     1        MINEXP,MX,NEGEP,NGRD,NXRES
CS    REAL
CD    DOUBLE PRECISION
     1   A,B,BETA,BETAIN,BETAH,CONV,EPS,EPSNEG,ONE,T,TEMP,TEMPA,
     2   TEMP1,TWO,XMAX,XMIN,Y,Z,ZERO
C-----------------------------------------------------------------------
CS    CONV(I) = REAL(I)
CD    CONV(I) = DBLE(I)
      ONE = CONV(1)
      TWO = ONE + ONE
      ZERO = ONE - ONE
C-----------------------------------------------------------------------
C  Determine IBETA, BETA ala Malcolm.
C-----------------------------------------------------------------------
      A = ONE
   10 A = A + A
         TEMP = A+ONE
         TEMP1 = TEMP-A
         IF (TEMP1-ONE .EQ. ZERO) GO TO 10
      B = ONE
   20 B = B + B
         TEMP = A+B
         ITEMP = INT(TEMP-A)
         IF (ITEMP .EQ. 0) GO TO 20
      IBETA = ITEMP
      BETA = CONV(IBETA)
C-----------------------------------------------------------------------
C  Determine IT, IRND.
C-----------------------------------------------------------------------
      IT = 0
      B = ONE
  100 IT = IT + 1
         B = B * BETA
         TEMP = B+ONE
         TEMP1 = TEMP-B
         IF (TEMP1-ONE .EQ. ZERO) GO TO 100
      IRND = 0
      BETAH = BETA / TWO
      TEMP = A+BETAH
      IF (TEMP-A .NE. ZERO) IRND = 1
      TEMPA = A + BETA
      TEMP = TEMPA+BETAH
      IF ((IRND .EQ. 0) .AND. (TEMP-TEMPA .NE. ZERO)) IRND = 2
C-----------------------------------------------------------------------
C  Determine NEGEP, EPSNEG.
C-----------------------------------------------------------------------
      NEGEP = IT + 3
      BETAIN = ONE / BETA
      A = ONE
      DO 200 I = 1, NEGEP
         A = A * BETAIN
  200 CONTINUE
      B = A
  210 TEMP = ONE-A
         IF (TEMP-ONE .NE. ZERO) GO TO 220
         A = A * BETA
         NEGEP = NEGEP - 1
      GO TO 210
  220 NEGEP = -NEGEP
      EPSNEG = A
C-----------------------------------------------------------------------
C  Determine MACHEP, EPS.
C-----------------------------------------------------------------------
      MACHEP = -IT - 3
      A = B
  300 TEMP = ONE+A
         IF (TEMP-ONE .NE. ZERO) GO TO 320
         A = A * BETA
         MACHEP = MACHEP + 1
      GO TO 300
  320 EPS = A
C-----------------------------------------------------------------------
C  Determine NGRD.
C-----------------------------------------------------------------------
      NGRD = 0
      TEMP = ONE+EPS
      IF ((IRND .EQ. 0) .AND. (TEMP*ONE-ONE .NE. ZERO)) NGRD = 1
C-----------------------------------------------------------------------
C  Determine IEXP, MINEXP, XMIN.
C
C  Loop to determine largest I and K = 2**I such that
C         (1/BETA) ** (2**(I))
C  does not underflow.
C  Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------------
      I = 0
      K = 1
      Z = BETAIN
      T = ONE + EPS
      NXRES = 0
  400 Y = Z
         Z = Y * Y
C-----------------------------------------------------------------------
C  Check for underflow here.
C-----------------------------------------------------------------------
         A = Z * ONE
         TEMP = Z * T
         IF ((A+A .EQ. ZERO) .OR. (ABS(Z) .GE. Y)) GO TO 410
         TEMP1 = TEMP * BETAIN
         IF (TEMP1*BETA .EQ. Z) GO TO 410
         I = I + 1
         K = K + K
      GO TO 400
  410 IF (IBETA .EQ. 10) GO TO 420
      IEXP = I + 1
      MX = K + K
      GO TO 450
C-----------------------------------------------------------------------
C  This segment is for decimal machines only.
C-----------------------------------------------------------------------
  420 IEXP = 2
      IZ = IBETA
  430 IF (K .LT. IZ) GO TO 440
         IZ = IZ * IBETA
         IEXP = IEXP + 1
      GO TO 430
  440 MX = IZ + IZ - 1
C-----------------------------------------------------------------------
C  Loop to determine MINEXP, XMIN.
C  Exit from loop is signaled by an underflow.
C-----------------------------------------------------------------------
  450 XMIN = Y
         Y = Y * BETAIN
C-----------------------------------------------------------------------
C  Check for underflow here.
C-----------------------------------------------------------------------
         A = Y * ONE
         TEMP = Y * T
         IF (((A+A) .EQ. ZERO) .OR. (ABS(Y) .GE. XMIN)) GO TO 460
         K = K + 1
         TEMP1 = TEMP * BETAIN
         IF ((TEMP1*BETA .NE. Y) .OR. (TEMP .EQ. Y)) THEN
               GO TO 450
            ELSE
               NXRES = 3
               XMIN = Y
         END IF
  460 MINEXP = -K
C-----------------------------------------------------------------------
C  Determine MAXEXP, XMAX.
C-----------------------------------------------------------------------
      IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500
      MX = MX + MX
      IEXP = IEXP + 1
  500 MAXEXP = MX + MINEXP
C-----------------------------------------------------------------
C  Adjust IRND to reflect partial underflow.
C-----------------------------------------------------------------
      IRND = IRND + NXRES
C-----------------------------------------------------------------
C  Adjust for IEEE-style machines.
C-----------------------------------------------------------------
      IF (IRND .GE. 2) MAXEXP = MAXEXP - 2
C-----------------------------------------------------------------
C  Adjust for machines with implicit leading bit in binary
C  significand, and machines with radix point at extreme
C  right of significand.
C-----------------------------------------------------------------
      I = MAXEXP + MINEXP
      IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1
      IF (I .GT. 20) MAXEXP = MAXEXP - 1
      IF (A .NE. Y) MAXEXP = MAXEXP - 2
      XMAX = ONE - EPSNEG
      IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG
      XMAX = XMAX / (BETA * BETA * BETA * XMIN)
      I = MAXEXP + MINEXP + 3
      IF (I .LE. 0) GO TO 520
      DO 510 J = 1, I
          IF (IBETA .EQ. 2) XMAX = XMAX + XMAX
          IF (IBETA .NE. 2) XMAX = XMAX * BETA
  510 CONTINUE
  520 RETURN
C---------- LAST CARD OF MACHAR ----------
      END
      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
C  Latest revision - March 19, 1987
C
C  Author - W. J. Cody
C           Argonne National Laboratory
C
C---------------------------------------------------------------------
      INTEGER IY,J,K
CS    REAL             CONV,C1,C2,C3,ONE,REN
CD    DOUBLE PRECISION CONV,C1,C2,C3,ONE,REN
      DATA IY/100001/
CS    DATA ONE,C1,C2,C3/1.0E0,2796203.0E0,1.0E-6,1.0E-12/
CD    DATA ONE,C1,C2,C3/1.0D0,2796203.0D0,1.0D-6,1.0D-12/
C---------------------------------------------------------------------
C  Statement functions for conversion between integer and float
C---------------------------------------------------------------------
CS    CONV(J) = REAL(J)
CD    CONV(J) = DBLE(J)
C---------------------------------------------------------------------
      J = K
      IY = IY * 125
      IY = IY - (IY/2796203) * 2796203
      REN = CONV(IY) / C1 * (ONE + C2 + C3)
      RETURN
C---------- Last card of REN ----------
      END
      SUBROUTINE TABLAT
C----------------------------------------------------------------------
C  Subroutine to tabulate the results of tests for complex
C  functions.
C
C  Generic Fortran subprograms required:
C
C         ABS, AIMAG (or DIMAG), REAL (or DBLE)
C
C  Latest revision - October 11, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    COMPLEX*16
     1       CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,CZ,CZZ
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
CS    REAL
CD    DOUBLE PRECISION
     1   CONVR,CONVI,ONE,W,Z,ZERO,ZZ
CS    COMPLEX
CD    COMPLEX*16
     1   CW,CZERO,ZDUM
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
CS    DATA ONE,ZERO/1.0E0,0.0E0/
CS    DATA CZERO/(0.0E0,0.0E0)/
CD    DATA ONE,ZERO/1.0D0,0.0D0/
CD    DATA CZERO/(0.0D0,0.0D0)/
C----------------------------------------------------------------------
C  Statement functions for conversion between complex and float
C----------------------------------------------------------------------
CS    CONVR(ZDUM) = REAL(ZDUM)
CS    CONVI(ZDUM) = AIMAG(ZDUM)
CD    CONVR(ZDUM) = DBLE(ZDUM)
CD    CONVI(ZDUM) = DIMAG(ZDUM)
C-----------------------------------------------------------------
C Handle vector error
C-----------------------------------------------------------------
      IF (CZ .NE. CZERO) THEN
            CW = (CZ-CZZ)/CZ
         ELSE
            CW = ONE
      END IF
      W = ABS(CW)
      IF (W .LT. ZERO) THEN
            K1 = K1 + 1
         ELSE IF (W .GT. ZERO) THEN
            K3 = K3 + 1
      END IF
      IF (W .GT. R6) THEN
         R6 = W
         CX1 = CX
         IF (PFLAG) CY1 = CV
         CF1 = CZ
      END IF
      R7 = R7 + W*W
C-----------------------------------------------------------------
C Handle error in real component
C-----------------------------------------------------------------
      IF (RFLAG) THEN
         Z = CONVR(CZ)
         ZZ = CONVR(CZZ)
         IF (Z .NE. ZERO) THEN
               W = (Z-ZZ)/Z
            ELSE
               W = ONE
         END IF
         IF (W .LT. ZERO) THEN
               KX1 = KX1 + 1
            ELSE IF (W .GT. ZERO) THEN
               KX3 = KX3 + 1
         END IF
         W = ABS(W)
         IF (W .GT. RX6) THEN
            RX6 = W
            CX2 = CX
            IF (PFLAG) CY2 = CV
            CF2 = CZ
         END IF
         RX7 = RX7 + W*W
      END IF
C-----------------------------------------------------------------
C Handle error in imaginary component
C-----------------------------------------------------------------
      IF (IFLAG) THEN
         Z = CONVI(CZ)
         ZZ = CONVI(CZZ)
         IF (Z .NE. ZERO) THEN
               W = (Z-ZZ)/Z
            ELSE
               W = ONE
         END IF
         IF (W .LT. ZERO) THEN
               KY1 = KY1 + 1
            ELSE IF (W .GT. ZERO) THEN
               KY3 = KY3 + 1
         END IF
         W = ABS(W)
         IF (W .GT. RY6) THEN
            RY6 = W
            CX3 = CX
            IF (PFLAG) CY3 = CV
            CF3 = CZ
         END IF
         RY7 = RY7 + W*W
      END IF
      RETURN
C---------- Last line of TABLAT ----------
      END
*
      SUBROUTINE REPORT(IOUT)
C----------------------------------------------------------------------
C  Subroutine to print the results of tests for complex functions.
C
C  Generic Fortran subprograms required:
C
C         ABS, LOG, MAX, SQRT
C
C  Latest revision - October 11, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    COMPLEX*16
     1       CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,CZ,CZZ
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
      INTEGER IOUT
CS    REAL
CD    DOUBLE PRECISION
     1   ALL9,W,ZERO
C----------------------------------------------------------------------
C  Mathematical constants
C----------------------------------------------------------------------
CS    DATA ALL9,ZERO/-999.0E0,0.0E0/
CD    DATA ALL9,ZERO/-999.0D0,0.0D0/
C----------------------------------------------------------------------
         WRITE (IOUT,1000) N,AX,BX,AY,BY
C-----------------------------------------------------------------
C Report vector error
C-----------------------------------------------------------------
         WRITE (IOUT,1001)
         K2 = N - K3 - K1
         R7 = SQRT(R7/XN)
         WRITE (IOUT,1011) K1,K2,K3
         WRITE (IOUT,1020) IT,IBETA
         IF (R6 .NE. ZERO) THEN
               W = LOG(ABS(R6))/ALBETA
            ELSE
               W = ALL9
         END IF
         IF (.NOT. PFLAG) THEN
               WRITE (IOUT,1021) R6,IBETA,W,CX1,CF1
            ELSE
               WRITE (IOUT,1022) R6,IBETA,W,CX1,CY1,CF1
         END IF
         W = MAX(AIT+W,ZERO)
         WRITE (IOUT,1023) IBETA,W
         IF (R7 .NE. ZERO) THEN
               W = LOG(ABS(R7))/ALBETA
            ELSE
               W = ALL9
         END IF
         WRITE (IOUT,1024) R7,IBETA,W
         W = MAX(AIT+W,ZERO)
         WRITE (IOUT,1023) IBETA,W
C-----------------------------------------------------------------
C Report error in real component
C-----------------------------------------------------------------
         IF (RFLAG) THEN
            WRITE (IOUT,1002)
            KX2 = N - KX3 - KX1
            RX7 = SQRT(RX7/XN)
            WRITE (IOUT,1011) KX1,KX2,KX3
            WRITE (IOUT,1020) IT,IBETA
            IF (RX6 .NE. ZERO) THEN
                  W = LOG(ABS(RX6))/ALBETA
               ELSE
                  W = ALL9
            END IF
            IF (.NOT. PFLAG) THEN
                  WRITE (IOUT,1021) RX6,IBETA,W,CX2,CF2
               ELSE
                  WRITE (IOUT,1022) RX6,IBETA,W,CX2,CY2,CF2
            END IF
            W = MAX(AIT+W,ZERO)
            WRITE (IOUT,1023) IBETA,W
            IF (RX7 .NE. ZERO) THEN
                  W = LOG(ABS(RX7))/ALBETA
               ELSE
                  W = ALL9
            END IF
            WRITE (IOUT,1024) RX7,IBETA,W
            W = MAX(AIT+W,ZERO)
            WRITE (IOUT,1023) IBETA,W
         END IF
C-----------------------------------------------------------------
C Report error in imaginary component
C-----------------------------------------------------------------
         IF (IFLAG) THEN
            WRITE (IOUT,1003)
            KY2 = N - KY3 - KY1
            RY7 = SQRT(RY7/XN)
            WRITE (IOUT,1011) KY1,KY2,KY3
            WRITE (IOUT,1020) IT,IBETA
            IF (RY6 .NE. ZERO) THEN
                  W = LOG(ABS(RY6))/ALBETA
               ELSE
                  W = ALL9
            END IF
            IF (.NOT. PFLAG) THEN
                  WRITE (IOUT,1021) RY6,IBETA,W,CX3,CF3
               ELSE
                  WRITE (IOUT,1022) RY6,IBETA,W,CX3,CY3,CF3
            END IF
            W = MAX(AIT+W,ZERO)
            WRITE (IOUT,1023) IBETA,W
            IF (RY7 .NE. ZERO) THEN
                  W = LOG(ABS(RY7))/ALBETA
               ELSE
                  W = ALL9
            END IF
            WRITE (IOUT,1024) RY7,IBETA,W
            W = MAX(AIT+W,ZERO)
            WRITE (IOUT,1023) IBETA,W
         END IF
      RETURN
C-----------------------------------------------------------------
 1000 FORMAT(I7,' Random arguments were tested from the box'/
     1 6X,'(',E15.4,',',E15.4,')'/6X,'(',E15.4,',',E15.4,1H)//)
 1001 FORMAT(' First, report the vector error'//)
 1002 FORMAT(' Then, report the error in the real component.'//)
 1003 FORMAT(' Then, report the error in the imaginary component.'//)
 1011 FORMAT(' Test function was larger',I6,' times,'/
     1     18X,' agreed',I6,' times, and'/
     2   14X,'was smaller',I6,' times.'//)
 1020 FORMAT(' There are',I4,' base',I4,
     1    ' significant digits in a floating-point number.'//)
 1021 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **',
     1  F7.2/4X,'occurred for Z =',(E17.6,E17.6)/
     2       4X,'with a value of ',(E17.6,E17.6))
 1022 FORMAT(' The maximum relative error of',E15.4,' = ',I4,' **',
     1  F7.2/4X,'occurred for Z =',(E17.6,E17.6)/
     2       13X,'and W =',(E17.6,E17.6)/
     3       4X,'with a value of ',(E17.6,E17.6))
 1023 FORMAT(' The estimated loss of base',I4,
     1  ' significant digits is',F7.2//)
 1024 FORMAT(' The root mean square relative error was',E15.4,
     1    ' = ',I4,' **',F7.2)
C---------- Last line of REPORT ----------
      END
*
      SUBROUTINE RESET
C----------------------------------------------------------------------
C  Subroutine to reset the paramaters for tests for complex
C  functions.
C
C  Latest revision - October 15, 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
CD    DOUBLE PRECISION
     1       AIT,ALBETA,AX,AY,BX,BY,RX6,RX7,RY6,RY7,R6,R7,XN
CS    COMPLEX
CD    COMPLEX*16
     1       CF1,CF2,CF3,CV,CX,CX1,CX2,CX3,CY1,CY2,CY3,CZ,CZZ
C----------------------------------------------------------------------
C  Declare local variables
C----------------------------------------------------------------------
CS    REAL
CD    DOUBLE PRECISION
     1   ZERO
CS    COMPLEX
CD    COMPLEX*16
     1   CZERO
C-----------------------------------------------------------------
CS    DATA CZERO/(0.0E0,0.0E0)/
CS    DATA ZERO/0.0E0/
CD    DATA CZERO/(0.0D0,0.0D0)/
CD    DATA ZERO/0.0D0/
C-----------------------------------------------------------------
      K1 = 0
      K3 = 0
      KX1 = 0
      KX3 = 0
      KY1 = 0
      KY3 = 0
      CX1 = CZERO
      CX2 = CZERO
      CX3 = CZERO
      CY1 = CZERO
      CY2 = CZERO
      CY3 = CZERO
      CF1 = CZERO
      CF2 = CZERO
      CF3 = CZERO
      R6 = ZERO
      R7 = ZERO
      RX6 = ZERO
      RX7 = ZERO
      RY6 = ZERO
      RY7 = ZERO
      RETURN
C---------- Last line of RESET ----------
      END
