C     TEST --- Subroutine to test the double precision POWER function
C
C     Eugene Spafford
C     Software Tools Subsystem Math Library Test Routine
C     School of Information and Computer Science
C     Georgia Institute of Technology
C     Atlanta, Georgia 30332
C
C     Adapted from:
C     "Software Manual for the Elementary Functions"
C       by William J. Cody, Jr.  &  William Waite
C      Prentice-Hall, Englewood Cliffs, NJ  1980
C
C     Coded April 1983 by Eugene Spafford
C
C     ----------------------------------------------------------
C
      SUBROUTINE TEST
C
      DOUBLE PRECISION ALXMAX,SQBETA,ONEP5,SCALE,DELY,Y1,Y2,XSQ
      DOUBLE PRECISION RNAME
$INSERT MACHAR.F.I
      DATA RNAME /'** D66'/
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 = -DMAX1(ALXMAX,-XALOG(XMIN))/XALOG(100D0)
      DELY = -C-C
      Y1 = ZERO
C
C     START WITH RANDOM ACCURACY TESTS
C
      DO 110 J = 1,4
        K1 = 0
        K3 = 0
        X1 = ZERO
        R6 = ZERO
        R7 = ZERO
        DEL = (B-A)/XN
        XL = A
C
        DO 70 I = 1,N
          X = DEL*RANDX(I1)+XL
          IF (J .NE. 1) GO TO 20
          ZZ = X**ONE
          Z = X
          GO TO 40
   20     CONTINUE
          W = SCALE*X
          X = (X+W)-W
          XSQ = X*X
          IF (J .EQ. 4) GO TO 30
          ZZ = XSQ**ONEP5
          Z = X*XSQ
          GO TO 40
   30     CONTINUE
          Y = DELY*RANDX(I1)+C
          Y2 = (Y/TWO+Y)-Y
          Y = Y2+Y2
          Z = X**Y
          ZZ = XSQ**Y2
   40     CONTINUE
          W = ONE
          IF (Z .NE. ZERO) W = (Z-ZZ)/Z
          IF (W .GT. ZERO) K1 = K1+1
          IF (W .GE. ZERO) GO TO 50
          K3 = K3+1
          W = -W
   50     CONTINUE
          IF (W .LE. R6) GO TO 60
          R6 = W
          X1 = X
          IF (J .EQ. 4) Y1 = Y
   60     CONTINUE
          R7 = R7+W*W
          XL = XL+DEL
   70   CONTINUE
C
        K2 = N-K1-K3
        R7 = XSQRT(R7/XN)
        IF (J .GT. 1) GO TO 80
        PRINT 130
        PRINT 160, N,A,B
        PRINT 170, RNAME,K1,K2,K3
        GO TO 100
   80   CONTINUE
        IF (J .EQ. 4) GO TO 90
        PRINT 140
        PRINT 160, N,A,B
        PRINT 180, RNAME,K1,K2,K3
        GO TO 100
   90   CONTINUE
        PRINT 150
        W = C+DELY
        PRINT 200, N,A,B,C,W
        PRINT 190, RNAME,K1,K2,K3
  100   CONTINUE
        PRINT 210, IT,IBETA
        W = -999.0D0
        IF (R6 .NE. ZERO) W = XALOG(DABS(R6))/ALBETA
        IF (J .NE. 4) PRINT 220, R6,IBETA,W,X1
        IF (J .EQ. 4) PRINT 250, R6,IBETA,W,X1,Y1
        W = DMAX1(AIT+W,ZERO)
        PRINT 230, IBETA,W
        W = -999.0D0
        IF (R7 .NE. ZERO) W = XALOG(DABS(R7))/ALBETA
        PRINT 240, R7,IBETA,W
        W = DMAX1(AIT+W,ZERO)
        PRINT 230, IBETA,W
        IF (J .EQ. 1) GO TO 110
        B = TEN
        A = 0.01D0
        IF (J .EQ. 3) GO TO 110
        A = ONE
        B = XEXP(ALXMAX/THREE)
  110 CONTINUE
C
C     SPECIAL TESTS
C
      PRINT 260
      PRINT 270
      B = TEN
C
      DO 120 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 320, X,Y,W
  120 CONTINUE
C
C     TEST OF ERROR RETURNS
C
      PRINT 280
      X = BETA
      Y = MINEXP
      PRINT 290, X,Y
      Z = X**Y
      PRINT 310, Z
      Y = MAXEXP-1
      PRINT 290, X,Y
      Z = X**Y
      PRINT 310, Z
      X = ZERO
      Y = TWO
      PRINT 290, X,Y
      Z = X**Y
      PRINT 310, Z
      X = -Y
      Y = ZERO
      PRINT 300, X,Y
      Z = X**Y
      PRINT 310, Z
      Y = TWO
      PRINT 300, X,Y
      Z = X**Y
      PRINT 310, Z
      X = ZERO
      Y = ZERO
      PRINT 300, X,Y
      Z = X**Y
      PRINT 310, Z
      PRINT 330
      RETURN
C
  130 FORMAT (/'Test of X**1.0 vs. X'//)
  140 FORMAT (/'Test of XSQ**1.5 vs. XSQ*X'//)
  150 FORMAT (/'Test of X**Y vs. XSQ**(Y/2)'//)
  160 FORMAT (I7,' Random arguments were tested in the interval '/6X,'('
     &,E12.4,',',E12.4,')'//)
  170 FORMAT (1X,A4,'(X**1.0) was larger',I6,' times,'/12X,' agreed',I6,
     &' times, and'/8X,'was smaller',I6,' times.'//)
  180 FORMAT (1X,A4,'(X**1.5) was larger',I6,' times,'/12X,' agreed',I6,
     &' times, and'/8X,'was smaller',I6,' times.'//)
  190 FORMAT (1X,A4,'(X**Y) was larger',I6,' times,'/12X,' agreed',I6,
     &' times, and'/8X,'was smaller',I6,' times.'//)
  200 FORMAT (I7,' Random arguments were tested from the region'/6X,
     &'X in (',E13.4,',',E13.4,') Y in (',E13.4,',',E13.4,')'//)
  210 FORMAT (' There are ',I3,' base',I3,
     &' significant digits in a floating point number.'//)
  220 FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2
     &/4X,'occurred for X =',E17.6)
  230 FORMAT (' The estimated loss of base',I3,' significant digits is',
     &F7.2)
  240 FORMAT (' The root mean square relative error was',E15.4,' = ',I3,
     &' **',F7.2)
  250 FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2
     &/4X,'occurred for X =',E17.6,' Y =',E17.6)
  260 FORMAT (//'Special Tests'//)
  270 FORMAT (' The identity X**Y = (1/X)**(-Y) will be tested.'//8X,'X'
     &,14X,'Y',9X,'(X**Y-(1/X)**(-Y))/X**Y'/)
  280 FORMAT (///'Test of Error Returns'//)
  290 FORMAT (' (',E14.7,' ) ** (',E14.7,' ) will be computed.'/
     &' This should not trigger an error message.'//)
  300 FORMAT (' (',E14.7,' ) ** (',E14.7,' ) will be computed.'/
     &' This should trigger an error message.'//)
  310 FORMAT (' The value returned is',E15.4///)
  320 FORMAT (2E15.7,6X,E15.7)
  330 FORMAT (10X,'***** This concludes the tests. *****'//)
      END