CALL TDSQRT(6) END SUBROUTINE TDSQRT (IOUT) * Need this to allow private dsqrt() testing EXTERNAL DSQRT DOUBLE PRECISION DSQRT EXTERNAL DRAN DOUBLE PRECISION DRAN C PROGRAM TO TEST DSQRT 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 SIX 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 EPS - THE SMALLEST POSITIVE FLOATING-POINT C NUMBER SUCH THAT 1.0+EPS .NE. 1.0 C EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT C NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0 C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C POWER OF THE RADIX C XMAX - THE LARGEST FINITE FLOATING-POINT NO. C C DRANDL(X) - A FUNCTION SUBPROGRAM RETURNING LOGARITHMICALLY C DISTRIBUTED RANDOM DOUBLE PRECISION NUMBERS. C IN PARTICULAR, C A * DRANDL(DLOG(B/A)) C IS LOGARITHMICALLY DISTRIBUTED OVER (A,B) C C DRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM DOUBLE C PRECISION NUMBERS UNIFORMLY DISTRIBUTED OVER C (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DFLOAT, DSQRT C C C LATEST REVISION - AUGUST 2, 1979 C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,K1,K2,K3,MACHEP, 1 MAXEXP,MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,DRANDL,R6, 1 R7,SQBETA,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO C C IOUT = 6 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, 1 MAXEXP,EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) SQBETA = DSQRT(BETA) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0D+00 ZERO = 0.0D+00 A = ONE / SQBETA B = ONE N = 2000 XN = DBLE(FLOAT(N)) C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 300 J = 1, 2 C = DLOG(B/A) K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO C DO 200 I = 1, N X = A * DRANDL(C) Y = X * X Z = DSQRT(Y) W = (Z - X) / X IF (W .GT. ZERO) K1 = K1 + 1 IF (W .LT. ZERO) K3 = K3 + 1 W = DABS(W) IF (W .LE. R6) GO TO 120 R6 = W X1 = X 120 R7 = R7 + W * W 200 CONTINUE C K2 = N - K1 - K3 R7 = DSQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (6,1020) IT,IBETA W = -999.0D+00 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 W = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W W = -999.0D+00 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W W = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,W A = ONE B = SQBETA 300 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1040) X = XMIN Y = DSQRT(X) WRITE (IOUT,1041) X,Y X = ONE - EPSNEG Y = DSQRT(X) WRITE (IOUT,1042) EPSNEG,Y X = ONE Y = DSQRT(X) WRITE (IOUT,1043) X,Y X = ONE + EPS Y = DSQRT(X) WRITE (IOUT,1044) EPS,Y X = XMAX Y = DSQRT(X) WRITE (IOUT,1045) X,Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = ZERO WRITE (IOUT,1051) X Y = DSQRT(X) WRITE (IOUT,1055) Y X = -ONE WRITE (IOUT,1052) X Y = DSQRT(X) WRITE (IOUT,1055) Y WRITE (IOUT,1100) RETURN 1000 FORMAT(24H1TEST OF DSQRT(X*X) - X //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / 1 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(20H DSQRT(X) WAS LARGER,I6,7H TIMES, / 1 13X,7H AGREED,I6,11H TIMES, AND / 2 9X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, 1 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, 1 F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, 1 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, 1 3H = ,I4,3H **,F7.2) 1040 FORMAT(26H1TEST OF SPECIAL ARGUMENTS//) 1041 FORMAT(21H DSQRT(XMIN) = DSQRT(,D25.17,4H) = ,D25.17//) 1042 FORMAT(27H DSQRT(1-EPSNEG) = DSQRT(1-,D15.7,4H) = ,D25.17//) 1043 FORMAT(20H DSQRT(1.0) = DSQRT(,D25.17,4H) = ,D25.17//) 1044 FORMAT(24H DSQRT(1+EPS) = DSQRT(1+,D15.7,4H) = ,D25.17//) 1045 FORMAT(21H DSQRT(XMAX) = DSQRT(,D25.17,4H) = ,D25.17//) 1050 FORMAT(22H1TEST OF ERROR RETURNS//) 1051 FORMAT(39H DSQRT WILL BE CALLED WITH THE ARGUMENT,D15.4/ 1 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(39H0DSQRT WILL BE CALLED WITH THE ARGUMENT,D15.4/ 1 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H DSQRT RETURNED THE VALUE,D15.4///) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DSQRT TEST PROGRAM ---------- END