PROGRAM STEST C C DRIVER FOR REAL MATH FUNCTION TESTS C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C INTEGER I,J C CALL EMAS3('DEFINE','10,STESTOUT,,C',J) CALL EMAS3('DEFINE','9,STESTINPUT',J) CALL SSETMESAGES WRITE(6,10000) 10000 FORMAT(' TESTS ARE NUMBERED AS FOLLOWS:',/ 1 ,' 0 - H/W PARMS 4 - ''**'' 8 - ATAN/ ATAN2',/ 2 ,' 1 - SQRT 5 - SIN/ COS 9 - SINH/ COSH',/ 3 ,' 2 - ALOG/ALOG10 6 - TAN/ COTAN 10 - TANH',/ 4 ,' 3 - EXP 7 - ASIN/ ACOS 11 - GAMMA,ALGAMA' 5 ,', ERF, ERFC',/ 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 STEST0 GOTO 99 C 105 CALL STEST5 GOTO 99 C 102 CALL STEST2 GOTO 99 C 103 CALL STEST3 GOTO 99 C 101 CALL STEST1 GOTO 99 C 108 CALL STEST8 GOTO 99 C 104 CALL STEST4 GOTO 99 C 106 CALL STEST6 GOTO 99 C 107 CALL STEST7 GOTO 99 C 109 CALL STEST9 GOTO 99 C 110 CALL STEST10 GOTO 99 C 111 CALL STEST11 GOTO 99 C 300 CALL STEST0 CALL STEST1 CALL STEST2 CALL STEST3 CALL STEST4 CALL STEST5 CALL STEST6 CALL STEST7 CALL STEST8 CALL STEST9 CALL STEST10 CALL STEST11 STOP END C SUBROUTINE SPRINTITLE(TEXT) CHARACTER*30 TEXT C WRITE(6,10000) TEXT 10000 FORMAT(//,1X,A30 1 ,/,21X,'J MRE RMS MICRO-SECS') RETURN END C SUBROUTINE SSUMARY(J,MRE,RMS,MSECS) INTEGER J REAL 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 STIMERI(J,N,IFUNC,ARGS,MRE,RMS,IOUT) INTEGER J,N,IOUT,IFUNC REAL ARGS(10000),ARGS2(10000),MRE,RMS C INTEGER I,K,KK REAL Y DOUBLE PRECISION T0,T1,T2,T3 C KK= 40000/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 STIMERI2(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= SQRT(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= ALOG(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= EXP(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= SIN(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= TAN(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= ASIN(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= ATAN(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= SINH(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= TANH(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= ALOG10(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= COS(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= COTAN(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= ACOS(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= ATAN2(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= COSH(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= GAMMA(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= ALGAMA(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= ERF(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= ERFC(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 SSUMARY(J,MRE,RMS,T3) RETURN END C SUBROUTINE SSETMESAGES C C CHARACTER*55 MESS(45) C C INTEGER ALOGSMALL, SQRTNEG, EXPLARGE, EXPSMALL, SINLARGE 1 , ASINLARGE, COSLARGE, ACOSLARGE, TANLARGE, COTANSMALL 2 ,ARCSINLARGE,ARCCOSLARGE,ARCTAN2ZERO, SINHLARGE, COSHLARGE 3 ,POWERNEG,POWERZERO,POWERBIG,POWERSMALL 4 ,GAMLARGE,GAMNEGINT,LGAMNEG,LGAMLARGE C ALOGSMALL= 1 ALOGSMALL = 21 SQRTNEG = 18 EXPLARGE = 19 EXPSMALL = 20 SINLARGE = 26 ASINLARGE = 29 COSLARGE = 27 ACOSLARGE = 30 TANLARGE = 28 COTANSMALL = 41 ARCTAN2ZERO = 31 SINHLARGE = 32 COSHLARGE = 33 POWERNEG = 23 POWERZERO = 22 POWERBIG = 24 POWERSMALL = 25 GAMLARGE = 42 GAMNEGINT = 43 LGAMNEG = 44 LGAMLARGE = 45 MESS(ALOGSMALL )= ' ALOG ARG NEGATIVE OR ZERO' MESS(SQRTNEG )= ' SQRT ARG NEGATIVE' MESS(EXPLARGE )= ' EXP ARG GREATER THAN 174.6731' MESS(EXPSMALL )= ' EXP ARG LESS THAN -180.2182, NOT USED' MESS(SINLARGE )= ' SIN ARG MAGNITUDE GREATER THAN 5.274E7' MESS(ASINLARGE )= ' ASIN ARG MAGNITUDE GREATER THAN 1.0' MESS(COSLARGE )= ' COS ARG MAGNITUDE GREATER THAN 5.274E7' MESS(ACOSLARGE )= ' ACOS ARG MAGNITUDE GREATER THAN 1.0' MESS(TANLARGE )= 1 ' TAN/COTAN ARG MAGNITUDE GREATER THAN 2.633E7' MESS(COTANSMALL )= ' COTAN ARG SMALL IN MAGNITUDE' MESS(ARCTAN2ZERO)= ' ATAN2 ARGS ZERO' MESS(SINHLARGE )= ' SINH ARG MAGNITUDE GREATER THAN 175.3662 ' MESS(COSHLARGE )= ' COSH ARG MAGNITUDE GREATER THAN 175.3662 ' MESS(POWERNEG )= 1 ' NEGATIVE S.P. VALUE RAISED TO A NON-INTEGER POWER' MESS(POWERZERO )= ' S.P. ZERO RAISED TO NON-POSITIVE POWER' MESS(POWERBIG )= ' S.P. VALUE RAISED TO S.P. POWER TOO BIG' MESS(POWERSMALL )= 1 ' S.P. VALUE TO S.P. POWER UNDERFLOWS. NOT USED' MESS(GAMLARGE )= ' GAMMA ARG MAGNITUDE GREATER THAN 57.0' MESS(GAMNEGINT )= ' GAMMA ARG NEAR ZERO OR NEGATIVE INTEGER' MESS(LGAMNEG )= ' ALGAMA ARG IS NEGATIVE' MESS(LGAMLARGE )= ' ALGAMA ARG IS GREATER THAN 4.29E73' RETURN C ENTRY SERROR(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 SMACHAR(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 REAL A,B,BETA,BETAIN,BETAM1,EPS,EPSNEG,ONE,XMAX,XMIN, *Y,Z,ZERO CD REAL 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 REAL 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.0E0 CD ZERO = 0.0E0 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 = NINT((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. ( ABS(Z) .GE. Y)) GO TO 410 CD IF ((A+A .EQ. ZERO) .OR. ( ABS(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. ( ABS(Y) .GE. XMIN)) GO TO 460 CD IF (((A+A) .EQ. ZERO) .OR. ( ABS(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 SMACHAR ---------- END C C******************************************************************* C REAL FUNCTION SRAN(I) DOUBLE PRECISION LOCRAN 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 LOCRAN= FLOAT(IY)/2796203.0E0 SRAN= LOCRAN C RETURN END C C******************************************************************* C SUBROUTINE STEST0 C PROGRAM TO MONITOR SUBROUTINE SMACHAR, REAL C C DATA REQUIRED C C NONE C C INTEGER IBETA, IEXP, IRND, IT, MACHEP, MAXEXP, MINEXP,NEGEP, NGRD INTEGER IOUT REAL EPS, EPSNEG, XMIN, XMAX C C OUTPUT UNIT SET TO 6 C IOUT = 10 CALL SMACHAR(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 (30H OUTPUT OF VALUES FROM SMACHAR/ * 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 STEST1 C PROGRAM TO TEST SQRT C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, SQRT 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 REAL A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,RANDL,R6,R7, *SQBETA,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO, SQRT,MRE,RMS REAL SRAN C C THE FOLLOWING STATEMENT FUNCTION HAS BEEN ADDED TO C PROVIDE THE REQUIRED FUNCTION RANDL C RANDL(X) = EXP(SRAN(0)*X) C CALL SPRINTITLE(' SQRT ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) SQBETA = SQRT(BETA) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0E0 ZERO = 0.0E0 A = ONE / SQBETA B = ONE N = 2000 XN = DBLE(FLOAT(N)) C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 300 J = 1, 2 C = ALOG(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 = SQRT(Y) ARGS(I)= Y W = (Z - X) / X IF (W .GT. ZERO) K1 = K1 + 1 IF (W .LT. ZERO) K3 = K3 + 1 W = ABS(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 = SQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0E0 IF (R6 .NE. ZERO) W= ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = ONE B = SQBETA C C CPU TIME CALCULATION AND SUMMARY PRINT C CALL STIMERI(J,N,1,ARGS,MRE,RMS,IOUT) 300 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1040) X = XMIN Y = SQRT(X) WRITE (IOUT,1041) X,Y X = ONE - EPSNEG Y = SQRT(X) WRITE (IOUT,1042) EPSNEG,Y X = ONE Y = SQRT(X) WRITE (IOUT,1043) X,Y X = ONE + EPS Y = SQRT(X) WRITE (IOUT,1044) EPS,Y X = XMAX Y = SQRT(X) WRITE (IOUT,1045) X,Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = ZERO WRITE (IOUT,1051) X Y = SQRT(X) WRITE (IOUT,1055) Y X = -ONE WRITE (IOUT,1052) X Y = SQRT(X) WRITE (IOUT,1055) Y WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,' TEST OF SQRT(X*X) - X '//) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(20H SQRT(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1040 FORMAT(26H TEST OF SPECIAL ARGUMENTS//) 1041 FORMAT(21H SQRT(XMIN) = SQRT(,E15.7,4H) = ,E15.7//) 1042 FORMAT(27H SQRT(1-EPSNEG) = SQRT(1-,E15.7,4H) = ,E15.7//) 1043 FORMAT(20H SQRT(1.0) = SQRT(,E15.7,4H) = ,E15.7//) 1044 FORMAT(24H SQRT(1+EPS) = SQRT(1+,E15.7,4H) = ,E15.7//) 1045 FORMAT(21H SQRT(XMAX) = SQRT(,E15.7,4H) = ,E15.7//) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(39H SQRT WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(39H0 SQRT WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H SQRT RETURNED THE VALUE,E15.4///) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF SQRT TEST PROGRAM ---------- END C C*********************************************************** C SUBROUTINE STEST5 C PROGRAM TO TEST SIN/ COS C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, COS, SIN, SQRT 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 REAL 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 REAL SRAN REAL SIN, COS C C=================================================================== C CALL SPRINTITLE(' SIN (J=1,2), COS (J=3) ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0E0 ZERO = 0.0E0 THREE = 3.0E0 A = ZERO B = 1.570796327E0 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 * SRAN(I1) + XL Y = X / THREE Y = (X + Y) - X X = THREE * Y ARGS(I)= Y IF (J .EQ. 3) GO TO 100 Z = SIN(X) ZZ = SIN(Y) W = ONE IF (Z .NE. ZERO) W = (Z - ZZ*(THREE-4.0E0*ZZ*ZZ)) / Z GO TO 110 100 Z = COS(X) ZZ = COS(Y) W = ONE IF (Z .NE. ZERO) W = (Z + ZZ*(THREE-4.0E0*ZZ*ZZ)) / Z 110 IF (W .GT. ZERO) K1 = K1 + 1 IF (W .LT. ZERO) K3 = K3 + 1 W = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = 18.84955592E0 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 STIMERI(J,N,5,ARGS,MRE,RMS,IOUT) ELSE CALL STIMERI(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 = ( SIN(A+C) - SIN(A-C)) / (C + C) WRITE (IOUT,1026) Z C WRITE (IOUT,1030) C DO 320 I = 1, 5 X = SRAN(I1) * A Z = SIN(X) + SIN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = SRAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - SIN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) C DO 340 I = 1, 5 X = SRAN(I1) * A Z = COS(X) - COS(-X) WRITE (IOUT,1060) X, Z 340 CONTINUE C WRITE (IOUT,1035) EXPON = DBLE(FLOAT(MINEXP)) * 0.75E0 X = BETA ** EXPON Y = SIN(X) WRITE (IOUT,1061) X, Y WRITE (IOUT,1040) Z = SQRT(BETAP) X = Z * (ONE - EPSNEG) Y = SIN(X) WRITE (IOUT,1061) X, Y Y = SIN(Z) WRITE (IOUT,1061) Z, Y X = Z * (ONE + EPS) Y = SIN(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = BETAP WRITE (IOUT,1052) X Y = SIN(X) WRITE (IOUT,1055) Y WRITE(IOUT,1072) X Y= COS(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,' TEST OF SIN(X) VS 3* SIN(X/3)-4* SIN(X/3)**3' * ,//) 1005 FORMAT(' TEST OF COS(X) VS 4* COS(X/3)**3-3* COS(X/3)' * ,//) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(19H SIN(X) WAS LARGER,I6,7H TIMES, / * 12X,7H AGREED,I6,11H TIMES, AND / * 8X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(19H COS(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1026 FORMAT(4H IF ,E13.6,21H IS NOT ALMOST 1.0E0,, * 4X,26H SIN HAS THE WRONG PERIOD. //) 1030 FORMAT(53H THE IDENTITY SIN(-X) = - SIN(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(52H THE IDENTITY SIN(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(52H THE IDENTITY COS(-X) = COS(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 SIN WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H SIN RETURNED THE VALUE,E15.4///) 1060 FORMAT(2E15.7/) 1072 FORMAT(38H COS WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(24H COS RETURNED THE VALUE,E15.4///) 1061 FORMAT(/6X,6H SIN(,E13.6,3H) =,E13.6) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF SIN/ COS TEST PROGRAM ---------- END C C***************************************************************** C SUBROUTINE STEST2 C PROGRAM TO TEST ALOG C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, ALOG10, AMAX1, SIGN, SQRT 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 REAL A,AIT,ALBETA,B,BETA,C,EIGHT,EPS,EPSNEG,HALF, *ONE,R6,R7,TENTH,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ C REAL SRAN REAL MRE,RMS,ALOG,ALOG10,CC1,CC2 DOUBLE PRECISION YY,XX,XL,DEL C DATA CC1/Z3D518600/, CC2/Z3E18B605/ C CALL SPRINTITLE('ALOG/ALOG10 ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) J = IT / 3 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.0E0 HALF = 0.5E0 EIGHT = 8.0E0 TENTH = 0.1E0 ONE = 1.0E0 C = ONE C DO 50 I = 1, J C = C / BETA 50 CONTINUE C B = ONE + C A = ONE - C N = 2000 C READ(5,*) N 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 = DBLE(B - A) / DBLE(XN) XL = A C AINCR=0.01 C X=0.5 C DO 200 I = 1, N C X=X+AINCR XX= DEL*SRAN(I1) + XL X= XX ARGS(I)= X IF (J.NE.1) GOTO 100 NNCNT= 0 9997 CONTINUE IF (X.NE.ONE) GOTO 9996 X= (B - A)*SRAN(I1) + A NNCNT= NNCNT + 1 IF (NNCNT.LE.10000) GOTO 9997 C 9996 CONTINUE C X= X + 16.0E0 C X= X - 16.0E0 ARGS(I)= X Y = (X - HALF) - HALF C Y= Y - HALF ZZ = ALOG(X) C XX= X C XX= LOG(XX) C X= XX C YY= XX + (XX - X) C Z= YY Z = ONE / 3.0E0 Z = Y * (Z - Y / 4.0E0) Z = (Z - HALF) * Y * Y + Y C W= ZERO GOTO 150 C 100 IF (J.NE.2) GOTO 110 X= (X + EIGHT) X= X - EIGHT Y= X + X/16.0E0 Z= ALOG(X) C ZZ= ALOG(Y) - 7.7746816434842581E-5 ZZ= ALOG(Y) - CC1 ZZ= ZZ - 31.0E0/512.0E0 GOTO 150 C 110 IF (J.NE.3) GOTO 120 X= (X + EIGHT) X= X - EIGHT Y= X + X*TENTH Z= ALOG10(X) C XX= 3.7706015822504075E-4 C ZZ= XX ZZ= ALOG10(Y) - CC2 ZZ= ZZ - 21.0E0/512.0E0 GOTO 150 C 120 CONTINUE Z= ALOG(X*X) ZZ= ALOG(X) ZZ= ZZ + ZZ C 150 CONTINUE W = ONE 151 CONTINUE ZP = Z IF (Z .NE. ZERO) W = (Z - ZZ) / Z Z = SIGN(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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .GT. 1) GO TO 230 A = SQRT(HALF) B = 15.0E0 / 16.0E0 GO TO 300 230 IF (J .GT. 2) GO TO 240 A = SQRT(TENTH) B = 0.9E0 GO TO 300 240 A = 16.0E0 B = 240.0E0 300 CONTINUE C C CPU TIME CALCULATION C IF (J.NE.3) THEN CALL STIMERI(J,N,2,ARGS,MRE,RMS,IOUT) ELSE CALL STIMERI(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 = SRAN(I1) X = X + X + 15.0E0 Y = ONE / X Z = ALOG(X) + ALOG(Y) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1040) X = ONE Y = ALOG(X) WRITE (IOUT,1041) Y X = XMIN Y = ALOG(X) WRITE (IOUT,1042) X, Y X = XMAX Y = ALOG(X) WRITE (IOUT,1043) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = -2.0E0 WRITE (IOUT,1052) X Y = ALOG(X) WRITE (IOUT,1055) Y X = ZERO WRITE (IOUT,1052) X Y = ALOG(X) WRITE (IOUT,1055) Y WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,' TEST OF ALOG(X) VS T.S. EXPANSION OF ALOG(1+Y)'//) 1001 FORMAT(' TEST OF ALOG(X) VS ALOG(17X/16)-ALOG(17/16)' //) 1002 FORMAT(' TEST OF ALOG(X*X) VS 2 * ALOG(X)' //) 1005 FORMAT(' TEST OF ALOG10(X) VS ALOG10(11X/10)' 1 ,'-ALOG10(11/10)' //) 1009 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,26H(1-EPS,1+EPS), WHERE EPS =, E15.4//) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(' ALOG(X) WAS LARGER',I6,7H TIMES, / * 15X,7H AGREED,I6,11H TIMES, AND / * 11X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(' ALOG10(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(52H THE IDENTITY ALOG(X) = -ALOG(1/X) WILL BE TESTED.// * 8X,1HX,9X,13HF(X) + F(1/X)/) 1040 FORMAT(//26H TEST OF SPECIAL ARGUMENTS //) 1041 FORMAT(13H ALOG(1.0) = ,E15.7//) 1042 FORMAT(19H ALOG(XMIN) = ALOG(,E15.7,4H) = ,E15.7//) 1043 FORMAT(19H ALOG(XMAX) = ALOG(,E15.7,4H) = ,E15.7//) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1052 FORMAT(38H ALOG WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H ALOG RETURNED THE VALUE,E15.4///) 1060 FORMAT(2E15.7/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF ALOG/ALOG10 TEST PROGRAM ---------- END C C***************************************************************** C SUBROUTINE STEST3 C PROGRAM TO TEST EXP C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, EXP, SQRT 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 AINT IS USED. THIS IS NOT IN THE C 1966 STANDARD, AND SPECIAL ACTION WILL BE NECESSARY IF THE C COMPILER LIBRARY DOES NOT PROVIDE AINT. 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 REAL 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 REAL SRAN REAL EXP C EXTERNAL EXP REAL XX,YY,CC1,CC2 C DATA CC1/Z3EA041DD/, CC2/Z3FF82A02/ C CALL SPRINTITLE(' EXP ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0E0 TWO = 2.0E0 TEN = 10.0E0 ZERO = 0.0E0 V = 0.0625E0 A = TWO B = ALOG(A) * 0.5E0 A = -B + V C C D = ALOG(0.9E0*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 * SRAN(I1) + XL C--------------------------------------------------------------------- C PURIFY ARGUMENTS C--------------------------------------------------------------------- Y = X - V IF (Y .LT. ZERO) X = Y + V ARGS(I)= X Z = EXP(X) ZZ = EXP(Y) IF (J .EQ. 1) GO TO 100 IF (IBETA .EQ. 10) GOTO 302 C XX= 2.4453321046920570389E-3 C XXRND= XX Z = Z * .0625E0 -Z *CC1 GOTO 110 C 302 CONTINUE IF (IBETA .EQ. 10) Z = Z * 6.0E-2 +Z * * 5.466789530794296106E-5 GO TO 110 100 CONTINUE C XX= 6.058693718652421388E-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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .EQ. 2) GO TO 270 V = 45.0E0 / 16.0E0 A = -TEN * B B = 4.0E0 * XMIN * BETA ** IT B = ALOG(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 STIMERI(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 = SRAN(I1) * BETA Y = -X Z = EXP(X) * EXP(Y) - ONE WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1040) X = ZERO Y = EXP(X) - ONE WRITE (IOUT,1041) Y X = AINT(ALOG(XMIN)) Y = EXP(X) WRITE (IOUT,1042) X, Y X = AINT(ALOG(XMAX)) C WRITE(6,1)X 1 FORMAT(1X,'X= ',E15.7) Y = EXP(X) WRITE (IOUT,1042) X, Y X = X / TWO V = X / TWO Y = EXP(X) Z = EXP (V) Z = Z * Z WRITE (IOUT,1043) X, Y, V, Z C--------------------------------------------------------------------- C TEST OF ERROR RETURNS C--------------------------------------------------------------------- WRITE (IOUT,1050) X = -ONE / SQRT(XMIN) WRITE (IOUT,1052) X Y = EXP(X) WRITE (IOUT,1061) Y X = -X WRITE (IOUT,1052) X Y = EXP(X) WRITE (IOUT,1061) Y WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,16H TEST OF EXP(X-,F7.4,18H) VS EXP(X)/ EXP(,F7.4,1H) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(21H EXP(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(54H THE IDENTITY EXP(X)* EXP(-X) = 1.0 WILL BE TESTED.// * 8X,1HX,9X,14HF(X)*F(-X) - 1 /) 1040 FORMAT(//26H TEST OF SPECIAL ARGUMENTS //) 1041 FORMAT(21H EXP(0.0) - 1.0E0 = ,E15.7/) 1042 FORMAT(6H EXP(,E13.6,3H) =,E13.6/) 1043 FORMAT(9H0IF EXP(,E13.6,4H) = ,E13.6,13H IS NOT ABOUT / * 6H EXP(,E13.6,7H)**2 = ,E13.6,26H THERE IS AN ARG RED ERROR) 1050 FORMAT(22H TEST OF ERROR RETURNS //) 1052 FORMAT(38H0 EXP WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1060 FORMAT(2E15.7/) 1061 FORMAT(24H EXP RETURNED THE VALUE,E15.4///) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF EXP TEST PROGRAM ---------- END C C***************************************************************** C SUBROUTINE STEST4 C PROGRAM TO TEST POWER FUNCTION (**), REAL C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, EXP, SQRT 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 REAL 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 REAL SRAN,POWQI CALL SPRINTITLE('X**Y, X**N ') IOUT= 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) C C ALXMAX = ALOG(XMAX) C ALXMAX = ALOG(XMAX*0.1) ZERO = 0.0E0 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 = -AMAX1(ALXMAX,-ALOG(XMIN))/ALOG(100.0E0) C C= -10 DELY = -C - C N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 Y1 = ZERO LNINTB= LOG(2.0E0**31) C----------------------------------------------------------------- C RANDOM ARGUMENT ACCURACY TESTS C----------------------------------------------------------------- DO 301 J = 1,5 K1 = 0 K3 = 0 X1 = ZERO R6 = ZERO R7 = ZERO DEL = (B - A) / XN XL = A C DO 200 I = 1, N X = DEL * SRAN(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.0E0 + 1E-14 Y = DELY * SRAN(I1) + C IF (J.GT.4) GOTO 5110 Y2 = (Y/TWO + Y) - Y Y = Y2 + Y2 ARGS2(I)= Y C 2.0E0**5 - 1.0E0 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.0E0*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.7E0*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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(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 = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS GOTO (300,302,303,304,305,300), J 302 CONTINUE A = ONE B = EXP(ALXMAX/3.0E0) GOTO 300 C 303 CONTINUE B = 10.0E0 A = 0.01E0 GOTO 300 C 304 CONTINUE A= -100.0 B= 100.0E0 C= -100.0E0 DELY= -C - C GOTO 300 C 305 CONTINUE A= -100.0 B= 100.0 C= ZERO DELY= 30.0E0 300 CONTINUE C C CPU TIME CALCULATION AND SUMMARY PRINT IF (J.LE.4) THEN NN= 4 ELSE NN= 19 + J ENDIF CALL STIMERI2(J,N,NN,ARGS,ARGS2,MRE,RMS,IOUT) C 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) B = 10.0E0 C DO 320 I = 1, 5 X = SRAN(I1) * B + ONE Y = SRAN(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 C J= 0 C I= 0 C WRITE(IOUT,1054) J,I C K= J**I C WRITE(IOUT,1056) K CC C J= 2 C I= -1 C WRITE(IOUT,1054) J,I C K= J**I C WRITE(IOUT,1056) K CC CC C J= 2 C I= 33 C WRITE(IOUT,1054) J,I C K= J**I C WRITE(IOUT,1056) K CC 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(,E15.4,1H,,E15.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 (,E15.4,1H,,E15.4,9H), Y IN (,E15.4,1H,,E15.4,1H)//) 5014 FORMAT(I7,45H RANDOM ARGUMENTS WERE TESTED FROM THE REGION / * 6X,6HX IN (,E15.4,1H,,E15.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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1024 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.6,4H Y =,E17.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 (,E14.7,7H ) ** (,E14.7,20H ) WILL BE COMPUTED.,/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(2H (,E14.7,7H ) ** (,E14.7,20H ) WILL BE COMPUTED.,/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1053 FORMAT(2H (,E14.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,E15.4///) 1056 FORMAT(22H THE VALUE RETURNED IS,I12///) 1060 FORMAT(2E15.7,6X,E15.7/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF POWER TEST PROGRAM (D.P.) ---------- END C****************************************************************************** SUBROUTINE STEST8 C PROGRAM TO TEST ATAN, ATAN2 C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, ATAN, ATAN2, SQRT 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 REAL 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, ATAN, ATAN2 * ,MRE,RMS C REAL SRAN C EXTERNAL ATAN, ATAN2 C C CALL SPRINTITLE(' ATAN (AND ATAN2) ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) ONE = 1.0E0 HALF = 0.5E0 TWO = 2.0E0 ZERO = 0.0E0 A = -0.0625E0 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 * SRAN(I1) + XL IF (J .EQ. 2) X = ((1.0E0+X*A)-ONE)*16.0E0 ARGS(I)= X Z = ATAN(X) IF (J .NE. 1) GO TO 100 XSQ = X * X EM = 17.0E0 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 - .0625E0 Y = Y / (ONE + X*A) ZZ = ( ATAN(Y) - 8.1190004042651526021E-5) + OB32 ZZ = ZZ + OB32 GO TO 110 105 Z = Z + Z Y = X / ((HALF + X * HALF)*((HALF - X) + HALF)) ZZ = ATAN(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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = B IF (J .EQ. 1) B = TWO - SQRT(3.0E0) IF (J .EQ. 2) B = SQRT(TWO) - ONE IF (J .EQ. 3) B = ONE 300 CONTINUE C C CPU TIME CALCULATION AND SUMMARY PRINT. C CALL STIMERI(J,N,8,ARGS,MRE,RMS,IOUT) 301 CONTINUE C----------------------------------------------------------------- C SPECIAL TESTS C----------------------------------------------------------------- WRITE (IOUT,1025) WRITE (IOUT,1030) A = 5.0E0 C DO 320 I = 1, 5 X = SRAN(I1) * A Z = ATAN(X) + ATAN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = SRAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - ATAN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) A = -TWO B = 4.0E0 C DO 340 I = 1, 5 X = SRAN(I1) * B + A Y = SRAN(I1) W = -Y Z = ATAN(X/Y) - ATAN2(X,Y) ZZ = ATAN(X/W) - ATAN2(X,W) WRITE (IOUT,1059) X, Y, Z, ZZ 340 CONTINUE C WRITE (IOUT,1035) EXPON = DBLE(FLOAT(MINEXP)) * 0.75E0 X = BETA ** EXPON Y = ATAN(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) WRITE (IOUT,1051) XMAX Z = ATAN(XMAX) WRITE (IOUT,1061) XMAX, Z X = ONE Y = ZERO WRITE (IOUT,1053) X, Y Z = ATAN2(X,Y) WRITE (IOUT,1062) X, Y, Z WRITE (IOUT,1053) XMIN, XMAX Z = ATAN2(XMIN,XMAX) WRITE (IOUT,1062) XMIN, XMAX, Z WRITE (IOUT,1053) XMAX, XMIN Z = ATAN2(XMAX,XMIN) WRITE (IOUT,1062) XMAX, XMIN, Z X = ZERO WRITE (IOUT,1054) X, Y Z = ATAN2(X,Y) WRITE (IOUT,1062) X, Y, Z WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,44H TEST OF ATAN(X) VS TRUNCATED TAYLOR SERIES //) 1001 FORMAT(21H TEST OF ATAN(X) VS , * 42H ATAN(1/16) + ATAN((X-1/16)/(1+X/16)) //) 1002 FORMAT(42H TEST OF 2* ATAN(X) VS ATAN(2X/(1-X*X)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(20H ATAN(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(55H THE IDENTITY ATAN(-X) = - ATAN(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY ATAN(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(53H THE IDENTITY ATAN(X/Y) = ATAN2(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 ATAN WILL BE CALLED WITH THE ARGUMENT,E15.7/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE/) 1053 FORMAT(41H ATAN2 WILL BE CALLED WITH THE ARGUMENTS/2E15.7/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE/) 1054 FORMAT(41H ATAN2 WILL BE CALLED WITH THE ARGUMENTS/2E15.7/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1059 FORMAT(4(E15.7,4X)/) 1060 FORMAT(2E15.7/) 1061 FORMAT(6X,7H ATAN(,E13.6,3H) =,E13.6/) 1062 FORMAT(6X,8H ATAN2(,2E13.6,3H) =,E13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF ATAN/ ATAN2 TEST PROGRAM ---------- END C C******************************************************************* C SUBROUTINE STEST10 C PROGRAM TO TEST TANH C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C THIS REFERS TO THE FUNCTION TANH 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 SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, SQRT 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 REAL 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 * , TANH,MRE,RMS C REAL SRAN C C EXTERNAL TANH C CALL SPRINTITLE(' TANH ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) AIT = DBLE(FLOAT(IT)) ZERO = 0.0E0 ONE = 1.0E0 HALF = 0.5E0 A = 0.125E0 B = ALOG(3.0E0) * HALF C = 1.2435300177159620805E-1 D = ALOG(2.0E0) + (AIT+ONE) * ALOG(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 * SRAN(I1) + XL ARGS(I)= X Z = TANH(X) Y = X - 0.125E0 ZZ = TANH(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 = ABS(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 = SQRT(R7/XN) WRITE (IOUT,1000) WRITE (IOUT,1010) N,A,B WRITE (IOUT,1011) K1,K2,K3 WRITE (IOUT,1020) IT,IBETA W = -999.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS A = B + A B = D CALL STIMERI(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 = SRAN(I1) Z = TANH(X) + TANH(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = SRAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - TANH(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) X = D B = 4.0E0 C DO 340 I = 1, 5 Z = ( TANH(X) - HALF) - HALF WRITE (IOUT,1060) X, Z X = X + B 340 CONTINUE C WRITE (IOUT,1035) EXPON = DBLE(FLOAT(MINEXP)) * 0.75E0 X = BETA ** EXPON Z = TANH(X) WRITE (IOUT,1061) X, Z WRITE (IOUT,1040) XMAX Z = TANH(XMAX) WRITE (IOUT,1061) XMAX, Z WRITE (IOUT,1040) XMIN Z = TANH(XMIN) WRITE (IOUT,1061) XMIN, Z X = ZERO WRITE (IOUT,1040) X Z = TANH(X) WRITE (IOUT,1061) X, Z WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,21H TEST OF TANH(X) VS , * 54H( TANH(X-1/8)+ TANH(1/8))/(1+ TANH(X-1/8) TANH(1/8)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(20H TANH(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(54H THE IDENTITY TANH(-X) = - TANH(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY TANH(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(53H THE IDENTITY TANH(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 TANH WILL BE CALLED WITH THE ARGUMENT, * E15.7) 1060 FORMAT(2E15.7/) 1061 FORMAT(6X,7H TANH(,E13.6,3H) =,E13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF TANH TEST PROGRAM ---------- END C************************************************************************ C********************************************************************* C SUBROUTINE STEST6 C PROGRAM TO TEST TAN/ COTAN C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C THIS REFERS TO FUNCTIONS TAN AND COTAN. NEITHER WERE C INCLUDED IN THE 1966 STANDARD, AND COTAN 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 TAN IS AVAILABLE BUT NOT COTAN, THEN THE TEST CAN BE RUN C PROVIDED THAT : C COTAN IS DECLARED REAL C A DUMMY FUNCTION COTAN IS ADDED C C C DATA REQUIRED C C NONE C C SUBPROGRAMS REQUIRED FROM THIS PACKAGE C C SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, SQRT 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 REAL 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 REAL TAN, COTAN,MRE,RMS REAL SRAN C C EXTERNAL TAN, COTAN C C THE FOLLOWING STATEMENT FUNCTION SHOULD BE ACTIVATED IF C THIS TEST IS TO BE RUN BUT NO FUNCTION COTAN (UNDER ANY C NAME) IS AVAILABLE. C C COTAN(Z) = 0.0E0 C CALL SPRINTITLE (' TAN (J=1,2,3), COTAN (J=4) ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) ZERO = 0.0E0 HALF = 0.5E0 AIT = DBLE(FLOAT(IT)) PI = 3.14159265E0 A = ZERO B = PI * 0.25E0 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 TAN AND COTAN C J = 1,3 IF TESTING TAN 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 * SRAN(I1) + XL Y = X * HALF X = Y + Y ARGS(I)= X IF (J .EQ. 4) GO TO 80 Z = TAN(X) ZZ = TAN(Y) W = 1.0E0 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 = COTAN(X) ZZ = COTAN(Y) W = 1.0E0 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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .NE. 1) GO TO 250 A = PI * 0.875E0 B = PI * 1.125E0 GO TO 300 250 A = PI * 6.0E0 B = A + PI * 0.25E0 300 CONTINUE IF (J.NE.4) THEN CALL STIMERI(J,N,6,ARGS,MRE,RMS,IOUT) ELSE CALL STIMERI(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 = SRAN(I1) * A Z = TAN(X) + TAN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = SRAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - TAN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1035) X = BETA ** (DBLE(FLOAT(MINEXP))*0.75E0) Y = TAN(X) WRITE (IOUT,1061) X, Y C1 = -225.0E0 C2 = -.950846454195142026E0 X = 11.0E0 Y = TAN(X) W = ((C1-Y)+C2)/(C1+C2) Z = ALOG( ABS(W))/ALBETA WRITE (IOUT,1040) W, IBETA, Z WRITE (IOUT,1061) X, Y W = AMAX1(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 = TAN(X) WRITE (IOUT,1055) Y X = BETAP WRITE (IOUT,1052) X Y = TAN(X) WRITE (IOUT,1055) Y X= 0.0E0 WRITE(IOUT,1072) X Y= COTAN(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,48H TEST OF TAN(X) VS 2* TAN(X/2)/(1- TAN(X/2)**2) //) 1005 FORMAT(40H TEST OF COTAN(X) VS ( COTAN(X/2)**2-1), * 16H/(2* COTAN(X/2)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(19H TAN(X) WAS LARGER,I6,7H TIMES, / * 12X,7H AGREED,I6,11H TIMES, AND / * 8X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(21H COTAN(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(51H THE IDENTITY TAN(-X) = - TAN(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(52H THE IDENTITY TAN(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 TAN(11) IS ,E15.7,3H = , * I4,3H **,F7.2,6H WHERE /) 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(38H TAN WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(38H0 TAN WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H TAN RETURNED THE VALUE,E15.4///) 1060 FORMAT(2E15.7/) 1061 FORMAT(6X,6H TAN(,E13.6,3H) =,E13.6/) 1071 FORMAT(' COTAN WILL BE CALLED WITH THE ARGUMENT',E15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1072 FORMAT('0 COTAN WILL BE CALLED WITH THE ARGUMENT',E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(' COTAN RETURNED THE VALUE',E15.4///) 1081 FORMAT(6X,' COTAN(',E13.6,3H) =,E13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF TAN/ COTAN TEST PROGRAM ---------- END C C******************************************************************** C SUBROUTINE STEST7 C PROGRAM TO TEST ASIN/ ACOS C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C THIS REFERS TO FUNCTIONS ASIN AND ACOS. 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 SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, ALOG10, AMAX1, NINT, SQRT 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 REAL 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 * , ASIN, ACOS,MRE,RMS C REAL SRAN C C EXTERNAL ASIN, ACOS C CALL SPRINTITLE(' ASIN (J=1,3), ACOS (J=2,4,5)') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) ZERO = 0.0E0 HALF = 0.5E0 AIT = DBLE(FLOAT(IT)) K = NINT(ALOG10(BETA**IT)) + 1 C WRITE(6,*) K C READ(5,*) K IF (IBETA .NE. 10) GO TO 20 C1 = 1.57E0 C2 = 7.9632679489661923132E-4 GO TO 30 20 C1 = 201.0E0/128.0E0 C2 = 4.8382679489661923132E-4 30 A = -0.125E0 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*SRAN(I1) + XL IF (J .LE. 2) GO TO 40 YSQ = HALF - HALF* ABS(X) X = (HALF - (YSQ+YSQ)) + HALF IF (J .EQ. 5) X = -X Y = SQRT(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 = ASIN(X) IF (L .LT. 0) Z = ACOS(X) C DO 60 M = 1, K SUM = YSQ*(SUM + 1.0E0/XM) XM = XM - 2.0E0 SUM = SUM*(XM/(XM+1.0E0)) 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.0E0 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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,RMS IF (J .NE. 2) GO TO 250 A = 0.75E0 B = 1.0E0 250 IF (J .NE. 4) GO TO 300 B = -A A = -1.0E0 C1 = C1 + C1 C2 = C2 + C2 L = -L 300 CONTINUE C IF (L.GT.0) CALL STIMERI(J,N,7,ARGS,MRE,RMS,IOUT) C IF (L.LT.0) CALL STIMERI(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 = SRAN(I1)*A Z = ASIN(X) + ASIN(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = SRAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - ASIN(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1035) X = BETA ** (DBLE(FLOAT(MINEXP))*0.75E0) Y = ASIN(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = 1.2E0 WRITE (IOUT,1052) X Y = ASIN(X) WRITE (IOUT,1055) Y WRITE(IOUT,1072) X Y= ACOS(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,34H TEST OF ASIN(X) VS TAYLOR SERIES //) 1005 FORMAT(34H TEST OF ACOS(X) VS TAYLOR SERIES //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(20H ASIN(X) WAS LARGER,I6,7H TIMES, / * 14X,7H AGREED,I6,11H TIMES, AND / * 10X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(20H ACOS(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(53H THE IDENTITY ASIN(-X) = - ASIN(X) WILL BE TESTED.// * ,8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY ASIN(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 ASIN WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H ASIN RETURNED THE VALUE,E15.4///) 1060 FORMAT(2E15.7/) 1072 FORMAT(39H ACOS WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(25H ACOS RETURNED THE VALUE,E15.4///) 1061 FORMAT(6X,7H ASIN(,E13.6,3H) =,E13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF ASIN/ ACOS TEST PROGRAM ---------- END C************************************************************************ C SUBROUTINE STEST9 C PROGRAM TO TEST SINH/ COSH C C C COMMON /SPRMVLS/ ARGS,ARGS2 REAL ARGS(10000),ARGS2(10000) C C THIS REFERS TO FUNCTIONS SINH AND COSH. 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 SMACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING C INFORMATION ON THE FLOATING-POINT ARITHMETIC C SYSTEM. NOTE THAT THE CALL TO SMACHAR 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 SRAN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1) C C C STANDARD FORTRAN SUBPROGRAMS REQUIRED C C ABS, ALOG, AMAX1, NINT, SQRT 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 REAL 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, COSH, SINH,MRE,RMS C REAL SRAN C C EXTERNAL COSH, SINH C CALL SPRINTITLE(' SINH (J=1,3), COSH (J=2,4) ') IOUT = 10 CALL SMACHAR(IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP, *EPS,EPSNEG,XMIN,XMAX) BETA = DBLE(FLOAT(IBETA)) ALBETA = ALOG(BETA) ALXMAX = ALOG(XMAX) AIT = DBLE(FLOAT(IT)) ZERO = 0.0E0 ONE = 1.0E0 THREE = 3.0E0 FIVE = 5.0E0 C0 = FIVE/16.0E0 + 1.152713683194269979E-2 A = ZERO B = 0.5E0 C = (AIT + ONE) * 0.35E0 IF (IBETA .EQ. 10) C = C * THREE N = 2000 XN = DBLE(FLOAT(N)) I1 = 0 I2 = 2 NIT = 2 - (NINT(ALOG(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 * SRAN(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.0E0 40 CONTINUE C ARGS(I)= X IF (J .EQ. 2) GO TO 50 W = X*XSQ*ZZ/6.0E0 ZZ = X + W Z = SINH(X) IF (IRND .NE. 0) GO TO 110 W = (X - ZZ) + W ZZ = ZZ + (W + W) GO TO 110 50 Z = COSH(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 = SINH(X) ZZ = ( SINH(Y) + SINH(W)) * C0 C WRITE(6,12121) Z,ZZ C12121 FORMAT(1X,2(Z16,1X)) GO TO 110 100 Z = COSH(X) ZZ = ( COSH(Y) + COSH(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 = ABS(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 = SQRT(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.0E0 IF (R6 .NE. ZERO) W = ALOG( ABS(R6))/ALBETA WRITE (IOUT,1021) R6,IBETA,W,X1 MRE = AMAX1(AIT+W,ZERO) WRITE (IOUT,1022) IBETA,MRE W = -999.0E0 IF (R7 .NE. ZERO) W = ALOG( ABS(R7))/ALBETA WRITE (IOUT,1023) R7,IBETA,W RMS = AMAX1(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 STIMERI(J,N,9,ARGS,MRE,RMS,IOUT) IF ((J.EQ.2).OR.(J.EQ.4)) * CALL STIMERI(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 = SRAN(I1) * A Z = SINH(X) + SINH(-X) WRITE (IOUT,1060) X, Z 320 CONTINUE C WRITE (IOUT,1031) BETAP = BETA ** IT X = SRAN(I1) / BETAP C DO 330 I = 1, 5 Z = X - SINH(X) WRITE (IOUT,1060) X, Z X = X / BETA 330 CONTINUE C WRITE (IOUT,1032) C DO 340 I = 1, 5 X = SRAN(I1) * A Z = COSH(X) - COSH(-X) WRITE (IOUT,1060) X, Z 340 CONTINUE C WRITE (IOUT,1035) X = BETA ** (DBLE(FLOAT(MINEXP))*0.75E0) Y = SINH(X) WRITE (IOUT,1061) X, Y C----------------------------------------------------------------- C TEST OF ERROR RETURNS C----------------------------------------------------------------- WRITE (IOUT,1050) X = ALXMAX + 0.125E0 WRITE (IOUT,1051) X Y = SINH(X) WRITE (IOUT,1055) Y X = BETAP WRITE (IOUT,1052) X Y = SINH(X) WRITE (IOUT,1055) Y WRITE(IOUT,1072) X Y= COSH(X) WRITE(IOUT,1075) Y C WRITE (IOUT,1100) RETURN 1000 FORMAT(1X,80('*') 1 ,/,47H TEST OF SINH(X) VS T.S. EXPANSION OF SINH(X) //) 1001 FORMAT(46H TEST OF SINH(X) VS C*( SINH(X+1)+ SINH(X-1)) //) 1005 FORMAT(47H TEST OF COSH(X) VS T.S. EXPANSION OF COSH(X) //) 1006 FORMAT(46H TEST OF COSH(X) VS C*( COSH(X+1)+ COSH(X-1)) //) 1010 FORMAT(I7,47H RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL / * 6X,1H(,E15.4,1H,,E15.4,1H)//) 1011 FORMAT(20H SINH(X) WAS LARGER,I6,7H TIMES, / * 13X,7H AGREED,I6,11H TIMES, AND / * 9X,11HWAS SMALLER,I6,7H TIMES.//) 1012 FORMAT(20H COSH(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,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 1025 FORMAT(14H SPECIAL TESTS//) 1030 FORMAT(53H THE IDENTITY SINH(-X) = - SINH(X) WILL BE TESTED.// * 8X,1HX,9X,12HF(X) + F(-X)/) 1031 FORMAT(53H THE IDENTITY SINH(X) = X , X SMALL, WILL BE TESTED.// * 8X,1HX,9X,8HX - F(X)/) 1032 FORMAT(52H THE IDENTITY COSH(-X) = COSH(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 SINH WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(39H0 SINH WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H SINH RETURNED THE VALUE,E15.4///) 1060 FORMAT(2E15.7/) 1072 FORMAT(39H COSH WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1075 FORMAT(25H COSH RETURNED THE VALUE,E15.4///) 1061 FORMAT(6X,7H SINH(,E13.6,3H) =,E13.6/) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) C ---------- LAST CARD OF SINH/ COSH TEST PROGRAM ---------- END C****************************************************************************** SUBROUTINE STEST11 C C TEST GAMMA,ALGAMA, ERF, ERFC C COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL EPS,EPSNEG,XMAX,XMIN CALL SMACHAR(IBETA, IT, IRND, NGRD, MACHEP, NEGEP, IEXP,MINEXP, * MAXEXP, EPS, EPSNEG, XMIN, XMAX) C CALL SSETMESAGES REWIND 9 CALL STEST11A CALL STEST11B CALL STEST11C CALL STEST11D RETURN END C SUBROUTINE STEST11A C GAMMA TEST PROGRAM TEXT COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. REAL ALOG10, GAMMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SRANTST C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL SPRINTITLE(' GAMMA ') 50000 FORMAT(1X,80('*'),/,' TEST OF GAMMA USING RESULTS FROM FILE.') C C DO 40 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 SRANTSTA(NIN, NOUT,JJ) C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= 56.9E0 WRITE(NOUT,1051) X Y= GAMMA(X) WRITE(NOUT,1055) Y X= 57.1E0 WRITE(NOUT,1052) X Y= GAMMA(X) WRITE(NOUT,1055) Y C X= -56.9E0 WRITE(NOUT,1051) X Y= GAMMA(X) WRITE(NOUT,1055) Y X= -57.1E0 WRITE(NOUT,1052) X Y= GAMMA(X) WRITE(NOUT,1055) Y C X= -1.0E-74 WRITE(NOUT,1051) X Y= GAMMA(X) WRITE(NOUT,1055) Y X= 1.0E-76 WRITE(NOUT,1052) X Y= GAMMA(X) WRITE(NOUT,1055) Y C X= -9.99999E-1 WRITE(NOUT,1051) X Y= GAMMA(X) WRITE(NOUT,1055) Y X= -1.000E0 WRITE(NOUT,1052) X Y= GAMMA(X) WRITE(NOUT,1055) Y C WRITE(NOUT,1100) RETURN 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(40H GAMMA WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(40H0 GAMMA WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(26H GAMMA RETURNED THE VALUE,E15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE SRANTSTA(NIN, NOUT,JJ) COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL 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 /SPRMVLS/ ARGS,ARG2 REAL ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. REAL SQRT , GAMMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SREADT 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= ALOG(16.0E0) IBETA= 16 X1= 0.0E0 AIT= 6.0E0 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.0E0 ER = 0.0E0 EM = 0.0E0 KT = 0 DO 40 J=1,NX CALL SREADT(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SSSAAF/ IFAIL = 1 Y = GAMMA(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL SERRCAL(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 /SSSAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. ABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL SPRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF ( ABS(ERR).GT. ABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL SPRNT(NOUT, X, Y, 1.0E0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = SQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0E0 IF (EM .NE. 0.0E0) W= ALOG( ABS(EM))/ALBETA WRITE (NOUT,1021) EM,IBETA,W,X1 MRE = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,MRE W = -999.0E0 IF (ER .NE. 0.0E0) W = ALOG( ABS(ER))/ALBETA WRITE (NOUT,1023) ER,IBETA,W RMS = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C J= 1 CALL STIMERI(J,NX,20,ARGS,MRE,RMS,NOUT) 1021 FORMAT(21H THE MAXIMUM ERROR OF,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PE11.2, 5X, E11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PE11.2, * 5H AND , E11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PE11.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 =, 1PE11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END SUBROUTINE SREADT(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 /SSSAAF/ C C .. SCALAR ARGUMENTS .. REAL X INTEGER NIN C .. C .. SCALARS IN COMMON .. REAL YEXP C .. COMMON /SSSAAF/ YEXP READ (NIN,99999) X, YEXP RETURN 99999 FORMAT (D20.6, D40.20) END SUBROUTINE SERRCAL(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 /SSSAAF/ C .. SCALAR ARGUMENTS .. REAL ERR, Y INTEGER IET, IFAIL C .. C .. SCALARS IN COMMON .. REAL YEXP C .. C .. LOCAL SCALARS .. REAL DERR, DY C .. COMMON /SSSAAF/ YEXP DY = Y DERR = YEXP - DY GO TO (20, 40), IET 20 ERR = DERR IFAIL = 0 RETURN 40 IF (Y.EQ.0.0E0) GO TO 60 ERR = DERR/Y IFAIL = 0 RETURN 60 ERR = 0.0E0 IFAIL = 1 RETURN END SUBROUTINE SPRNT(NOUT, X, Y, ERR) C PRINT SUBROUTINE(MACHINE DEPENDENT) C YEXP IS PASSED VIA COMMON /SSSAAF/ AND IS IN EXTRA PRECISION C PRINTING IS ON DEVICE NO. NOUT C .. SCALAR ARGUMENTS .. REAL ERR, X, Y INTEGER NOUT C .. C .. SCALARS IN COMMON .. REAL YEXP C .. COMMON /SSSAAF/ YEXP WRITE (NOUT,99999) X, Y, YEXP, ERR RETURN 99999 FORMAT (1X, 1PD30.11, D30.11, D30.12, 0PF15.2) END C****************************************************************************** SUBROUTINE STEST11B C ALGAMA TEST PROGRAM TEXT COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. REAL ALOG10, ALGAMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SRANTST C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL SPRINTITLE('ALGAMA ') 50000 FORMAT(1X,80('*'),/,' TEST OF ALGAMA 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 SRANTSTB(NIN, NOUT,J) 55 CONTINUE C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= 1.0D70 WRITE(NOUT,1051) X Y= ALGAMA(X) WRITE(NOUT,1055) Y X= 1.0D75 WRITE(NOUT,1052) X Y= ALGAMA(X) WRITE(NOUT,1055) Y C X= -56.9E0 WRITE(NOUT,1052) X Y= ALGAMA(X) WRITE(NOUT,1055) Y C C WRITE(NOUT,1100) RETURN 1050 FORMAT(22H TEST OF ERROR RETURNS//) 1051 FORMAT(40H ALGAMA WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 41H THIS SHOULD NOT TRIGGER AN ERROR MESSAGE//) 1052 FORMAT(40H0ALGAMA WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(26H ALGAMA RETURNED THE VALUE,E15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE SRANTSTB(NIN, NOUT,JJ) COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL 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 /SPRMVLS/ ARGS,ARG2 REAL ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. REAL SQRT ,ALGAMA C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SREADT 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= ALOG(16.0E0) IBETA= 16 X1= 0.0E0 AIT= 6.0E0 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.0E0 ER = 0.0E0 EM = 0.0E0 KT = 0 DO 40 J=1,NX CALL SREADT(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SSSAAF/ IFAIL = 1 Y = ALGAMA(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL SERRCAL(Y, IET, ERR, IFAIL) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SSSAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. ABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL SPRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF ( ABS(ERR).GT. ABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL SPRNT(NOUT, X, Y, 1.0E0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = SQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0E0 IF (EM .NE. 0.0E0) W= ALOG( ABS(EM))/ALBETA WRITE (NOUT,1021) EM,IBETA,W,X1 MRE = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,MRE W = -999.0E0 IF (ER .NE. 0.0E0) W = ALOG( ABS(ER))/ALBETA WRITE (NOUT,1023) ER,IBETA,W RMS = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C CALL STIMERI(JJ,NX,21,ARGS,MRE,RMS,NOUT) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PE11.2, 5X, E11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PE11.2, * 5H AND , E11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PE11.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 =, 1PE11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END C****************************************************************************** SUBROUTINE STEST11C C ERF TEST PROGRAM TEXT COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. REAL ALOG10, ERF C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT REAL S15AEF CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SRANTSTC C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL SPRINTITLE(' ERF ') 50000 FORMAT(1X,80('*'),/,' TEST OF ERF 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 SRANTSTC(NIN, NOUT) C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= -10.0E0 WRITE(NOUT,1051) X I= 0 C Y= S15AEF(X,I) Y= ERF(X) C Y= ERF(X) WRITE(NOUT,1055) Y X= 10.0E0 WRITE(NOUT,1051) X Y= ERF(X) WRITE(NOUT,1055) Y C WRITE(NOUT,1100) RETURN 1050 FORMAT(23H TEST OF EXTREME VALUES//) 1051 FORMAT(38H ERF WILL BE CALLED WITH THE ARGUMENT,E15.4) 1052 FORMAT(38H0 ERF WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(24H ERF RETURNED THE VALUE,E15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE SRANTSTC(NIN, NOUT) COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL 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 REAL S15AEF COMMON /SPRMVLS/ ARGS,ARG2 REAL ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. REAL SQRT , ERF C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SREADTC 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= ALOG(16.0E0) IBETA= 16 X1= 0.0E0 AIT= 6.0E0 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.0E0 ER = 0.0E0 EM = 0.0E0 KT = 0 DO 40 J=1,NX CALL SREADTC(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SSSAAF/ IFAIL = 1 JJ= 0 C Y = S15AEF(X,JJ) Y= ERF(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL SERRCAL(Y, IET, ERR, IFAIL) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SSSAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. ABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL SPRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF ( ABS(ERR).GT. ABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL SPRNT(NOUT, X, Y, 1.0E0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = SQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0E0 IF (EM .NE. 0.0E0) W= ALOG( ABS(EM))/ALBETA WRITE (NOUT,1021) EM,IBETA,W,X1 MRE = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,MRE W = -999.0E0 IF (ER .NE. 0.0E0) W = ALOG( ABS(ER))/ALBETA WRITE (NOUT,1023) ER,IBETA,W RMS = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C J= 1 CALL STIMERI(J,NX,22,ARGS,MRE,RMS,NOUT) 1021 FORMAT(30H THE MAXIMUM RELATIVE ERROR OF,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PE11.2, 5X, E11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PE11.2, * 5H AND , E11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PE11.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 =, 1PE11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END SUBROUTINE SREADTC(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 /SSSAAF/ C C .. SCALAR ARGUMENTS .. REAL X INTEGER NIN C .. C .. SCALARS IN COMMON .. REAL YEXP C .. COMMON /SSSAAF/ YEXP READ (NIN,99999) X, YEXP RETURN C99999 FORMAT (D20.6, D40.20) 99999 FORMAT(2(1X,D30.19)) END C****************************************************************************** SUBROUTINE STEST11D C ERFC TEST PROGRAM TEXT COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL DUMMY, ERR, REPER, X, Y INTEGER IA, IB, IFAIL, IHI, ILOW, ISKIP, ITEST, J, JFAIL,NIN, * NOUT, NUMTST C .. FUNCTION REFERENCES .. REAL ALOG10, ERFC C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT REAL S15AEF CHARACTER*72 TEXT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SRANTSTD C .. DATA NIN, NOUT /9,10/, NUMTST /4/ WRITE(NOUT,50000) CALL SPRINTITLE(' ERFC ') 50000 FORMAT(1X,80('*'),/,' TEST OF ERFC 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 SRANTSTD(J,NIN, NOUT) C C ERROR RETURN TESTS C WRITE(NOUT,1050) X= -10.0E0 WRITE(NOUT,1051) X I= 0 C Y= S15AEF(X,I) Y= ERFC(X) C Y= ERFC(X) WRITE(NOUT,1055) Y X= 20.0E0 WRITE(NOUT,1051) X Y= ERFC(X) WRITE(NOUT,1055) Y C WRITE(NOUT,1100) RETURN 1050 FORMAT(23H TEST OF EXTREME VALUES//) 1051 FORMAT(39H ERFC WILL BE CALLED WITH THE ARGUMENT,E15.4) 1052 FORMAT(39H0 ERFC WILL BE CALLED WITH THE ARGUMENT,E15.4/ * 37H THIS SHOULD TRIGGER AN ERROR MESSAGE//) 1055 FORMAT(25H ERFC RETURNED THE VALUE,E15.4//) 1100 FORMAT(25H THIS CONCLUDES THE TESTS ) RETURN END SUBROUTINE SRANTSTD(JJ,NIN, NOUT) COMMON /SMCHPRM/ EPS,EPSNEG,XMAX,XMIN * ,IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD INTEGER IBETA,IEXP,IRND,IT,MACHEP,MAXEXP,MINEXP,NEGEP,NGRD REAL 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 .. REAL 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 REAL S15AEF COMMON /SPRMVLS/ ARGS,ARG2 REAL ARGS(10000),ARG2(10000) C .. LOCAL ARRAYS .. INTEGER ETP(6), SAM(6) C .. FUNCTION REFERENCES .. REAL SQRT , ERFC C *** IMPLEMENTATION DEPENDENT DECLARATION *** C REAL DFLOAT C .. SUBROUTINE REFERENCES .. C SERRCAL, SPRNT, SREADTC 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= ALOG(16.0E0) IBETA= 16 X1= 0.0E0 AIT= 6.0E0 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.0E0 ER = 0.0E0 EM = 0.0E0 KT = 0 DO 40 J=1,NX CALL SREADTC(NIN, X) C READS ON INPUT DEVICE NO. NIN,RETURNING ARGUMENT X AND PUTTING C THE EXPECTED VALUE YEXP INTO COMMON BLOCK /SSSAAF/ IFAIL = 1 Y= ERFC(X) ARGS(J)= X IFAIL= 0 IF (IFAIL.NE.0) GO TO 20 CALL SERRCAL(Y, IET, ERR, IFAIL) IF (JJ.EQ.2) ERR= 0.5E0*ERR/(X*X) C CALCULATES THE ERROR AND RETURNS IT V1A ERR C YEXP IS OBTAINED FROM COMMON BLOCK /SSSAAF/ C IFAIL.NE.0 IF YEXP.EQ.0 AND REL.ERR.(IET=2) REQUESTED IF (J.LE.1 .OR. ABS(ERR).GE.PRTER .OR. IFAIL.NE.0) * CALL SPRNT(NOUT, X, Y, ERR/REPER) IF (IFAIL.NE.0) GO TO 40 EA = EA + ERR ER = ER + ERR*ERR IF ( ABS(ERR).GT. ABS(EM)) THEN EM = ERR X1= X ENDIF KT = KT + 1 GO TO 40 20 WRITE (NOUT,99996) IFAIL CALL SPRNT(NOUT, X, Y, 1.0E0) 40 CONTINUE EA = EA/DFLOAT(KT) ER = SQRT(ER/DFLOAT(KT)) FEA = EA/REPER FER = ER/REPER FEM = EM/REPER WRITE (NOUT,99995) EA, FEA, ER, FER, EM, FEM W = -999.0E0 IF (EM .NE. 0.0E0) W= ALOG( ABS(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 = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,MRE W = -999.0E0 IF (ER .NE. 0.0E0) W = ALOG( ABS(ER))/ALBETA IF (JJ.EQ.1) WRITE (NOUT,1023) ER,IBETA,W IF (JJ.EQ.2) WRITE (NOUT,2023) ER,IBETA,W RMS = AMAX1(AIT+W,0.0E0) WRITE (NOUT,1022) IBETA,RMS C C CPU AND SUMMARY C CALL STIMERI(JJ,NX,22,ARGS,MRE,RMS,NOUT) 1021 FORMAT(21H THE MAXIMUM ERROR OF,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.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,E15.4, * 3H = ,I4,3H **,F7.2) 2021 FORMAT(32H THE MAXIMUM SCALED-REL-ERROR OF,E15.4,3H = ,I4,3H **, * F7.2/4X,16HOCCURRED FOR X =,E17.6) 2023 FORMAT(42H THE ROOT MEAN SQUARE SACLED-REL-ERROR WAS,E15.4, * 3H = ,I4,3H **,F7.2) RETURN 99999 FORMAT (I4, 30X, 1PE11.2, 5X, E11.2/20X, I1/18X, I1/18X, I4) 99998 FORMAT (1H0, I4, 30H ARGUMENTS IN SAMPLE, BETWEEN , 1PE11.2, * 5H AND , E11.2/21H SAMPLE DISTRIBUTION , 3A1) 99997 FORMAT (15H ERROR MEASURE , 3A1, 25H, REPRESENTATION ERROR = , * 1PE11.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 = , * 1PE11.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 =, 1PE11.2, 2H (, 0PF11.2, 9H)*REP.ERR, * 1H./13H RMS ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR./ * 13H MAX ERROR =, 1PE11.2, 2H (, 0PF11.2, 10H)*REP.ERR.) END