home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C!T ROM 2
/
ctrom_ii_b.zip
/
ctrom_ii_b
/
PROGRAM
/
PASCAL
/
TPL60N19
/
TESTPRGS
/
DSQRT.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1993-02-14
|
8KB
|
266 lines
PROGRAM DSqrt; { 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 DSQRT
C
C DATA REQUIRED
C
C NONE
C
C SUBPROGRAMS REQUIRED FROM THIS PACKAGE
C
C MACHAR - AN ENVIRONMENTAL INQUIRY PROGRAM PROVIDING
C INFORMATION ON THE FLOATING-POINT ARITHMETIC
C SYSTEM. NOTE THAT THE CALL TO MACHAR CAN
C BE DELETED PROVIDED THE FOLLOWING SIX
C PARAMETERS ARE ASSIGNED THE VALUES INDICATED
C
C IBETA - THE RADIX OF THE FLOATING-POINT SYSTEM
C IT - THE NUMBER OF BASE-IBETA DIGITS IN THE
C SIGNIFICAND OF A FLOATING-POINT NUMBER
C EPS - THE SMALLEST POSITIVE FLOATING-POINT
C NUMBER SUCH THAT 1.0+EPS .NE. 1.0
C EPSNEG - THE SMALLEST POSITIVE FLOATING-POINT
C NUMBER SUCH THAT 1.0-EPSNEG .NE. 1.0
C XMIN - THE SMALLEST NON-VANISHING FLOATING-POINT
C POWER OF THE RADIX
C XMAX - THE LARGEST FINITE FLOATING-POINT NO.
C
C RANDL(X) - A FUNCTION SUBPROGRAM RETURNING LOGARITHMICALLY
C DISTRIBUTED RANDOM REAL NUMBERS. IN PARTICULAR,
C A * RANDL(DLOG(B/A))
C IS LOGARITHMICALLY DISTRIBUTED OVER (A,B)
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, DMAX1, DFLOAT, DSQRT
C
C
C LATEST REVISION - AUGUST 2, 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;
FUNCTION RANDL(X: REAL): REAL;
{
C
C RETURNS PSEUDO RANDOM NUMBERS LOGARITHMICALLY DISTRIBUTED
C OVER (1,EXP(X)). THUS A*RANDL(LN(B/A)) IS LOGARITHMICALLY
C DISTRIBUTED IN (A,B).
C
C OTHER SUBROUTINES REQUIRED
C
C EXP(X) - THE EXPONENTIAL ROUTINE
C
C REN(K) - A FUNCTION PROGRAM RETURNING RANDOM REAL
C NUMBERS UNIFORMLY DISTRIBUTED OVER (0,1).
C THE ARGUMENT K IS A DUMMY.
C
C
}
CONST K:LONGINT=1;
BEGIN
RANDL := EXP (X*REN(K));
END;
VAR I,IBETA,IEXP,IOUT,IRND,IT,J,K1,K2,
K3,MACHEP,MAXEXP,MINEXP,N,NEGEP,NGRD: LONGINT;
A,AIT,ALBETA,B,BETA,C,EPS,EPSNEG,ONE,
R6,R7,SQBETA,W,X,XMAX,XMIN,XN,X1,Y,Z,ZERO: REAL;
LABEL 100, 110, 120, 150, 160, 210, 220, 230, 240, 300;
BEGIN
N := 1000000; { 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;
SQBETA:= SQRT (BETA);
ALBETA:= LN (BETA);
AIT := (IT);
ONE := 1;
ZERO := 0;
A := ONE / SQBETA;
B := ONE;
XN := N;
{-----------------------------------------------------------------}
{ RANDOM ARGUMENT ACCURACY TESTS }
{-----------------------------------------------------------------}
FOR J := 1 TO 2 DO BEGIN
C := LN (B/A);
K1 := 0;
K3 := 0;
X1 := ZERO;
R6 := ZERO;
R7 := ZERO;
FOR I := 1 TO N DO BEGIN
X := A * RANDL(C);
Y := X * X;
Z := SQRT (Y);
IF X <> ZERO THEN
W := (Z - X) / X
ELSE IF Z <> 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;
END;
K2 := N - K1 - K3;
R7 := SQRT (R7/XN);
WRITELN;
WRITELN;
WRITELN ('TEST OF SQRT(X*X) - X ');
WRITELN;
WRITELN (N, ' RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL');
WRITELN ('(', A, ',', B, ')');
WRITELN;
WRITELN ('SQRT (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.0;
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;
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 := ONE;
B := SQBETA;
END;
{-----------------------------------------------------------------}
{ SPECIAL TESTS }
{-----------------------------------------------------------------}
WRITELN;
WRITELN;
WRITELN ('TEST OF SPECIAL ARGUMENTS');
WRITELN;
X := XMIN;
Y := SQRT (X);
WRITELN ('SQRT (XMIN) = SQRT ( ', X:18, ') = ', Y:18);
WRITELN;
X := ONE - EPSNEG;
Y := SQRT(X);
WRITELN ('SQRT(1-EPSNEG) = SQRT (1-', EPSNEG:18, ') = ', Y:18);
WRITELN;
X := ONE;
Y := SQRT(X);
WRITELN ('SQRT (1.0) = SQRT ( ', X:18, ') = ', Y:18);
WRITELN;
X := ONE + EPS;
Y := SQRT(X);
WRITELN ('SQRT (1+EPS) = SQRT (1+', EPS:18, ') = ', Y:18);
WRITELN;
X := XMAX;
Y := SQRT(X);
WRITELN ('SQRT (XMAX) = SQRT ( ', X:18, ') = ', Y:18);
WRITELN;
{-----------------------------------------------------------------}
{ TEST OF ERROR RETURNS }
{-----------------------------------------------------------------}
WRITELN;
WRITELN;
WRITELN ('TEST OF ERROR RETURNS');
WRITELN;
X := ZERO;
WRITELN ('SQRT WILL BE CALLED WITH THE ARGUMENT ', X:15);
WRITELN ('THIS SHOULD NOT TRIGGER AN ERROR MESSAGE');
Y := SQRT(X);
WRITELN ('SQRT RETURNED THE VALUE ', Y:15);
X := -ONE;
WRITELN ('SQRT WILL BE CALLED WITH THE ARGUMENT ', X:15);
WRITELN ('THIS SHOULD TRIGGER AN ERROR MESSAGE');
Y := SQRT(X);
WRITELN ('SQRT RETURNED THE VALUE ', Y:15);
WRITELN;
WRITELN ('THIS CONCLUDES THE TESTS');
END. { DSqrt }