C TEST --- Subroutine to test the single precision LOG & LOG10 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 REAL U1,U2,ALOG,ALOG10 DOUBLE PRECISION RNAME,RNAME2 $INSERT MACHAR.F.I DATA RNAME,RNAME2/'ALOG66','ALOG10'/ C C BETA = IBETA ALBETA = XALOG(BETA) AIT = IT J = IT/3 C = ONE C DO 10 I = 1,J C = C/BETA 10 CONTINUE C B = ONE+C A = ONE-C C DO 170 J = 1,4 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B-A)/XN XL = A C DO 80 I = 1,N X = DEL*RANDX(I1)+XL GO TO (20,30,40,50), J GO TO 60 C 20 CONTINUE Y = (X-HALF)-HALF ZZ = ALOG(X) Z = ONE/THREE Z = Y*(Z-Y/FOUR) Z = (Z-HALF)*Y*Y+Y GO TO 60 C 30 CONTINUE U1 = X+EIGHT U2 = ZZ+ZERO X = U1-EIGHT Y = X+X/16.0 Z = ALOG(X) ZZ = ALOG(Y)-7.7746816434843E-5 ZZ = ZZ-(31.0/512.0) GO TO 60 C 40 CONTINUE U1 = X+EIGHT U2 = ZZ+ZERO X = U1-EIGHT Y = X+X*TENTH Z = ALOG10(X) ZZ = ALOG10(Y)-3.7706015822504E-4 ZZ = ZZ-(21.0/512.0) GO TO 60 C 50 CONTINUE Z = ALOG(X*X) ZZ = ALOG(X) ZZ = ZZ+ZZ C 60 CONTINUE W = ONE IF (Z .NE. ZERO) W = (Z-ZZ)/Z Z = SIGN(W,Z) IF (Z .GT. ZERO) K1 = K1+1 IF (Z .LT. ZERO) K3 = K3+1 W = ABS(W) IF (W .LE. R6) GO TO 70 R6 = W X1 = X 70 CONTINUE R7 = R7+W*W XL = XL+DEL 80 CONTINUE C K2 = N-K3-K1 GO TO (90,100,110,120), J GO TO 130 90 CONTINUE PRINT 190 PRINT 230, N,C PRINT 250, RNAME,K1,K2,K3 GO TO 130 100 CONTINUE PRINT 200 PRINT 240, N,A,B PRINT 250, RNAME,K1,K2,K3 GO TO 130 110 CONTINUE PRINT 220 PRINT 240, N,A,B PRINT 250, RNAME2,K1,K2,K3 GO TO 130 120 CONTINUE PRINT 210 PRINT 240, N,A,B PRINT 250, RNAME,K1,K2,K3 130 CONTINUE R7 = XSQRT(R7/XN) PRINT 260, IT,IBETA W = -999.0 IF (R6 .NE. ZERO) W = XALOG(ABS(R6))/ALBETA PRINT 270, R6,IBETA,W,X1 W = AMAX1(AIT+W,ZERO) PRINT 280, IBETA,W W = -999.0 IF (R7 .NE. ZERO) W = XALOG(ABS(R7))/ALBETA PRINT 290, R7,IBETA,W W = AMAX1(AIT+W,ZERO) PRINT 280, IBETA,W IF (J .NE. 1) GO TO 140 A = XSQRT(HALF) B = 15.0/16.0 GO TO 160 140 CONTINUE IF (J .NE. 2) GO TO 150 A = XSQRT(TENTH) B = 0.9 GO TO 160 150 CONTINUE A = 16.0 B = 240.0 160 CONTINUE 170 CONTINUE C C SPECIAL TESTS C PRINT 300 PRINT 350 C DO 180 I = 1,5 X = RANDX(I1) X = X+X+15.0 Y = ONE/X Z = ALOG(X)+ALOG(Y) PRINT 330, X,Z 180 CONTINUE C PRINT 310 X = ONE Y = ALOG(X) PRINT 360, Y X = XMIN Y = ALOG(X) PRINT 370, X,Y X = XMAX Y = ALOG(X) PRINT 380, X,Y C C TEST OF ERROR RETURNS C PRINT 320 X = -TWO PRINT 390, X Y = ALOG(X) PRINT 400, Y X = ZERO PRINT 390, X Y = ALOG(X) PRINT 400, Y PRINT 340 RETURN C C C ***** END OF ROUTINE TO TEST THE LOG FUNCTIONS ***** C 190 FORMAT (/// &'Test of ALOG(X) vs. Taylor Series Expansion of ALOG(1+Y)'//) 200 FORMAT (///'Test of ALOG(X) vs. ALOG(17X/16)-ALOG(17/16)'//) 210 FORMAT (///'Test of ALOG(X*X) vs. 2 * ALOG(X)'//) 220 FORMAT (///'Test of ALOG10(X) vs. ALOG10(11X/10) - ALOG10(11/10)'/ &/) 230 FORMAT (I7,' random arguments were tested from the interval'/6X, &'(1-EPS, 1+EPS), where EPS =',E15.4//) 240 FORMAT (I6,' Random arguments were tested in the interval '/6X,'(' &,E12.4,',',E12.4,')'//) 250 FORMAT (1X,A6,'(X) was larger',I6,' times,'/14X,' agreed',I6, &' times, and'/10X,'was smaller',I6,' times.'//) 260 FORMAT (' There are ',I3,' base',I3, &' significant digits in a floating point number.'//) 270 FORMAT (' The maximum relative error of',E12.4,' = ',I3,' **',F7.2 &/4X,'occurred for X =',E17.6) 280 FORMAT (' The estimated loss of base',I3,' significant digits is', &F7.2) 290 FORMAT (' The root mean square relative error was',E15.4,' = ',I3, &' **',F7.2) 300 FORMAT (//'Special Tests'//) 310 FORMAT (///'Test of Special Arguments'//) 320 FORMAT (///'Test of Error Returns'//) 330 FORMAT (2E15.7/) 340 FORMAT (10X,'***** This concludes the tests. *****'//) 350 FORMAT (' The identity ALOG(X) = - ALOG(1/X) will be tested.'//8X, &'X',9X,'F(X) + F(1/X)'/) 360 FORMAT (' ALOG(1.0) = ',E15.7//) 370 FORMAT (' ALOG(XMIN) = ALOG(',E15.7,') =',E15.7//) 380 FORMAT (' ALOG(XMAX) = ALOG(',E15.7,') =',E15.7//) 390 FORMAT (' ALOG will be called with the argument',E15.4/ &' This should cause an error.'//) 400 FORMAT (' ALOG returned the value',E15.4//) END