home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C!T ROM 2
/
ctrom_ii_b.zip
/
ctrom_ii_b
/
PROGRAM
/
PASCAL
/
TPL60N19
/
TESTPRGS
/
MAINVARS.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1992-02-28
|
18KB
|
581 lines
{$a+,n-,x-,s-,i-,r-,b-,v-}
unit mainvars;
interface
{ (C) Apr 19 1983 in BASIC version by:
Professor W M Kahan,
567 Evans Hall.
Electrical Engineering & Computer Science Dept.
University of California
Berkeley, California 94720
USA
converted to Pascal by:
B A Wichmann
National Physical Laboratory
Teddington Middx
TW11 OLW
UK
further massaging by dmg =
David M. Gay
AT&T Bell Labs
600 Mountain Avenue
Murray Hill, NJ 07974
and a couple of bug fixes from dgh = sun!dhough (29 May 1986)
See the article by Richard Karpinski in the February 1985 issue
of BYTE Magazine.
You may copy this program freely if you acknowledge its source.
Comments on the Pascal version to NPL or dmg, please. }
const
{integer constants}
NoTrials = 20;
{Number of tests for commutativity. }
type
Guard = (Yes, No);
Rounding = (Chopped, Rounded, Other);
Message = packed array [1..40] of char;
WhichOp = packed array [1..14] of char;
Class = (Flaw, Defect, SeriousDefect, Failure);
var
{input: text;}
{Small floating point constants.}
Zero, { 0.0; }
Half, { 0.5; }
One, { 1.0; }
Two, { 2.0; }
Three, { 3.0; }
Four, { 4.0; }
Five, { 5.0; }
Eight, { 8.0; }
Nine, { 9.0; }
TwentySeven, { 27.0; }
ThirtyTwo, { 32.0; }
TwoForty, { 240.0; }
MinusOne, { -1.0; }
OneAndHalf: { 1.5; } real;
MyZero: integer;
NoTimes, Index: integer;
ch: char;
AInverse, A1: real;
Radix, BInverse, RadixD2, BMinusU2: real;
C, CInverse: real;
D, FourD: real;
E0, E1, Exp2, MinSqrtError: real;
SqrtError, MaxSqrtError, E9: real;
Third: real;
F6, F9: real;
H, HInverse: real;
I: integer;
StickyBit, J: real;
M, N, N1: real;
Precision: real;
Q, Q9: real;
R, R9: real;
T, Underflow, S: real;
OneUlp, UnderflowThreshold, U1, U2: real;
V, V0, V9: real;
W: real;
X, X1, X2, X8, RandomNumber1: real;
Y, Y1, Y2, RandomNumber2: real;
Z, PseudoZero, Z1, Z2, Z9: real;
NoErrors: array [Class] of integer;
Milestone: integer;
PageNo: integer;
GMult, GDiv, GAddSub: Guard;
RMult, RDiv, RAddSub, RSqrt: Rounding;
Continue, Break, Done, NotMonot, Monot, AnomolousArithmetic, IEEE,
SquareRootWrong, UnderflowNotGradual: Boolean;
{ Computed constants. }
{U1 gap below 1.0, i.e, 1.0-U1 is next number below 1.0 }
{U2 gap above 1.0, i.e, 1.0+U2 is next number above 1.0 }
procedure Page;
function Int (X: real): real;
function Sign (X: real): real;
procedure Pause;
procedure Instructions;
procedure Heading;
procedure Characteristics;
procedure History;
procedure notify(T: WhichOp);
procedure TestCondition (K: Class; Valid: Boolean; T: Message);
function Random: real;
procedure SqrtXMinX (ErrorKind: Class);
procedure NewD;
procedure SubRout3750;
function Power (X, Y: real): real;
procedure DoesYequalX;
procedure SubRout3980;
procedure PrintIfNPositive;
procedure TestPartialUnderflow;
implementation
procedure Page;
begin
(* write(#$C) *) {FF in TURBO Pascal} writeln; writeln;
end;
function Int (X: real): real;
{ simulates BASIC INT-function, which is defined as:
INT(X) is the greatest integer value less than or
equal to X. }
function LargeTrunc (X: real): real;
var
start, acc, y, p: real;
trunced: integer; (* dgh *)
begin (* LargeTrunc *)
if abs (X) < maxint then begin
trunced := trunc(X);
LargeTrunc := trunced;
end
else
begin
start := abs (X);
acc := 0.0;
repeat
y := start;
p := 1.0;
while y > maxint - 1.0 do
begin
y := y / Radix;
p := p * Radix;
end;
trunced := trunc(y);
acc := acc + trunced * p;
start := start - trunced * p;
until start < 1.0;
if X < 0.0 then
LargeTrunc := - acc
else
LargeTrunc := acc
end;
end (* LargeTrunc *);
begin (* Int *)
if X > 0.0 then
Int := LargeTrunc (X)
else if LargeTrunc (X - 0.5) = X then
Int := X
else
Int := LargeTrunc (X) - 1;
end (* Int *);
function Sign (X: real): real;
begin (* Sign *)
if X < 0.0 then
Sign := - 1.0
else
Sign := + 1.0;
end (* Sign *);
procedure Pause;
var
ch: char;
begin (* Pause *)
writeln ('To continue, press any key and newline:');
readln (input);
while not eoln (input) do
read (input, ch);
Page;
write ('Diagnosis resumes after milestone no ', Milestone);
writeln (' Page: ', PageNo);
writeln;
Milestone := Milestone + 1;
PageNo := PageNo + 1;
end (* Pause *);
procedure Instructions;
begin (* Instructions *)
writeln ('Lest this program stop prematurely, ',
'i.e. before displaying');
writeln (' "END OF TEST",');
writeln ('try to persuade the computer NOT to',
' terminate execution whenever an');
writeln ('error like Over/Underflow or Division by Zero occurs,',
' but rather');
writeln ('to persevere with a surrogate value after, ',
' perhaps, displaying some');
writeln ('warning. If persuasion avails naught, don''t despair'
, ' but run this');
writeln ('program anyway to see how many milestones it passes,',
' and then');
writeln ('amend it to make further progress.');
writeln ('Answer questions with Y, y, N or n',
' (unless otherwise indicated).');
writeln;
end (* Instructions *);
procedure Heading;
begin (* Heading *)
writeln ('Users are invited to help debug and augment',
' this program so it will');
writeln ('cope with unanticipated and newly uncovered',
' arithmetic pathologies.');
writeln ('Please send suggestions and interesting results to');
writeln(' Richard Karpinski');
writeln(' Computer Center U-76');
writeln(' University of California');
writeln(' San Francisco, CA 94143-0704, USA');
writeln;
writeln('In doing so, please include the following information:');
writeln(' Version: 10 February 1989');
writeln(' Computer:'); writeln;
writeln(' Compiler:'); writeln;
writeln(' Optimization level:'); writeln;
writeln(' Other relevant compiler options:'); writeln;
end (* Heading *);
procedure Characteristics;
begin (* Characteristics *)
writeln (
'Running this program should reveal these characteristics');
writeln (' Radix = 1, 2, 4, 8, 10, 16, 100, 256, or ...');
writeln (' Precision = number of significant digits carried.');
writeln (' U2 = Radix/Radix^Precision = One Ulp (OneUlpnit in the');
writeln (' Last Place) of 1.000xxx .');
writeln (' U1 = 1/Radix^Precision = One Ulp of numbers',
' a little less than 1.0 .');
writeln (' Adequacy of guard digits for Mult., Div., and Subt.');
writeln (' Whether arithmetic is chopped, correctly rounded, ',
'or something else');
writeln (' for Mult., Div., Add/Subt. and Sqrt.');
writeln (' Whether a Sticky Bit is used correctly for rounding.');
writeln (' UnderflowThreshold = an Underflow Threshold.');
writeln (' E0 and PseudoZero tell whether underflow is abrupt,',
' gradual, or fuzzy.');
writeln (' V = an overflow threshold, roughly.');
writeln (' V0 tells, roughly, whether Infinity is represented.');
writeln (' Comparisions are checked for consistency with',
' subtraction');
writeln (' and for contamination with pseudo-zeros.');
writeln (' Sqrt is tested. Y^X is not tested.');
writeln (' Extra-precise subexpressions are revealed',
' but NOT YET tested.');
writeln (' Decimal-Binary conversion is NOT YET tested',
' for accuracy.');
end (* Characteristics *);
procedure History;
begin (* History *)
writeln ('The program attempts to discriminate among');
writeln (' FLAWs, like lack of a sticky bit,');
writeln (' SERIOUS DEFECTs, like lack of a guard digit, and');
writeln (' FAILUREs, like 2+2 = 5 .');
writeln ('Failures may confound subsequent diagnoses.');
writeln;
writeln ('The diagnostic capabilities of this program go beyond',
' an earlier');
writeln ('program called "MACHAR", which can be found at the',
' end of the');
writeln ('book "Software Manual for the Elementary Functions',
'" (1980) by');
writeln ('W. J. Cody and W. Waite. Although both programs',
' try to discover');
writeln ('the Radix, Precision and range (over/underflow',
' thresholds)');
writeln ('of the arithmetic, this program tries to cope',
' with a wider variety');
writeln ('of pathologies, and to say how well the',
' arithmetic is implemented.');
writeln ('BASIC version of this program (C) 1983 by Professor',
' W. M. Kahan;');
writeln ('see source comments for more history.');
writeln;
writeln ('The program is based upon a conventional',
' radix representation for');
writeln ('floating-point numbers, but also allows',
' logarithmic encoding');
writeln ('as used by certain early WANG machines.');
writeln;
end (* History *);
procedure notify(T: WhichOp);
begin
writeln('The ', T, ' test is inconsistent;');
writeln(' PLEASE NOTIFY KARPINSKI !');
end;
procedure TestCondition (K: Class; Valid: Boolean; T: Message);
begin (* TestCondition *)
if not Valid then
begin
NoErrors [K] := NoErrors [K] + 1;
case K of
Flaw:
write ('FLAW');
Defect:
write ('DEFECT');
SeriousDefect:
write ('SERIOUS DEFECT');
Failure:
write ('FAILURE');
end;
writeln (': ', T);
end;
end (* TestCondition *);
function Random: real;
var
X, Y: real;
begin (* Random *)
X := RandomNumber1 + R9;
Y := X * X;
Y := Y * Y;
X := X * Y;
Y := X - Int (X);
RandomNumber1 := Y + X * 0.000005;
Random := RandomNumber1;
end (* Random *);
procedure SqrtXMinX (ErrorKind: Class);
begin (* SqrtXMinX *)
SqrtError := ((sqrt (X * X) - X * BInverse) - (X - X * BInverse))
/ OneUlp;
if SqrtError <> 0 then
begin
if SqrtError < MinSqrtError then
MinSqrtError := SqrtError;
if SqrtError > MaxSqrtError then
MaxSqrtError := SqrtError;
J := J + 1;
writeln;
if ErrorKind = SeriousDefect then
write ('SERIOUS ');
writeln ('DEFECT: sqrt( ', X * X, ' - ', X, ') = ');
writeln (OneUlp * SqrtError, ' instead of correct value 0 .');
end;
end (* SqrtXMinX *);
procedure NewD;
begin (* NewD *)
X := Z1 * Q;
X := Int (Half - X / Radix) * Radix + X;
Q := (Q - X * Z) / Radix + X * X * (D / Radix);
Z := Z - Two * X * D;
if Z <= Zero then
begin
Z := - Z;
Z1 := - Z1;
end;
D := Radix * D;
end (* NewD *);
procedure SubRout3750;
begin (* SubRout3750 *)
if not ((X - Radix < Z2 - Radix) or (X - Z2 > W - Z2)) then
begin
I := I + 1;
X2 := sqrt (X * D);
Y2 := (X2 - Z2) - (Y - Z2);
X2 := X8 / (Y - Half);
X2 := X2 - Half * X2 * X2;
SqrtError := (Y2 + Half) + (Half - X2);
if SqrtError < MinSqrtError then
MinSqrtError := SqrtError;
SqrtError := Y2 - X2;
if SqrtError > MaxSqrtError then
MaxSqrtError := SqrtError;
end;
end (* SubRout3750 *);
{ Sun version of Power (from dgh)...
function pow(x, y : real ): real; external c;
function Power(X, Y: real): real;
var xx, yy: real;
begin
xx := X;
yy := Y;
Power := pow(xx,yy);
end;
}
function Power (X, Y: real): real;
var z:real;
iy:integer;
reciproc:boolean;
(* Crude power function: see Cody & Waite for a better one. *)
begin (* Power *)
if Y = 0.0 then
Power := 1.0
else if (X = 0.0) and (Y > 0.0) then
Power := 0.0
else if (Y = Int (Y)) AND (Abs(Y) <= MAXINT) THEN BEGIN
IY := TRUNC (Y);
reciproc := iy < 0;
iy := abs (iy);
z := 1.0;
while iy > 0 do begin
if odd(iy) then begin
z := z*x;
end else begin
end; {endif}
iy := iy div 2;
if iy >0 then begin
x := sqr (x);
end else begin
end; {endif}
end; {endwhile}
if reciproc then begin
power := 1 / z;
end else begin
power := z;
end {endif}
end
else if (X < 0.0) and (Y = Int(Y)) then
if odd(trunc(Y)) then Power := -exp (Y * ln (-X))
else Power := exp (Y * ln (-X))
else
Power := exp (Y * ln (X));
end (* Power *);
procedure DoesYequalX;
begin (* DoesYequalX *)
if Y <> X then
begin
if N <= 0 then
begin
if (Z = Zero) and (Q <= Zero) then write('WARNING') (* dgh: Y --> Z *)
else begin
NoErrors [Defect] := NoErrors [Defect] + 1;
write('DEFECT');
end;
writeln (': computed (', Z, ') ^ (', Q, ') = ');
writeln (Y, ', which compares unequal to correct ', X, ';'); (* dgh: V --> Y *)
writeln (' they differ by ', Y - X);
end;
N := N + 1;
{ ... count discrepancies. }
end;
end (* DoesYequalX *);
procedure SubRout3980;
begin (* SubRout3980 *)
repeat
Q := I;
Y := Power (Z, Q);
DoesYequalX;
I := I + 1;
if I <= M then
X := Z * X;
until (X >= W) or (I > M);
end (* SubRout3980 *);
procedure PrintIfNPositive;
begin (* PrintIfNPositive *)
if N > 0 then
writeln ('Similar discrepancies have occurred ', N, ' times.');
end (* PrintIfNPositive *);
procedure TestPartialUnderflow;
begin (* TestPartialUnderflow *)
N := 0;
if Z <> 0 then
begin
writeln ('Since comparison denies Z = 0, evaluating');
writeln ('(Z + Z) / Z should be safe.');
write ('What the machine gets for (Z + Z) / Z is: ');
Q9 := (Z + Z) / Z;
writeln (Q9);
if (abs (Q9 - Two) < Radix * U2) then
begin
write ('This is O.K., provided Over/Underflow');
writeln (' has NOT just been signaled.');
end
else if (Q9 < One) or (Q9 > Two) then
begin
N := 1;
NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
writeln ('This is a VERY SERIOUS DEFECT!');
end
else
begin
N := 1;
NoErrors [Defect] := NoErrors [Defect] + 1;
writeln ('This is a DEFECT.');
end;
V9 := Z * One;
RandomNumber1 := V9;
V9 := One * Z;
RandomNumber2 := V9;
V9 := Z / One;
if (Z = RandomNumber1) and (Z = RandomNumber2)
and (Z = V9) then
begin
if N > 0 then
Pause
end
else
begin
N := 1;
NoErrors [Defect] := NoErrors [Defect] + 1;
writeln ('DEFECT: What prints as Z = ', Z, 'compares');
write ('different from ');
if not (Z = RandomNumber1) then
writeln ('Z * 1 = ', RandomNumber1);
if not ((Z = RandomNumber2)
or (RandomNumber2 = RandomNumber1)) then
writeln ('1 * Z = ', RandomNumber2);
if not (Z = V9) then
writeln ('Z / 1 = ', V9);
if RandomNumber2 <> RandomNumber1 then
begin
NoErrors [Defect] := NoErrors [Defect] + 1;
writeln ('DEFECT Multiplication does not commute;');
writeln ('comparison alleges that 1 * Z = ',
RandomNumber2);
writeln ('differs from Z * 1 = ', RandomNumber1);
end;
if N > 0 then
Pause;
end;
end;
end (* TestPartialUnderflow *);
end.