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
      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,CONVR,C1,C2,DELX,DELY,EPS,EPSNEG,FOUR,HALF,
     2   ONE,REN,SCALE,TEN,W,X,XL,XMAX,XMIN,Y,ZERO
CS    COMPLEX
      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/
      DATA ZERO,HALF,ONE,FOUR,TEN/0.0D0,0.5D0,1.0D0,4.0D0,10.0D0/
      DATA CONE,CTWO/(1.0D0,0.0D0),(2.0D0,0.0D0)/
      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)
      CONV(NDUM) = DBLE(NDUM)
      CONVC(X,Y) = DCMPLX(X,Y)
      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
