C     TESTPOWER --- TEST THE POWER FUNCTION
C
C        Eugene Spafford
C        SWT Math Library Test Support
C        School of Information and Computer Science
C        Georgia Institute of Technology
C        Atlanta, GA  30332
C
C
      SUBROUTINE TEST
C
      REAL ALXMAX,SQBETA,ONEP5,SCALE,DELY,Y1,Y2,XSQ
      DOUBLE PRECISION RNAME
$INSERT MACHAR.F.I
      DATA RNAME /'** 66 '/
C
C
      BETA = IBETA
      SQBETA = XSQRT(BETA)
      ALBETA = XALOG(BETA)
      AIT = IT
      ALXMAX = XALOG(XMAX)
      ONEP5 = (TWO+ONE)/TWO
      SCALE = ONE
      J = (IT+1)/2
C
      DO 10 I = 1,J
        SCALE = SCALE*BETA
   10 CONTINUE
C
      A = ONE/BETA
      B = ONE
      C = -AMAX1(ALXMAX,-XALOG(XMIN))/XALOG(100E0)
      DELY = -C-C
      Y1 = ZERO
C
C     START WITH RANDOM ACCURACY TESTS
C
      DO 300  J = 1,4
        K1 = 0
        K3 = 0
        X1 = ZERO
        R6 = ZERO
        R7 = ZERO
        DEL = (B-A)/XN
        XL = A
C
        DO 20 I = 1,N
          X = DEL*RANDX(I1)+XL
          IF (J .NE. 1) GO TO 50
              ZZ = X**ONE
              Z = X
              GO TO 110
50        CONTINUE
              W = SCALE*X
              Z = ONE
              X = X+W
              Z = TWO
              X = X-W
              XSQ = X*X
              IF (J .EQ. 4) GO TO 70
                  ZZ = XSQ**ONEP5
                  Z = X*XSQ
                  GO TO 110
70            CONTINUE
                  Y = DELY*RANDX(I1)+C
                  Y2 = (Y/TWO+Y)-Y
                  Y = Y2+Y2
                  Z = X**Y
                  ZZ = XSQ**Y2
110   CONTINUE
          W = ONE
          IF (Z .NE. ZERO) W = (Z-ZZ)/Z
          IF (W .GT. ZERO) K1 = K1+1
          IF (W .GE. ZERO) GO TO 119
              K3 = K3+1
              W = -W
119   CONTINUE
          IF (W .LE. R6) GO TO 120
              R6 = W
              X1 = X
              IF (J .EQ. 4) Y1 = Y
120   CONTINUE
          R7 = R7+W*W
          XL = XL+DEL
   20   CONTINUE
C
      K2 = N-K1-K3
      R7 = XSQRT(R7/XN)
      IF (J .GT. 1) GO TO 210
        PRINT 1000
        PRINT 1010, N,A,B
        PRINT 1011, RNAME,K1,K2,K3
        GO TO 220
210   CONTINUE
      IF (J .EQ. 4) GO TO 215
        PRINT 1001
        PRINT 1010, N,A,B
        PRINT 1012, RNAME,K1,K2,K3
        GO TO 220
215   CONTINUE
        PRINT 1002
        W = C+DELY
        PRINT 1014, N,A,B,C,W
        PRINT 1013, RNAME, K1,K2,K3
220   CONTINUE
      PRINT 1020, IT,IBETA
      W = -999.0
      IF (R6 .NE. ZERO) W = XALOG(ABS(R6))/ALBETA
      IF (J .NE. 4) PRINT 1021, R6,IBETA,W,X1
      IF (J .EQ. 4) PRINT 1024, R6,IBETA,W,X1,Y1
      W = AMAX1(AIT+W,ZERO)
      PRINT 1022, IBETA,W
      W = -999.0
      IF (R7 .NE. ZERO) W = XALOG(ABS(R7))/ALBETA
      PRINT 1023, R7,IBETA,W
      W = AMAX1(AIT+W,ZERO)
      PRINT 1022, IBETA,W
      IF (J .EQ. 1) GO TO 300
        B = TEN
        A = 0.01E0
        IF (J .EQ. 3) GO TO 300
          A = ONE
          B = XEXP(ALXMAX/THREE)
300   CONTINUE
C
C     SPECIAL TESTS
C
      PRINT 1025
      PRINT 1030
      B = TEN
C
      DO 320 I =1,5
      X = RANDX(I1)*B+ONE
      Y = RANDX(I1)*B+ONE
      Z = X**Y
      ZZ = (ONE/X)**(-Y)
      W = (Z-ZZ)/Z
      PRINT 1060, X,Y,W
320   CONTINUE
C
C     TEST OF ERROR RETURNS
C
      PRINT 1050
      X = BETA
      Y = MINEXP
      PRINT 1051, X, Y
      Z = X**Y
      PRINT 1055, Z
      Y = MAXEXP-1
      PRINT 1051, X, Y
      Z = X**Y
      PRINT 1055, Z
      X = ZERO
      Y = TWO
      PRINT 1051, X, Y
      Z = X**Y
      PRINT 1055, Z
      X = -Y
      Y = ZERO
      PRINT 1052, X, Y
      Z = X**Y
      PRINT 1055, Z
      Y = TWO
      PRINT 1052, X, Y
      Z = X**Y
      PRINT 1055, Z
      X = ZERO
      Y = ZERO
      PRINT 1052, X, Y
      Z = X**Y
      PRINT 1055, Z
      PRINT 1100
      RETURN
C
1000  FORMAT (/'Test of X**1.0 vs. X'//)
1001  FORMAT (/'Test of XSQ**1.5 vs. XSQ*X'//)
1002  FORMAT (/'Test of X**Y vs. XSQ**(Y/2)'//)
1010  FORMAT (I7,' Random arguments were tested in the interval '/6X,'('
     $,E12.4,',',E12.4,')'//)
1011  FORMAT (1X,A4,'(X**1.0) was larger',I6,' times,'/12X,' agreed',I6,
     $' times, and'/8X,'was smaller',I6,' times.'//)
1012  FORMAT (1X,A4,'(X**1.5) was larger',I6,' times,'/12X,' agreed',I6,
     $' times, and'/8X,'was smaller',I6,' times.'//)
1013  FORMAT (1X,A4,'(X**Y) was larger',I6,' times,'/12X,' agreed',I6,
     $' times, and'/8X,'was smaller',I6,' times.'//)
1014  FORMAT (I7,' Random arguments were tested from the region'/
     &6X,'X in (',E13.4,',',E13.4,') Y in (',E13.4,',',E13.4,')'//)
1020  FORMAT (' There are ',I3,' base',I3,
     $' significant digits in a floating point number.'//)
1021  FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2
     $/4X,'occurred for X =',E17.6)
1022  FORMAT (' The estimated loss of base',I3,' significant digits is',
     $F7.2)
1023  FORMAT (' The root mean square relative error was',E15.4,' = ',I3,
     $' **',F7.2)
1024  FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2
     $/4X,'occurred for X =',E17.6,' Y =',E17.6)
1025  FORMAT (//'Special Tests'//)
1030  FORMAT (' The identity X**Y = (1/X)**(-Y) will be tested.'//
     &8X,'X',14X,'Y',9X,'(X**Y-(1/X)**(-Y))/X**Y'/)
1040  FORMAT (///'Test of Special Arguments'//)
1050  FORMAT (///'Test of Error Returns'//)
1051  FORMAT (' (',E14.7,' ) ** (',E14.7,' ) will be computed.'/
     &' This should not trigger an error message.'//)
1052  FORMAT (' (',E14.7,' ) ** (',E14.7,' ) will be computed.'/
     &' This should trigger an error message.'//)
1055  FORMAT (' The value returned is',E15.4///)
1060  FORMAT (2E15.7,6X,E15.7)
1100  FORMAT (10X,'***** This concludes the tests. *****'//)
      END