home *** CD-ROM | disk | FTP | other *** search
- {$a+,n-,x-,s-,i-,r-,b-,v-}
-
- unit Unit1;
- interface
- uses mainvars;
- procedure start;
- procedure mile2060;
-
- implementation
- procedure start;
-
-
- begin (* PARA *)
-
- {First two assignments use integer right-hand sides.}
- Zero := 0;
- One := 1;
- Two := One + One;
- Three := Two + One;
- Four := Three + One;
- Five := Four + One;
- Eight := Four + Four;
- Nine := Three * Three;
- TwentySeven := Nine * Three;
- ThirtyTwo := Four * Eight;
- TwoForty := Four * Five * Three * Four;
- MinusOne := -One;
- Half := One / Two;
- OneAndHalf := One + Half;
-
- NoErrors [Failure] := 0;
- NoErrors [SeriousDefect] := 0;
- NoErrors [Defect] := 0;
- NoErrors [Flaw] := 0;
- PageNo := 0;
- {=============================================}
- Milestone := 0;
- {=============================================}
- writeln ('Type any character to start the program.');
- { assign(input,'con:');} { for TURBO Pascal version 2 }
- { reset (input); } { for old Cray Pascal }
- while not eoln (input) do
- read (input, ch);
- Instructions;
- Pause;
- Heading;
- Pause;
- Characteristics;
- Pause;
- History;
- {=============================================}
- Milestone := 7;
- {=============================================}
- Pause;
- writeln ('Program is now RUNNING tests on small integers:');
- TestCondition (Failure, (Zero + Zero = Zero) and (One - One = Zero)
- and (One > Zero)
- and (One + One = Two), ' 0+0<>0 or 1-1<>0 or 1<=0 or 1+1<>2 '
- );
- Z := - Zero;
- if Z <> 0.0 then
- begin
- NoErrors [Failure] := NoErrors [Failure] + 1;
- writeln ('Comparison alleges that -0.0 is Non-zero!');
- U2 := 0.001;
- Radix := 1;
- TestPartialUnderflow;
- end;
- TestCondition (Failure, (Three = Two + One) and (Four = Three + One)
- and (Four + Two * (- Two) = Zero)
- and (Four - Three - One = Zero),
- ' 3<>2+1, 4<>3+1, 4+2*(-2)<>0 or 4-3-1<>0');
- TestCondition (Failure, (MinusOne = - One)
- and (MinusOne + One = Zero ) and (One + MinusOne = Zero)
- and (MinusOne + abs (One) = Zero)
- and (MinusOne + MinusOne * MinusOne = Zero),
- '-1+1<>0, -1+abs(1)<>0 or -1+(-1)*(-1)<>0');
- TestCondition (Failure, Half + MinusOne + Half = Zero,
- ' 1/2 + (-1) + 1/2 <> 0 ');
- {=============================================}
- Milestone := 10;
- {=============================================}
- TestCondition (Failure, (Nine = Three * Three)
- and (TwentySeven = Nine * Three) and (Eight = Four + Four)
- and (ThirtyTwo = Eight * Four)
- and (ThirtyTwo - TwentySeven - Four - One = Zero),
- '9<>3*3, 27<>9*3, 32<>8*4 or 32-27-4-1<>0');
- TestCondition (Failure, (Five = Four + One)
- and (TwoForty = Four * Five * Three * Four)
- and (TwoForty / Three - Four * Four * Five = Zero)
- and ( TwoForty / Four - Five * Three * Four = Zero)
- and ( TwoForty / Five - Four * Three * Four = Zero),
- '5<>4+1,240/3<>80,240/4<>60, or 240/5<>48');
- if NoErrors [Failure] = 0 then
- begin
- writeln (' -1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.');
- writeln
- end;
- writeln ('Searching for Radix and Precision.');
- W := One;
- repeat
- W := W + W;
- Y := W + One;
- Z := Y - W;
- Y := Z - One;
- until (MinusOne + abs (Y) >= Zero);
- {.. now W is just big enough that |((W+1)-W)-1| >= 1 ...}
- Precision := 0;
- Y := One;
- repeat
- Radix := W + Y;
- Y := Y + Y;
- Radix := Radix - W;
- until (Radix <> Zero);
- if Radix < Two then
- Radix := One;
- writeln ('Radix = ', Radix);
- if Radix <> 1 then
- begin
- W := One;
- repeat
- Precision := Precision + One;
- W := W * Radix;
- Y := W + One;
- until (Y - W) <> One;
- {... now W = Radix^Precision is barely too big to satisfy (W+1)-W = 1
- ...}
- end;
- U1 := One / W;
- U2 := Radix * U1;
- writeln ('Closest relative separation found is U1 = ', U1);
- writeln;
- writeln ('Recalculating radix and precision');
- E0 := Radix;
- E1 := U1;
- E9 := U2;
- {save old values}
- X := Four / Three;
- Third := X - One;
- F6 := Half - Third;
- X := F6 + F6;
- X := abs (X - Third);
- if X < U2 then
- X := U2;
- {... now X = (unknown no.) ulps of 1+...}
- repeat
- U2 := X;
- Y := Half * U2 + ThirtyTwo * U2 * U2;
- Y := One + Y;
- X := Y - One;
- until (U2 <= X) or (X <= Zero);
- {... now U2 = 1 ulp of 1 + ... }
- X := Two / Three;
- F6 := X - Half;
- Third := F6 + F6;
- X := Third - Half;
- X := abs (X + F6);
- if X < U1 then
- X := U1;
- {... now X = (unknown no.) ulps of 1 -... }
- repeat
- U1 := X;
- Y := Half * U1 + ThirtyTwo * U1 * U1;
- Y := Half - Y;
- X := Half + Y;
- Y := Half - X;
- X := Half + Y;
- until (U1 <= X) or (X <= Zero);
- {... now U1 = 1 ulp of 1 - ... }
- if U1 = E1 then
- writeln (' confirms closest relative separation U1 .')
- else
- writeln (' gets better closest relative separation U1 = ', U1);
- W := One / U1;
- F9 := (Half - U1) + Half;
- Radix := Int (0.01 + U2 / U1);
- if Radix = E0 then
- writeln ('Radix confirmed.')
- else
- writeln ('MYSTERY: recalculated Radix = ', Radix);
- TestCondition (Defect, Radix <= Eight + Eight,
- 'Radix is too big: roundoff problems ');
- TestCondition (Flaw, (Radix = Two) or (Radix = 10)
- or (Radix = One), 'Radix is not as good as 2 or 10. ');
- end (*start*);
-
- procedure mile2060;
- begin
-
- {=============================================}
- Milestone := 20;
- {=============================================}
- TestCondition (Failure, F9 - Half < Half,
- ' (1-U1)-1/2 < 1/2 is FALSE, prog. fails?');
- X := F9;
- I := 1;
- Y := X - Half;
- Z := Y - Half;
- TestCondition (Failure, (X <> One)
- or (Z = Zero), 'Comparison is fuzzy,X=1 but X-1/2-1/2<>1');
- X := One + U2;
- I := 0;
- {=============================================}
- Milestone := 25;
- {=============================================}
- BMinusU2 := Radix - One;
- BMinusU2 := (BMinusU2 - U2) + One;
- if Radix <> One then
- begin {... BMinusU2 = nextafter(Radix, 0) }
- X := - TwoForty * ln (U1) / ln (Radix);
- Y := Int (Half + X);
- if abs (X - Y) * Four < One then
- X := Y;
- Precision := X / TwoForty;
- Y := Int (Half + Precision);
- if abs (Precision - Y) * TwoForty < Half then
- Precision := Y;
- { Purify integers }
- end;
- if (Precision <> Int (Precision)) or (Radix = One) then
- begin
- writeln ('Precision cannot be characterized by an integer',
- ' number of sig. digits,');
- writeln ('but, by itself, this is a minor flaw.');
- end;
- if Radix = One then
- writeln ('logarithmic encoding has precision characterized',
- 'solely by U1.')
- else
- writeln ('The number of significant digits of the Radix is ',
- Precision);
- TestCondition (SeriousDefect, U2 * Nine * Nine * TwoForty < One,
- ' Precision worse than 5 decimal figures ');
- {=============================================}
- Milestone := 30;
- {=============================================}
- { Test for extra-precise subepressions }
- X := abs (((Four / Three - One) - One / Four) * Three - One / Four);
- repeat
- Z2 := X;
- X := (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
- until (Z2 <= X) or (X <= Zero);
- Y := abs ((Three / Four - Two / Three) * Three - One / Four);
- Z := Y;
- X := Y;
- repeat
- Z1 := Z;
- Z := (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
- + One / Two)) + One / Two;
- until (Z1 <= Z) or (Z <= Zero);
- repeat
- repeat
- Y1 := Y;
- Y := (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
- )) + Half;
- until (Y1 <= Y) or (Y <= Zero);
- X1 := X;
- X := ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
- until (X1 <= X) or (X <= Zero);
- if (X1 <> Y1) or (X1 <> Z1) then
- begin
- NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
- writeln ('SERIOUS DEFECT: Disagreements among the values X1, Y1, Z1');
- writeln ('resp. ', X1, Y1, Z1);
- writeln ('are symptoms of inconsistencies introduced');
- writeln ('by extra-precise evaluation of allegedly');
- writeln ('"optimized" arithmetic subexpressions.');
- writeln ('Possibly some part of this test is inconsistent.');
- if (X1 = U1) or (Y1 = U1) or (Z1 = U1) then
- writeln ('That feature is not tested further by this program.');
- end
- else if (Z1 <> U1) or (Z2 <> U2) then
- begin
- if (Z1 >= U1) or (Z2 >= U2) then
- begin
- NoErrors [Failure] := NoErrors [Failure] + 1;
- writeln ('FAILURE: Precision ', Precision);
- writeln ('U1 = ', U1, ' Z1 - U1 = ', Z1 - U1);
- writeln ('U2 = ', U2, ' Z2 - U2 = ', Z2 - U2);
- end
- else begin
- if (Z1 <= Zero) or (Z2 <= Zero) then begin
- writeln ('Because of unusual Radix = ', Radix);
- writeln (' or exact rational arithmetic a result');
- writeln (' Z1 = ', Z1, ' or Z2 = ', Z2);
- writeln (' of an extra precision test is inconsistent.');
- if Z1 = Z2 then
- end;
- if (Z1 <> Z2) or (Z1 > Zero) then begin
- X := Z1 / U1;
- Y := Z2 / U2;
- if Y > X then X := Y;
- Q := - ln (X);
- writeln ('Some subexpressions appear to be calculated');
- writeln ('extra precisely with about ');
- writeln (Q / ln (Radix), 'extra B-digits i.e. ');
- writeln ('roughly ', Q / ln (10),
- ' extra significant decimals.');
- end;
- writeln ('That feature is not tested further by this program.')
- end
- end;
- Pause;
- {=============================================}
- Milestone := 35;
- {=============================================}
- if Radix >= Two then
- begin
- X := W / (Radix * Radix);
- Y := X + One;
- Z := Y - X;
- T := Z + U2;
- X := T - Z;
- TestCondition (Failure, X = U2,
- 'Subtraction is not normlzd X=Y,X+Z<>Y+Z!');
- if X = U2 then
- writeln ('Subtraction appears to be normalized,',
- ' as it should be.');
- end;
- writeln;
- writeln ('Checking for guard digit on *, /, and -.');
- Y := F9 * One;
- Z := One * F9;
- X := F9 - Half;
- Y := (Y - Half) - X;
- Z := (Z - Half) - X;
- X := One + U2;
- T := X * Radix;
- R := Radix * X;
- X := T - Radix;
- X := X - Radix * U2;
- T := R - Radix;
- T := T - Radix * U2;
- X := X * (Radix - One);
- T := T * (Radix - One);
- if (X = Zero) and (Y = Zero) and (Z = Zero) and (T = Zero) then
- GMult := Yes
- else
- begin
- GMult := No;
- TestCondition (SeriousDefect, false,
- ' * lacks guard digit, 1*X <> X ');
- end;
- Z := Radix * U2;
- X := One + Z;
- Y := abs ((X + Z) - X * X) - U2;
- X := One - U2;
- Z := abs ((X - U2) - X * X) - U1;
- TestCondition (Failure, (Y <= Zero)
- and (Z <= Zero), ' * gets too many final digits wrong. ');
- Y := One - U2;
- X := One + U2;
- Z := One / Y;
- Y := Z - X;
- X := One / Three;
- Z := Three / Nine;
- X := X - Z;
- T := Nine / TwentySeven;
- Z := Z - T;
- TestCondition (Defect, (X = Zero) and (Y = Zero)
- and (Z = Zero), 'Division error > ulp, 1/3 <> 3/9 <> 9/27');
- Y := F9 / One;
- X := F9 - Half;
- Y := (Y - Half) - X;
- X := One + U2;
- T := X / One;
- X := T - X;
- if (X = Zero) and (Y = Zero) and (Z = Zero) then
- GDiv := Yes
- else
- begin
- GDiv := No;
- TestCondition (SeriousDefect, false,
- ' Division lacks guard digit so X/1 <> X');
- end;
- X := One / (One + U2);
- Y := X - Half - Half;
- TestCondition (SeriousDefect, Y < Zero,
- ' Computed value of 1/1.000..1 >= 1. ');
- X := One - U2;
- Y := One + Radix * U2;
- Z := X * Radix;
- T := Y * Radix;
- R := Z / Radix;
- StickyBit := T / Radix;
- X := R - X;
- Y := StickyBit - Y;
- TestCondition (Failure, (X = Zero) and (Y = Zero),
- ' * &or / gets too many last digits wrong');
- Y := One - U1;
- X := One - F9;
- Y := One - Y;
- T := Radix - U2;
- Z := Radix - BMinusU2;
- T := Radix - T;
- if (X = U1) and (Y = U1) and (Z = U2) and (T = U2) then
- GAddSub := Yes
- else
- begin
- GAddSub := No;
- TestCondition (SeriousDefect, false,
- '- lacks guard dig.,cancellation obscured');
- end;
-
- if (F9 <> One) and (F9 - One >= Zero) then begin
- TestCondition (SeriousDefect, false,
- 'comparison alleges (1-U1) < 1 although');
- writeln(' subtration yields (1-U1) - 1 = 0 , thereby vitiating');
- writeln(' such precautions against division by zero as');
- writeln(' ... if (X=1.0) then ..... else .../(X-1.0)...');
- end;
- if (GMult = Yes) and (GDiv = Yes) and (GAddSub = Yes) then
- writeln (' *, /, and - have guard digits, as they should.');
- {=============================================}
- Milestone := 40;
- {=============================================}
- Pause;
- writeln ('Checking rounding on multiply, divide and add/subtract.');
- RMult := Other;
- RDiv := Other;
- RAddSub := Other;
- RadixD2 := Radix / Two;
- A1 := Two;
- Done := false;
- repeat
- AInverse := Radix;
- repeat
- X := AInverse;
- AInverse := AInverse / A1;
- until Int (AInverse) <> AInverse;
- Done := (X = One) or (A1 > Three);
- if not Done then
- A1 := Nine + One;
- until Done;
- if X = One then
- A1 := Radix;
- AInverse := One / A1;
- X := A1;
- Y := AInverse;
- Done := false;
- repeat
- Z := X * Y - Half;
- TestCondition (Failure, Z = Half,
- ' X * (1/X) differs from 1. ');
- Done := X = Radix;
- X := Radix;
- Y := One / X;
- until Done;
- Y2 := One + U2;
- Y1 := One - U2;
- X := OneAndHalf - U2;
- Y := OneAndHalf + U2;
- Z := (X - U2) * Y2;
- T := Y * Y1;
- Z := Z - X;
- T := T - X;
- X := X * Y2;
- Y := (Y + U2) * Y1;
- X := X - OneAndHalf;
- Y := Y - OneAndHalf;
- if (X = Zero) and (Y = Zero) and (Z = Zero) and (T <= Zero) then
- begin
- X := (OneAndHalf + U2) * Y2;
- Y := OneAndHalf - U2 - U2;
- Z := OneAndHalf + U2 + U2;
- T := (OneAndHalf - U2) * Y1;
- X := X - (Z + U2);
- StickyBit := Y * Y1;
- S := Z * Y2;
- T := T - Y;
- Y := (U2 - Y) + StickyBit;
- Z := S - (Z + U2 + U2);
- StickyBit := (Y2 + U2) * Y1;
- Y1 := Y2 * Y1;
- StickyBit := StickyBit - Y2;
- Y1 := Y1 - Half;
- if (X = Zero) and (Y = Zero) and (Z = Zero) and (T = Zero)
- and ( StickyBit = Zero) and (Y1 = Half) then
- begin
- RMult := Rounded;
- writeln ('Multiplication appears to round correctly.');
- end
- else if (X + U2 = Zero) and (Y < Zero) and (Z + U2 = Zero)
- and (T < Zero) and (StickyBit + U2 = Zero)
- and (Y1 < Half) then
- begin
- RMult := Chopped;
- writeln ('Multiplication appears to chop.');
- end
- else
- writeln ('* is neither chopped nor correctly rounded.');
- if (RMult = Rounded) and (GMult = No) then
- notify('multiplication');
- end
- else
- writeln ('* is neither chopped nor correctly rounded.');
- {=============================================}
- Milestone := 45;
- {=============================================}
- Y2 := One + U2;
- Y1 := One - U2;
- Z := OneAndHalf + U2 + U2;
- X := Z / Y2;
- T := OneAndHalf - U2 - U2;
- Y := (T - U2) / Y1;
- Z := (Z + U2) / Y2;
- X := X - OneAndHalf;
- Y := Y - T;
- T := T / Y1;
- Z := Z - (OneAndHalf + U2);
- T := (U2 - OneAndHalf) + T;
- if not ((X > Zero) or (Y > Zero) or (Z > Zero) or (T > Zero)) then
- begin
- X := OneAndHalf / Y2;
- Y := OneAndHalf - U2;
- Z := OneAndHalf + U2;
- X := X - Y;
- T := OneAndHalf / Y1;
- Y := Y / Y1;
- T := T - (Z + U2);
- Y := Y - Z;
- Z := Z / Y2;
- Y1 := (Y2 + U2) / Y2;
- Z := Z - OneAndHalf;
- Y2 := Y1 - Y2;
- Y1 := (F9 - U1) / F9;
- if (X = Zero) and (Y = Zero) and (Z = Zero) and (T = Zero)
- and (Y2 = Zero) and (Y2 = Zero)
- and (Y1 - Half = F9 - Half ) then
- begin
- RDiv := Rounded;
- writeln ('Division appears to round correctly.');
- if GDiv = No then notify(' division ');
- end
- else if (X < Zero) and (Y < Zero) and (Z < Zero) and (T < Zero)
- and (Y2 < Zero) and (Y1 - Half < F9 - Half) then
- begin
- RDiv := Chopped;
- writeln ('Division appears to chop.');
- end;
- end;
- if RDiv = Other then
- writeln ('/ is neither chopped nor correctly rounded.');
- BInverse := One / Radix;
- TestCondition (Failure, (BInverse * Radix - Half = Half),
- ' Radix * ( 1 / Radix ) differs from 1. ');
- {=============================================}
- Milestone := 50;
- {=============================================}
- TestCondition (Failure, ((F9 + U1) - Half = Half)
- and ((BMinusU2 + U2 ) - One = Radix - One),
- 'Incomplete carry-propagation in Addition');
- X := One - U1 * U1;
- Y := One + U2 * (One - U2);
- Z := F9 - Half;
- X := (X - Half) - Z;
- Y := Y - One;
- if (X = Zero) and (Y = Zero) then
- begin
- RAddSub := Chopped;
- writeln ('Add/Subtract appears to be chopped.');
- end;
- if GAddSub = Yes then
- begin
- X := (Half + U2) * U2;
- Y := (Half - U2) * U2;
- X := One + X;
- Y := One + Y;
- X := (One + U2) - X;
- Y := One - Y;
- if (X = Zero) and (Y = Zero) then
- begin
- X := (Half + U2) * U1;
- Y := (Half - U2) * U1;
- X := One - X;
- Y := One - Y;
- X := F9 - X;
- Y := One - Y;
- if (X = Zero) and (Y = Zero) then
- begin
- RAddSub := Rounded;
- writeln ('Addition/Subtraction appears to round correctly.');
- if GAddSub = No then notify(' add/subtract ');
- end
- else
- writeln ('Addition/Subtraction neither rounds nor chops.');
- end
- else
- writeln ('Addition/Subtraction neither rounds nor chops.');
- end
- else
- writeln ('Addition/Subtraction neither rounds nor chops.');
- S := One;
- X := One + Half * (One + Half);
- Y := (One + U2) * Half;
- Z := X - Y;
- T := Y - X;
- StickyBit := Z + T;
- if StickyBit <> 0 then
- begin
- S := 0;
- NoErrors [Flaw] := NoErrors [Flaw] + 1;
- write('FLAW: (X - Y) + (Y - X) is non zero!');
- end;
- StickyBit := Zero;
- if (GMult = Yes) and (GDiv = Yes) and (GAddSub = Yes)
- and (RMult = Rounded) and (RDiv = Rounded)
- and (RAddSub = Rounded) and (Int (RadixD2) = RadixD2) then
- begin
- writeln (' Checking for sticky bit.');
- X := (Half + U1) * U2;
- Y := Half * U2;
- Z := One + Y;
- T := One + X;
- if (Z - One <= Zero) and (T - One >= U2) then
- begin
- Z := T + Y;
- Y := Z - X;
- if (Z - T >= U2) and (Y - T = Zero) then
- begin
- X := (Half + U1) * U1;
- Y := Half * U1;
- Z := One - Y;
- T := One - X;
- if (Z - One = Zero) and (T - F9 = Zero) then
- begin
- inline($fa/$fb);
- Z := (Half - U1) * U1;
- T := F9 - Z;
- Q := F9 - Y;
- if (T - F9 = Zero) and (F9 - U1 - Q = Zero) then
- begin
- Z := (One + U2) * OneAndHalf;
- T := (OneAndHalf + U2) - Z + U2;
- X := One + Half / Radix;
- Y := One + Radix * U2;
- Z := X * Y;
- if (T = Zero) and (X + Radix * U2 - Z = Zero) then
- begin
- if Radix <> Two then
- begin
- X := Two + U2;
- Y := X / Two;
- if (Y - One = Zero) then
- StickyBit := S;
- end
- else StickyBit := S;
- end;
- end;
- end;
- end;
- end;
- end;
- if StickyBit = One then
- writeln ('Sticky bit apparently used correctly.')
- else writeln ('Sticky bit used incorrectly or not at all.');
- if (GMult = No) or (GDiv = No) or (GAddSub = No) or (RMult = Other)
- or (RDiv = Other) or (RAddSub = Other) then begin
- TestCondition (Flaw, false,
- 'lack(s) of guard digits or failure(s) to');
- writeln('correctly round or chop (noted above) count as one');
- writeln('flaw in the final tally below.')
- end;
-
-
- {=============================================}
- Milestone := 60;
- {=============================================}
- writeln;
- writeln ('Does Multiplication commute? Testing on ', NoTrials,
- ' random pairs.');
- R9 := sqrt (3.0);
- RandomNumber1 := Third;
- I := 1;
- repeat
- X := Random;
- Y := Random;
- Z9 := Y * X;
- Z := X * Y;
- Z9 := Z - Z9;
- I := I + 1;
- until (I > NoTrials) or (Z9 <> Zero);
- if I = NoTrials then
- begin
- RandomNumber1 := One + Half / Three;
- RandomNumber2 := (U2 + U1) + One;
- Z := RandomNumber1 * RandomNumber2;
- Y := RandomNumber2 * RandomNumber1;
- Z9 := (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
- Three) * ((U2 + U1) + One);
- end;
- if not ((I = NoTrials) or (Z9 = Zero)) then
- begin
- NoErrors [Defect] := NoErrors [Defect] + 1;
- writeln ('DEFECT: X * Y = Y * X trail fails.');
- end
- else
- writeln ('No failures found in ', NoTrials, ' integer pairs.');
- end;
- end.