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
      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,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
      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)/
      DATA ZERO,ONE16,HALF,ONE/0.0D0,0.0625D0,0.5D0,1.0D0/
      DATA ONEP6,THREE,TEN,SIXTEN/1.625D0,3.0D0,10.0D0,16.0D0/
      DATA EXPDEL/(6.2416044877018563681D-2,6.6487597751003112768D-2)/
      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)
      CONV(NDUM) = DBLE(NDUM)
      CONVR(ZDUM) = DBLE(ZDUM)
      CONVI(ZDUM) = DIMAG(ZDUM)
      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
