PROGRAM DTEST C C DRIVER FOR DOUBLE PRECISION MATH FUNCTION TESTS C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) C INTEGER I,J C CALL EMAS3('DEFINE','10,DTESTOUT,,C',J) CALL EMAS3('DEFINE','9,DTESTINPUT',J) CALL SETMESSAGES WRITE(6,10000) 10000 FORMAT(' TESTS ARE NUMBERED AS FOLLOWS:',/ 1 ,' 0 - H/W PARMS 4 - ''**'' 8 - DATAN/DATAN2',/ 2 ,' 1 - DSQRT 5 - DSIN/DCOS 9 - DSINH/DCOSH',/ 3 ,' 2 - DLOG/DLOG10 6 - DTAN/DCOTAN 10 - DTANH',/ 4 ,' 3 - DEXP 7 - DASIN/DACOS 11 - DGAMMA,DLGAMA' 5 ,',DERF,DERFC',/ 6 ,' GIVE TEST NUMBER 0 TO 11, -1 TO STOP, 100 FOR ALL') 99 CONTINUE READ(5,*) I IF (I.EQ.0) GOTO 100 IF (I.EQ.100) GOTO 300 IF (I.EQ.-1) STOP GOTO (101,102,103,104,105,106,107,108,109,110,111), I C 100 CALL DTEST0 GOTO 99 C 105 CALL DTEST5 GOTO 99 C 102 CALL DTEST2 GOTO 99 C 103 CALL DTEST3 GOTO 99 C 101 CALL DTEST1 GOTO 99 C 108 CALL DTEST8 GOTO 99 C 104 CALL DTEST4 GOTO 99 C 106 CALL DTEST6 GOTO 99 C 107 CALL DTEST7 GOTO 99 C 109 CALL DTEST9 GOTO 99 C 110 CALL DTEST10 GOTO 99 C 111 CALL DTEST11 GOTO 99 C 300 CALL DTEST0 CALL DTEST1 CALL DTEST2 CALL DTEST3 CALL DTEST4 CALL DTEST5 CALL DTEST6 CALL DTEST7 CALL DTEST8 CALL DTEST9 CALL DTEST10 CALL DTEST11 STOP END C SUBROUTINE PRINTTITLE(TEXT) CHARACTER*30 TEXT C WRITE(6,10000) TEXT 10000 FORMAT(//,1X,A30 1 ,/,21X,'J MRE RMS MICRO-SECS') RETURN END C SUBROUTINE SUMMARY(J,MRE,RMS,MSECS) INTEGER J DOUBLE PRECISION MRE,RMS,MSECS C WRITE(6,10000) J,MRE,RMS,MSECS 10000 FORMAT(20X,I2,5X,F6.2,5X,F6.2,10X,F6.2) RETURN END C SUBROUTINE TIMERI(J,N,IFUNC,ARGS,MRE,RMS,IOUT) INTEGER J,N,IOUT,IFUNC DOUBLE PRECISION ARGS(10000),ARGS2(10000),MRE,RMS C INTEGER I,K,KK DOUBLE PRECISION T0,T1,T2,T3,Y C KK= 20000/N C CALL EMAS3CPUTIME(T0) DO 1 K= 1,KK DO 400 I= 1,N Y= ARGS(I) 400 CONTINUE 1 CONTINUE C CALL EMAS3CPUTIME(T1) GOTO 200 C ENTRY TIMERI2(J,N,IFUNC,ARGS,ARGS2,MRE,RMS,IOUT) C KK= 20000/N CALL EMAS3CPUTIME(T0) DO 2 K= 1,KK DO 500 I= 1,N Y= ARGS(I) Y= ARGS2(I) 500 CONTINUE 2 CONTINUE C CALL EMAS3CPUTIME(T1) C 200 CONTINUE C GOTO (101,102,103,104,105,106,107,108,109,110,111,112,113,114 * ,115,116,117,118,119,120,121,122,123,124,125), IFUNC C 101 CONTINUE CALL EMAS3CPUTIME(T2) DO 51 K= 1,KK DO 501 I= 1,N Y= DSQRT(ARGS(I)) 501 CONTINUE 51 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 102 CONTINUE CALL EMAS3CPUTIME(T2) DO 52 K= 1,KK DO 502 I= 1,N Y= DLOG(ARGS(I)) 502 CONTINUE 52 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 103 CONTINUE CALL EMAS3CPUTIME(T2) DO 53 K= 1,KK DO 503 I= 1,N Y= DEXP(ARGS(I)) 503 CONTINUE 53 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 104 CONTINUE CALL EMAS3CPUTIME(T2) DO 54 K= 1,KK DO 504 I= 1,N Y= ARGS(I)**ARGS2(I) 504 CONTINUE 54 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 105 CONTINUE CALL EMAS3CPUTIME(T2) DO 55 K= 1,KK DO 505 I= 1,N Y= DSIN(ARGS(I)) 505 CONTINUE 55 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 106 CONTINUE CALL EMAS3CPUTIME(T2) DO 56 K= 1,KK DO 506 I= 1,N Y= DTAN(ARGS(I)) 506 CONTINUE 56 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 107 CONTINUE CALL EMAS3CPUTIME(T2) DO 57 K= 1,KK DO 507 I= 1,N Y= DASIN(ARGS(I)) 507 CONTINUE 57 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 108 CONTINUE CALL EMAS3CPUTIME(T2) DO 58 K= 1,KK DO 508 I= 1,N Y= DATAN(ARGS(I)) 508 CONTINUE 58 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 109 CONTINUE CALL EMAS3CPUTIME(T2) DO 59 K= 1,KK DO 509 I= 1,N Y= DSINH(ARGS(I)) 509 CONTINUE 59 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 110 CONTINUE CALL EMAS3CPUTIME(T2) DO 60 K= 1,KK DO 510 I= 1,N Y= DTANH(ARGS(I)) 510 CONTINUE 60 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 111 CONTINUE CALL EMAS3CPUTIME(T2) DO 61 K= 1,KK C SHOULD NOT GET HERE DO 511 I= 1,N C Y= FUNC(ARGS(I)) 511 CONTINUE 61 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 112 CONTINUE CALL EMAS3CPUTIME(T2) DO 62 K= 1,KK DO 512 I= 1,N Y= DLOG10(ARGS(I)) 512 CONTINUE 62 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 113 CONTINUE CALL EMAS3CPUTIME(T2) DO 63 K= 1,KK DO 513 I= 1,N C SHOULD NOT GET HERE C Y= FUNC(ARGS(I)) 513 CONTINUE 63 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 114 CONTINUE CALL EMAS3CPUTIME(T2) DO 64 K= 1,KK DO 514 I= 1,N C SHOULD NOT GET HERE C Y= FUNC(ARGS(I)) 514 CONTINUE 64 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 115 CONTINUE CALL EMAS3CPUTIME(T2) DO 65 K= 1,KK DO 515 I= 1,N Y= DCOS(ARGS(I)) 515 CONTINUE 65 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 116 CONTINUE CALL EMAS3CPUTIME(T2) DO 66 K= 1,KK DO 516 I= 1,N Y= DCOTAN(ARGS(I)) 516 CONTINUE 66 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 117 CONTINUE CALL EMAS3CPUTIME(T2) DO 67 K= 1,KK DO 517 I= 1,N Y= DACOS(ARGS(I)) 517 CONTINUE 67 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 118 CONTINUE CALL EMAS3CPUTIME(T2) DO 68 K= 1,KK DO 518 I= 1,N Y= DATAN2(ARGS(I),ARGS2(I)) 518 CONTINUE 68 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 119 CONTINUE CALL EMAS3CPUTIME(T2) DO 69 K= 1,KK DO 519 I= 1,N Y= DCOSH(ARGS(I)) 519 CONTINUE 69 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 120 CONTINUE CALL EMAS3CPUTIME(T2) DO 70 K= 1,KK DO 520 I= 1,N Y= DGAMMA(ARGS(I)) 520 CONTINUE 70 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 121 CONTINUE CALL EMAS3CPUTIME(T2) DO 71 K= 1,KK DO 521 I= 1,N Y= DLGAMA(ARGS(I)) 521 CONTINUE 71 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 122 CONTINUE CALL EMAS3CPUTIME(T2) DO 72 K= 1,KK DO 522 I= 1,N Y= DERF(ARGS(I)) 522 CONTINUE 72 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 123 CONTINUE CALL EMAS3CPUTIME(T2) DO 73 K= 1,KK DO 523 I= 1,N Y= DERFC(ARGS(I)) 523 CONTINUE 73 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 124 CONTINUE CALL EMAS3CPUTIME(T2) DO 74 K= 1,KK DO 524 I= 1,N Y= ARGS(I)**NINT(ARGS2(I)) 524 CONTINUE 74 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 125 CONTINUE CALL EMAS3CPUTIME(T2) DO 75 K= 1,KK DO 525 I= 1,N Y= NINT(ARGS(I))**NINT(ARGS2(I)) 525 CONTINUE 75 CONTINUE CALL EMAS3CPUTIME(T3) GOTO 100 C 100 CONTINUE T3= 1.0D6*(T3 - T2 - (T1 - T0))/DBLE(KK*N) C WRITE(IOUT,8000) T3 8000 FORMAT(' AVERAGE EVALUATION TIME: ',G12.5,' MICRO-SECONDS.',//) C CALL SUMMARY(J,MRE,RMS,T3) RETURN END C SUBROUTINE SETMESSAGES C C CHARACTER*55 MESS(40) C C INTEGER DLOGSMALL,DSQRTNEG,DEXPLARGE,DEXPSMALL,DSINLARGE 1 ,DASINLARGE,DCOSLARGE,DACOSLARGE,DTANLARGE,DCOTANSMALL 2 ,DARCSINLARGE,DARCCOSLARGE,DARCTAN2ZERO,DSINHLARGE,DCOSHLARGE 3 ,DPOWERNEG,DPOWERZERO 4 ,DGAMLARGE,DGAMNEGINT,DLGAMNEG,DLGAMLARGE 5 ,IPOWERZERO,IPOWEXPNEG,IPOWLARGE C DLOGSMALL= 1 DSQRTNEG= 2 DEXPLARGE= 3 DEXPSMALL= 4 DSINLARGE= 5 DASINLARGE= 6 DCOSLARGE= 7 DACOSLARGE= 8 DTANLARGE= 9 DCOTANSMALL= 10 DARCSINLARGE= 11 DARCCOSLARGE= 12 DARCTAN2ZERO= 13 DSINHLARGE= 14 DCOSHLARGE= 15 DPOWERNEG= 16 DPOWERZERO= 17 C DGAMLARGE= 34 DGAMNEGINT= 35 DLGAMNEG= 36 DLGAMLARGE= 37 IPOWEXPNEG= 38 IPOWERZERO= 39 IPOWLARGE= 40 MESS(DSQRTNEG)= 'DSQRT ARG NEGATIVE' MESS(DLOGSMALL)= 'DLOG ARG NEGATIVE OR ZERO.' MESS(DEXPLARGE)= 'DEXP ARG GREATER THAN 174.6731' MESS(DEXPSMALL)= 'DEXP ARG LESS THAN -180.2182' MESS(DPOWERNEG)= 1 'NEGATIVE D.P. VALUE RAISED TO A NON-INTEGER POWER' MESS(DPOWERZERO)= 'D.P. ZERO RAISED TO NON-POSITIVE POWER' C C NOTE THIS LAST ERROR IS NOT USED. DEFAULT RESULT ZERO IS ASSIGNED. C MESS(DSINLARGE)= 'DSIN ARG MAGNITUDE GREATER THAN 7.205D16' MESS(DCOSLARGE)= 'DCOS ARG MAGNITUDE GREATER THAN 7.205D16' MESS(DTANLARGE)= * 'DTAN/DCOTAN ARG MAGNITUDE GREATER THAN 3.521D15' MESS(DCOTANSMALL)= 'DCOTAN ARG SMALL IN MAGNITUDE' MESS(DASINLARGE)= 'DASIN ARG MAGNITUDE GREATER THAN 1.0' MESS(DACOSLARGE)= 'DACOS ARG MAGNITUDE GREATER THAN 1.0' MESS(DARCTAN2ZERO)= 'DATAN2 ARGS ZERO.' MESS(DSINHLARGE)= 'DSINH ARG MAGNITUDE GREATER THAN 175.3662' MESS(DCOSHLARGE)= 'DCOSH ARG MAGNITUDE GREATER THAN 175.3662' C MESS(DGAMLARGE)= 'DGAMMA ARG MAGNITUDE GREATER THAN 57.0' MESS(DGAMNEGINT)= 'DGAMMA ARG NEAR ZERO OR NEGATIVE INTEGER' MESS(DLGAMNEG)= 'DLGAMA ARG IS NEGATIVE' MESS(DLGAMLARGE)= 'DLGAMA ARG IS GREATER THAN 4.29D73' MESS(IPOWERZERO)= 'INTEGER ZERO RAISED TO ZERO INTEGER POWER' MESS(IPOWEXPNEG)= 'INTEGER RAISED TO NEGATIVE INTEGER POWER' MESS(IPOWLARGE)= 'INTEGER TO INTEGER POWER OVERFLOWS' RETURN C ENTRY TERROR(I) IOUT= 10 WRITE(IOUT,10000) I,MESS(I) IF (IOUT.NE.6) WRITE(6,10001) I,MESS(I) 10000 FORMAT(' ** ERROR ',I2,': ',A55,/) 10001 FORMAT(/' ** ERROR ',I2,': ',A55) RETURN END C C******************************************************************* C SUBROUTINE MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP, *MAXEXP,EPS,EPSNEG,XMIN,XMAX) C INTEGER I,IBETA,IEXP,IRND,IT,IZ,J,K,MACHEP,MAXEXP,MINEXP,MX,NEGEP, *NGRD DOUBLE PRECISION A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN, *Y,Z,ZERO CD DOUBLE PRECISION A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX, CD 1 XMIN,Y,Z,ZERO C C THIS SUBROUTINE IS INTENDED TO DETERMINE THE CHARACTERISTICS C OF THE FLOATING-POINT ARITHMETIC SYSTEM THAT ARE SPECIFIED C BELOW. THE FIRST THREE ARE DETERMINED ACCORDING TO AN C ALGORITHM DUE TO M. MALCOLM, CACM 15 (1972), PP. 949-951, C INCORPORATING SOME, BUT NOT ALL, OF THE IMPROVEMENTS C SUGGESTED BY M. GENTLEMAN AND S. MAROVICH, CACM 17 (1974), C PP. 276-277. THE VERSION GIVEN HERE IS FOR SINGLE PRECISION. C CARDS CONTAINING CD IN COLUMNS 1 AND 2 CAN BE USED TO C CONVERT THE SUBROUTINE TO DOUBLE PRECISION BY REPLACING C EXISTING CARDS IN THE OBVIOUS MANNER. C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C C IBETA - THE RADIX OF THE FLOATING-POINT REPRESENTATION C IT - THE NUMBER OF BASE IBETA DIGITS IN THE FLOATING-POINT C SIGNIFICAND C IRND - 0 IF FLOATING-POINT ADDITION CHOPS, C 1 IF FLOATING-POINT ADDITION ROUNDS C NGRD - THE NUMBER OF GUARD DIGITS FOR MULTIPLICATION. IT IS C 0 IF IRND=1, OR IF IRND=0 AND ONLY IT BASE IBETA C DIGITS PARTICIPATE IN THE POST NORMALIZATION SHIFT C OF THE FLOATING-POINT SIGNIFICAND IN MULTIPLICATION C 1 IF IRND=0 AND MORE THAN IT BASE IBETA DIGITS C PARTICIPATE IN THE POST NORMALIZATION SHIFT OF THE C FLOATING-POINT SIGNIFICAND IN MULTIPLICATION C MACHEP - THE LARGEST NEGATIVE INTEGER SUCH THAT C 1.0+FLOAT(IBETA)**MACHEP .NE. 1.0, EXCEPT THAT C MACHEP IS BOUNDED BELOW BY -(IT+3) C NEGEPS - THE LARGEST NEGATIVE INTEGER SUCH THAT C 1.0-FLOAT(IBETA)**NEGEPS .NE. 1.0, EXCEPT THAT C NEGEPS IS BOUNDED BELOW BY -(IT+3) C IEXP - THE NUMBER OF BITS (DECIMAL PLACES IF IBETA = 10) C RESERVED FOR THE REPRESENTATION OF THE EXPONENT C (INCLUDING THE BIAS OR SIGN) OF A FLOATING-POINT C NUMBER C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE INTEGER SUCH THAT C FLOAT(IBETA)**MINEXP IS A POSITIVE FLOATING-POINT C NUMBER C MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT FOR A FINITE C FLOATING-POINT NUMBER C EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH C THAT 1.0+EPS .NE. 1.0. IN PARTICULAR, IF EITHER C IBETA = 2 OR IRND = 0, EPS = FLOAT(IBETA)**MACHEP. C OTHERWISE, EPS = (FLOAT(IBETA)**MACHEP)/2 C EPSNEG - A SMALL POSITIVE FLOATING-POINT NUMBER SUCH THAT C 1.0-EPSNEG .NE. 1.0. IN PARTICULAR, IF IBETA = 2 C OR IRND = 0, EPSNEG = FLOAT(IBETA)**NEGEPS. C OTHERWISE, EPSNEG = (IBETA**NEGEPS)/2. BECAUSE C NEGEPS IS BOUNDED BELOW BY -(IT+3), EPSNEG MAY NOT C BE THE SMALLEST NUMBER WHICH CAN ALTER 1.0 BY C SUBTRACTION. C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT POWER OF THE C RADIX. IN PARTICULAR, XMIN = FLOAT(IBETA)**MINEXP C XMAX - THE LARGEST FINITE FLOATING-POINT NUMBER. IN C PARTICULAR XMAX = (1.0-EPSNEG)*FLOAT(IBETA)**MAXEXP C NOTE - ON SOME MACHINES XMAX WILL BE ONLY THE C SECOND, OR PERHAPS THIRD, LARGEST NUMBER, BEING C TOO SMALL BY 1 OR 2 UNITS IN THE LAST DIGIT OF C THE SIGNIFICAND. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C----------------------------------------------------------------- C CALL ICL9HEMASK(64,IMPH) ONE = DBLE(FLOAT(1)) CD ONE = DBLE(FLOAT(1)) ZERO = 0.0D0 CD ZERO = 0.0D0 C----------------------------------------------------------------- C DETERMINE IBETA,BETA ALA MALCOLM C----------------------------------------------------------------- A = ONE 10 A = A + A IF (((A+ONE)-A)-ONE .EQ. ZERO) GO TO 10 B = ONE 20 B = B + B IF ((A+B)-A .EQ. ZERO) GO TO 20 IBETA = IDINT((A+B)-A) C IBETA=16 CD IBETA = INT(SNGL((A + B) - A)) BETA = DBLE(FLOAT(IBETA)) CD BETA = DBLE(FLOAT(IBETA)) C----------------------------------------------------------------- C DETERMINE IT, IRND C----------------------------------------------------------------- IT = 0 B = ONE 100 IT = IT + 1 B = B * BETA IF (((B+ONE)-B)-ONE .EQ. ZERO) GO TO 100 IRND = 0 BETAM1 = BETA - ONE IF ((A+BETAM1)-A .NE. ZERO) IRND = 1 C IT=28 C IRND=0 C NGRD=1 C----------------------------------------------------------------- C DETERMINE NEGEP, EPSNEG C----------------------------------------------------------------- NEGEP = IT + 3 BETAIN = ONE / BETA A = ONE C DO 200 I = 1, NEGEP A = A * BETAIN 200 CONTINUE C B = A 210 IF ((ONE-A)-ONE .NE. ZERO) GO TO 220 A = A * BETA NEGEP = NEGEP - 1 GO TO 210 220 NEGEP = -NEGEP EPSNEG = A IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 300 A = (A*(ONE+A)) / (ONE+ONE) IF ((ONE-A)-ONE .NE. ZERO) EPSNEG = A C----------------------------------------------------------------- C DETERMINE MACHEP, EPS C----------------------------------------------------------------- 300 MACHEP = -IT - 3 A = B 310 IF((ONE+A)-ONE .NE. ZERO) GO TO 320 A = A * BETA MACHEP = MACHEP + 1 GO TO 310 320 EPS = A IF ((IBETA .EQ. 2) .OR. (IRND .EQ. 0)) GO TO 350 A = (A*(ONE+A)) / (ONE+ONE) IF ((ONE+A)-ONE .NE. ZERO) EPS = A C----------------------------------------------------------------- C DETERMINE NGRD C----------------------------------------------------------------- 350 NGRD = 0 IF ((IRND .EQ. 0) .AND. ((ONE+EPS)*ONE-ONE) .NE. ZERO) NGRD = 1 C----------------------------------------------------------------- C DETERMINE IEXP, MINEXP, XMIN C C LOOP TO DETERMINE LARGEST I AND K = 2**I SUCH THAT C (1/BETA) ** (2**(I)) C DOES NOT UNDERFLOW C EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW. C----------------------------------------------------------------- I = 0 K = 1 Z = BETAIN 400 Y = Z Z = Y * Y C----------------------------------------------------------------- C CHECK FOR UNDERFLOW HERE C----------------------------------------------------------------- A = Z * ONE IF ((A+A .EQ. ZERO) .OR. (DABS(Z) .GE. Y)) GO TO 410 CD IF ((A+A .EQ. ZERO) .OR. (DABS(Z) .GE. Y)) GO TO 410 I = I + 1 K = K + K GO TO 400 410 IF (IBETA .EQ. 10) GO TO 420 IEXP = I + 1 MX = K + K GO TO 450 C----------------------------------------------------------------- C FOR DECIMAL MACHINES ONLY C----------------------------------------------------------------- 420 IEXP = 2 IZ = IBETA 430 IF (K .LT. IZ) GO TO 440 IZ = IZ * IBETA IEXP = IEXP + 1 GO TO 430 440 MX = IZ + IZ - 1 C----------------------------------------------------------------- C LOOP TO DETERMINE MINEXP, XMIN C EXIT FROM LOOP IS SIGNALED BY AN UNDERFLOW. C----------------------------------------------------------------- 450 XMIN = Y Y = Y * BETAIN C----------------------------------------------------------------- C CHECK FOR UNDERFLOW HERE C----------------------------------------------------------------- A = Y * ONE IF (((A+A) .EQ. ZERO) .OR. (DABS(Y) .GE. XMIN)) GO TO 460 CD IF (((A+A) .EQ. ZERO) .OR. (DABS(Y) .GE. XMIN)) GO TO 460 K = K + 1 GO TO 450 460 MINEXP = -K C----------------------------------------------------------------- C DETERMINE MAXEXP, XMAX C----------------------------------------------------------------- IF ((MX .GT. K+K-3) .OR. (IBETA .EQ. 10)) GO TO 500 MX = MX + MX IEXP = IEXP + 1 500 MAXEXP = MX + MINEXP C----------------------------------------------------------------- C ADJUST FOR MACHINES WITH IMPLICIT LEADING C BIT IN BINARY SIGNIFICAND AND MACHINES WITH C RADIX POINT AT EXTREME RIGHT OF SIGNIFICAND C----------------------------------------------------------------- I = MAXEXP + MINEXP IF ((IBETA .EQ. 2) .AND. (I .EQ. 0)) MAXEXP = MAXEXP - 1 IF (I .GT. 20) MAXEXP = MAXEXP - 1 IF (A .NE. Y) MAXEXP = MAXEXP - 2 XMAX = ONE - EPSNEG IF (XMAX*ONE .NE. XMAX) XMAX = ONE - BETA * EPSNEG XMAX = XMAX / (BETA * BETA * BETA * XMIN) I = MAXEXP + MINEXP + 3 IF (I .LE. 0) GO TO 520 C DO 510 J = 1, I IF (IBETA .EQ. 2) XMAX = XMAX + XMAX IF (IBETA .NE. 2) XMAX = XMAX * BETA 510 CONTINUE C 520 RETURN C ---------- LAST CARD OF MACHAR ---------- END C C******************************************************************* C DOUBLE PRECISION FUNCTION RAN(I) INTEGER I C THIS FUNCTION MUST RETURN A RANDOM NUMBER BETWEEN 0 AND 1 C INTEGER IY,J DATA IY/100001/ C J= I IY= IY*125 IY= IY - (IY/2796203)*2796203 RAN= FLOAT(IY)/2796203.0D0 C RETURN END C C******************************************************************* C SUBROUTINE DTEST0 C PROGRAM TO MONITOR SUBROUTINE MACHAR, DOUBLE PRECISION C C DATA REQUIRED C C NONE C C INTEGER IBETA, IEXP, IRND, IT, MACHEP, MAXEXP, MINEXP,NEGEP, NGRD INTEGER IOUT DOUBLE PRECISION EPS, EPSNEG, XMIN, XMAX C C OUTPUT UNIT SET TO 6 C IOUT = 10 CALL MACHAR(IBETA, IT, IRND, NGRD, MACHEP, NEGEP, IEXP,MINEXP, * MAXEXP, EPS, EPSNEG, XMIN, XMAX) WRITE(6,2000) IT,IBETA 2000 FORMAT(/,' THE TESTS USE ',I3,' BASE ',I2,' DIGITS.',/) WRITE (IOUT, 1000) WRITE (IOUT, 1001) IBETA, IT, IRND, NGRD, MACHEP, NEGEP,IEXP, * MINEXP, MAXEXP WRITE (IOUT, 1002) EPS, EPSNEG, XMIN, XMAX WRITE (IOUT, 1003) EPS, EPSNEG, XMIN, XMAX C C THE OUTPUT FOR REALS IS REPEATED SO THAT, IF REQUIRED, C VALUES IN BINARY/HEXADECIMAL CAN ALSO BE OBTAINED C RETURN C C IF REQUIRED, MODIFY THIS FORMAT STATEMENT TO INCLUDE C DETAILS OF THE MACHINE YOU ARE TESTING C 1000 FORMAT (29H OUTPUT OF VALUES FROM MACHAR/ * 52H FIRST LINE IS IBETA, IT, IRND, NGRD, MACHEP, NEGEP, *, 22H ,IEXP, MINEXP, MAXEXP/ * 39H SECOND LINE IS EPS, EPSNEG, XMIN, XMAX/ * 24H THIRD LINE IS AS SECOND/) 1001 FORMAT (1X, 9I8) 1002 FORMAT (1X, 4D20.13) 1003 FORMAT (1X,4(1X,Z16,3X)) C C IF REQUIRED MODIFY THIS FORMAT C END C C************************************************************** C SUBROUTINE DTEST1 C PROGRAM TO TEST DSQRT C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) 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 RANDL(X) - A FUNCTION SUBPROGRAM RETURNING LOGARITHMICALLY C DISTRIBUTED RANDOM REAL NUMBERS. IN PARTICULAR, C A * RANDL(ALOG(B/A)) C IS LOGARITHMICALLY DISTRIBUTED OVER (A,B) C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,J,K1,K2,K3,MACHEP,MAXEXP,MINEXP, *N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,RANDL,R6,R7, *SQBETA,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO,DSQRT,MRE,RMS DOUBLE PRECISION RAN C C THE FOLLOWING STATEMENT FUNCTION HAS BEEN ADDED TO C PROVIDE THE REQUIRED FUNCTION RANDL C RANDL(X) = DEXP(RAN(0)*X) C CALL PRINTTITLE('DSQRT ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) SQBETA = DSQRT(BETA) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0D0 ZERO = 0.0D0 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 * RANDL(C) Y = X * X Z = DSQRT(Y) ARGS(I)= 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 (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W= DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = ONE B = SQBETA C C CPU TIME CALCULATION AND SUMMARY PRINT C CALL TIMERI(J,N,1,ARGS,MRE,RMS,IOUT) 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(1X,80('*') 1 ,/,' TEST OF DSQRT(X*X) - X '//) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(20H DSQRT(X) WAS LARGER,I6,7H TIMES, / * 13X,7H AGREED,I6,11H TIMES, AND / * 9X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1040 FORMAT(26H TEST OF SPECIAL ARGUMENTS//) 1041 FORMAT(21H DSQRT(XMIN) = DSQRT(,D15.7,4H) = ,D15.7//) 1042 FORMAT(27H DSQRT(1-EPSNEG) = DSQRT(1-,D15.7,4H) = ,D15.7//) 1043 FORMAT(20H DSQRT(1.0) = DSQRT(,D15.7,4H) = ,D15.7//) 1044 FORMAT(24H DSQRT(1+EPS) = DSQRT(1+,D15.7,4H) = ,D15.7//) 1045 FORMAT(21H DSQRT(XMAX) = DSQRT(,D15.7,4H) = ,D15.7//) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(39H DSQRT WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(39H0DSQRT WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 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 C C*********************************************************** C SUBROUTINE DTEST5 C PROGRAM TO TEST DSIN/DCOS C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) 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 FIVE 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 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE 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 C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DCOS, DSIN, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,BETAP,C,DEL,EPS,EPSNEG,EXPON *,ONE,R6,R7,THREE,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ * ,MRE,RMS C DOUBLE PRECISION RAN DOUBLE PRECISION DSIN,DCOS C C=================================================================== C CALL PRINTTITLE('DSIN (J=1,2), DCOS (J=3) ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0D0 ZERO = 0.0D0 THREE = 3.0D0 A = ZERO B = 1.570796327D0 C = B N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 300 J = 1, 3 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL Y = X / THREE Y = (X + Y) - X X = THREE * Y ARGS(I)= Y IF (J .EQ. 3) GO TO 100 Z = DSIN(X) ZZ = DSIN(Y) W = ONE IF (Z .NE. ZERO) W = (Z - ZZ*(THREE-4.0D0*ZZ*ZZ)) / Z GO TO 110 100 Z = DCOS(X) ZZ = DCOS(Y) W = ONE IF (Z .NE. ZERO) W = (Z + ZZ*(THREE-4.0D0*ZZ*ZZ)) / Z 110 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 XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (J .EQ. 3) GO TO 210 WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 GO TO 220 210 WRITE (IOUT,1005) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1012) K1,K2,K3 220 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = 18.84955592D0 IF (J .EQ. 2) A = B + C B = A + C C C CPU TIME CALCULATION AND SUMMARY PRINT C IF (J.NE.3) THEN CALL TIMERI(J,N,5,ARGS,MRE,RMS,IOUT) ELSE CALL TIMERI(J,N,15,ARGS,MRE,RMS,IOUT) ENDIF C C 300 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) C = ONE / BETA ** (IT/2) Z = (DSIN(A+C) - DSIN(A-C)) / (C + C) WRITE (IOUT,1026) Z C WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1) * A Z = DSIN(X) + DSIN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = RAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - DSIN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) C DO 340 I = 1, 5 X = RAN(I1) * A Z = DCOS(X) - DCOS(-X) WRITE (IOUT,1060) X, Z 340 CONTINUE C WRITE (IOUT,1035) EXPON = DBLE(FLOAT(MINEXP)) * 0.75D0 X = BETA ** EXPON Y = DSIN(X) WRITE (IOUT,1061) X, Y WRITE (IOUT,1040) Z = DSQRT(BETAP) X = Z * (ONE - EPSNEG) Y = DSIN(X) WRITE (IOUT,1061) X, Y Y = DSIN(Z) WRITE (IOUT,1061) Z, Y X = Z * (ONE + EPS) Y = DSIN(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = BETAP WRITE (IOUT,1052) X Y = DSIN(X) WRITE (IOUT,1055) Y WRITE(IOUT,1072) X Y= DCOS(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,' TEST OF DSIN(X) VS 3*DSIN(X/3)-4*DSIN(X/3)**3' * ,//) 1005 FORMAT(' TEST OF DCOS(X) VS 4*DCOS(X/3)**3-3*DCOS(X/3)' * ,//) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(19H DSIN(X) WAS LARGER,I6,7H TIMES, / * 12X,7H AGREED,I6,11H TIMES, AND / * 8X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(19H DCOS(X) WAS LARGER,I6,7H TIMES, / * 12X,7H AGREED,I6,11H TIMES, AND / * 8X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1026 FORMAT(4H IF ,D13.6,21H IS NOT ALMOST 1.0D0,, * 4X,26HDSIN HAS THE WRONG PERIOD. //) 1030 FORMAT(53H THE IDENTITY DSIN(-X) = -DSIN(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(52H THE IDENTITY DSIN(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(52H THE IDENTITY DCOS(-X) = DCOS(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) - F(-X)/) 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT.) 1040 FORMAT(49H THE FOLLOWING THREE LINES ILLUSTRATE THE LOSS IN, * 13H SIGNIFICANCE/36H FOR LARGE ARGUMENTS. THE ARGUMENTS, * 17H ARE CONSECUTIVE.) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1052 FORMAT(38H DSIN WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H DSIN RETURNED THE VALUE,D15.4///) 1060 FORMAT(2D15.7/) 1072 FORMAT(38H DCOS WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(24H DCOS RETURNED THE VALUE,D15.4///) 1061 FORMAT(/6X,6H DSIN(,D13.6,3H) =,D13.6) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DSIN/DCOS TEST PROGRAM ---------- END C C***************************************************************** C SUBROUTINE DTEST2 C PROGRAM TO TEST DLOG C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) 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 FOUR 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 XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C POWER OF THE RADIX C XMAX - THE LARGEST FINITE FLOATING-POINT NO. C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DLOG10, DMAX1, DSIGN, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,C,DEL,EIGHT,EPS,EPSNEG,HALF, *ONE,R6,R7,TENTH,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ C DOUBLE PRECISION RAN DOUBLE PRECISION MRE,RMS,DLOG,DLOG10 DOUBLE PRECISION YY,XX,CC1,CC2 C DATA CC1/Z3D5186008B15330C/, CC2/Z3E18B6050C56E8DC/ C CALL PRINTTITLE('DLOG/DLOG10 ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) J = IT / 2 C C NOTE THAT TEXT AND TEST PROGRAM DISSAGREE IN CODY AND WAITE. C TEXT SAYS B**(-T/2) AND CODE SAYS J= IT/3 . C THE TEST RESULTS FOR /3 CANNOT AGREE WITH TEST RESULTS IN CODY AND WAITE. C ZERO = 0.0D0 HALF = 0.5D0 EIGHT = 8.0D0 TENTH = 0.1D0 ONE = 1.0D0 C = ONE C DO 50 I = 1, J C = C / BETA 50 CONTINUE C B = ONE + C A = ONE - C N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 301 J = 1, 4 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C AINCR=0.01 C X=0.5 C DO 200 I = 1, N C X=X+AINCR X= DEL*RAN(I1) + XL ARGS(I)= X IF (J.NE.1) GOTO 100 Y = (X - HALF) - HALF C Y= Y - HALF ZZ = DLOG(X) C XX= X C XX= LOG(XX) C X= XX C YY= XX + (XX - X) C Z= YY Z = ONE / 3.0D0 Z = Y * (Z - Y / 4.0D0) Z = (Z - HALF) * Y * Y + Y GOTO 150 C 100 IF (J.NE.2) GOTO 110 X= (X + EIGHT) X= X - EIGHT Y= X + X/16.0D0 Z= DLOG(X) C ZZ= DLOG(Y) - 7.7746816434842581D-5 ZZ= DLOG(Y) - CC1 ZZ= ZZ - 31.0D0/512.0D0 GOTO 150 C 110 IF (J.NE.3) GOTO 120 X= (X + EIGHT) X= X - EIGHT Y= X + X*TENTH Z= DLOG10(X) C XX= 3.7706015822504075D-4 C ZZ= XX ZZ= DLOG10(Y) - CC2 ZZ= ZZ - 21.0D0/512.0D0 GOTO 150 C 120 CONTINUE Z= DLOG(X*X) ZZ= DLOG(X) ZZ= ZZ + ZZ C 150 CONTINUE W = ONE ZP = Z IF (Z .NE. ZERO) W = (Z - ZZ) / Z Z = DSIGN(W,Z) C WRITE(6,1)X,Y,ZP,ZZ,W 1 FORMAT(1X,'X= ',E15.7,'Y= ',E15.7,'ZP= ',E15.7,'ZZ= ',E15.7, *'W= ',E15.7) IF (Z .GT. ZERO) K1 = K1 + 1 IF (Z .LT. ZERO) K3 = K3 + 1 W = DABS(W) IF (W .LE. R6) GO TO 160 R6 = W X1 = X 160 R7 = R7 + W*W XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (J .EQ. 1) WRITE (IOUT,1000) IF (J .EQ. 2) WRITE (IOUT,1001) IF (J .EQ. 3) WRITE (IOUT,1005) IF (J .EQ. 4) WRITE (IOUT,1002) IF (J .EQ. 1) WRITE (IOUT,1009) N,C IF (J .NE. 1) WRITE (IOUT,1010) N,A,B IF (J .NE. 3) WRITE (IOUT,1011) K1,K2,K3 IF (J .EQ. 3) WRITE (IOUT,1012) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .GT. 1) GO TO 230 A = DSQRT(HALF) B = 15.0D0 / 16.0D0 GO TO 300 230 IF (J .GT. 2) GO TO 240 A = DSQRT(TENTH) B = 0.9D0 GO TO 300 240 A = 16.0D0 B = 240.0D0 300 CONTINUE C C CPU TIME CALCULATION C IF (J.NE.3) THEN CALL TIMERI(J,N,2,ARGS,MRE,RMS,IOUT) ELSE CALL TIMERI(J,N,12,ARGS,MRE,RMS,IOUT) ENDIF 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1) X = X + X + 15.0D0 Y = ONE / X Z = DLOG(X) + DLOG(Y) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1040) X = ONE Y = DLOG(X) WRITE (IOUT,1041) Y X = XMIN Y = DLOG(X) WRITE (IOUT,1042) X, Y X = XMAX Y = DLOG(X) WRITE (IOUT,1043) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = -2.0D0 WRITE (IOUT,1052) X Y = DLOG(X) WRITE (IOUT,1055) Y X = ZERO WRITE (IOUT,1052) X Y = DLOG(X) WRITE (IOUT,1055) Y WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,' TEST OF DLOG(X) VS T.S. EXPANSION OF DLOG(1+Y)'//) 1001 FORMAT(' TEST OF DLOG(X) VS DLOG(17X/16)-DLOG(17/16)' //) 1002 FORMAT(' TEST OF DLOG(X*X) VS 2 * DLOG(X)' //) 1005 FORMAT(' TEST OF DLOG10(X) VS DLOG10(11X/10)' 1 ,'-DLOG10(11/10)' //) 1009 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,26H(1-EPS,1+EPS), WHERE EPS =, D15.4//) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(' DLOG(X) WAS LARGER',I6,7H TIMES, / * 15X,7H AGREED,I6,11H TIMES, AND / * 11X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(' DLOG10(X) WAS LARGER',I6,7H TIMES, / * 17X,7H AGREED,I6,11H TIMES, AND / * 13X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(52H THE IDENTITY DLOG(X) = -DLOG(1/X) WILL BE TESTED.// * 8X,1HX,9X,13HF(X) + F(1/X)/) 1040 FORMAT(//26H TEST OF SPECIAL ARGUMENTS //) 1041 FORMAT(13H DLOG(1.0) = ,D15.7//) 1042 FORMAT(19H DLOG(XMIN) = DLOG(,D15.7,4H) = ,D15.7//) 1043 FORMAT(19H DLOG(XMAX) = DLOG(,D15.7,4H) = ,D15.7//) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1052 FORMAT(38H DLOG WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H DLOG RETURNED THE VALUE,D15.4///) 1060 FORMAT(2D15.7/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DLOG/DLOG10 TEST PROGRAM ---------- END C C***************************************************************** C SUBROUTINE DTEST3 C PROGRAM TO TEST DEXP C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) 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 FOUR 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 XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C POWER OF THE RADIX C XMAX - THE LARGEST FINITE FLOATING-POINT NO. C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DEXP, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C NOTE THAT THE FUNCTION DINT IS USED. THIS IS NOT IN THE C 1966 STANDARD, AND SPECIAL ACTION WILL BE NECESSARY IF THE C COMPILER LIBRARY DOES NOT PROVIDE DINT. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,D,DEL,EPS,EPSNEG,ONE,R6, *R7,TWO,TEN,V,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ,XXRND * ,MRE,RMS C DOUBLE PRECISION RAN DOUBLE PRECISION DEXP C EXTERNAL DEXP DOUBLE PRECISION XX,YY,CC1,CC2 C DATA CC1/Z3EA041DD6AE20424/, CC2/Z3FF82A021C7EAE19/ C CALL PRINTTITLE('DEXP ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0D0 TWO = 2.0D0 TEN = 10.0D0 ZERO = 0.0D0 V = 0.0625D0 A = TWO B = DLOG(A) * 0.5D0 A = -B + V C C D = DLOG(0.9D0*XMAX) C D = ALOG(0.1E0*XMAX) N = 2000 C READ (5,*)N XN = DBLE(FLOAT(N)) I1 = 0 C--------------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C--------------------------------------------------------------------- DO 301 J = 1, 3 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL C--------------------------------------------------------------------- C PURIFY ARGUMENTS C--------------------------------------------------------------------- Y = X - V IF (Y .LT. ZERO) X = Y + V ARGS(I)= X Z = DEXP(X) ZZ = DEXP(Y) IF (J .EQ. 1) GO TO 100 IF (IBETA .EQ. 10) GOTO 302 C XX= 2.4453321046920570389D-3 C XXRND= XX Z = Z * .0625D0 -Z *CC1 GOTO 110 C 302 CONTINUE IF (IBETA .EQ. 10) Z = Z * 6.0D-2 +Z * * 5.466789530794296106D-5 GO TO 110 100 CONTINUE C XX= 6.058693718652421388D-2 C XXRND= XX Z = Z - Z * CC2 110 W = ONE IF (ZZ .NE. ZERO) W = (Z - ZZ) / ZZ IF (W .LT. ZERO) K1 = K1 + 1 IF (W .GT. ZERO) K3 = K3 + 1 W = DABS(W) IF (W .LE. R6) GO TO 120 R6 = W X1 = X 120 R7 = R7 + W*W XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) WRITE (IOUT,1000) V, V WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .EQ. 2) GO TO 270 V = 45.0D0 / 16.0D0 A = -TEN * B B = 4.0D0 * XMIN * BETA ** IT B = DLOG(B) GO TO 300 270 A = -TWO * A B = TEN * A IF (B .LT. D) B = D 300 CONTINUE C C CPU TIME CALCULATION C CALL TIMERI(J,N,3,ARGS,MRE,RMS,IOUT) C 301 CONTINUE C--------------------------------------------------------------------- C SPECIAL TESTS C--------------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1) * BETA Y = -X Z = DEXP(X) * DEXP(Y) - ONE WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1040) X = ZERO Y = DEXP(X) - ONE WRITE (IOUT,1041) Y X = DINT(DLOG(XMIN)) Y = DEXP(X) WRITE (IOUT,1042) X, Y X = DINT(DLOG(XMAX)) C WRITE(6,1)X 1 FORMAT(1X,'X= ',E15.7) Y = DEXP(X) WRITE (IOUT,1042) X, Y X = X / TWO V = X / TWO Y = DEXP(X) Z = DEXP (V) Z = Z * Z WRITE (IOUT,1043) X, Y, V, Z C--------------------------------------------------------------------- C TEST OF ERROR RETURNS C--------------------------------------------------------------------- WRITE (IOUT,1050) X = -ONE / DSQRT(XMIN) WRITE (IOUT,1052) X Y = DEXP(X) WRITE (IOUT,1061) Y X = -X WRITE (IOUT,1052) X Y = DEXP(X) WRITE (IOUT,1061) Y WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,16H TEST OF DEXP(X-,F7.4,18H) VS DEXP(X)/DEXP(,F7.4,1H) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(21H DEXP(X-V) WAS LARGER,I6,7H TIMES, / * 14X,7H AGREED,I6,11H TIMES, AND / * 10X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(54H THE IDENTITY DEXP(X)*DEXP(-X) = 1.0 WILL BE TESTED.// * 8X,1HX,9X,14HF(X)*F(-X) - 1 /) 1040 FORMAT(//26H TEST OF SPECIAL ARGUMENTS //) 1041 FORMAT(21H DEXP(0.0) - 1.0D0 = ,D15.7/) 1042 FORMAT(6H DEXP(,D13.6,3H) =,D13.6/) 1043 FORMAT(9H0IF DEXP(,D13.6,4H) = ,D13.6,13H IS NOT ABOUT / * 6H DEXP(,D13.6,7H)**2 = ,D13.6,26H THERE IS AN ARG RED ERROR) 1050 FORMAT(22H TEST OF ERROR RETURNS //) 1052 FORMAT(38H0DEXP WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1060 FORMAT(2D15.7/) 1061 FORMAT(24H DEXP RETURNED THE VALUE,D15.4///) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DEXP TEST PROGRAM ---------- END C C***************************************************************** C SUBROUTINE DTEST4 C PROGRAM TO TEST POWER FUNCTION (**), DOUBLE PRECISION C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) 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 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C MAXEXP - THE LARGEST POSITIVE INTEGER EXPONENT C FOR A FINITE FLOATING-POINT NUMBER C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C POWER OF THE RADIX C XMAX - THE LARGEST FINITE FLOATING-POINT C NUMBER C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DEXP, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,ALXMAX,B,BETA,C,DEL,DELY,EPS,EPSNEG, *ONE,ONEP5,R6,R7,SCALE,TWO,W,X,XL,XMAX,XMIN,XN,XSQ,X1,Y,Y1,Y2, *Z,ZERO,ZZ,LNINTB * ,MRE,RMS C DOUBLE PRECISION RAN,POWQI CALL PRINTTITLE('X**Y, X**N, I**N ') IOUT= 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) C C ALXMAX = DLOG(XMAX) C ALXMAX = ALOG(XMAX*0.1) ZERO = 0.0D0 ONE = DBLE(FLOAT(1)) TWO = ONE + ONE ONEP5 = (TWO + ONE) / TWO SCALE = ONE J = (IT+1) / 2 C DO 20 I = 1, J SCALE = SCALE * BETA 20 CONTINUE C A = ONE / BETA B = ONE C = -DMAX1(ALXMAX,-DLOG(XMIN))/DLOG(100.0D0) C C= -10 DELY = -C - C N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 Y1 = ZERO LNINTB= LOG(2.0D0**31) C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 301 J = 1,6 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL ARGS(I)= X IF (J .NE. 1) GO TO 50 ARGS2(I)= ONE C ZZ = X ** ONE ZZ= (X)**(ONE) Z = X GO TO 110 50 W = SCALE * X IF (J.EQ.5) W = W *BETA*BETA X = (X + W) X= X - W IF (X.EQ.ZERO) X= ONE XSQ = X * X IF (J .GE. 4) GO TO 70 ARGS(I)= XSQ ARGS2(I)= ONEP5 C ZZ = XSQ ** ONEP5 ZZ= (XSQ)**(ONEP5) Z = X * XSQ GO TO 110 70 CONTINUE ARGS(I)= X C 1.0D0 + 1D-14 Y = DELY * RAN(I1) + C IF (J.GT.4) GOTO 5110 Y2 = (Y/TWO + Y) - Y Y = Y2 + Y2 ARGS2(I)= Y C 2.0D0**5 - 1.0D0 Z= (X)**(Y) ZZ= (XSQ)**(Y2) GOTO 110 C C TEST OF D**I C 5110 CONTINUE IF (J.GT.5) GOTO 6110 IF (ABS(X).NE.ONE) Y= 170.0D0*Y/(-C*LOG(ABS(X))) IY= NINT(Y) ARGS2(I)= FLOAT(IY) Z= (X)**IY ZZ= ABS(X)**FLOAT(IY) IF ((X.LT.0).AND.(2*(IY/2).NE.IY)) ZZ= -ZZ C ZZ= (X*XSQ)**IY3 GOTO 110 C 6110 CONTINUE IF (X.NE.ONE) Y= 0.7D0*LNINTB*Y/(DELY*LOG(ABS(X))) IX= NINT(X) IF (IX.EQ.0) IX= 1 IY= NINT(Y) ARGS2(I)= FLOAT(IY) ARGS(I)= FLOAT(IX) C WRITE(6,*) X,Y,IX,IY Z= DBLE(IX**IY) ZZ= ABS(DBLE(IX))**DBLE(IY) IF ((IX.LT.0).AND.(2*(IY/2).NE.IY)) ZZ= -ZZ C WRITE(6,*) Z,ZZ C ZZ= (X*XSQ)**IY3 110 W = ONE IF (Z .NE. ZERO) W = (Z - ZZ) / Z 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 IF (J .GE. 4) Y1 = Y 120 R7 = R7 + W * W XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (J .GT. 1) GO TO 210 WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 GO TO 220 210 IF (J .GE. 4) GO TO 215 WRITE (IOUT,1001) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1012) K1,K2,K3 GO TO 220 215 IF (J.EQ.4) WRITE (IOUT,1002) IF (J.EQ.5) WRITE(IOUT,5002) IF (J.EQ.6) WRITE(IOUT,6002) W = C + DELY IF (J.GT.4) GOTO 5215 WRITE (IOUT,1014) N,A,B,C,W GOTO 6215 C 5215 CONTINUE IC= INT(C) IW= INT(W) IF (J.GT.5) GOTO 5216 WRITE(IOUT,5014) N,A,B,IC,IW GOTO 6215 C 5216 IIA= INT(A) IIB= INT(B) WRITE(IOUT,6014) N,IIA,IIB,IC,IW C 6215 CONTINUE IF (J.EQ.4) WRITE (IOUT,1013) K1,K2,K3 IF (J.EQ.5) WRITE(IOUT,5013) K1,K2,K3 IF (J.EQ.6) WRITE(IOUT,6013) K1,K2,K3 220 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA IF (J .LT. 4) WRITE (IOUT,1021) R6,IBETA,W,X1 IF (J .GE. 4) WRITE (IOUT,1024) R6,IBETA,W,X1,Y1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS GOTO (300,302,303,304,305,300), J 302 CONTINUE A = ONE B = DEXP(ALXMAX/3.0D0) GOTO 300 C 303 CONTINUE B = 10.0D0 A = 0.01D0 GOTO 300 C 304 CONTINUE A= -100.0 B= 100.0D0 C= -100.0D0 DELY= -C - C GOTO 300 C 305 CONTINUE A= -100.0 B= 100.0 C= ZERO DELY= 30.0D0 300 CONTINUE C C CPU TIME CALCULATION AND SUMMARY PRINT IF (J.LE.4) THEN NN= 4 ELSE NN= 19 + J ENDIF CALL TIMERI2(J,N,NN,ARGS,ARGS2,MRE,RMS,IOUT) C 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) B = 10.0D0 C DO 320 I = 1, 5 X = RAN(I1) * B + ONE Y = RAN(I1) * B + ONE C Z = X ** Y Z= (X)**(Y) C ZZ = (ONE/X) ** (-Y) ZZ= (ONE/X)**(-Y) W = (Z - ZZ) / Z WRITE (IOUT,1060) X, Y, W 320 CONTINUE C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = BETA Y = DBLE(FLOAT(MINEXP)) Y=-64.0 WRITE (IOUT,1051) X, Y C Z = X ** Y Z= (X)**(Y) WRITE (IOUT,1055) Z Y = DBLE(FLOAT(MAXEXP-1)) WRITE (IOUT,1051) X, Y C Z = X ** Y Z= (X)**(Y) WRITE (IOUT,1055) Z X = ZERO Y = TWO WRITE (IOUT,1051) X, Y C Z = X ** Y Z= (X)**(Y) WRITE (IOUT,1055) Z X = -Y Y = ZERO WRITE (IOUT,1052) X, Y C Z = X ** Y Z= (X)**(Y) WRITE (IOUT,1055) Z Y = TWO WRITE (IOUT,1052) X, Y C Z = X ** Y Z= (X)**(Y) WRITE (IOUT,1055) Z X = ZERO Y = ZERO WRITE (IOUT,1052) X, Y C Z = X ** Y Z= (X)**(Y) WRITE (IOUT,1055) Z C C ERROR EXIT IN R**I C X= ZERO I= -1 WRITE(IOUT,1053) X,I Z= X**I WRITE(IOUT,1055) Z C C ERROR EXITS IN J**I C J= 0 I= 0 WRITE(IOUT,1054) J,I K= J**I WRITE(IOUT,1056) K C J= 2 I= -1 WRITE(IOUT,1054) J,I K= J**I WRITE(IOUT,1056) K C C J= 2 I= 33 WRITE(IOUT,1054) J,I K= J**I WRITE(IOUT,1056) K C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,20H TEST OF X**1.0 VS X //) 1001 FORMAT(26H TEST OF XSQ**1.5 VS XSQ*X //) 1002 FORMAT(27H TEST OF X**Y VS XSQ**(Y/2) //) 5002 FORMAT(' TEST OF X**N VS X**FLOAT(N)' //) 6002 FORMAT(' TEST OF FLOAT(I**N) VS FLOAT(I)**FLOAT(N)' //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(18H X**1.0 WAS LARGER,I6,7H TIMES, / * 11X,7H AGREED,I6,11H TIMES, AND / * 7X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(18H X**1.5 WAS LARGER,I6,7H TIMES, / * 11X,7H AGREED,I6,11H TIMES, AND / * 7X,11HWAS SMALLER,I6,7H TIMES.//) 1013 FORMAT(18H X**Y WAS LARGER,I6,7H TIMES, / * 11X,7H AGREED,I6,11H TIMES, AND / * 7X,11HWAS SMALLER,I6,7H TIMES.//) 5013 FORMAT(18H X**N WAS LARGER,I6,7H TIMES, / * 11X,7H AGREED,I6,11H TIMES, AND / * 7X,11HWAS SMALLER,I6,7H TIMES.//) 6013 FORMAT(18H I**N WAS LARGER,I6,7H TIMES, / * 11X,7H AGREED,I6,11H TIMES, AND / * 7X,11HWAS SMALLER,I6,7H TIMES.//) 1014 FORMAT(I7,45H RANDOM ARGUMENTS WERE TESTED FROM THE REGION / * 6X,6HX IN (,D15.4,1H,,D15.4,9H), Y IN (,D15.4,1H,,D15.4,1H)//) 5014 FORMAT(I7,45H RANDOM ARGUMENTS WERE TESTED FROM THE REGION / * 6X,6HX IN (,D15.4,1H,,D15.4,9H), J IN (,I12,1H,,I12,1H)//) 6014 FORMAT(I7,45H RANDOM ARGUMENTS WERE TESTED FROM THE REGION / * 6X,6HI IN (,I12,1H,,I12,9H), J IN (,I12,1H,,I12,1H)//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1024 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6,4H Y =,D17.6) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(54H THE IDENTITY X ** Y = (1/X) ** (-Y) WILL BE TESTED. * //8X,1HX,14X,1HY,9X,24H(X**Y-(1/X)**(-Y) / X**Y /) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(2H (,D14.7,7H ) ** (,D14.7,20H ) WILL BE COMPUTED.,/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(2H (,D14.7,7H ) ** (,D14.7,20H ) WILL BE COMPUTED.,/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1053 FORMAT(2H (,D14.7,7H ) ** (I5,') WILL BE COMPUTED.',/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1054 FORMAT(2H (,I5,7H ) ** (,I5,20H ) WILL BE COMPUTED.,/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(22H THE VALUE RETURNED IS,D15.4///) 1056 FORMAT(22H THE VALUE RETURNED IS,I12///) 1060 FORMAT(2D15.7,6X,D15.7/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF POWER TEST PROGRAM (D.P.) ---------- END C****************************************************************************** SUBROUTINE DTEST8 C PROGRAM TO TEST DATAN, DATAN2 C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) 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 IRND - 0 IF FLOATING-POINT ADDITION CHOPS, C 1 IF FLOATING-POINT ADDITION ROUNDS C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C POWER OF THE RADIX C XMAX - THE LARGEST FINITE FLOATING-POINT NO. C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DATAN, DATAN2, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,II,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,BETAP,DEL,EM,EPS,EPSNEG, *EXPON,HALF,OB32,ONE,R6,R7,SUM,TWO,W,X,XL,XMAX,XMIN,XN,XSQ,X1, *Y,Z,ZERO,ZZ,DATAN,DATAN2 * ,MRE,RMS C DOUBLE PRECISION RAN C EXTERNAL DATAN,DATAN2 C C CALL PRINTTITLE('DATAN (AND DATAN2) ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0D0 HALF = 0.5D0 TWO = 2.0D0 ZERO = 0.0D0 A = -0.0625D0 B = -A OB32 = B * HALF N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 301 J = 1, 4 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL IF (J .EQ. 2) X = ((1.0D0+X*A)-ONE)*16.0D0 ARGS(I)= X Z = DATAN(X) IF (J .NE. 1) GO TO 100 XSQ = X * X EM = 17.0D0 SUM = XSQ / EM C DO 80 II = 1, 7 EM = EM - TWO SUM = (ONE/EM - SUM) * XSQ 80 CONTINUE C SUM = -X * SUM ZZ = X + SUM SUM = (X - ZZ) + SUM IF (IRND .EQ. 0) ZZ = ZZ + (SUM + SUM) GO TO 110 100 IF (J .NE. 2) GO TO 105 Y = X - .0625D0 Y = Y / (ONE + X*A) ZZ = (DATAN(Y) - 8.1190004042651526021D-5) + OB32 ZZ = ZZ + OB32 GO TO 110 105 Z = Z + Z Y = X / ((HALF + X * HALF)*((HALF - X) + HALF)) ZZ = DATAN(Y) 110 W = ONE IF (Z .NE. ZERO) W = (Z - ZZ) / Z 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 XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (J .EQ. 1) WRITE (IOUT,1000) IF (J .EQ. 2) WRITE (IOUT,1001) IF (J .GT. 2) WRITE (IOUT,1002) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = B IF (J .EQ. 1) B = TWO - DSQRT(3.0D0) IF (J .EQ. 2) B = DSQRT(TWO) - ONE IF (J .EQ. 3) B = ONE 300 CONTINUE C C CPU TIME CALCULATION AND SUMMARY PRINT. C CALL TIMERI(J,N,8,ARGS,MRE,RMS,IOUT) 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) A = 5.0D0 C DO 320 I = 1, 5 X = RAN(I1) * A Z = DATAN(X) + DATAN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = RAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - DATAN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) A = -TWO B = 4.0D0 C DO 340 I = 1, 5 X = RAN(I1) * B + A Y = RAN(I1) W = -Y Z = DATAN(X/Y) - DATAN2(X,Y) ZZ = DATAN(X/W) - DATAN2(X,W) WRITE (IOUT,1059) X, Y, Z, ZZ 340 CONTINUE C WRITE (IOUT,1035) EXPON = DBLE(FLOAT(MINEXP)) * 0.75D0 X = BETA ** EXPON Y = DATAN(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) WRITE (IOUT,1051) XMAX Z = DATAN(XMAX) WRITE (IOUT,1061) XMAX, Z X = ONE Y = ZERO WRITE (IOUT,1053) X, Y Z = DATAN2(X,Y) WRITE (IOUT,1062) X, Y, Z WRITE (IOUT,1053) XMIN, XMAX Z = DATAN2(XMIN,XMAX) WRITE (IOUT,1062) XMIN, XMAX, Z WRITE (IOUT,1053) XMAX, XMIN Z = DATAN2(XMAX,XMIN) WRITE (IOUT,1062) XMAX, XMIN, Z X = ZERO WRITE (IOUT,1054) X, Y Z = DATAN2(X,Y) WRITE (IOUT,1062) X, Y, Z WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,44H TEST OF DATAN(X) VS TRUNCATED TAYLOR SERIES //) 1001 FORMAT(21H TEST OF DATAN(X) VS , * 42HDATAN(1/16) + DATAN((X-1/16)/(1+X/16)) //) 1002 FORMAT(42H TEST OF 2*DATAN(X) VS DATAN(2X/(1-X*X)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(20H DATAN(X) WAS LARGER,I6,7H TIMES,/ * 13X,7H AGREED,I6,11H TIMES, AND / * 9X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(55H THE IDENTITY DATAN(-X) = -DATAN(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY DATAN(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(53H THE IDENTITY DATAN(X/Y) = DATAN2(X,Y) WILL BE TESTED / * 57H THE FIRST COLUMN OF RESULTS SHOULD BE 0, THE SECOND +-PI// * 8X,1HX,17X,1HY,7X,17HF1(X/Y)-F2(X,Y) * ,21HF1(X/(-Y))-F2(X,(-Y)) /) 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT. /) 1050 FORMAT(22H TEST OF ERROR RETURNS //) 1051 FORMAT(39H DATAN WILL BE CALLED WITH THE ARGUMENT,D15.7/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE/) 1053 FORMAT(41H DATAN2 WILL BE CALLED WITH THE ARGUMENTS/2D15.7/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE/) 1054 FORMAT(41H DATAN2 WILL BE CALLED WITH THE ARGUMENTS/2D15.7/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1059 FORMAT(4(D15.7,4X)/) 1060 FORMAT(2D15.7/) 1061 FORMAT(6X,7H DATAN(,D13.6,3H) =,D13.6/) 1062 FORMAT(6X,8H DATAN2(,2D13.6,3H) =,D13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DATAN/DATAN2 TEST PROGRAM ---------- END C C******************************************************************* C SUBROUTINE DTEST10 C PROGRAM TO TEST DTANH C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) C C THIS REFERS TO THE FUNCTION DTANH WHICH IS NOT INCLUDED C IN THE 1966 STANDARD. C IF THIS IS NOT AVAILABLE THE TEST SHOULD BE OMITTED. C C IT IS ASSUMED THAT NO TYPE DECLARATION IS NEEDED. 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 FIVE 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 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT C POWER OF THE RADIX C XMAX - THE LARGEST FINITE FLOATING-POINT NO. C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,BETAP,C,D,DEL,EPS,EPSNEG, *EXPON,HALF,ONE,R6,R7,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ * ,DTANH,MRE,RMS C DOUBLE PRECISION RAN C C EXTERNAL DTANH C CALL PRINTTITLE('DTANH ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) AIT = DBLE(FLOAT(IT)) ZERO = 0.0D0 ONE = 1.0D0 HALF = 0.5D0 A = 0.125D0 B = DLOG(3.0D0) * HALF C = 1.2435300177159620805D-1 D = DLOG(2.0D0) + (AIT+ONE) * DLOG(BETA) * HALF N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 300 J = 1, 2 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL ARGS(I)= X Z = DTANH(X) Y = X - 0.125D0 ZZ = DTANH(Y) ZZ = (ZZ + C) / (ONE + C*ZZ) W = ONE IF (Z .NE. ZERO) W = (Z - ZZ) / Z 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 XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = B + A B = D CALL TIMERI(J,N,10,ARGS,MRE,RMS,IOUT) 300 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1) Z = DTANH(X) + DTANH(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = RAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - DTANH(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) X = D B = 4.0D0 C DO 340 I = 1, 5 Z = (DTANH(X) - HALF) - HALF WRITE (IOUT,1060) X, Z X = X + B 340 CONTINUE C WRITE (IOUT,1035) EXPON = DBLE(FLOAT(MINEXP)) * 0.75D0 X = BETA ** EXPON Z = DTANH(X) WRITE (IOUT,1061) X, Z WRITE (IOUT,1040) XMAX Z = DTANH(XMAX) WRITE (IOUT,1061) XMAX, Z WRITE (IOUT,1040) XMIN Z = DTANH(XMIN) WRITE (IOUT,1061) XMIN, Z X = ZERO WRITE (IOUT,1040) X Z = DTANH(X) WRITE (IOUT,1061) X, Z WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,21H TEST OF DTANH(X) VS , * 54H(DTANH(X-1/8)+DTANH(1/8))/(1+DTANH(X-1/8)DTANH(1/8)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(20H DTANH(X) WAS LARGER,I6,7H TIMES, / * 13X,7H AGREED,I6,11H TIMES, AND / * 9X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(54H THE IDENTITY DTANH(-X) = -DTANH(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY DTANH(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(53H THE IDENTITY DTANH(X) = 1 , X LARGE, WILL BE TESTED.// * 8X,1HX,9X,8H1 - F(X)/) 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT. /) 1040 FORMAT(52H THE FUNCTION DTANH WILL BE CALLED WITH THE ARGUMENT, * D15.7) 1060 FORMAT(2D15.7/) 1061 FORMAT(6X,7H DTANH(,D13.6,3H) =,D13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DTANH TEST PROGRAM ---------- END C************************************************************************ C********************************************************************* C SUBROUTINE DTEST6 C PROGRAM TO TEST DTAN/DCOTAN C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) C C THIS REFERS TO FUNCTIONS DTAN AND DCOTAN. NEITHER WERE C INCLUDED IN THE 1966 STANDARD, AND DCOTAN IS NOT IN C FORTRAN 77. C IF THESE ARE NOT AVAILABLE THE TEST SHOULD BE OMITTED. C C IT IS ASSUMED THAT NO TYPE DECLARATION IS NEEDED. C C IF DTAN IS AVAILABLE BUT NOT DCOTAN, THEN THE TEST CAN BE RUN C PROVIDED THAT : C DCOTAN IS DECLARED DOUBLE PRECISION C A DUMMY FUNCTION DCOTAN IS ADDED C 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 THREE 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 MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,K3,MACHEP,MAXEXP, *MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,BETAP,C1,C2,DEL,EPS,EPSNEG, *HALF,PI,R6,R7,W,X,XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ DOUBLE PRECISION DTAN,DCOTAN,MRE,RMS DOUBLE PRECISION RAN C C EXTERNAL DTAN,DCOTAN C C THE FOLLOWING STATEMENT FUNCTION SHOULD BE ACTIVATED IF C THIS TEST IS TO BE RUN BUT NO FUNCTION DCOTAN (UNDER ANY C NAME) IS AVAILABLE. C C DCOTAN(Z) = 0.0D0 C CALL PRINTTITLE ('DTAN (J=1,2,3), DCOTAN (J=4) ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) ZERO = 0.0D0 HALF = 0.5D0 AIT = DBLE(FLOAT(IT)) PI = 3.14159265D0 A = ZERO B = PI * 0.25D0 N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- C C THE FOLLOWING LINE SHOULD USE J = 1,4 IF TESTING DTAN AND DCOTAN C J = 1,3 IF TESTING DTAN ONLY C DO 301 J = 1, 4 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL Y = X * HALF X = Y + Y ARGS(I)= X IF (J .EQ. 4) GO TO 80 Z = DTAN(X) ZZ = DTAN(Y) W = 1.0D0 IF (Z .EQ. ZERO) GO TO 110 W = ((HALF-ZZ)+HALF)*((HALF+ZZ)+HALF) W = (Z - (ZZ+ZZ)/W) / Z GO TO 110 80 Z = DCOTAN(X) ZZ = DCOTAN(Y) W = 1.0D0 IF (Z .EQ. ZERO) GO TO 110 W = ((HALF-ZZ)+HALF)*((HALF+ZZ)+HALF) W = (Z+W/(ZZ+ZZ))/Z 110 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 XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (J .NE. 4) WRITE (IOUT,1000) IF (J .EQ. 4) WRITE (IOUT,1005) WRITE (IOUT,1010) N,A,B IF (J .NE. 4) WRITE (IOUT,1011) K1,K2,K3 IF (J .EQ. 4) WRITE (IOUT,1012) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .NE. 1) GO TO 250 A = PI * 0.875D0 B = PI * 1.125D0 GO TO 300 250 A = PI * 6.0D0 B = A + PI * 0.25D0 300 CONTINUE IF (J.NE.4) THEN CALL TIMERI(J,N,6,ARGS,MRE,RMS,IOUT) ELSE CALL TIMERI(J,N,16,ARGS,MRE,RMS,IOUT) ENDIF 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1) * A Z = DTAN(X) + DTAN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = RAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - DTAN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1035) X = BETA ** (DBLE(FLOAT(MINEXP))*0.75D0) Y = DTAN(X) WRITE (IOUT,1061) X, Y C1 = -225.0D0 C2 = -.950846454195142026D0 X = 11.0D0 Y = DTAN(X) W = ((C1-Y)+C2)/(C1+C2) Z = DLOG(DABS(W))/ALBETA WRITE (IOUT,1040) W, IBETA, Z WRITE (IOUT,1061) X, Y W = DMAX1(AIT+Z,ZERO) WRITE (IOUT,1022) IBETA, W C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = BETA ** (IT/2) WRITE (IOUT,1051) X Y = DTAN(X) WRITE (IOUT,1055) Y X = BETAP WRITE (IOUT,1052) X Y = DTAN(X) WRITE (IOUT,1055) Y X= 0.0D0 WRITE(IOUT,1072) X Y= DCOTAN(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,48H TEST OF DTAN(X) VS 2*DTAN(X/2)/(1-DTAN(X/2)**2) //) 1005 FORMAT(40H TEST OF DCOTAN(X) VS (DCOTAN(X/2)**2-1), * 16H/(2*DCOTAN(X/2)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(19H DTAN(X) WAS LARGER,I6,7H TIMES, / * 12X,7H AGREED,I6,11H TIMES, AND / * 8X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(21H DCOTAN(X) WAS LARGER,I6,7H TIMES, / * 14X,7H AGREED,I6,11H TIMES, AND / * 10X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(51H THE IDENTITY DTAN(-X) = -DTAN(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(52H THE IDENTITY DTAN(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT. /) 1040 FORMAT(34H THE RELATIVE ERROR IN DTAN(11) IS ,D15.7,3H = , * I4,3H **,F7.2,6H WHERE /) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(38H DTAN WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(38H0DTAN WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H DTAN RETURNED THE VALUE,D15.4///) 1060 FORMAT(2D15.7/) 1061 FORMAT(6X,6H DTAN(,D13.6,3H) =,D13.6/) 1071 FORMAT(' DCOTAN WILL BE CALLED WITH THE ARGUMENT',D15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1072 FORMAT('0DCOTAN WILL BE CALLED WITH THE ARGUMENT',D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(' DCOTAN RETURNED THE VALUE',D15.4///) 1081 FORMAT(6X,' DCOTAN(',D13.6,3H) =,D13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DTAN/DCOTAN TEST PROGRAM ---------- END C C******************************************************************** C SUBROUTINE DTEST7 C PROGRAM TO TEST DASIN/DACOS C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) C C THIS REFERS TO FUNCTIONS DASIN AND DACOS. C THESE WERE NOT INCLUDED IN THE 1966 STANDARD. C IF THESE ARE NOT AVAILABLE THE TEST SHOULD BE OMITTED. C IT IS POSSIBLE THAT THE LIBRARY FUNCTIONS ARE CALLED C DARSIN AND DARCOS, IN WHICH CASE EDIT THE TEST FIRST. C C IT IS ASSUMED THAT NO TYPE DECLARATION IS NEEDED. 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 FOUR 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 IRND - 0 IF FLOATING-POINT ADDITION CHOPS, C 1 IF FLOATING-POINT ADDITION ROUNDS C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DLOG10, DMAX1, IDINT, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K,K1,K2,K3,L,M,MACHEP, *MAXEXP,MINEXP,N,NEGEP,NGRD DOUBLE PRECISION A,AIT,ALBETA,B,BETA,BETAP,C1,C2,DEL,EPS,EPSNEG, *HALF,R6,R7,S,SUM,W,X,XL,XM,XMAX,XMIN,XN,X1,Y,YSQ,Z,ZERO,ZZ * ,DASIN,DACOS,MRE,RMS C DOUBLE PRECISION RAN C C EXTERNAL DASIN,DACOS C CALL PRINTTITLE('DASIN (J=1,3), DACOS (J=2,4,5)') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) ZERO = 0.0D0 HALF = 0.5D0 AIT = DBLE(FLOAT(IT)) K = IDINT(DLOG10(BETA**IT)) + 1 IF (IBETA .NE. 10) GO TO 20 C1 = 1.57D0 C2 = 7.9632679489661923132D-4 GO TO 30 20 C1 = 201.0D0/128.0D0 C2 = 4.8382679489661923132D-4 30 A = -0.125D0 B = -A N = 2000 C READ (5,*) N XN = DBLE(FLOAT(N)) I1 = 0 L = -1 C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 301 J = 1, 5 K1 = 0 K3 = 0 L = -L X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL*RAN(I1) + XL IF (J .LE. 2) GO TO 40 YSQ = HALF - HALF*DABS(X) X = (HALF - (YSQ+YSQ)) + HALF IF (J .EQ. 5) X = -X Y = DSQRT(YSQ) Y = Y + Y GO TO 50 40 Y = X YSQ = Y*Y 50 SUM = ZERO XM = DBLE(FLOAT(K+K+1)) ARGS(I)= X IF (L .GT. 0) Z =DASIN(X) IF (L .LT. 0) Z =DACOS(X) C DO 60 M = 1, K SUM = YSQ*(SUM + 1.0D0/XM) XM = XM - 2.0D0 SUM = SUM*(XM/(XM+1.0D0)) 60 CONTINUE C SUM = SUM*Y IF ((J .NE. 1) .AND. (J .NE. 4)) GO TO 70 ZZ = Y + SUM SUM = (Y - ZZ) + SUM IF (IRND .NE. 1) ZZ = ZZ + (SUM+SUM) GO TO 110 70 S = C1 + C2 SUM = ((C1 - S) + C2) - SUM ZZ = S + SUM SUM = ((S - ZZ) + SUM) - Y S = ZZ ZZ = S + SUM SUM = (S - ZZ) + SUM IF (IRND .NE. 1) ZZ = ZZ + (SUM+SUM) 110 W = 1.0D0 C WRITE(6,* ) Z,ZZ IF (Z .NE. ZERO) W = (Z-ZZ)/Z 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 XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) IF (L .LT. 0) GO TO 210 WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 GO TO 220 210 WRITE (IOUT,1005) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1012) K1,K2,K3 220 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .NE. 2) GO TO 250 A = 0.75D0 B = 1.0D0 250 IF (J .NE. 4) GO TO 300 B = -A A = -1.0D0 C1 = C1 + C1 C2 = C2 + C2 L = -L 300 CONTINUE C IF (L.GT.0) CALL TIMERI(J,N,7,ARGS,MRE,RMS,IOUT) C IF (L.LT.0) CALL TIMERI(J,N,17,ARGS,MRE,RMS,IOUT) C 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1)*A Z = DASIN(X) + DASIN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = RAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - DASIN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1035) X = BETA ** (DBLE(FLOAT(MINEXP))*0.75D0) Y = DASIN(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = 1.2D0 WRITE (IOUT,1052) X Y = DASIN(X) WRITE (IOUT,1055) Y WRITE(IOUT,1072) X Y= DACOS(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,34H TEST OF DASIN(X) VS TAYLOR SERIES //) 1005 FORMAT(34H TEST OF DACOS(X) VS TAYLOR SERIES //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(20H DASIN(X) WAS LARGER,I6,7H TIMES, / * 14X,7H AGREED,I6,11H TIMES, AND / * 10X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(20H DACOS(X) WAS LARGER,I6,7H TIMES, / * 14X,7H AGREED,I6,11H TIMES, AND / * 10X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(53H THE IDENTITY DASIN(-X) = -DASIN(X) WILL BE TESTED.// * ,8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY DASIN(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT. /) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1052 FORMAT(39H DASIN WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H DASIN RETURNED THE VALUE,D15.4///) 1060 FORMAT(2D15.7/) 1072 FORMAT(39H DACOS WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(25H DACOS RETURNED THE VALUE,D15.4///) 1061 FORMAT(6X,7H DASIN(,D13.6,3H) =,D13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DASIN/DACOS TEST PROGRAM ---------- END C************************************************************************ C SUBROUTINE DTEST9 C PROGRAM TO TEST DSINH/DCOSH C C C COMMON /PRMVLS/ ARGS,ARGS2 DOUBLE PRECISION ARGS(10000),ARGS2(10000) C C THIS REFERS TO FUNCTIONS DSINH AND DCOSH. C THESE WERE NOT INCLUDED IN THE 1966 STANDARD. C IF THESE ARE NOT AVAILABLE THE TEST SHOULD BE OMITTED. C C IT IS ASSUMED THAT NO TYPE DECLARATION IS NEEDED. 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 IRND - 0 IF FLOATING-POINT ADDITION CHOPS, C 1 IF FLOATING-POINT ADDITION ROUNDS C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE C INTEGER SUCH THAT FLOAT(IBETA)**MINEXP C IS A POSITIVE FLOATING-POINT NUMBER C EPS - THE SMALLEST POSITIVE FLOATING-POINT C NUMBER SUCH THAT 1.0+EPS .NE. 1.0 C XMAX - THE LARGEST FINITE FLOATING-POINT NO. C C RAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C DABS, DLOG, DMAX1, IDINT, DSQRT C C NOTE THAT THE CONSTRUCT DBLE(FLOAT( )) IS USED. THIS SHOULD C BE SATISFACTORY AS ALL INTEGERS CONCERNED ARE SMALL. IF C DIFFICULTIES ARISE THE USE OF DFLOAT MIGHT BE NECESSARY. C C AUTHOR - W. J. CODY C ARGONNE NATIONAL LABORATORY C C C INTEGER I,IBETA,IEXP,II,IOUT,IRND,IT,I1,I2,J,K1,K2,K3,MACHEP, *MAXEXP,MINEXP,N,NEGEP,NGRD,NIT DOUBLE PRECISION A,AIND,AIT,ALBETA,ALXMAX,B,BETA,BETAP,C,C0,DEL, *DEN,EPS,EPSNEG,FIVE,ONE,R6,R7,THREE,W,X,XL,XMAX,XMIN,XN,X1, *XSQ,Y,Z,ZERO,ZZ, DCOSH, DSINH,MRE,RMS C DOUBLE PRECISION RAN C C EXTERNAL DCOSH, DSINH C CALL PRINTTITLE('DSINH (J=1,3), DCOSH (J=2,4) ') IOUT = 10 CALL MACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = DLOG(BETA) ALXMAX = DLOG(XMAX) AIT = DBLE(FLOAT(IT)) ZERO = 0.0D0 ONE = 1.0D0 THREE = 3.0D0 FIVE = 5.0D0 C0 = FIVE/16.0D0 + 1.152713683194269979D-2 A = ZERO B = 0.5D0 C = (AIT + ONE) * 0.35D0 IF (IBETA .EQ. 10) C = C * THREE N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 I2 = 2 NIT = 2 - (IDINT(DLOG(EPS)*THREE))/20 AIND = DBLE(FLOAT(NIT+NIT+1)) C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 301 J = 1, 4 IF (J .NE. 2) GO TO 30 AIND = AIND - ONE I2 = 1 30 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * RAN(I1) + XL IF (J .GT. 2) GO TO 80 XSQ = X * X ZZ = ONE DEN = AIND C DO 40 II = I2, NIT W = ZZ * XSQ/(DEN*(DEN-ONE)) ZZ = W + ONE DEN = DEN - 2.0D0 40 CONTINUE C ARGS(I)= X IF (J .EQ. 2) GO TO 50 W = X*XSQ*ZZ/6.0D0 ZZ = X + W Z = DSINH(X) IF (IRND .NE. 0) GO TO 110 W = (X - ZZ) + W ZZ = ZZ + (W + W) GO TO 110 50 Z = DCOSH(X) IF (IRND .NE. 0) GO TO 110 W = (ONE - ZZ) + W ZZ = ZZ + (W + W) GO TO 110 80 Y = X X = Y - ONE ARGS(I)= X W = X - ONE IF (J .EQ. 4) GO TO 100 Z = DSINH(X) ZZ = ( DSINH(Y) + DSINH(W)) * C0 C WRITE(6,12121) Z,ZZ C12121 FORMAT(1X,2(Z16,1X)) GO TO 110 100 Z = DCOSH(X) ZZ = (DCOSH(Y) + DCOSH(W)) * C0 110 W = ONE IF (Z .NE. ZERO) W = (Z - ZZ)/Z 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 XL = XL + DEL 200 CONTINUE C K2 = N - K3 - K1 R7 = DSQRT(R7/XN) I = (J/2) * 2 IF (J .EQ. 1) WRITE (IOUT,1000) IF (J .EQ. 2) WRITE (IOUT,1005) IF (J .EQ. 3) WRITE (IOUT,1001) IF (J .EQ. 4) WRITE (IOUT,1006) WRITE (IOUT,1010) N,A,B IF (I .NE. J) WRITE (IOUT,1011) K1,K2,K3 IF (I .EQ. J) WRITE (IOUT,1012) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0D0 IF (R6 .NE. ZERO) W = DLOG(DABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0D0 IF (R7 .NE. ZERO) W = DLOG(DABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = DMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .NE. 2) GO TO 300 B = ALXMAX A = THREE 300 CONTINUE IF ((J.EQ.1).OR.(J.EQ.3)) * CALL TIMERI(J,N,9,ARGS,MRE,RMS,IOUT) IF ((J.EQ.2).OR.(J.EQ.4)) * CALL TIMERI(J,N,19,ARGS,MRE,RMS,IOUT) 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) C DO 320 I = 1, 5 X = RAN(I1) * A Z = DSINH(X) + DSINH(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = RAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - DSINH(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) C DO 340 I = 1, 5 X = RAN(I1) * A Z = DCOSH(X) - DCOSH(-X) WRITE (IOUT,1060) X, Z 340 CONTINUE C WRITE (IOUT,1035) X = BETA ** (DBLE(FLOAT(MINEXP))*0.75D0) Y = DSINH(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = ALXMAX + 0.125D0 WRITE (IOUT,1051) X Y = DSINH(X) WRITE (IOUT,1055) Y X = BETAP WRITE (IOUT,1052) X Y = DSINH(X) WRITE (IOUT,1055) Y WRITE(IOUT,1072) X Y= DCOSH(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,47H TEST OF DSINH(X) VS T.S. EXPANSION OF DSINH(X) //) 1001 FORMAT(46H TEST OF DSINH(X) VS C*(DSINH(X+1)+DSINH(X-1)) //) 1005 FORMAT(47H TEST OF DCOSH(X) VS T.S. EXPANSION OF DCOSH(X) //) 1006 FORMAT(46H TEST OF DCOSH(X) VS C*(DCOSH(X+1)+DCOSH(X-1)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,D15.4,1H,,D15.4,1H)//) 1011 FORMAT(20H DSINH(X) WAS LARGER,I6,7H TIMES, / * 13X,7H AGREED,I6,11H TIMES, AND / * 9X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(20H DCOSH(X) WAS LARGER,I6,7H TIMES, / * 13X,7H AGREED,I6,11H TIMES, AND / * 9X,11HWAS SMALLER,I6,7H TIMES.//) 1020 FORMAT(10H THERE ARE,I4,5H BASE,I4, * 46H SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER //) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(53H THE IDENTITY DSINH(-X) = -DSINH(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY DSINH(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(52H THE IDENTITY DCOSH(-X) = DCOSH(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) - F(-X)/) 1035 FORMAT(43H TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT. /) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(39H DSINH WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(39H0DSINH WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H DSINH RETURNED THE VALUE,D15.4///) 1060 FORMAT(2D15.7/) 1072 FORMAT(39H DCOSH WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(25H DCOSH RETURNED THE VALUE,D15.4///) 1061 FORMAT(6X,7H DSINH(,D13.6,3H) =,D13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF DSINH/DCOSH TEST PROGRAM ---------- END C****************************************************************************** SUBROUTINE DTEST11 C C TEST DGAMMA,DLGAMA,DERF,DERFC C COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN CALL MACHAR(IBETA, IT, IRND, NGRD, MACHEP, NEGEP, IEXP,MINEXP, * MAXEXP, EPS, EPSNEG, XMIN, XMAX) C CALL SETMESSAGES REWIND 9 CALL DTEST11A CALL DTEST11B CALL DTEST11C CALL DTEST11D RETURN END C SUBROUTINE DTEST11A C DGAMMA TEST PROGRAM TEXT COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C ACTUAL TEST APPEAR AS SEPARATE BLOCKS OF CODE,THE DATA C CONTROLS THE SWITCHING TO WHICH EVER BLOCK IS REQUIRED C C .. LOCAL SCALARS .. DOUBLE PRECISION DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. DOUBLE PRECISION DLOG10, DGAMMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, RANTST C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL PRINTTITLE('DGAMMA ') 50000 FORMAT(1X,80('*'),/,' TEST OF DGAMMA USING RESULTS FROM FILE.') C C DO 55 JJ= 1,2 READ(NIN,50001) TEXT WRITE(NOUT,50002) TEXT 50001 FORMAT(A72) 50002 FORMAT(1X,A72) C C TEST TYPE 1 40 CALL RANTSTA(NIN, NOUT,JJ) 55 CONTINUE C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= 56.9D0 WRITE(NOUT,1051) X Y= DGAMMA(X) WRITE(NOUT,1055) Y X= 57.1D0 WRITE(NOUT,1052) X Y= DGAMMA(X) WRITE(NOUT,1055) Y C X= -56.9D0 WRITE(NOUT,1051) X Y= DGAMMA(X) WRITE(NOUT,1055) Y X= -57.1D0 WRITE(NOUT,1052) X Y= DGAMMA(X) WRITE(NOUT,1055) Y C X= -1.0D-74 WRITE(NOUT,1051) X Y= DGAMMA(X) WRITE(NOUT,1055) Y X= 1.0D-76 WRITE(NOUT,1052) X Y= DGAMMA(X) WRITE(NOUT,1055) Y C X= -9.99999D-1 WRITE(NOUT,1051) X Y= DGAMMA(X) WRITE(NOUT,1055) Y X= -1.000D0 WRITE(NOUT,1052) X Y= DGAMMA(X) WRITE(NOUT,1055) Y C WRITE(NOUT,1100) RETURN 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(40H DGAMMA WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(40H0DGAMMA WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(26H DGAMMA RETURNED THE VALUE,D15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE RANTSTA(NIN, NOUT,JJ) COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C TEST TYPE 1 C INPUT ON CR NO. NIN C OUTPUT ON LP NO.NOUT C .. SCALAR ARGUMENTS .. INTEGER NIN, NOUT C .. FUNCTION ARGUMENTS .. C .. C .. LOCAL SCALARS .. DOUBLE PRECISION EA, EM, ER, ERR, FEA, FEM, FER, PRTER, REPER, X, * XL, XU,Y,ALBETA,AIT,W,MRE,RMSX1 INTEGER I, IET, IFAIL, ISAM, I1, I2, J, KT, MULT, NX COMMON /PRMVLS/ ARGS,ARG2 DOUBLE PRECISION ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. DOUBLE PRECISION DSQRT ,DGAMMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, READT C .. DATA ETP(1), ETP(2), ETP(3), ETP(4), ETP(5), ETP(6) /1HA,1HB,1HS, *1HR,1HE,1HL/ DATA SAM(1), SAM(2), SAM(3), SAM(4), SAM(5), SAM(6) /1HU,1HN,1HI, *1HL,1HO,1HG/ ALBETA= DLOG(16.0D0) IBETA= 16 X1= 0.0D0 AIT= 14.0D0 READ (NIN,99999) NX, XL, XU, ISAM, IET, MULT I2 = ISAM*3 I1 = I2 - 2 WRITE (NOUT,99998) NX, XL, XU, (SAM(I),I=I1,I2) REPER = EPS PRTER = DFLOAT(MULT)*REPER I2 = IET*3 I1 = I2 - 2 IF (JJ.EQ.1) 1 WRITE (NOUT,99997) (ETP(I),I=I1,I2), REPER, MULT IF (JJ.EQ.2) 1 WRITE (NOUT,99987) (ETP(I),I=I1,I2), REPER, MULT EA = 0.0D0 ER = 0.0D0 EM = 0.0D0 KT = 0 DO 40 J=1,NX CALL READT(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SDDAAF/ IFAIL = 1 Y = DGAMMA(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL ERRCAL(Y, IET, ERR, IFAIL) IF (JJ.EQ.2) ERR= ERR/(X) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SDDAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. DABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL PRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF (DABS(ERR).GT.DABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL PRNT(NOUT, X, Y, 1.0D0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = DSQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0D0 IF (EM .NE. 0.0D0) W= DLOG(DABS(EM))/ALBETA WRITE (NOUT,1021) EM,IBETA,W,X1 MRE = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,MRE W = -999.0D0 IF (ER .NE. 0.0D0) W = DLOG(DABS(ER))/ALBETA WRITE (NOUT,1023) ER,IBETA,W RMS = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C J= 1 CALL TIMERI(J,NX,20,ARGS,MRE,RMS,NOUT) 1021 FORMAT(21H THE MAXIMUM ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(31H THE ROOT MEAN SQUARE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PD11.2, 5X, D11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PD11.2, * 5H AND , D11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PD11.2/46H PRINT IN FULL FIRST POINT, ALSO PRINT IN FULL/ * 14H IF(ERROR.GE. , I4, 10H*REP.ERR.)/1H0, 29X, 1HX, 29X, 1HY, * 26X, 4HYEXP, 4X, 11HERR/REP.ERR) 99987 FORMAT (15H ERROR MEASURE , 3A1,'/X',25H, REPRESENTATION ERROR = , * 1PD11.2/46H PRINT IN FULL FIRST POINT, ALSO PRINT IN FULL/ * 14H IF(ERROR.GE. , I4, 10H*REP.ERR.)/1H0, 29X, 1HX, 29X, 1HY, * 26X, 4HYEXP, 4X, 11HERR/REP.ERR) 99996 FORMAT (12H FAILURE NO., I2, 10H CAUSED BY) 99995 FORMAT (13H0MEAN ERROR =, 1PD11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END SUBROUTINE READT(NIN, X) C THIS ROUTINE MAY BE MACHINE DEPENDENT C IT READS FROM DEVICE NO. NIN A ARGUMENT X (STD PRECISION) C AND THE EXPECTED VALUE YEXP (EXTRA PRECISION) C X IS RETURNED VIA THE ARGUMENT LIST C YEXP VIA COMMON BLOCK /SDDAAF/ C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X INTEGER NIN C .. C .. SCALARS IN COMMON .. DOUBLE PRECISION YEXP C .. COMMON /SDDAAF/ YEXP READ (NIN,99999) X, YEXP RETURN 99999 FORMAT (D20.6, D40.20) END SUBROUTINE ERRCAL(Y, IET, ERR, IFAIL) C TAKES VALUE OF FUNCTION Y (STD PRECISION) AND CALCULATES C THE VALUE OF THE ERROR (TYPE IET,1 ABS,2 REL) C YEXP IS PASSED VIA COMMON BLOCK /SDDAAF/ C .. SCALAR ARGUMENTS .. DOUBLE PRECISION ERR, Y INTEGER IET, IFAIL C .. C .. SCALARS IN COMMON .. DOUBLE PRECISION YEXP C .. C .. LOCAL SCALARS .. DOUBLE PRECISION DERR, DY C .. COMMON /SDDAAF/ YEXP DY = Y C WRITE(6,123) Y,YEXP 123 FORMAT(2(1X,Z16)) DERR = YEXP - DY GO TO (20, 40), IET 20 ERR = DERR IFAIL = 0 RETURN 40 IF (Y.EQ.0.0D0) GO TO 60 ERR = DERR/Y IFAIL = 0 RETURN 60 ERR = 0.0D0 IFAIL = 1 RETURN END SUBROUTINE PRNT(NOUT, X, Y, ERR) C PRINT SUBROUTINE(MACHINE DEPENDENT) C YEXP IS PASSED VIA COMMON /SDDAAF/ AND IS IN EXTRA PRECISION C PRINTING IS ON DEVICE NO. NOUT C .. SCALAR ARGUMENTS .. DOUBLE PRECISION ERR, X, Y INTEGER NOUT C .. C .. SCALARS IN COMMON .. DOUBLE PRECISION YEXP C .. COMMON /SDDAAF/ YEXP WRITE (NOUT,99999) X, Y, YEXP, ERR RETURN 99999 FORMAT (1X, 1PD30.11, D30.11, D30.12, 0PF15.2) END C****************************************************************************** SUBROUTINE DTEST11B C DLGAMA TEST PROGRAM TEXT COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C ACTUAL TEST APPEAR AS SEPARATE BLOCKS OF CODE,THE DATA C CONTROLS THE SWITCHING TO WHICH EVER BLOCK IS REQUIRED C C .. LOCAL SCALARS .. DOUBLE PRECISION DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. DOUBLE PRECISION DLOG10, DLGAMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, RANTST C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL PRINTTITLE('DLGAMA ') 50000 FORMAT(1X,80('*'),/,' TEST OF DLGAMA USING RESULTS FROM FILE.') C C DO 55 J= 1,2 READ(NIN,50001) TEXT WRITE(NOUT,50002) TEXT 50001 FORMAT(A72) 50002 FORMAT(1X,A72) C C TEST TYPE 1 40 CALL RANTSTB(NIN, NOUT,J) 55 CONTINUE C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= 1.0D70 WRITE(NOUT,1051) X Y= DLGAMA(X) WRITE(NOUT,1055) Y X= 1.0D75 WRITE(NOUT,1052) X Y= DLGAMA(X) WRITE(NOUT,1055) Y C X= -56.9D0 WRITE(NOUT,1052) X Y= DLGAMA(X) WRITE(NOUT,1055) Y C C WRITE(NOUT,1100) RETURN 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(40H DLGAMA WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(40H0DLGAMA WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(26H DLGAMA RETURNED THE VALUE,D15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE RANTSTB(NIN, NOUT,JJ) COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C TEST TYPE 1 C INPUT ON CR NO. NIN C OUTPUT ON LP NO.NOUT C .. SCALAR ARGUMENTS .. INTEGER NIN, NOUT C .. FUNCTION ARGUMENTS .. C .. C .. LOCAL SCALARS .. DOUBLE PRECISION EA, EM, ER, ERR, FEA, FEM, FER, PRTER, REPER, X, * XL, XU,Y,ALBETA,AIT,W,MRE,RMSX1 INTEGER I, IET, IFAIL, ISAM, I1, I2, J, KT, MULT, NX COMMON /PRMVLS/ ARGS,ARG2 DOUBLE PRECISION ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. DOUBLE PRECISION DSQRT ,DLGAMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, READT C .. DATA ETP(1), ETP(2), ETP(3), ETP(4), ETP(5), ETP(6) /1HA,1HB,1HS, *1HR,1HE,1HL/ DATA SAM(1), SAM(2), SAM(3), SAM(4), SAM(5), SAM(6) /1HU,1HN,1HI, *1HL,1HO,1HG/ ALBETA= DLOG(16.0D0) IBETA= 16 X1= 0.0D0 AIT= 14.0D0 READ (NIN,99999) NX, XL, XU, ISAM, IET, MULT I2 = ISAM*3 I1 = I2 - 2 WRITE (NOUT,99998) NX, XL, XU, (SAM(I),I=I1,I2) REPER = EPS PRTER = DFLOAT(MULT)*REPER I2 = IET*3 I1 = I2 - 2 WRITE (NOUT,99997) (ETP(I),I=I1,I2), REPER, MULT EA = 0.0D0 ER = 0.0D0 EM = 0.0D0 KT = 0 DO 40 J=1,NX CALL READT(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SDDAAF/ IFAIL = 1 Y = DLGAMA(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL ERRCAL(Y, IET, ERR, IFAIL) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SDDAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. DABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL PRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF (DABS(ERR).GT.DABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL PRNT(NOUT, X, Y, 1.0D0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = DSQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0D0 IF (EM .NE. 0.0D0) W= DLOG(DABS(EM))/ALBETA WRITE (NOUT,1021) EM,IBETA,W,X1 MRE = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,MRE W = -999.0D0 IF (ER .NE. 0.0D0) W = DLOG(DABS(ER))/ALBETA WRITE (NOUT,1023) ER,IBETA,W RMS = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C CALL TIMERI(JJ,NX,21,ARGS,MRE,RMS,NOUT) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PD11.2, 5X, D11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PD11.2, * 5H AND , D11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PD11.2/46H PRINT IN FULL FIRST POINT, ALSO PRINT IN FULL/ * 14H IF(ERROR.GE. , I4, 10H*REP.ERR.)/1H0, 29X, 1HX, 29X, 1HY, * 26X, 4HYEXP, 4X, 11HERR/REP.ERR) 99996 FORMAT (12H FAILURE NO., I2, 10H CAUSED BY) 99995 FORMAT (13H0MEAN ERROR =, 1PD11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END C****************************************************************************** SUBROUTINE DTEST11C C DERF TEST PROGRAM TEXT COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C ACTUAL TEST APPEAR AS SEPARATE BLOCKS OF CODE,THE DATA C CONTROLS THE SWITCHING TO WHICH EVER BLOCK IS REQUIRED C C .. LOCAL SCALARS .. DOUBLE PRECISION DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. DOUBLE PRECISION DLOG10, DERF C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT DOUBLE PRECISION S15AEF CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, RANTSTC C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL PRINTTITLE('DERF ') 50000 FORMAT(1X,80('*'),/,' TEST OF DERF USING RESULTS FROM FILE.') C C READ(NIN,50001) TEXT WRITE(NOUT,50002) TEXT 50001 FORMAT(A72) 50002 FORMAT(1X,A72) C C TEST TYPE 1 40 CALL RANTSTC(NIN, NOUT) C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= -10.0D0 WRITE(NOUT,1051) X I= 0 C Y= S15AEF(X,I) Y= DERF(X) C Y= DERF(X) WRITE(NOUT,1055) Y X= 10.0D0 WRITE(NOUT,1051) X Y= DERF(X) WRITE(NOUT,1055) Y C WRITE(NOUT,1100) RETURN 1050 FORMAT(23H TEST OF EXTREME VALUES//) 1051 FORMAT(38H DERF WILL BE CALLED WITH THE ARGUMENT,D15.4) 1052 FORMAT(38H0DERF WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H DERF RETURNED THE VALUE,D15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE RANTSTC(NIN, NOUT) COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C TEST TYPE 1 C INPUT ON CR NO. NIN C OUTPUT ON LP NO.NOUT C .. SCALAR ARGUMENTS .. INTEGER NIN, NOUT C .. FUNCTION ARGUMENTS .. C .. C .. LOCAL SCALARS .. DOUBLE PRECISION EA, EM, ER, ERR, FEA, FEM, FER, PRTER, REPER, X, * XL, XU,Y,ALBETA,AIT,W,MRE,RMSX1 INTEGER I, IET, IFAIL, ISAM, I1, I2, J, KT, MULT, NX DOUBLE PRECISION S15AEF COMMON /PRMVLS/ ARGS,ARG2 DOUBLE PRECISION ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. DOUBLE PRECISION DSQRT ,DERF C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, READTC C .. DATA ETP(1), ETP(2), ETP(3), ETP(4), ETP(5), ETP(6) /1HA,1HB,1HS, *1HR,1HE,1HL/ DATA SAM(1), SAM(2), SAM(3), SAM(4), SAM(5), SAM(6) /1HU,1HN,1HI, *1HL,1HO,1HG/ ALBETA= DLOG(16.0D0) IBETA= 16 X1= 0.0D0 AIT= 14.0D0 READ (NIN,99999) NX, XL, XU, ISAM, IET, MULT I2 = ISAM*3 I1 = I2 - 2 WRITE (NOUT,99998) NX, XL, XU, (SAM(I),I=I1,I2) REPER = EPS PRTER = DFLOAT(MULT)*REPER I2 = IET*3 I1 = I2 - 2 WRITE (NOUT,99997) (ETP(I),I=I1,I2), REPER, MULT EA = 0.0D0 ER = 0.0D0 EM = 0.0D0 KT = 0 DO 40 J=1,NX CALL READTC(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SDDAAF/ IFAIL = 1 JJ= 0 Y= DERF(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL ERRCAL(Y, IET, ERR, IFAIL) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SDDAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. DABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL PRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF (DABS(ERR).GT.DABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL PRNT(NOUT, X, Y, 1.0D0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = DSQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0D0 IF (EM .NE. 0.0D0) W= DLOG(DABS(EM))/ALBETA WRITE (NOUT,1021) EM,IBETA,W,X1 MRE = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,MRE W = -999.0D0 IF (ER .NE. 0.0D0) W = DLOG(DABS(ER))/ALBETA WRITE (NOUT,1023) ER,IBETA,W RMS = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C J= 1 CALL TIMERI(J,NX,22,ARGS,MRE,RMS,NOUT) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(40H THE ROOT MEAN SQUARE RELATIVE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PD11.2, 5X, D11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PD11.2, * 5H AND , D11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PD11.2/46H PRINT IN FULL FIRST POINT, ALSO PRINT IN FULL/ * 14H IF(ERROR.GE. , I4, 10H*REP.ERR.)/1H0, 29X, 1HX, 29X, 1HY, * 26X, 4HYEXP, 4X, 11HERR/REP.ERR) 99996 FORMAT (12H FAILURE NO., I2, 10H CAUSED BY) 99995 FORMAT (13H0MEAN ERROR =, 1PD11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END SUBROUTINE READTC(NIN, X) C THIS ROUTINE MAY BE MACHINE DEPENDENT C IT READS FROM DEVICE NO. NIN A ARGUMENT X (STD PRECISION) C AND THE EXPECTED VALUE YEXP (EXTRA PRECISION) C X IS RETURNED VIA THE ARGUMENT LIST C YEXP VIA COMMON BLOCK /SDDAAF/ C C .. SCALAR ARGUMENTS .. DOUBLE PRECISION X INTEGER NIN C .. C .. SCALARS IN COMMON .. DOUBLE PRECISION YEXP C .. COMMON /SDDAAF/ YEXP READ (NIN,99999) X, YEXP RETURN C99999 FORMAT (D20.6, D40.20) 99999 FORMAT(2(1X,D30.19)) END C****************************************************************************** SUBROUTINE DTEST11D C DERFC TEST PROGRAM TEXT COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C ACTUAL TEST APPEAR AS SEPARATE BLOCKS OF CODE,THE DATA C CONTROLS THE SWITCHING TO WHICH EVER BLOCK IS REQUIRED C C .. LOCAL SCALARS .. DOUBLE PRECISION DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. DOUBLE PRECISION DLOG10, DERFC C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT DOUBLE PRECISION S15AEF CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, RANTSTD C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL PRINTTITLE('DERFC ') 50000 FORMAT(1X,80('*'),/,' TEST OF DERFC USING RESULTS FROM FILE.') C C READ(NIN,50001) TEXT WRITE(NOUT,50002) TEXT 50001 FORMAT(A72) 50002 FORMAT(1X,A72) C C TEST TYPE 1 DO 40 J= 1,2 40 CALL RANTSTD(NIN, NOUT,J) C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= -10.0D0 WRITE(NOUT,1051) X I= 0 C Y= S15AEF(X,I) Y= DERFC(X) C Y= DERFC(X) WRITE(NOUT,1055) Y X= 20.0D0 WRITE(NOUT,1051) X Y= DERFC(X) WRITE(NOUT,1055) Y C WRITE(NOUT,1100) RETURN 1050 FORMAT(23H TEST OF EXTREME VALUES//) 1051 FORMAT(39H DERFC WILL BE CALLED WITH THE ARGUMENT,D15.4) 1052 FORMAT(39H0DERFC WILL BE CALLED WITH THE ARGUMENT,D15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H DERFC RETURNED THE VALUE,D15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE RANTSTD(NIN, NOUT,JJ) COMMON /MCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD DOUBLE PRECISION EPS,EPSNEG,XMAX,XMIN C TEST TYPE 1 C INPUT ON CR NO. NIN C OUTPUT ON LP NO.NOUT C .. SCALAR ARGUMENTS .. INTEGER NIN, NOUT C .. FUNCTION ARGUMENTS .. C .. C .. LOCAL SCALARS .. DOUBLE PRECISION EA, EM, ER, ERR, FEA, FEM, FER, PRTER, REPER, X, * XL, XU,Y,ALBETA,AIT,W,MRE,RMSX1 INTEGER I, IET, IFAIL, ISAM, I1, I2, J, KT, MULT, NX DOUBLE PRECISION S15AEF COMMON /PRMVLS/ ARGS,ARG2 DOUBLE PRECISION ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. DOUBLE PRECISION DSQRT ,DERFC C *** IMPLEMENTATION DEPENDENT DECLARATION *** C DOUBLE PRECISION DFLOAT C .. SUBROUTINE REFERENCES .. C ERRCAL, PRNT, READTC C .. DATA ETP(1), ETP(2), ETP(3), ETP(4), ETP(5), ETP(6) /1HA,1HB,1HS, *1HR,1HE,1HL/ DATA SAM(1), SAM(2), SAM(3), SAM(4), SAM(5), SAM(6) /1HU,1HN,1HI, *1HL,1HO,1HG/ ALBETA= DLOG(16.0D0) IBETA= 16 X1= 0.0D0 AIT= 14.0D0 READ (NIN,99999) NX, XL, XU, ISAM, IET, MULT I2 = ISAM*3 I1 = I2 - 2 WRITE (NOUT,99998) NX, XL, XU, (SAM(I),I=I1,I2) REPER = EPS PRTER = DFLOAT(MULT)*REPER I2 = IET*3 I1 = I2 - 2 IF (JJ.EQ.1) 1 WRITE (NOUT,99997) (ETP(I),I=I1,I2), REPER, MULT IF (JJ.EQ.2) 1 WRITE (NOUT,89997) (ETP(I),I=I1,I2), REPER, MULT EA = 0.0D0 ER = 0.0D0 EM = 0.0D0 KT = 0 DO 40 J=1,NX CALL READTC(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SDDAAF/ IFAIL = 1 Y= DERFC(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL ERRCAL(Y, IET, ERR, IFAIL) IF (JJ.EQ.2) ERR= 0.5D0*ERR/(X*X) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SDDAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. DABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL PRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF (DABS(ERR).GT.DABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL PRNT(NOUT, X, Y, 1.0D0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = DSQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0D0 IF (EM .NE. 0.0D0) W= DLOG(DABS(EM))/ALBETA IF (JJ.EQ.1) WRITE (NOUT,1021) EM,IBETA,W,X1 IF (JJ.EQ.2) WRITE (NOUT,2021) ER,IBETA,W,X1 MRE = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,MRE W = -999.0D0 IF (ER .NE. 0.0D0) W = DLOG(DABS(ER))/ALBETA IF (JJ.EQ.1) WRITE (NOUT,1023) ER,IBETA,W IF (JJ.EQ.2) WRITE (NOUT,2023) ER,IBETA,W RMS = DMAX1(AIT+W,0.0D0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C CALL TIMERI(JJ,NX,22,ARGS,MRE,RMS,NOUT) 1021 FORMAT(21H THE MAXIMUM ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 1022 FORMAT(27H THE ESTIMATED LOSS OF BASE,I4, * 22H SIGNIFICANT DIGITS IS,F7.2//) 1023 FORMAT(31H THE ROOT MEAN SQUARE ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) 2021 FORMAT(32H THE MAXIMUM SCALED-REL-ERROR OF,D15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,D17.6) 2023 FORMAT(42H THE ROOT MEAN SQUARE SACLED-REL-ERROR WAS,D15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PD11.2, 5X, D11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PD11.2, * 5H AND , D11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PD11.2/46H PRINT IN FULL FIRST POINT, ALSO PRINT IN FULL/ * 14H IF(ERROR.GE. , I4, 10H*REP.ERR.)/1H0, 29X, 1HX, 29X, 1HY, * 26X, 4HYEXP, 4X, 11HERR/REP.ERR) 89997 FORMAT (15H ERROR MEASURE , 3A1,'/2*X**2' 1 , 25H, REPRESENTATION ERROR = , * 1PD11.2/46H PRINT IN FULL FIRST POINT, ALSO PRINT IN FULL/ * 14H IF(ERROR.GE. , I4, 10H*REP.ERR.)/1H0, 29X, 1HX, 29X, 1HY, * 26X, 4HYEXP, 4X, 11HERR/REP.ERR) 99996 FORMAT (12H FAILURE NO., I2, 10H CAUSED BY) 99995 FORMAT (13H0MEAN ERROR =, 1PD11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PD11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END