home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
C!T ROM 2
/
ctrom_ii_b.zip
/
ctrom_ii_b
/
PROGRAM
/
PASCAL
/
TPL60N19
/
TESTPRGS
/
UNIT2.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1992-02-28
|
25KB
|
822 lines
{$a+,n-,x-,s-,i-,r-,b-,v-}
unit Unit2;
interface
uses mainvars;
procedure mile70170;
implementation
procedure mile70170;
begin
{=============================================}
Milestone := 70;
{=============================================}
writeln;
writeln ('Running test of square root(x).');
TestCondition (Failure, (Zero = sqrt (Zero))
and (- Zero = sqrt (- Zero))
and (One = sqrt (One)), ' Square root of 0.0, -0.0 or 1.0 wrong '
);
MinSqrtError := Zero;
MaxSqrtError := Zero;
J := 0;
X := Radix;
OneUlp := U2;
SqrtXMinX (SeriousDefect);
X := BInverse;
OneUlp := BInverse * U1;
SqrtXMinX (SeriousDefect);
X := U1;
OneUlp := U1 * U1;
SqrtXMinX (SeriousDefect);
if J <> 0 then
begin
NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
Pause;
end;
writeln ('Testing if sqrt(X * X) = X for ', NoTrials, ' integers X.');
J := 0;
X := Two;
Y := Radix;
if (Radix <> One) then
repeat
X := Y;
Y := Radix * Y;
until (Y - X >= NoTrials);
OneUlp := X * U2;
I := 1;
Continue := true;
while (I <= NoTrials) and Continue do (* dgh: 10 --> NoTrials *)
begin
X := X + One;
SqrtXMinX (Defect);
if J > 0 then
begin
Continue := false;
NoErrors [Defect] := NoErrors [Defect] + 1;
end;
I := I + 1;
end;
writeln ('Test for Sqrt Monotonicity.');
I := - 1;
X := BMinusU2;
Y := Radix;
Z := Radix + Radix * U2;
NotMonot := false;
Monot := false;
while not (NotMonot or Monot) do
begin
I := I + 1;
X := sqrt (X);
Q := sqrt (Y);
Z := sqrt (Z);
if (X > Q) or (Q > Z) then
NotMonot := true
else
begin
Q := Int (Q + Half);
if (I > 0) or (Radix = Q * Q) then
Monot := true
else if I > 0 then
begin
if I > 1 then
Monot := true
else
begin
Y := Y * BInverse;
X := Y - U1;
Z := Y + U1;
end
end
else
begin
Y := Q;
X := Y - U2;
Z := Y + U2;
end
end
end;
if Monot then
writeln ('Sqrt has passed a test for Monotonicity.')
else
begin
NoErrors [Defect] := NoErrors [Defect] + 1;
writeln('DEFECT: Sqrt(X) is non-monotonic for X near ', Y);
end;
{=============================================}
Milestone := 80;
{=============================================}
MinSqrtError := MinSqrtError + Half;
MaxSqrtError := MaxSqrtError - Half;
Y := (sqrt (One + U2) - One) / U2;
SqrtError := (Y - One) + U2 / Eight;
if SqrtError > MaxSqrtError then
MaxSqrtError := SqrtError;
SqrtError := Y + U2 / Eight;
if SqrtError < MinSqrtError then
MinSqrtError := SqrtError;
Y := ((sqrt (F9) - U2) - (One - U2)) / U1;
SqrtError := Y + U1 / Eight;
if SqrtError > MaxSqrtError then
MaxSqrtError := SqrtError;
SqrtError := (Y + One) + U1 / Eight;
if SqrtError < MinSqrtError then
MinSqrtError := SqrtError;
OneUlp := U2;
X := OneUlp;
for Index := 1 to 3 do
begin
Y := sqrt ((X + U1 + X) + F9);
Y := ((Y - U2) - ((One - U2) + X)) / OneUlp;
Z := ((U1 - X) + F9) * Half * X * X / OneUlp;
SqrtError := (Y + Half) + Z;
if SqrtError < MinSqrtError then
MinSqrtError := SqrtError;
SqrtError := (Y - Half) + Z;
if SqrtError > MaxSqrtError then
MaxSqrtError := SqrtError;
if ((Index = 1) or (Index = 3)) then
X := OneUlp * Sign (X) * Int (Eight / (Nine * sqrt (OneUlp)))
else
begin
OneUlp := U1;
X := - OneUlp;
end;
end;
{=============================================}
Milestone := 85;
{=============================================}
SquareRootWrong := false;
AnomolousArithmetic := false;
RSqrt := Other; (* ~dgh *)
if Radix <> One then
begin
writeln ('Testing whether sqrt is rounded or chopped: ');
D := Int (Half + Power (Radix, One + Precision - Int (Precision)))
;
{ ... = Radix^(1 + fract) if Precision = integer + fract. }
X := D / Radix;
Y := D / A1;
if (X <> Int (X)) or (Y <> Int (Y)) then
begin
AnomolousArithmetic := true;
end
else
begin
X := Zero;
Z2 := X;
Y := One;
Y2 := Y;
Z1 := Radix - One;
FourD := Four * D;
repeat
if Y2 > Z2 then
begin
Q := Radix;
Y1 := Y;
repeat
X1 := abs (Q + Int (Half - Q / Y1) * Y1);
Q := Y1;
Y1 := X1;
until X1 <= Zero;
if Q <= One then
begin
Z2 := Y2;
Z := Y;
end;
end;
Y := Y + Two;
X := X + Eight;
Y2 := Y2 + X;
if Y2 >= FourD then
Y2 := Y2 - FourD;
until Y >= D;
X8 := FourD - Z2;
Q := (X8 + Z * Z) / FourD;
X8 := X8 / Eight;
if Q <> Int (Q) then
AnomolousArithmetic := true
else
begin
Break := false;
repeat
X := Z1 * Z;
X := X - Int (X / Radix) * Radix;
if X = One then
Break := true
else
Z1 := Z1 - One;
until Break or (Z1 <= 0);
if (Z1 <= 0) and (not Break) then
AnomolousArithmetic := true
else
begin
if Z1 > RadixD2 then
Z1 := Z1 - Radix;
repeat
NewD;
until U2 * D >= F9;
if D * Radix - D <> W - D then
AnomolousArithmetic := true
else
begin
Z2 := D;
I := 0;
Y := D + (One + Z) * Half;
X := D + Z + Q;
SubRout3750;
Y := D + (One - Z) * Half + D;
X := D - Z + D;
X := X + Q + X;
SubRout3750;
NewD;
if D - Z2 <> W - Z2 then
AnomolousArithmetic := true
else
begin
Y := (D - Z2) + (Z2 + (One - Z) * Half);
X := (D - Z2) + (Z2 - Z + Q);
SubRout3750;
Y := (One + Z) * Half;
X := Q;
SubRout3750;
if I = 0 then
AnomolousArithmetic := true;
end
end
end
end
end;
if (I = 0) or AnomolousArithmetic then
begin
NoErrors [Failure] := NoErrors [Failure] + 1;
writeln ('FAILURE: Anomolous arithmetic with ',
'integer < Radix^Precision = ');
writeln (W, ' fails test whether sqrt rounds or chops.');
SquareRootWrong := true;
end
end;
if not AnomolousArithmetic then
begin
if not ((MinSqrtError < 0) or (MaxSqrtError > 0)) then
begin
RSqrt := Rounded;
writeln ('Square root appears to be correctly rounded.');
end
else if (MaxSqrtError + U2 > U2 - Half) or (MinSqrtError > Half)
or (MinSqrtError + Radix < Half) then
SquareRootWrong := true
else
begin
RSqrt := Chopped;
writeln ('Square root appears to be chopped.');
end;
end;
if SquareRootWrong then
begin
writeln ('Square root is neither chopped nor correctly rounded.');
writeln ('Observed errors run from ', MinSqrtError - Half);
writeln ('to ', Half + MaxSqrtError, ' ulps.');
TestCondition (SeriousDefect, MaxSqrtError - MinSqrtError < Radix
* Radix, 'sqrt gets too many last digits wrong. ');
end;
{=============================================}
Milestone := 90;
{=============================================}
Pause;
(* fix from dgh: Wichman had effectively commented out here to Milestone 110 *)
writeln ('Testing powers Z^i for small integers Z and i.');
N := 0;
{ ... test power of zero. }
I := 0;
Z := - Zero;
M := 3;
Break := false;
repeat
X := One;
SubRout3980;
if I <= 10 then
begin
I := 1023;
SubRout3980;
end;
if Z = MinusOne then
Break := true
else
begin
Z := MinusOne;
{ .. if(-1)^N is invalid, replace MinusOne by One. }
I := - 4;
end;
until Break;
PrintIfNPositive;
N1 := N;
N := 0;
Z := A1;
M := Int (Two * ln (W) / ln (A1));
Break := false;
repeat
X := Z;
I := 1;
SubRout3980;
if Z = AInverse then
Break := true
else
Z := AInverse;
until Break;
{=============================================}
Milestone := 100;
{=============================================}
{ Power of Radix have been tested, }
{ next try a few primes }
M := NoTrials;
Z := Three;
repeat
X := Z;
I := 1;
SubRout3980;
repeat
Z := Z + Two;
until (Three * Int (Z / Three) <> Z);
until (Z >= Eight * Three);
if N > 0 then
begin
writeln ('Error like this may invalidate financial ');
writeln ('calculations involving interest rates.');
end;
PrintIfNPositive;
N := N + N1;
if N = 0 then
writeln ('... no discrepancies found.');
writeln;
if N > 0 then
Pause;
{=============================================}
Milestone := 110;
{=============================================}
writeln ('Seeking Underflow thresholds UnderflowThreshold and E0');
D := U1;
if (Precision <> Int (Precision)) then
begin
D := BInverse;
X := Precision;
repeat
D := D * BInverse;
X := X - One;
until X <= Zero;
end;
Y := One;
Z := D;
{ ... D is power of 1/Radix < 1. }
repeat
C := Y;
Y := Z;
Z := Y * Y;
until not ((Y > Z) and (Z + Z > Z));
Y := C;
Z := Y * D;
repeat
C := Y;
Y := Z;
Z := Y * D;
until not ((Y > Z) and (Z + Z > Z));
if Radix < Two then
HInverse := Two
else
HInverse := Radix;
H := One / HInverse;
{ ... 1/HInverse = H = Min(1/Radix, 1/2) }
CInverse := One / C;
E0 := C;
Z := E0 * H;
{ ...1/Radix^(BIG integer) << 1 << CInverse = 1/C }
repeat
Y := E0;
E0 := Z;
Z := E0 * H;
until not ((E0 > Z) and (Z + Z > Z));
UnderflowThreshold := E0;
E1 := Zero;
Q := Zero;
E9 := U2;
S := One + E9;
D := C * S;
if D <= C then
begin
E9 := Radix * U2;
S := One + E9;
D := C * S;
if D <= C then
begin
writeln ('FAILURE: multiplication gets too many last digits wrong.');
NoErrors [Failure] := NoErrors [Failure] + 1;
Underflow := E0;
Y1 := Zero;
PseudoZero := Z;
Pause;
end
end
else
begin
Underflow := D;
PseudoZero := Underflow * H;
UnderflowThreshold := Zero;
repeat
Y1 := Underflow;
Underflow := PseudoZero;
if E1 + E1 <= E1 then
begin
Y2 := Underflow * HInverse;
E1 := abs (Y1 - Y2);
Q := Y1;
if (UnderflowThreshold = Zero) and (Y1 <> Y2) then
UnderflowThreshold := Y1;
end;
PseudoZero := PseudoZero * H;
until not ((Underflow > PseudoZero)
and (PseudoZero + PseudoZero > PseudoZero));
end;
{ Comment line 4530 .. 4560 }
if PseudoZero <> Zero then
begin
writeln;
Z := PseudoZero;
{ ... Test PseudoZero for "phoney-zero" violating }
{ ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
... }
if PseudoZero <= 0 then
begin
NoErrors [Failure] := NoErrors [Failure] + 1;
writeln ('FAILURE: Positive expressions can underflow to an ');
writeln ('allegedly negative value PseudoZero that prints out as:');
writeln (PseudoZero);
X := - PseudoZero;
if X <= 0 then
begin
writeln ('But -PseudoZero, which should then be positive, isn''t;');
writeln (' it prints out as: ', X);
end
end
else
begin
NoErrors [Flaw] := NoErrors [Flaw] + 1;
writeln ('FLAW: Underflow can stick at an allegedly positive');
writeln ('value PseudoZero that prints out as: ', PseudoZero);
end;
TestPartialUnderflow;
end;
{=============================================}
Milestone := 120;
{=============================================}
if (CInverse * Y > CInverse * Y1) then
begin
S := H * S;
E0 := Underflow;
end;
if not ((E1 = 0) or (E1 = E0)) then
begin
NoErrors [Defect] := NoErrors [Defect] + 1;
write ('DEFECT: ');
if E1 < E0 then
begin
writeln ('Products underflow at a higher');
writeln ('threshold than differences.');
if PseudoZero = Zero then
E0 := E1;
end
else
begin
writeln ('Difference underflows at a higher');
writeln ('threshold than products.');
end;
end;
writeln ('Smallest strictly positive number found is E0 =', E0);
Z := E0;
TestPartialUnderflow;
Underflow := E0;
if N = 1 then
Underflow := Y;
I := 4;
if E1 = Zero then
I := 3;
if UnderflowThreshold = Zero then
I := I - 2;
UnderflowNotGradual := true;
case I of
1:
begin
UnderflowThreshold := Underflow;
if (CInverse * Q) <> ((CInverse * Y) * S) then
begin
NoErrors [Failure] := NoErrors [Failure] + 1;
UnderflowThreshold := Y;
writeln ('FAILURE: Either accuracy deteriorates as numbers');
writeln ('approach a threshold UnderflowThreshold = ',
UnderflowThreshold);
writeln ('coming down from ', C, ' or else multiplication');
writeln ('gets too many last digits wrong.');
end;
Pause;
end;
2:
begin
NoErrors [Failure] := NoErrors [Failure] + 1;
writeln ('FAILURE: Underflow confuses Comparison, which alleges that');
writeln ('Q = Y while denying that |Q - Y| = 0 ; these values');
write ('print out as Q = ', Q, ' Y = ', Y2, ' |Q - Y| = ');
writeln (abs (Q - Y2));
UnderflowThreshold := Q;
end;
3:
begin
X := X;
{FOR PRETTYPRINTER, IS DUMMY}
end;
4:
if (Q = UnderflowThreshold) and (E1 = E0)
and (abs ( UnderflowThreshold - E1 / E9) <= E1) then
begin
UnderflowNotGradual := false;
writeln ('Underflow is gradual; it incurs Absolute Error =')
;
writeln ('(roundoff in UnderflowThreshold) < E0.');
Y := E0 * CInverse;
Y := Y * (OneAndHalf + U2);
X := CInverse * (One + U2);
Y := Y / X;
IEEE := (Y = E0);
end;
end;
if UnderflowNotGradual then
begin
writeln;
R := sqrt (Underflow / UnderflowThreshold);
if R <= H then
begin
Z := R * UnderflowThreshold;
X := Z * (One + R * H * (One + H));
end
else
begin
Z := UnderflowThreshold;
X := Z * (One + H * H * (One + H));
end;
if not ((X = Z) or (X - Z <> 0)) then
begin
NoErrors [Flaw] := NoErrors [Flaw] + 1;
writeln ('FLAW: X = ', X, ' is unequal to Z = ', Z);
write ('yet X - Z yields ');
Z9 := X - Z;
writeln (Z9, '. Should this NOT signal Underflow,');
writeln ('this is a SERIOUS DEFECT that causes');
writeln ('confusion when innocent statements like');
write (' if X = Z then ... else');
writeln (' ... (f(X) - f(Z)) / (X - Z) ...');
writeln ('encounter Division by Zero although actually');
writeln ('X / Z = 1 + ', (X / Z - Half) - Half);
end;
end;
writeln ('The Underflow threshold is ', UnderflowThreshold,
' below which');
writeln ('calculation may suffer larger Relative error then',
' merely roundoff.');
Y2 := U1 * U1;
Y := Y2 * Y2;
Y2 := Y * U1;
if Y2 <= UnderflowThreshold then
begin
if Y > E0 then
begin
NoErrors [Defect] := NoErrors [Defect] + 1;
I := 5;
end
else
begin
NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
write ('SERIOUS ');
I := 4;
end;
writeln ('DEFECT: Range is too narrow; U1^', I, ' Underflows.');
end;
{=============================================}
Milestone := 130;
{=============================================}
{ SKIP
Y := - Int (Half - TwoForty * ln (UnderflowThreshold) / ln (HInverse)
) / TwoForty;
Y2 := Y + Y;
writeln ('Since underflow occurs below the threshold');
writeln ('UnderflowThreshold = ( ', HInverse, ' ) ^ ( ', Y, ') only u
nderflow');
writeln ('should afflict the expression ( ', HInverse, ' ) ^ ( ', Y2,
' )');
write ('actually calculating yields: ');
V9 := Power (HInverse, Y2);
writeln (V9);
if not ((V9 >= 0) and (V9 <= (Radix + Radix + E9) * UnderflowThreshol
d)) then
begin
NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
writeln ('SERIOUS DEFECT: this is not between 0 and');
writeln (' UnderflowThreshold = ', UnderflowThreshold);
end
else if not (V9 > UnderflowThreshold * (One + E9)) then
writeln ('this computed value is O.K.')
else
begin
NoErrors [Defect] := NoErrors [Defect] + 1;
writeln ('DEFECT: This is not between 0 and UnderflowThreshold = ',
UnderflowThreshold);
end;
end SKIP }
{=============================================}
Milestone := 140;
{=============================================}
writeln;
{ ...calculate Exp2 = exp(2) = 7.389056099... }
X := 0;
I := 2;
Y := Two * Three;
Q := 0;
N := 0;
repeat
Z := X;
I := I + 1;
Y := Y / (I + I);
R := Y + Q;
X := Z + R;
Q := (Z - X) + R;
until X <= Z;
Z := (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
X := Z * Z;
Exp2 := X * X;
X := F9;
Y := X - U1;
write ('Testing X^((X + 1) / (X - 1)) vs. exp(2) = ');
writeln (Exp2, ' as X -> 1.');
Break := false;
I := 1;
while (not Break) do
begin
Z := X - BInverse;
Z := (X + One) / (Z - (One - BInverse));
Q := Power (X, Z) - Exp2;
if abs (Q) > TwoForty * U2 then
begin
Break := true;
N := 1;
NoErrors [Defect] := NoErrors [Defect] + 1;
writeln ('DEFECT: Calculated ', Power(X,Z), ' for');
writeln (' (1 + (', (X - BInverse) - (One - BInverse),
')) ^ (', Z, ')');
writeln (' which differs from the correct value by Q =', Q);
writeln (' This much error may spoil financial',
' calculations involving');
writeln (' tiny interest rates.');
end
else
begin
Z := (Y - X) * Two + Y;
X := Y;
Y := Z;
Z := One + (X - F9) * (X - F9);
if ((Z > One) and (I < NoTrials)) then
I := I + 1
else if X > One then
begin
if N = 0 then
writeln ('Accuracy seems adequate.');
Break := true;
end
else
begin
X := One + U2;
Y := U2 + U2 + X;
I := 1;
end
end
end;
{=============================================}
Milestone := 150;
{=============================================}
writeln ('Testing powers Z^Q at four nearly extreme values.');
N := 0;
Z := A1;
Q := Int (Half - ln (C) / ln (A1));
Break := false;
repeat
X := CInverse;
Y := Power (Z, Q);
DoesYequalX;
Q := - Q;
X := C;
Y := Power (Z, Q);
DoesYequalX;
if Z < One then
Break := true
else
Z := AInverse;
until Break;
PrintIfNPositive;
if N = 0 then
writeln (' ... no discrepancies found.');
writeln;
if N > 0 then
Pause;
{=============================================}
Milestone := 160;
{=============================================}
Pause;
writeln ('Searching for Overflow threshold: ',
'This may generate an error.');
writeln ('Try a few values for N, starting with a large one,');
writeln ('and take the one that just does not stop the machine.');
writeln ('Did you find the correct value for N yet?');
writeln ('Hint: the number for Turbo-Pascal is 45');
readln (input);
while not eoln (input) do
read (input, ch);
Break := false;
repeat
writeln ('N = ');
readln (input);
while not eoln (input) do
read (input, NoTimes);
Y := - CInverse;
V9 := HInverse * Y;
Index := 1;
I := 0;
repeat
V := Y;
Y := V9;
V9 := HInverse * Y;
if V9 >= Y then begin I := 1; ch := 'Y'; Index := NoTimes end;
Index := Index + 1;
until Index > NoTimes;
if I = 0 then Y := V9;
if (ch = 'N') or (ch = 'n') then
writeln ('N seems not large enough, try again.')
else
begin
writeln ('O.K.');
Break := true;
end;
until Break;
Z := V9;
writeln ('Can "Z = -Y" overflow? Trying it on Y = ', Y);
V9 := - Y;
V0 := V9;
if (V - Y = V + V0) then
writeln ('Seems O.K.')
else
begin
NoErrors [Flaw] := NoErrors [Flaw] + 1;
writeln ('Finds a FLAW: -(-Y) differs from Y.');
end;
if Z <> Y then
begin
NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
writeln ('SERIOUS DEFECT:');
writeln (' overflow past ', Y, ' shrinks to ', Z);
end;
if (I = 1) then
begin
Y := V * (HInverse * U2 - HInverse);
Z := Y + ((One - HInverse) * U2) * V;
if Z < V0 then
Y := Z;
if Y < V0 then
V := Y;
if V0 - V < V0 then
V := V0;
end
else begin
V := Y * (HInverse * U2 - HInverse);
V := V + ((One - HInverse) * U2) * Y;
end;
writeln ('Overflow threshold is V = ', V);
if (I = 1) then writeln ('Overflow saturates at V0 = ', V0)
else begin writeln('There is no saturation value because');
writeln('the system traps on overflow.') end;
write ('No Overflow should get signaled for V * 1 = ');
V9 := V * One;
writeln (V9);
write (' nor for V / 1 = ');
V9 := V / One;
writeln (V9);
writeln ('Any overflow signal separating this * from the one');
writeln ('above is a DEFECT.');
{=============================================}
Milestone := 170;
{=============================================}
TestCondition (Failure, (- V < V) and (- V0 < V0)
and (- UnderflowThreshold < V)
and (UnderflowThreshold < V),
'Comparisons are confused by Overflow. ');
end;
end.