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
$INSERT MACHAR.F77.I
C
DOUBLE PRECISION ALXMAX,SQBETA,ONEP5,SCALE,DELY,Y1,Y2,XSQ,U1
DOUBLE PRECISION ZPOWR
EXTERNAL ZPOWR
CHARACTER*6 RNAME
C
C
IF (WHICH) THEN
RNAME = 'DPWR$M'
ELSE
RNAME = '** D77'
ENDIF
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 30 J = 1,4
K1 = 0
K3 = 0
X1 = ZERO
R6 = ZERO
R7 = ZERO
DEL = (B-A)/XN
XL = A
U1 = XPOWER(XMAX,ONE/THREE)
C
DO 20 I = 1,N
X = DEL*RANDX(I1)+XL
IF (J .EQ. 1) THEN
ZZ = ZPOWR(X,ONE)
Z = X
C
ELSE
W = SCALE*X
X = (X+W)-W
IF (X .GT. U1) X = X+X-U1
XSQ = X*X
IF (J .NE. 4) THEN
ZZ = ZPOWR(XSQ,ONEP5)
Z = X*XSQ
C
ELSE
Y = DELY*RANDX(I1)+C
Y2 = (Y/TWO+Y)-Y
Y = Y2+Y2
Z = ZPOWR(X,Y)
ZZ = ZPOWR(XSQ,Y2)
ENDIF
ENDIF
W = ONE
IF (Z .NE. ZERO) W = (Z-ZZ)/Z
IF (W .GT. ZERO) K1 = K1+1
IF (W .LT. ZERO) THEN
K3 = K3+1
W = -W
ENDIF
IF (W .GT. R6) THEN
R6 = W
X1 = X
IF (J .EQ. 4) Y1 = Y
ENDIF
R7 = R7+W*W
XL = XL+DEL
20 CONTINUE
C
K2 = N-K1-K3
R7 = XSQRT(R7/XN)
IF (J .LE. 1) THEN
PRINT 50
PRINT 80, N,A,B
PRINT 90, RNAME,K1,K2,K3
C
ELSEIF (J .NE. 4) THEN
PRINT 60
PRINT 80, N,A,B
PRINT 100, RNAME,K1,K2,K3
C
ELSE
PRINT 70
W = C+DELY
PRINT 120, N,A,B,C,W
PRINT 110, RNAME,K1,K2,K3
ENDIF
PRINT 130, IT,IBETA
W = -999.0D0
IF (R6 .NE. ZERO) W = XALOG(DABS(R6))/ALBETA
IF (J .NE. 4) THEN
PRINT 140, R6,IBETA,W,X1
C
ELSE
PRINT 170, R6,IBETA,W,X1,Y1
ENDIF
W = DMAX1(AIT+W,ZERO)
PRINT 150, IBETA,W
W = -999.0D0
IF (R7 .NE. ZERO) W = XALOG(DABS(R7))/ALBETA
PRINT 160, R7,IBETA,W
W = DMAX1(AIT+W,ZERO)
PRINT 150, IBETA,W
IF (J .NE. 1) THEN
B = TEN
A = 0.01D0
IF (J .NE. 3) THEN
A = ONE
B = XEXP(ALXMAX/THREE)
ENDIF
ENDIF
30 CONTINUE
C
C SPECIAL TESTS
C
PRINT 180
PRINT 190
B = TEN
C
DO 40 I = 1,5
X = RANDX(I1)*B+ONE
Y = RANDX(I1)*B+ONE
Z = ZPOWR(X,Y)
ZZ = ZPOWR((ONE/X),-Y)
W = (Z-ZZ)/Z
PRINT 240, X,Y,W
40 CONTINUE
C
C TEST OF ERROR RETURNS
C
PRINT 200
X = BETA
Y = MINEXP
PRINT 210, X,Y
Z = ZPOWR(X,Y)
PRINT 230, Z
Y = MAXEXP-1
PRINT 210, X,Y
Z = ZPOWR(X,Y)
PRINT 230, Z
X = ZERO
Y = TWO
PRINT 210, X,Y
Z = ZPOWR(X,Y)
PRINT 230, Z
X = -Y
Y = ZERO
PRINT 220, X,Y
Z = ZPOWR(X,Y)
PRINT 230, Z
Y = TWO
PRINT 220, X,Y
Z = ZPOWR(X,Y)
PRINT 230, Z
X = ZERO
Y = ZERO
PRINT 220, X,Y
Z = ZPOWR(X,Y)
PRINT 230, Z
PRINT 250
RETURN
C
50 FORMAT (/'Test of X**1.0 vs. X'//)
60 FORMAT (/'Test of XSQ**1.5 vs. XSQ*X'//)
70 FORMAT (/'Test of X**Y vs. XSQ**(Y/2)'//)
80 FORMAT (I7,' Random arguments were tested in the interval '/6X,'('
&,E12.4,',',E12.4,')'//)
90 FORMAT (1X,A6,'(X**1.0) was larger',I6,' times,'/14X,' agreed',I6,
&' times, and'/10X,'was smaller',I6,' times.'//)
100 FORMAT (1X,A4,'(X**1.5) was larger',I6,' times,'/12X,' agreed',I6,
&' times, and'/8X,'was smaller',I6,' times.'//)
110 FORMAT (1X,A4,'(X**Y) was larger',I6,' times,'/12X,' agreed',I6,
&' times, and'/8X,'was smaller',I6,' times.'//)
120 FORMAT (I7,' Random arguments were tested from the region'/6X,
&'X in (',E13.4,',',E13.4,') Y in (',E13.4,',',E13.4,')'//)
130 FORMAT (' There are ',I3,' base',I3,
&' significant digits in a floating point number.'//)
140 FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2
&/4X,'occurred for X =',E17.6)
150 FORMAT (' The estimated loss of base',I3,' significant digits is',
&F7.2)
160 FORMAT (' The root mean square relative error was',E15.4,' = ',I3,
&' **',F7.2)
170 FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2
&/4X,'occurred for X =',E17.6,' Y =',E17.6)
180 FORMAT (//'Special Tests'//)
190 FORMAT (' The identity X**Y = (1/X)**(-Y) will be tested.'//8X,'X'
&,14X,'Y',9X,'(X**Y-(1/X)**(-Y))/X**Y'/)
200 FORMAT (///'Test of Error Returns'//)
210 FORMAT (' (',E14.7,' ) ** (',E14.7,' ) will be computed.'/
&' This should not trigger an error message.'//)
220 FORMAT (' (',E14.7,' ) ** (',E14.7,' ) will be computed.'/
&' This should trigger an error message.'//)
230 FORMAT (' The value returned is',E15.4///)
240 FORMAT (2E15.7,6X,E15.7)
250 FORMAT (10X,'***** This concludes the tests. *****'//)
END