home *** CD-ROM | disk | FTP | other *** search
- PROGRAM DLog; { ported from Fortran original 05-01-92 Norbert Juffa }
-
- {$A+,B-,D-,E+,F-,G-,I-,L-,N-,O-,R-,S-,V-,X-}
-
- USES MachArit;
-
- {
- C PROGRAM TO TEST DLOG
- C
- C DATA REQUIRED
- C
- C NONE
- C
- C SUBPROGRAMS REQUIRED FROM THIS PACKAGE
- C
- C MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
- C INFORMATION ON THE FLOATING-POINT ARITHMETIC
- C SYSTEM. NOTE THAT THE CALL TO MACHAR CAN
- C BE DELETED PROVIDED THE FOLLOWING FOUR
- C PARAMETERS ARE ASSIGNED THE VALUES INDICATED
- C
- C IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM
- C IT - THE NUMBER OF BASE-IBETA DIGITS IN THE
- C SIGNIFICAND OF A FLOATING-POINT NUMBER
- C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT
- C POWER OF THE RADIX
- C XMAX - THE LARGEST FINITE FLOATING-POINT NO.
- C
- C REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
- C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
- C
- C
- C STANDARD FORTRAN SUBPROGRAMS REQUIRED
- C
- C DABS, DLOG, DLOG10, DMAX1, DFLOAT, DSIGN, DSQRT
- C
- C
- C LATEST REVISION - DECEMBER 6, 1979
- C
- C AUTHOR - W. J. CODY
- C ARGONNE NATIONAL LABORATORY
- C
- C
- }
-
-
-
- FUNCTION REN (K: LONGINT): REAL;
-
- {
- DOUBLE PRECISION FUNCTION REN(K)
- C
- C RANDOM NUMBER GENERATOR - BASED ON ALGORITHM 266 BY PIKE AND
- C HILL (MODIFIED BY HANSSON), COMMUNICATIONS OF THE ACM,
- C VOL. 8, NO. 10, OCTOBER 1965.
- C
- C THIS SUBPROGRAM IS INTENDED FOR USE ON COMPUTERS WITH
- C FIXED POINT WORDLENGTH OF AT LEAST 29 BITS. IT IS
- C BEST IF THE FLOATING POINT SIGNIFICAND HAS AT MOST
- C 29 BITS.
- C
- }
-
- VAR J: LONGINT;
- CONST IY: LONGINT = 100001;
-
- BEGIN
- J := K;
- IY := IY * 125;
- IY := IY - (IY DIV 2796203) * 2796203;
- REN:= 1.0 * (IY) / 2796203.0e0 * (1.0e0 + 1.0e-6 + 1.0e-12);
- END;
-
-
- FUNCTION LOG (X: REAL): REAL;
- BEGIN
- LOG := LN (X) * 0.43429448190325182765112891891660508;
- END;
-
-
- FUNCTION MAX1 (A, B:REAL): REAL;
- BEGIN
- IF A > B THEN
- MAX1 := A
- ELSE
- MAX1 := B;
- END;
-
-
-
- VAR I,IBETA,IEXP,IOUT,IRND,IT,I1,J,K1,K2,
- K3,MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
-
- A,AIT,ALBETA,B,BETA,C,DEL,EIGHT,EPS,
- EPSNEG,HALF,ONE,T, R6,R7,TENTH,W,X,
- XL,XMAX,XMIN,XN,X1,Y,Z,ZERO,ZZ,FOUR,
- TWO,THREE,NINETENTH,FIFTEEN,SIXTEEN,
- TWENTYONE,THIRTYONE,TWOHUNDREDFORTY,
- FIVEHUNDREDTWELVE: REAL;
-
- LABEL 100, 110, 120, 150, 160, 220, 230, 240, 300;
-
- BEGIN
-
- N := 10000; { number of trials }
-
- MACHAR (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
- EPS,EPSNEG,XMIN,XMAX);
- PRINTPARAM (IBETA,IT,IRND,NGRD,MACHEP,NEGEP,IEXP,MINEXP,MAXEXP,
- EPS,EPSNEG,XMIN,XMAX);
- BETA := IBETA;
- ALBETA := LN (BETA);
- AIT := IT;
- J := IT DIV 3;
- ZERO := 0;
- ONE := 1;
- TWO := 2;
- THREE := 3;
- FOUR := 4;
- EIGHT := 8;
- FIFTEEN := 15;
- SIXTEEN := 16;
- THIRTYONE := 31;
- TWENTYONE := 21;
- TENTH := 0.1;
- HALF := 0.5;
- NINETENTH := 0.9;
- TWOHUNDREDFORTY := 240;
- FIVEHUNDREDTWELVE:= 512;
- C := ONE;
-
- FOR I := 1 TO J DO BEGIN
- C := C / BETA;
- END;
-
- B := ONE + C;
- A := ONE - C;
- XN := N;
- I1 := 0;
-
- {-----------------------------------------------------------------}
- { RANDOM ARGUMENT ACCURACY TESTS }
- {-----------------------------------------------------------------}
-
- FOR J := 1 TO 4 DO BEGIN
- K1 := 0;
- K3 := 0;
- X1 := ZERO;
- R6 := ZERO;
- R7 := ZERO;
- DEL:= (B - A) / XN;
- XL := A;
-
- FOR I := 1 TO N DO BEGIN
- X := DEL * REN(I1) + XL;
- IF J <> 1 THEN
- GOTO 100;
- Y := X - HALF;
- Y := Y - HALF;
- ZZ:= LN (X);
- Z := (Y * (ONE / THREE - Y / FOUR) - HALF) * Y * Y + Y;
- GOTO 150;
- 100: IF J <> 2 THEN
- GOTO 110;
- X := X + EIGHT;
- X := X - EIGHT;
- Y := X / SIXTEEN;
- Y := X + Y;
- Z := LN (X);
- ZZ:= LN (Y);
- ZZ:= ZZ - 7.7746816434842581e-5; { Ln (17/16) - 31/512) }
- ZZ:= ZZ - THIRTYONE/FIVEHUNDREDTWELVE;
- GOTO 150;
- 110: IF J <> 3 THEN
- GOTO 120;
- X := X + EIGHT;
- X := X - EIGHT;
- T := X * TENTH;
- Y := X + T;
- Z := LOG (X);
- ZZ:= LOG (Y);
- ZZ:= ZZ - 3.7706015822504075e-4; { Log10 (11/10) - 21/512) }
- ZZ:= ZZ - TWENTYONE/FIVEHUNDREDTWELVE;
- GOTO 150;
- 120: T := X * X;
- Z := LN (T);
- ZZ:= LN (X);
- ZZ:= ZZ + ZZ;
-
- 150: IF Z <> ZERO THEN
- W := (Z - ZZ) / Z
- ELSE IF ZZ <> ZERO THEN
- W := ONE;
- IF W > ZERO THEN
- K1 := K1 + 1;
- IF W < ZERO THEN
- K3 := K3 + 1;
- W := ABS (W);
- IF W <= R6 THEN
- GOTO 160;
- R6 := W;
- X1 := X;
- 160: R7 := R7 + W * W;
- XL := XL + DEL;
- END;
-
- K2 := N - K3 - K1;
- R7 := SQRT (R7/XN);
-
- IF J = 1 THEN BEGIN
- WRITELN;
- WRITELN ;
- WRITELN ('TEST OF LN (X) VS T.S. EXPANSION OF LN(1+Y)');
- WRITELN;
- END;
- IF J = 2 THEN BEGIN
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF LN(X) VS LN(17X/16)-LN(17/16)');
- WRITELN;
- END;
- IF J = 3 THEN BEGIN
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF LOG10(X) VS LOG10(11X/10)-LOG10(11/10)');
- WRITELN;
- END;
- IF J = 4 THEN BEGIN
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF LN (X*X) VS 2*LN(X)');
- WRITELN;
- END;
- IF J = 1 THEN BEGIN
- WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
- WRITELN ('(1-EPS,1+EPS), WHERE EPS = ', C);
- WRITELN;
- END;
- IF J <> 1 THEN BEGIN
- WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
- WRITELN ('(', A, ',', B, ')');
- WRITELN;
- END;
- IF J <> 3 THEN BEGIN
- WRITELN ('LN (X) WAS LARGER', K1:6, ' TIMES');
- WRITELN (' AGREED', K2:6, ' TIMES');
- WRITELN (' AND WAS SMALLER', K3:6, ' TIMES');
- END;
- IF J = 3 THEN BEGIN
- WRITELN ('LOG (X) WAS LARGER', K1:6, ' TIMES');
- WRITELN (' AGREED', K2:6, ' TIMES');
- WRITELN (' AND WAS SMALLER', K3:6, ' TIMES');
- END;
- WRITELN;
- WRITELN ('THERE ARE ', IT, ' BASE ', IBETA,
- ' SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER');
- WRITELN;
- W := -999;
- IF R6 <> ZERO THEN
- W := LN (ABS(R6))/ALBETA;
- WRITELN ('THE MAXIMUM RELATIVE ERROR OF ', R6:12,
- ' = ', IBETA, ' **', W:7:2);
- WRITELN ('OCCURED FOR X = ', X1);
- W := MAX1 (AIT+W,ZERO);
- WRITELN;
- WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
- ' SIGNIFICANT DIGITS IS ', W:7:2);
- W := -999.0;
- IF R7 <> ZERO THEN
- W := LN (ABS(R7))/ALBETA;
- WRITELN;
- WRITELN ('THE ROOT MEAN SQUARE RELATIVE ERROR WAS', R7:12,
- ' = ', IBETA, ' **' , W:7:2);
- W := MAX1 (AIT+W,ZERO);
- WRITELN;
- WRITELN ('THE ESTIMATED LOSS OF BASE ', IBETA,
- ' SIGNIFICANT DIGITS IS ', W:7:2);
- IF J > 1 THEN
- GOTO 230;
- A := SQRT (HALF);
- B := FIFTEEN / SIXTEEN;
- GOTO 300;
- 230: IF J > 2 THEN
- GOTO 240;
- A := SQRT (TENTH);
- B := NINETENTH;
- GOTO 300;
- 240: A := SIXTEEN;
- B := TWOHUNDREDFORTY;
- 300:
- END;
-
- {-----------------------------------------------------------------}
- { SPECIAL TESTS }
- {-----------------------------------------------------------------}
-
- WRITELN;
- WRITELN;
- WRITELN ('SPECIAL TESTS');
- WRITELN;
- WRITELN ('THE IDENTITY LN (X) = - LN (1/X) WILL BE TESTED');
- WRITELN;
- WRITELN (' X F(X) + F(1/X)');
- WRITELN;
-
- FOR I := 1 TO 5 DO BEGIN
- X := REN(I1);
- T := X + X;
- X := T + FIFTEEN;
- Y := ONE / X;
- T := LN (X);
- Z := LN (Y);
- Z := Z + T;
- WRITELN (X:18, Z:18);
- END;
-
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF SPECIAL ARGUMENTS');
- WRITELN;
- X := ONE;
- Y := LN (X);
- WRITELN ('LN (1.0) = ', Y:15);
- X := XMIN;
- Y := LN (X);
- WRITELN ('LN (XMIN)= LN (', X:10, ') = ', Y:15);
- X := XMAX;
- Y := LN (X);
- WRITELN ('LN (XMAX)= LN (', X:10, ') = ', Y:15);
-
- {-----------------------------------------------------------------}
- { TEST OF ERROR RETURNS }
- {-----------------------------------------------------------------}
-
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF ERROR RETURNS');
- WRITELN;
- X := -TWO;
- WRITELN ('LN WILL BE CALLED WITH THE ARGUMENT ', X:15);
- WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE');
- Y := LN (X);
- WRITELN ('LN RETURNED THE VALUE ', Y:15);
- X := ZERO;
- WRITELN ('LN WILL BE CALLED WITH THE ARGUMENT ', X:15);
- WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE');
- Y := LN (X);
- WRITELN ('LN RETURNED THE VALUE ', Y:15);
- WRITELN;
- WRITELN ('THIS CONCLUDES THE TESTS');
- END. { DLog }