{$d- -- turn debug off to enable division by zero tests
#  // p5c - force c pre processor
}
(* Note: the statements "input: text;", "assign(input,'con:');",
 *    and "reset(input);" appear below as comments; some version of
 *    Pascal require you to activate one or more of these statements.
 *
 *  Some versions of TURBO Pascal (e.g. PC versions >= 4) require
 *  splitting the following source into several "units".  The goo
 *  between pairs of !! lines gives a way to do this.  If you have
 *  this file on a UNIX system, you can simply pipe it through
 *         sed /!!/d | /bin/sh
 *  to create files mainvars.pas, unit1.pas, unit2.pas, and par.pas;
 *  the first 3 are "units" needed in the fourth.  If using a UNIX
 *  system is inconvenient, you can do the splitting by hand:
 *  omit the lines that contain !! (that's what "sed /!!/d" does)
 *  and put the lines between each "cat >..." and the following
 *  "//GO.SYSIN DD" line into the file named on these lines.
 *)


program PARA (input, output);

{ (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 }

# if 0
   procedure Page;
     begin
       (* write(#$C) *) {FF in TURBO Pascal} writeln; writeln;
       end;
#endif

   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 *);

#if 0
   { alternative p5c version of Power ... }
   function Power(X, Y: real): real;
   function pow(x, y : real ): real; external;
   begin
      Power := pow(X,Y);
   end;

#else

   function Power (X, Y: real): real;

      (* 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 (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 *);
#endif

   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, ', relative error is ', (Y - X)/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 ('comparrison allegs that 1 * Z = ',
                     RandomNumber2);
               writeln ('differs from Z * 1 = ', RandomNumber1);
               end;
            if N > 0 then
               Pause;
            end;
         end;
      end (* TestPartialUnderflow *);

   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 and newline 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 = ', round(Radix):1);
   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.        ');


{=============================================}
   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 ',
               round(Precision):1);
   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;
(* For TURBO Pascal versions < 3, put everything after
   the next line into file "paran2.pas" . *)

{=============================================}
   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
               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.');

{=============================================}
   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 not ((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 underflow');
   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) * UnderflowThreshold)) 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.');
   write ('Did you find the correct value for N yet? ');
   ch := 'N';
   readln (input);
   while not eoln (input) do
      read (input, ch);
   Break := false;
   repeat
      write ('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.   ');

{=============================================}
   Milestone := 175;
{=============================================}
   writeln;
   for Index := 1 to 3 do
      begin
      case Index of
         1:
            Z := UnderflowThreshold;
         2:
            Z := E0;
         3:
            Z := PseudoZero;
         end;
      if Z <> 0 then
         begin
         V9 := sqrt (Z);
         Y := V9 * V9;
         if (Y / (One - Radix * E9) < Z)
               or (Y > (One + Radix * E9) * Z) then (* dgh: + E9 --> * E9 *)
            begin
            if V9 > U1 then
               begin
               NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
               write ('SERIOUS DEFECT:');
               end
            else
               begin
               NoErrors [Defect] := NoErrors [Defect] + 1;
               write ('DEFECT:');
               end;
            writeln ('  Comparison alleges that what prints as Z = ', Z);
            writeln ('is too far from sqrt(Z) ^ 2 = ', Y);
            end;
         end;
      end;

{=============================================}
   Milestone := 180;
{=============================================}
      for Index := 1 to 2 do
         begin
         if Index = 1 then
            Z := V
         else
            Z := V0;
         V9 := sqrt (Z);
         X := (One - Radix * E9) * V9;
         V9 := V9 * X;
         if ((V9 < (One - Two * Radix * E9) * Z) or (V9 > Z)) then
            begin
            Y := V9;
            if X < W then
               begin
               NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
               write ('SERIOUS ');
               end
            else
               NoErrors [Defect] := NoErrors [Defect] + 1;
            writeln ('DEFECT:  Comparison alleges that Z = ', Z);
            writeln ('is too far from sqrt(Z) ^ 2 is: ', Y);
            end;
         end;
{=============================================}
   Milestone := 190;
{=============================================}
   Pause;
   X := UnderflowThreshold * V;
   Y := Radix * Radix;
   if not ((X * Y >= One) and (X <= Y)) then
      begin
      if ((X * Y >= U1) and (X <= Y / U1)) then
         begin
         NoErrors [Flaw] := NoErrors [Flaw] + 1;
         write ('FLAW:');
         end
      else
         begin
         NoErrors [Defect] := NoErrors [Defect] + 1;
         write ('DEFECT: Badly');
         end;
      writeln (' unbalanced range; UnderflowThreshold * V = ');
      writeln (X, ' is too far from 1 .');
      end;
{=============================================}
   Milestone := 200;
{=============================================}
   for Index := 1 to 5 do
      begin
      X := F9;
      case Index of
         1:
            begin { Dummy Body }
            X := X;
            end;
         2:
            X := One + U2;
         3:
            X := V;
         4:
            X := UnderflowThreshold;
         5:
            X := Radix;
         end;
      Y := X;
      V9 := (Y / X - Half) - Half;
      if V9 <> 0 then
         begin
         if (V9 = - U1) and (Index < 5) then
            begin
            NoErrors [Flaw] := NoErrors [Flaw] + 1;
            write ('FLAW:');
            end
         else
            begin
            NoErrors [SeriousDefect] := NoErrors [SeriousDefect] + 1;
            write ('SERIOUS DEFECT:');
            end;
         writeln ('  X / X differs from 1 when X = ', X);
         writeln ('  instead, X / X - 1/2 - 1/2 = ', V9);
         writeln;
         end;
      end;
{=============================================}
   Milestone := 210;
{=============================================}
   MyZero := 0;
   writeln;
   writeln ('What message and/or values does Division by Zero produce?')
      ;
   writeln ('This can interupt your program. You can ',
         'skip this part if you wish.');
   write ('Do you wish to compute 1 / 0? ');
   readln (input);
   read (input, ch);
   if (ch = 'Y') or (ch = 'y') then
      writeln ('Trying to compute 1 / 0 produces: ', One / MyZero)
   else
      writeln ('O.K.');
   write ('Do you wish to compute 0 / 0? ');
   readln (input);
   read (input, ch);
   if (ch = 'Y') or (ch = 'y') then
      writeln ('Trying to compute 0 / 0 produces: ', MyZero / MyZero)
   else
      writeln ('O.K.');
{=============================================}
   Milestone := 220;
{=============================================}
   Pause;
   writeln;
   if NoErrors[Failure] > 0 then begin
      write ('The number of  FAILUREs  encountered =        ');
      writeln (NoErrors [Failure]);
      end;
   if NoErrors[SeriousDefect] > 0 then begin
      write ('The number of  SERIOUS DEFECTs  encountered = ');
      writeln (NoErrors [SeriousDefect]);
      end;
   if NoErrors[Defect] > 0 then begin
      write ('The number of  DEFECTs  encountered =         ');
      writeln (NoErrors [Defect]);
      end;
   if NoErrors[Flaw] > 0 then begin
      write ('The number of  FLAWs  encountered =           ');
      writeln (NoErrors [Flaw]);
      end;
   if (NoErrors [Failure] + NoErrors [SeriousDefect] + NoErrors [Defect]
         + NoErrors [Flaw]) > 0 then
      begin
      writeln;
      if (NoErrors [Failure] + NoErrors [SeriousDefect] + NoErrors [
            Defect] = 0) and (NoErrors [Flaw] > 0) then
         begin
         write ('The arithmetic diagnosed seems ');
         writeln ('Satisfactory though flawed.');
         end;
      if (NoErrors [Failure] + NoErrors [SeriousDefect] = 0)
            and ( NoErrors [Defect] > 0) then
         begin
         writeln ('The arithmetic diagnosed may be Acceptable');
         writeln ('despite inconvenient Defects.');
         end;
         (* dgh: Defect --> SeriousDefect in next line *)
      if (NoErrors [Failure] + NoErrors [SeriousDefect] > 0) then
         begin
         write ('The arithmetic diagnosed has ');
         writeln ('unacceptable Serious Defects.');
         end;
      if (NoErrors [Failure] > 0) then
         writeln ('Potentially fatal FAILURE may have spoiled this',
                  ' program''s subsequent diagnoses.');
      end
   else
      begin
      writeln ('No failures, defects nor flaws have been discovered.');
      if not ((RMult = Rounded) and (RDiv = Rounded)
            and (RAddSub = Rounded) and (RSqrt = Rounded)) then
         writeln ('The arithmetic diagnosed seems Satisfactory.')
      else begin
        if (StickyBit >= One)
            and ((Radix - Two) * (Radix - Nine - One) = 0) then begin
         write ('Rounding appears to conform to ');
         write ('the proposed IEEE standard P');
         if (Radix = Two)
               and ((Precision - Four * Three * Two) * ( Precision -
               TwentySeven - TwentySeven + One) = Zero) then
            write ('754')
         else
            write ('854');
         if IEEE then writeln('.')
          else begin
            writeln(',');
            writeln ('except possibly for Double Rounding',
               ' during Gradual Underflow.');
            end;
         end;
        writeln ('The arithmetic diagnosed appears to be Excellent!')
        end;
      end;
   writeln ('END OF TEST.');
   end (* PARA *).

