home *** CD-ROM | disk | FTP | other *** search
- PROGRAM DATan; { 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, Power;
-
- {
- C PROGRAM TO TEST DATAN, DATAN2
- C
- C DATA REQUIRED
- C
- C NONE
- C
- C SUBPROGRAMS REQUIRED FROM THIS PACKAGE
- C
- C MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
- C INFORMATION ON THE FLOATING-POINT ARITHMETIC
- C SYSTEM. NOTE THAT THE CALL TO MACHAR CAN
- C BE DELETED PROVIDED THE FOLLOWING SIX
- C PARAMETERS ARE ASSIGNED THE VALUES INDICATED
- C
- C IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM
- C IT - THE NUMBER OF BASE-IBETA DIGITS IN THE
- C SIGNIFICAND OF A FLOATING-POINT NUMBER
- C IRND - 0 IF FLOATING-POINT ADDITION CHOPS,
- C 1 IF FLOATING-POINT ADDITION ROUNDS
- C MINEXP - THE LARGEST IN MAGNITUDE NEGATIVE
- C INTEGER SUCH THAT DFLOAT(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 REN(K) - A FUNCTION SUBPROGRAM RETURNING RANDOM REAL
- C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1)
- C
- C STANDARD FORTRAN SUBPROGRAMS REQUIRED
- C
- C DABS, DLOG, DMAX1, DATAN, DATAN2, DFLOAT, 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 MAX1 (A, B:REAL): REAL;
- BEGIN
- IF A > B THEN
- MAX1 := A
- ELSE
- MAX1 := B;
- END;
-
-
-
- VAR I,IBETA,IEXP,IOUT,IRND,II,IT,I1,J,K1,
- K2,K3,MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
-
- 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,T,Z,
- THREE,FIVE,SIXTEEN,SEVENTEEN,ZERO,ZZ,
- THREEFORTH,SIXTEENTH: REAL;
-
- LABEL 100,105,110,120,200;
-
- 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;
- ZERO := 0;
- ONE := 1;
- TWO := 2;
- THREE := 3;
- FIVE := 5;
- SIXTEEN := 16;
- SEVENTEEN := 17;
- SIXTEENTH := 0.0625;
- HALF := 0.5;
- THREEFORTH:= 0.75;
- A := -SIXTEENTH;
- B := -A;
- OB32 := B * HALF;
- 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 = 2 THEN
- X := ((ONE+X*A)-ONE)*SIXTEEN;
- Z := ARCTAN (X);
- IF J <> 1 THEN
- GOTO 100;
- XSQ := X * X;
- EM := SEVENTEEN;
- SUM := XSQ / EM;
- FOR II := 1 TO 7 DO BEGIN
- EM := EM - TWO;
- SUM:= (ONE/EM - SUM) * XSQ;
- END;
- SUM := -X * SUM;
- ZZ := X + SUM;
- SUM := (X - ZZ) + SUM;
- IF IRND = 0 THEN
- ZZ := ZZ + (SUM + SUM);
- GOTO 110;
- 100: IF J <> 2 THEN
- GOTO 105;
- Y := (X - SIXTEENTH) / (ONE + X * A);
- ZZ:= ARCTAN (Y) - 8.1190004042651526021e-5; { arctan (1/16) - 1/16 }
- ZZ:= ZZ + OB32;
- ZZ:= ZZ + OB32;
- GOTO 110;
- 105: Z := Z + Z;
- Y := X / ((HALF + X * HALF)*((HALF - X) + HALF));
- ZZ:= ARCTAN (Y);
- 110: 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 120;
- R6 := W;
- X1 := X;
- 120: R7 := R7 + W * W;
- XL := XL + DEL;
- 200: END;
-
- K2 := N - K3 - K1;
- R7 := SQRT(R7/XN);
-
- IF J = 1 THEN BEGIN
- WRITELN;
- WRITELN ;
- WRITELN ('TEST OF ARCTAN(X) VS TRUNCATED TAYLOR SERIES');
- WRITELN;
- END;
- IF J = 2 THEN BEGIN
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF ARCTAN(X) VS ARCTAN(1/16) + ARCTAN((X-1/16)/(1+X/16))');
- WRITELN;
- END;
- IF J > 2 THEN BEGIN
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF 2*ARCTAN(X) VS ARCTAN(2X((1-X*X))');
- WRITELN;
- END;
-
- WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
- WRITELN ('(', A, ',', B, ')');
- WRITELN;
-
- WRITELN ('ARCTAN (X) WAS LARGER', K1:6, ' TIMES');
- WRITELN (' AGREED', K2:6, ' TIMES');
- WRITELN (' AND WAS SMALLER', K3:6, ' TIMES');
-
- 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);
- A := B;
- IF J = 1 THEN
- B := TWO - SQRT (THREE);
- IF J = 2 THEN
- B := SQRT (TWO) - ONE;
- IF J = 3 THEN
- B := ONE;
- END;
-
- {-----------------------------------------------------------------}
- { SPECIAL TESTS }
- {-----------------------------------------------------------------}
-
- WRITELN;
- WRITELN;
- WRITELN ('SPECIAL TESTS');
- WRITELN;
- WRITELN ('THE IDENTITY ARCTAN(-X) = -ARCTAN(X) WILL BE TESTED');
- WRITELN;
- WRITELN (' X F(X) + F(-X)');
- WRITELN;
-
- A := FIVE;
-
- FOR I := 1 TO 5 DO BEGIN
- X := REN(I1) * A;
- Z := ARCTAN(X) + ARCTAN(-X);
- WRITELN (X:18, Z:18);
- END;
-
- WRITELN;
- WRITELN;
- WRITELN ('THE IDENTITY ARCTAN(X) = X , X SMALL, WILL BE TESTED.');
- WRITELN;
- WRITELN (' X X-F(X)');
- WRITELN;
- BETAP := POW (BETA,IT);
- X := REN(I1) / BETAP;
-
- FOR I := 1 TO 5 DO BEGIN
- Z := X - ARCTAN (X);
- WRITELN (X:18, Z:18);
- X := X / BETA;
- END;
-
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT');
- WRITELN;
- EXPON := MINEXP * THREEFORTH;
- X := POW (BETA,EXPON);
- Y := ARCTAN(X);
- WRITELN ('ARCTAN (', X:18, ') = ', Y:18);
-
- {-----------------------------------------------------------------}
- { TEST OF ERROR RETURNS }
- {-----------------------------------------------------------------}
-
- WRITELN;
- WRITELN;
- WRITELN ('TEST OF ERROR RETURNS');
- WRITELN;
- WRITELN ('ARCTAN WILL BE CALLED WITH THE ARGUMENT ', XMAX:18);
- WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE');
- Z := ARCTAN(XMAX);
- WRITELN ('ARCTAN (', XMAX:18, ') = ', Z:18);
- WRITELN;
- WRITELN ('THIS CONCLUDES THE TESTS');
- END. { DATan }