src/tests/Paranoia.c

   1 /*      A C version of Kahan's Floating Point Test "Paranoia"
   2
   3                         Thos Sumner, UCSF, Feb. 1985
   4                         David Gay, BTL, Jan. 1986
   5
   6         This is a rewrite from the Pascal version by
   7
   8                         B. A. Wichmann, 18 Jan. 1985
   9
  10         (and does NOT exhibit good C programming style).
  11
  12         Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
  13         compile with -DKR_headers or insert
  14 #define KR_headers
  15         at the beginning if you have an old-style C compiler.
  16
  17 (C) Apr 19 1983 in BASIC version by:
  18         Professor W. M. Kahan,
  19         567 Evans Hall
  20         Electrical Engineering & Computer Science Dept.
  21         University of California
  22         Berkeley, California 94720
  23         USA
  24
  25 converted to Pascal by:
  26         B. A. Wichmann
  27         National Physical Laboratory
  28         Teddington Middx
  29         TW11 OLW
  30         UK
  31
  32 converted to C by:
  33
  34         David M. Gay            and     Thos Sumner
  35         AT&T Bell Labs                  Computer Center, Rm. U-76
  36         600 Mountain Avenue             University of California
  37         Murray Hill, NJ 07974           San Francisco, CA 94143
  38         USA                             USA
  39
  40 with simultaneous corrections to the Pascal source (reflected
  41 in the Pascal source available over netlib).
  42 [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
  43
  44 Reports of results on various systems from all the versions
  45 of Paranoia are being collected by Richard Karpinski at the
  46 same address as Thos Sumner.  This includes sample outputs,
  47 bug reports, and criticisms.
  48
  49 You may copy this program freely if you acknowledge its source.
  50 Comments on the Pascal version to NPL, please.
  51
  52
  53 The C version catches signals from floating-point exceptions.
  54 If signal(SIGFPE,...) is unavailable in your environment, you may
  55 #define NOSIGNAL to comment out the invocations of signal.
  56
  57 This source file is too big for some C compilers, but may be split
  58 into pieces.  Comments containing "SPLIT" suggest convenient places
  59 for this splitting.  At the end of these comments is an "ed script"
  60 (for the UNIX(tm) editor ed) that will do this splitting.
  61
  62 By #defining Single when you compile this source, you may obtain
  63 a single-precision C version of Paranoia.
  64
  65
  66 The following is from the introductory commentary from Wichmann's work:
  67
  68 The BASIC program of Kahan is written in Microsoft BASIC using many
  69 facilities which have no exact analogy in Pascal.  The Pascal
  70 version below cannot therefore be exactly the same.  Rather than be
  71 a minimal transcription of the BASIC program, the Pascal coding
  72 follows the conventional style of block-structured languages.  Hence
  73 the Pascal version could be useful in producing versions in other
  74 structured languages.
  75
  76 Rather than use identifiers of minimal length (which therefore have
  77 little mnemonic significance), the Pascal version uses meaningful
  78 identifiers as follows [Note: A few changes have been made for C]:
  79
  80
  81 BASIC   C               BASIC   C               BASIC   C
  82
  83    A                       J                       S    StickyBit
  84    A1   AInverse           J0   NoErrors           T
  85    B    Radix                    [Failure]         T0   Underflow
  86    B1   BInverse           J1   NoErrors           T2   ThirtyTwo
  87    B2   RadixD2                  [SeriousDefect]   T5   OneAndHalf
  88    B9   BMinusU2           J2   NoErrors           T7   TwentySeven
  89    C                             [Defect]          T8   TwoForty
  90    C1   CInverse           J3   NoErrors           U    OneUlp
  91    D                             [Flaw]            U0   UnderflowThreshold
  92    D4   FourD              K    PageNo             U1
  93    E0                      L    Milestone          U2
  94    E1                      M                       V
  95    E2   Exp2               N                       V0
  96    E3                      N1                      V8
  97    E5   MinSqEr            O    Zero               V9
  98    E6   SqEr               O1   One                W
  99    E7   MaxSqEr            O2   Two                X
 100    E8                      O3   Three              X1
 101    E9                      O4   Four               X8
 102    F1   MinusOne           O5   Five               X9   Random1
 103    F2   Half               O8   Eight              Y
 104    F3   Third              O9   Nine               Y1
 105    F6                      P    Precision          Y2
 106    F9                      Q                       Y9   Random2
 107    G1   GMult              Q8                      Z
 108    G2   GDiv               Q9                      Z0   PseudoZero
 109    G3   GAddSub            R                       Z1
 110    H                       R1   RMult              Z2
 111    H1   HInverse           R2   RDiv               Z9
 112    I                       R3   RAddSub
 113    IO   NoTrials           R4   RSqrt
 114    I3   IEEE               R9   Random9
 115
 116    SqRWrng
 117
 118 All the variables in BASIC are true variables and in consequence,
 119 the program is more difficult to follow since the "constants" must
 120 be determined (the glossary is very helpful).  The Pascal version
 121 uses Real constants, but checks are added to ensure that the values
 122 are correctly converted by the compiler.
 123
 124 The major textual change to the Pascal version apart from the
 125 identifiersis that named procedures are used, inserting parameters
 126 wherehelpful.  New procedures are also introduced.  The
 127 correspondence is as follows:
 128
 129
 130 BASIC       Pascal
 131 lines
 132
 133   90- 140   Pause
 134  170- 250   Instructions
 135  380- 460   Heading
 136  480- 670   Characteristics
 137  690- 870   History
 138 2940-2950   Random
 139 3710-3740   NewD
 140 4040-4080   DoesYequalX
 141 4090-4110   PrintIfNPositive
 142 4640-4850   TestPartialUnderflow
 143
 144 =*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=
 145
 146 Below is an "ed script" that splits para.c into 10 files
 147 of the form part[1-8].c, subs.c, and msgs.c, plus a header
 148 file, paranoia.h, that these files require.
 149
 150 r paranoia.c
 151 $
 152 ?SPLIT
 153  .d
 154 +d
 155 -,$w msgs.c
 156 -,$d
 157 ?SPLIT
 158  .d
 159 +d
 160 -,$w subs.c
 161 -,$d
 162 ?part8
 163 +d
 164 ?include
 165  .,$w part8.c
 166  .,$d
 167 -d
 168 ?part7
 169 +d
 170 ?include
 171  .,$w part7.c
 172  .,$d
 173 -d
 174 ?part6
 175 +d
 176 ?include
 177  .,$w part6.c
 178  .,$d
 179 -d
 180 ?part5
 181 +d
 182 ?include
 183  .,$w part5.c
 184  .,$d
 185 -d
 186 ?part4
 187 +d
 188 ?include
 189  .,$w part4.c
 190  .,$d
 191 -d
 192 ?part3
 193 +d
 194 ?include
 195  .,$w part3.c
 196  .,$d
 197 -d
 198 ?part2
 199 +d
 200 ?include
 201  .,$w part2.c
 202  .,$d
 203 ?SPLIT
 204  .d
 205 1,/^#include/-1d
 206 1,$w part1.c
 207 /Computed constants/,$d
 208 1,$s/^int/extern &/
 209 1,$s/^FLOAT/extern &/
 210 1,$s/^char/extern &/
 211 1,$s! = .*!;!
 212 /^Guard/,/^Round/s/^/extern /
 213 /^jmp_buf/s/^/extern /
 214 /^Sig_type/s/^/extern /
 215 s/$/\
 216 extern void sigfpe(INT);/
 217 w paranoia.h
 218 q
 219
 220 */
 221
 222 #include <stdio.h>
 223 #ifndef NOSIGNAL
 224 #include <signal.h>
 225 #endif
 226 #include <setjmp.h>
 227
 228 #ifdef Single
 229 #define FLOAT float
 230 #define FABS(x) (float)fabs((double)(x))
 231 #define FLOOR(x) (float)floor((double)(x))
 232 #define LOG(x) (float)log((double)(x))
 233 #define POW(x,y) (float)pow((double)(x),(double)(y))
 234 #define SQRT(x) (float)sqrt((double)(x))
 235 #else
 236 #define FLOAT double
 237 #define FABS(x) fabs(x)
 238 #define FLOOR(x) floor(x)
 239 #define LOG(x) log(x)
 240 #define POW(x,y) pow(x,y)
 241 #define SQRT(x) sqrt(x)
 242 #endif
 243
 244 jmp_buf ovfl_buf;
 245 #ifdef KR_headers
 246 #define VOID /* void */
 247 #define INT /* int */
 248 #define FP /* FLOAT */
 249 #define CHARP /* char * */
 250 #define CHARPP /* char ** */
 251 extern double fabs(), floor(), log(), pow(), sqrt();
 252 extern void exit();
 253 typedef void (*Sig_type)();
 254 FLOAT Sign(), Random();
 255 extern void BadCond();
 256 extern void SqXMinX();
 257 extern void TstCond();
 258 extern void notify();
 259 extern int read();
 260 #else
 261 #define VOID void
 262 #define INT int
 263 #define FP FLOAT
 264 #define CHARP char *
 265 #define CHARPP char **
 266 #ifdef __STDC__
 267 #include <stdlib.h>
 268 #include <math.h>
 269 #else
 270 #ifdef __cplusplus
 271 extern "C" {
 272 #endif
 273 extern double fabs(double), floor(double), log(double);
 274 extern double pow(double,double), sqrt(double);
 275 extern void exit(INT);
 276 #ifdef __cplusplus
 277         }
 278 #endif
 279 #endif
 280 typedef void (*Sig_type)(int);
 281 FLOAT Sign(FLOAT), Random(void);
 282 extern void BadCond(int, char*);
 283 extern void SqXMinX(int);
 284 extern void TstCond(int, int, char*);
 285 extern void notify(char*);
 286 extern int read(int, char*, int);
 287 #endif
 288 #undef V9
 289 extern void Characteristics(VOID);
 290 extern void Heading(VOID);
 291 extern void History(VOID);
 292 extern void Instructions(VOID);
 293 extern void IsYeqX(VOID);
 294 extern void NewD(VOID);
 295 extern void Pause(VOID);
 296 extern void PrintIfNPositive(VOID);
 297 extern void SR3750(VOID);
 298 extern void SR3980(VOID);
 299 extern void TstPtUf(VOID);
 300
 301 Sig_type sigsave;
 302
 303 #define KEYBOARD 0
 304
 305 FLOAT Radix, BInvrse, RadixD2, BMinusU2;
 306
 307 /*Small floating point constants.*/
 308 FLOAT Zero = 0.0;
 309 FLOAT Half = 0.5;
 310 FLOAT One = 1.0;
 311 FLOAT Two = 2.0;
 312 FLOAT Three = 3.0;
 313 FLOAT Four = 4.0;
 314 FLOAT Five = 5.0;
 315 FLOAT Eight = 8.0;
 316 FLOAT Nine = 9.0;
 317 FLOAT TwentySeven = 27.0;
 318 FLOAT ThirtyTwo = 32.0;
 319 FLOAT TwoForty = 240.0;
 320 FLOAT MinusOne = -1.0;
 321 FLOAT OneAndHalf = 1.5;
 322 /*Integer constants*/
 323 int NoTrials = 20; /*Number of tests for commutativity. */
 324 #define False 0
 325 #define True 1
 326
 327 /* Definitions for declared types
 328         Guard == (Yes, No);
 329         Rounding == (Chopped, Rounded, Other);
 330         Message == packed array [1..40] of char;
 331         Class == (Flaw, Defect, Serious, Failure);
 332           */
 333 #define Yes 1
 334 #define No  0
 335 #define Chopped 2
 336 #define Rounded 1
 337 #define Other   0
 338 #define Flaw    3
 339 #define Defect  2
 340 #define Serious 1
 341 #define Failure 0
 342 typedef int Guard, Rounding, Class;
 343 typedef char Message;
 344
 345 /* Declarations of Variables */
 346 int Indx;
 347 char ch[8];
 348 FLOAT AInvrse, A1;
 349 FLOAT C, CInvrse;
 350 FLOAT D, FourD;
 351 FLOAT E0, E1, Exp2, E3, MinSqEr;
 352 FLOAT SqEr, MaxSqEr, E9;
 353 FLOAT Third;
 354 FLOAT F6, F9;
 355 FLOAT H, HInvrse;
 356 int I;
 357 FLOAT StickyBit, J;
 358 FLOAT MyZero;
 359 FLOAT Precision;
 360 FLOAT Q, Q9;
 361 FLOAT R, Random9;
 362 FLOAT T, Underflow, S;
 363 FLOAT OneUlp, UfThold, U1, U2;
 364 FLOAT V, V0, V9;
 365 FLOAT W;
 366 FLOAT X, X1, X2, X8, Random1;
 367 FLOAT Y, Y1, Y2, Random2;
 368 FLOAT Z, PseudoZero, Z1, Z2, Z9;
 369 int ErrCnt[4];
 370 int fpecount;
 371 int Milestone;
 372 int PageNo;
 373 int M, N, N1;
 374 Guard GMult, GDiv, GAddSub;
 375 Rounding RMult, RDiv, RAddSub, RSqrt;
 376 int Break, Done, NotMonot, Monot, Anomaly, IEEE,
 377                 SqRWrng, UfNGrad;
 378 /* Computed constants. */
 379 /*U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */
 380 /*U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */
 381
 382 /* floating point exception receiver */
 383  void
 384 sigfpe(INT x)
 385 {
 386         fpecount++;
 387         printf("\n* * * FLOATING-POINT ERROR %d * * *\n", x);
 388         fflush(stdout);
 389         if (sigsave) {
 390 #ifndef NOSIGNAL
 391                 signal(SIGFPE, sigsave);
 392 #endif
 393                 sigsave = 0;
 394                 longjmp(ovfl_buf, 1);
 395                 }
 396         exit(1);
 397 }
 398
 399 main(VOID)
 400 {
 401         /* First two assignments use integer right-hand sides. */
 402         Zero = 0;
 403         One = 1;
 404         Two = One + One;
 405         Three = Two + One;
 406         Four = Three + One;
 407         Five = Four + One;
 408         Eight = Four + Four;
 409         Nine = Three * Three;
 410         TwentySeven = Nine * Three;
 411         ThirtyTwo = Four * Eight;
 412         TwoForty = Four * Five * Three * Four;
 413         MinusOne = -One;
 414         Half = One / Two;
 415         OneAndHalf = One + Half;
 416         ErrCnt[Failure] = 0;
 417         ErrCnt[Serious] = 0;
 418         ErrCnt[Defect] = 0;
 419         ErrCnt[Flaw] = 0;
 420         PageNo = 1;
 421         /*=============================================*/
 422         Milestone = 0;
 423         /*=============================================*/
 424 #ifndef NOSIGNAL
 425         signal(SIGFPE, sigfpe);
 426 #endif
 427         Instructions();
 428         Pause();
 429         Heading();
 430         Pause();
 431         Characteristics();
 432         Pause();
 433         History();
 434         Pause();
 435         /*=============================================*/
 436         Milestone = 7;
 437         /*=============================================*/
 438         printf("Program is now RUNNING tests on small integers:\n");
 439
 440         TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
 441                    && (One > Zero) && (One + One == Two),
 442                         "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
 443         Z = - Zero;
 444         if (Z != 0.0) {
 445                 ErrCnt[Failure] = ErrCnt[Failure] + 1;
 446                 printf("Comparison alleges that -0.0 is Non-zero!\n");
 447                 U2 = 0.001;
 448                 Radix = 1;
 449                 TstPtUf();
 450                 }
 451         TstCond (Failure, (Three == Two + One) && (Four == Three + One)
 452                    && (Four + Two * (- Two) == Zero)
 453                    && (Four - Three - One == Zero),
 454                    "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
 455         TstCond (Failure, (MinusOne == (0 - One))
 456                    && (MinusOne + One == Zero ) && (One + MinusOne == Zero)
 457                    && (MinusOne + FABS(One) == Zero)
 458                    && (MinusOne + MinusOne * MinusOne == Zero),
 459                    "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
 460         TstCond (Failure, Half + MinusOne + Half == Zero,
 461                   "1/2 + (-1) + 1/2 != 0");
 462         /*=============================================*/
 463         /*SPLIT
 464         {
 465                 extern void part2(VOID), part3(VOID), part4(VOID),
 466                         part5(VOID), part6(VOID), part7(VOID);
 467                 int part8(VOID);
 468
 469                 part2();
 470                 part3();
 471                 part4();
 472                 part5();
 473                 part6();
 474                 part7();
 475                 return part8();
 476                 }
 477         }
 478 #include "paranoia.h"
 479 void part2(VOID){
 480 */
 481         Milestone = 10;
 482         /*=============================================*/
 483         TstCond (Failure, (Nine == Three * Three)
 484                    && (TwentySeven == Nine * Three) && (Eight == Four + Four)
 485                    && (ThirtyTwo == Eight * Four)
 486                    && (ThirtyTwo - TwentySeven - Four - One == Zero),
 487                    "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
 488         TstCond (Failure, (Five == Four + One) &&
 489                         (TwoForty == Four * Five * Three * Four)
 490                    && (TwoForty / Three - Four * Four * Five == Zero)
 491                    && ( TwoForty / Four - Five * Three * Four == Zero)
 492                    && ( TwoForty / Five - Four * Three * Four == Zero),
 493                   "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
 494         if (ErrCnt[Failure] == 0) {
 495                 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
 496                 printf("\n");
 497                 }
 498         printf("Searching for Radix and Precision.\n");
 499         W = One;
 500         do  {
 501                 W = W + W;
 502                 Y = W + One;
 503                 Z = Y - W;
 504                 Y = Z - One;
 505                 } while (MinusOne + FABS(Y) < Zero);
 506         /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/
 507         Precision = Zero;
 508         Y = One;
 509         do  {
 510                 Radix = W + Y;
 511                 Y = Y + Y;
 512                 Radix = Radix - W;
 513                 } while ( Radix == Zero);
 514         if (Radix < Two) Radix = One;
 515         printf("Radix = %f .\n", Radix);
 516         if (Radix != 1) {
 517                 W = One;
 518                 do  {
 519                         Precision = Precision + One;
 520                         W = W * Radix;
 521                         Y = W + One;
 522                         } while ((Y - W) == One);
 523                 }
 524         /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1
 525                                                       ...*/
 526         U1 = One / W;
 527         U2 = Radix * U1;
 528         printf("Closest relative separation found is U1 = %.7e .\n\n", U1);
 529         printf("Recalculating radix and precision\n ");
 530
 531         /*save old values*/
 532         E0 = Radix;
 533         E1 = U1;
 534         E9 = U2;
 535         E3 = Precision;
 536
 537         X = Four / Three;
 538         Third = X - One;
 539         F6 = Half - Third;
 540         X = F6 + F6;
 541         X = FABS(X - Third);
 542         if (X < U2) X = U2;
 543
 544         /*... now X = (unknown no.) ulps of 1+...*/
 545         do  {
 546                 U2 = X;
 547                 Y = Half * U2 + ThirtyTwo * U2 * U2;
 548                 Y = One + Y;
 549                 X = Y - One;
 550                 } while ( ! ((U2 <= X) || (X <= Zero)));
 551
 552         /*... now U2 == 1 ulp of 1 + ... */
 553         X = Two / Three;
 554         F6 = X - Half;
 555         Third = F6 + F6;
 556         X = Third - Half;
 557         X = FABS(X + F6);
 558         if (X < U1) X = U1;
 559
 560         /*... now  X == (unknown no.) ulps of 1 -... */
 561         do  {
 562                 U1 = X;
 563                 Y = Half * U1 + ThirtyTwo * U1 * U1;
 564                 Y = Half - Y;
 565                 X = Half + Y;
 566                 Y = Half - X;
 567                 X = Half + Y;
 568                 } while ( ! ((U1 <= X) || (X <= Zero)));
 569         /*... now U1 == 1 ulp of 1 - ... */
 570         if (U1 == E1) printf("confirms closest relative separation U1 .\n");
 571         else printf("gets better closest relative separation U1 = %.7e .\n", U1);
 572         W = One / U1;
 573         F9 = (Half - U1) + Half;
 574         Radix = FLOOR(0.01 + U2 / U1);
 575         if (Radix == E0) printf("Radix confirmed.\n");
 576         else printf("MYSTERY: recalculated Radix = %.7e .\n", Radix);
 577         TstCond (Defect, Radix <= Eight + Eight,
 578                    "Radix is too big: roundoff problems");
 579         TstCond (Flaw, (Radix == Two) || (Radix == 10)
 580                    || (Radix == One), "Radix is not as good as 2 or 10");
 581         /*=============================================*/
 582         Milestone = 20;
 583         /*=============================================*/
 584         TstCond (Failure, F9 - Half < Half,
 585                    "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
 586         X = F9;
 587         I = 1;
 588         Y = X - Half;
 589         Z = Y - Half;
 590         TstCond (Failure, (X != One)
 591                    || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
 592         X = One + U2;
 593         I = 0;
 594         /*=============================================*/
 595         Milestone = 25;
 596         /*=============================================*/
 597         /*... BMinusU2 = nextafter(Radix, 0) */
 598         BMinusU2 = Radix - One;
 599         BMinusU2 = (BMinusU2 - U2) + One;
 600         /* Purify Integers */
 601         if (Radix != One)  {
 602                 X = - TwoForty * LOG(U1) / LOG(Radix);
 603                 Y = FLOOR(Half + X);
 604                 if (FABS(X - Y) * Four < One) X = Y;
 605                 Precision = X / TwoForty;
 606                 Y = FLOOR(Half + Precision);
 607                 if (FABS(Precision - Y) * TwoForty < Half) Precision = Y;
 608                 }
 609         if ((Precision != FLOOR(Precision)) || (Radix == One)) {
 610                 printf("Precision cannot be characterized by an Integer number\n");
 611                 printf("of significant digits but, by itself, this is a minor flaw.\n");
 612                 }
 613         if (Radix == One)
 614                 printf("logarithmic encoding has precision characterized solely by U1.\n");
 615         else printf("The number of significant digits of the Radix is %f .\n",
 616                         Precision);
 617         TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
 618                    "Precision worse than 5 decimal figures  ");
 619         /*=============================================*/
 620         Milestone = 30;
 621         /*=============================================*/
 622         /* Test for extra-precise subepressions */
 623         X = FABS(((Four / Three - One) - One / Four) * Three - One / Four);
 624         do  {
 625                 Z2 = X;
 626                 X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
 627                 } while ( ! ((Z2 <= X) || (X <= Zero)));
 628         X = Y = Z = FABS((Three / Four - Two / Three) * Three - One / Four);
 629         do  {
 630                 Z1 = Z;
 631                 Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
 632                         + One / Two)) + One / Two;
 633                 } while ( ! ((Z1 <= Z) || (Z <= Zero)));
 634         do  {
 635                 do  {
 636                         Y1 = Y;
 637                         Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
 638                                 )) + Half;
 639                         } while ( ! ((Y1 <= Y) || (Y <= Zero)));
 640                 X1 = X;
 641                 X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
 642                 } while ( ! ((X1 <= X) || (X <= Zero)));
 643         if ((X1 != Y1) || (X1 != Z1)) {
 644                 BadCond(Serious, "Disagreements among the values X1, Y1, Z1,\n");
 645                 printf("respectively  %.7e,  %.7e,  %.7e,\n", X1, Y1, Z1);
 646                 printf("are symptoms of inconsistencies introduced\n");
 647                 printf("by extra-precise evaluation of arithmetic subexpressions.\n");
 648                 notify("Possibly some part of this");
 649                 if ((X1 == U1) || (Y1 == U1) || (Z1 == U1))  printf(
 650                         "That feature is not tested further by this program.\n") ;
 651                 }
 652         else  {
 653                 if ((Z1 != U1) || (Z2 != U2)) {
 654                         if ((Z1 >= U1) || (Z2 >= U2)) {
 655                                 BadCond(Failure, "");
 656                                 notify("Precision");
 657                                 printf("\tU1 = %.7e, Z1 - U1 = %.7e\n",U1,Z1-U1);
 658                                 printf("\tU2 = %.7e, Z2 - U2 = %.7e\n",U2,Z2-U2);
 659                                 }
 660                         else {
 661                                 if ((Z1 <= Zero) || (Z2 <= Zero)) {
 662                                         printf("Because of unusual Radix = %f", Radix);
 663                                         printf(", or exact rational arithmetic a result\n");
 664                                         printf("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
 665                                         notify("of an\nextra-precision");
 666                                         }
 667                                 if (Z1 != Z2 || Z1 > Zero) {
 668                                         X = Z1 / U1;
 669                                         Y = Z2 / U2;
 670                                         if (Y > X) X = Y;
 671                                         Q = - LOG(X);
 672                                         printf("Some subexpressions appear to be calculated extra\n");
 673                                         printf("precisely with about %g extra B-digits, i.e.\n",
 674                                                 (Q / LOG(Radix)));
 675                                         printf("roughly %g extra significant decimals.\n",
 676                                                 Q / LOG(10.));
 677                                         }
 678                                 printf("That feature is not tested further by this program.\n");
 679                                 }
 680                         }
 681                 }
 682         Pause();
 683         /*=============================================*/
 684         /*SPLIT
 685         }
 686 #include "paranoia.h"
 687 void part3(VOID){
 688 */
 689         Milestone = 35;
 690         /*=============================================*/
 691         if (Radix >= Two) {
 692                 X = W / (Radix * Radix);
 693                 Y = X + One;
 694                 Z = Y - X;
 695                 T = Z + U2;
 696                 X = T - Z;
 697                 TstCond (Failure, X == U2,
 698                         "Subtraction is not normalized X=Y,X+Z != Y+Z!");
 699                 if (X == U2) printf(
 700                         "Subtraction appears to be normalized, as it should be.");
 701                 }
 702         printf("\nChecking for guard digit in *, /, and -.\n");
 703         Y = F9 * One;
 704         Z = One * F9;
 705         X = F9 - Half;
 706         Y = (Y - Half) - X;
 707         Z = (Z - Half) - X;
 708         X = One + U2;
 709         T = X * Radix;
 710         R = Radix * X;
 711         X = T - Radix;
 712         X = X - Radix * U2;
 713         T = R - Radix;
 714         T = T - Radix * U2;
 715         X = X * (Radix - One);
 716         T = T * (Radix - One);
 717         if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)) GMult = Yes;
 718         else {
 719                 GMult = No;
 720                 TstCond (Serious, False,
 721                         "* lacks a Guard Digit, so 1*X != X");
 722                 }
 723         Z = Radix * U2;
 724         X = One + Z;
 725         Y = FABS((X + Z) - X * X) - U2;
 726         X = One - U2;
 727         Z = FABS((X - U2) - X * X) - U1;
 728         TstCond (Failure, (Y <= Zero)
 729                    && (Z <= Zero), "* gets too many final digits wrong.\n");
 730         Y = One - U2;
 731         X = One + U2;
 732         Z = One / Y;
 733         Y = Z - X;
 734         X = One / Three;
 735         Z = Three / Nine;
 736         X = X - Z;
 737         T = Nine / TwentySeven;
 738         Z = Z - T;
 739         TstCond(Defect, X == Zero && Y == Zero && Z == Zero,
 740                 "Division lacks a Guard Digit, so error can exceed 1 ulp\n\
 741 or  1/3  and  3/9  and  9/27 may disagree");
 742         Y = F9 / One;
 743         X = F9 - Half;
 744         Y = (Y - Half) - X;
 745         X = One + U2;
 746         T = X / One;
 747         X = T - X;
 748         if ((X == Zero) && (Y == Zero) && (Z == Zero)) GDiv = Yes;
 749         else {
 750                 GDiv = No;
 751                 TstCond (Serious, False,
 752                         "Division lacks a Guard Digit, so X/1 != X");
 753                 }
 754         X = One / (One + U2);
 755         Y = X - Half - Half;
 756         TstCond (Serious, Y < Zero,
 757                    "Computed value of 1/1.000..1 >= 1");
 758         X = One - U2;
 759         Y = One + Radix * U2;
 760         Z = X * Radix;
 761         T = Y * Radix;
 762         R = Z / Radix;
 763         StickyBit = T / Radix;
 764         X = R - X;
 765         Y = StickyBit - Y;
 766         TstCond (Failure, X == Zero && Y == Zero,
 767                         "* and/or / gets too many last digits wrong");
 768         Y = One - U1;
 769         X = One - F9;
 770         Y = One - Y;
 771         T = Radix - U2;
 772         Z = Radix - BMinusU2;
 773         T = Radix - T;
 774         if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2)) GAddSub = Yes;
 775         else {
 776                 GAddSub = No;
 777                 TstCond (Serious, False,
 778                         "- lacks Guard Digit, so cancellation is obscured");
 779                 }
 780         if (F9 != One && F9 - One >= Zero) {
 781                 BadCond(Serious, "comparison alleges  (1-U1) < 1  although\n");
 782                 printf("  subtraction yields  (1-U1) - 1 = 0 , thereby vitiating\n");
 783                 printf("  such precautions against division by zero as\n");
 784                 printf("  ...  if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
 785                 }
 786         if (GMult == Yes && GDiv == Yes && GAddSub == Yes) printf(
 787                 "     *, /, and - appear to have guard digits, as they should.\n");
 788         /*=============================================*/
 789         Milestone = 40;
 790         /*=============================================*/
 791         Pause();
 792         printf("Checking rounding on multiply, divide and add/subtract.\n");
 793         RMult = Other;
 794         RDiv = Other;
 795         RAddSub = Other;
 796         RadixD2 = Radix / Two;
 797         A1 = Two;
 798         Done = False;
 799         do  {
 800                 AInvrse = Radix;
 801                 do  {
 802                         X = AInvrse;
 803                         AInvrse = AInvrse / A1;
 804                         } while ( ! (FLOOR(AInvrse) != AInvrse));
 805                 Done = (X == One) || (A1 > Three);
 806                 if (! Done) A1 = Nine + One;
 807                 } while ( ! (Done));
 808         if (X == One) A1 = Radix;
 809         AInvrse = One / A1;
 810         X = A1;
 811         Y = AInvrse;
 812         Done = False;
 813         do  {
 814                 Z = X * Y - Half;
 815                 TstCond (Failure, Z == Half,
 816                         "X * (1/X) differs from 1");
 817                 Done = X == Radix;
 818                 X = Radix;
 819                 Y = One / X;
 820                 } while ( ! (Done));
 821         Y2 = One + U2;
 822         Y1 = One - U2;
 823         X = OneAndHalf - U2;
 824         Y = OneAndHalf + U2;
 825         Z = (X - U2) * Y2;
 826         T = Y * Y1;
 827         Z = Z - X;
 828         T = T - X;
 829         X = X * Y2;
 830         Y = (Y + U2) * Y1;
 831         X = X - OneAndHalf;
 832         Y = Y - OneAndHalf;
 833         if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
 834                 X = (OneAndHalf + U2) * Y2;
 835                 Y = OneAndHalf - U2 - U2;
 836                 Z = OneAndHalf + U2 + U2;
 837                 T = (OneAndHalf - U2) * Y1;
 838                 X = X - (Z + U2);
 839                 StickyBit = Y * Y1;
 840                 S = Z * Y2;
 841                 T = T - Y;
 842                 Y = (U2 - Y) + StickyBit;
 843                 Z = S - (Z + U2 + U2);
 844                 StickyBit = (Y2 + U2) * Y1;
 845                 Y1 = Y2 * Y1;
 846                 StickyBit = StickyBit - Y2;
 847                 Y1 = Y1 - Half;
 848                 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
 849                         && ( StickyBit == Zero) && (Y1 == Half)) {
 850                         RMult = Rounded;
 851                         printf("Multiplication appears to round correctly.\n");
 852                         }
 853                 else    if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
 854                                 && (T < Zero) && (StickyBit + U2 == Zero)
 855                                 && (Y1 < Half)) {
 856                                 RMult = Chopped;
 857                                 printf("Multiplication appears to chop.\n");
 858                                 }
 859                         else printf("* is neither chopped nor correctly rounded.\n");
 860                 if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
 861                 }
 862         else printf("* is neither chopped nor correctly rounded.\n");
 863         /*=============================================*/
 864         Milestone = 45;
 865         /*=============================================*/
 866         Y2 = One + U2;
 867         Y1 = One - U2;
 868         Z = OneAndHalf + U2 + U2;
 869         X = Z / Y2;
 870         T = OneAndHalf - U2 - U2;
 871         Y = (T - U2) / Y1;
 872         Z = (Z + U2) / Y2;
 873         X = X - OneAndHalf;
 874         Y = Y - T;
 875         T = T / Y1;
 876         Z = Z - (OneAndHalf + U2);
 877         T = (U2 - OneAndHalf) + T;
 878         if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
 879                 X = OneAndHalf / Y2;
 880                 Y = OneAndHalf - U2;
 881                 Z = OneAndHalf + U2;
 882                 X = X - Y;
 883                 T = OneAndHalf / Y1;
 884                 Y = Y / Y1;
 885                 T = T - (Z + U2);
 886                 Y = Y - Z;
 887                 Z = Z / Y2;
 888                 Y1 = (Y2 + U2) / Y2;
 889                 Z = Z - OneAndHalf;
 890                 Y2 = Y1 - Y2;
 891                 Y1 = (F9 - U1) / F9;
 892                 if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
 893                         && (Y2 == Zero) && (Y2 == Zero)
 894                         && (Y1 - Half == F9 - Half )) {
 895                         RDiv = Rounded;
 896                         printf("Division appears to round correctly.\n");
 897                         if (GDiv == No) notify("Division");
 898                         }
 899                 else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
 900                         && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
 901                         RDiv = Chopped;
 902                         printf("Division appears to chop.\n");
 903                         }
 904                 }
 905         if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
 906         BInvrse = One / Radix;
 907         TstCond (Failure, (BInvrse * Radix - Half == Half),
 908                    "Radix * ( 1 / Radix ) differs from 1");
 909         /*=============================================*/
 910         /*SPLIT
 911         }
 912 #include "paranoia.h"
 913 void part4(VOID){
 914 */
 915         Milestone = 50;
 916         /*=============================================*/
 917         TstCond (Failure, ((F9 + U1) - Half == Half)
 918                    && ((BMinusU2 + U2 ) - One == Radix - One),
 919                    "Incomplete carry-propagation in Addition");
 920         X = One - U1 * U1;
 921         Y = One + U2 * (One - U2);
 922         Z = F9 - Half;
 923         X = (X - Half) - Z;
 924         Y = Y - One;
 925         if ((X == Zero) && (Y == Zero)) {
 926                 RAddSub = Chopped;
 927                 printf("Add/Subtract appears to be chopped.\n");
 928                 }
 929         if (GAddSub == Yes) {
 930                 X = (Half + U2) * U2;
 931                 Y = (Half - U2) * U2;
 932                 X = One + X;
 933                 Y = One + Y;
 934                 X = (One + U2) - X;
 935                 Y = One - Y;
 936                 if ((X == Zero) && (Y == Zero)) {
 937                         X = (Half + U2) * U1;
 938                         Y = (Half - U2) * U1;
 939                         X = One - X;
 940                         Y = One - Y;
 941                         X = F9 - X;
 942                         Y = One - Y;
 943                         if ((X == Zero) && (Y == Zero)) {
 944                                 RAddSub = Rounded;
 945                                 printf("Addition/Subtraction appears to round correctly.\n");
 946                                 if (GAddSub == No) notify("Add/Subtract");
 947                                 }
 948                         else printf("Addition/Subtraction neither rounds nor chops.\n");
 949                         }
 950                 else printf("Addition/Subtraction neither rounds nor chops.\n");
 951                 }
 952         else printf("Addition/Subtraction neither rounds nor chops.\n");
 953         S = One;
 954         X = One + Half * (One + Half);
 955         Y = (One + U2) * Half;
 956         Z = X - Y;
 957         T = Y - X;
 958         StickyBit = Z + T;
 959         if (StickyBit != Zero) {
 960                 S = Zero;
 961                 BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
 962                 }
 963         StickyBit = Zero;
 964         if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
 965                 && (RMult == Rounded) && (RDiv == Rounded)
 966                 && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
 967                 printf("Checking for sticky bit.\n");
 968                 X = (Half + U1) * U2;
 969                 Y = Half * U2;
 970                 Z = One + Y;
 971                 T = One + X;
 972                 if ((Z - One <= Zero) && (T - One >= U2)) {
 973                         Z = T + Y;
 974                         Y = Z - X;
 975                         if ((Z - T >= U2) && (Y - T == Zero)) {
 976                                 X = (Half + U1) * U1;
 977                                 Y = Half * U1;
 978                                 Z = One - Y;
 979                                 T = One - X;
 980                                 if ((Z - One == Zero) && (T - F9 == Zero)) {
 981                                         Z = (Half - U1) * U1;
 982                                         T = F9 - Z;
 983                                         Q = F9 - Y;
 984                                         if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
 985                                                 Z = (One + U2) * OneAndHalf;
 986                                                 T = (OneAndHalf + U2) - Z + U2;
 987                                                 X = One + Half / Radix;
 988                                                 Y = One + Radix * U2;
 989                                                 Z = X * Y;
 990                                                 if (T == Zero && X + Radix * U2 - Z == Zero) {
 991                                                         if (Radix != Two) {
 992                                                                 X = Two + U2;
 993                                                                 Y = X / Two;
 994                                                                 if ((Y - One == Zero)) StickyBit = S;
 995                                                                 }
 996                                                         else StickyBit = S;
 997                                                         }
 998                                                 }
 999                                         }
1000                                 }
1001                         }
1002                 }
1003         if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
1004         else printf("Sticky bit used incorrectly or not at all.\n");
1005         TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
1006                         RMult == Other || RDiv == Other || RAddSub == Other),
1007                 "lack(s) of guard digits or failure(s) to correctly round or chop\n\
1008 (noted above) count as one flaw in the final tally below");
1009         /*=============================================*/
1010         Milestone = 60;
1011         /*=============================================*/
1012         printf("\n");
1013         printf("Does Multiplication commute?  ");
1014         printf("Testing on %d random pairs.\n", NoTrials);
1015         Random9 = SQRT(3.0);
1016         Random1 = Third;
1017         I = 1;
1018         do  {
1019                 X = Random();
1020                 Y = Random();
1021                 Z9 = Y * X;
1022                 Z = X * Y;
1023                 Z9 = Z - Z9;
1024                 I = I + 1;
1025                 } while ( ! ((I > NoTrials) || (Z9 != Zero)));
1026         if (I == NoTrials) {
1027                 Random1 = One + Half / Three;
1028                 Random2 = (U2 + U1) + One;
1029                 Z = Random1 * Random2;
1030                 Y = Random2 * Random1;
1031                 Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
1032                         Three) * ((U2 + U1) + One);
1033                 }
1034         if (! ((I == NoTrials) || (Z9 == Zero)))
1035                 BadCond(Defect, "X * Y == Y * X trial fails.\n");
1036         else printf("     No failures found in %d integer pairs.\n", NoTrials);
1037         /*=============================================*/
1038         Milestone = 70;
1039         /*=============================================*/
1040         printf("\nRunning test of square root(x).\n");
1041         TstCond (Failure, (Zero == SQRT(Zero))
1042                    && (- Zero == SQRT(- Zero))
1043                    && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
1044         MinSqEr = Zero;
1045         MaxSqEr = Zero;
1046         J = Zero;
1047         X = Radix;
1048         OneUlp = U2;
1049         SqXMinX (Serious);
1050         X = BInvrse;
1051         OneUlp = BInvrse * U1;
1052         SqXMinX (Serious);
1053         X = U1;
1054         OneUlp = U1 * U1;
1055         SqXMinX (Serious);
1056         if (J != Zero) Pause();
1057         printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
1058         J = Zero;
1059         X = Two;
1060         Y = Radix;
1061         if ((Radix != One)) do  {
1062                 X = Y;
1063                 Y = Radix * Y;
1064                 } while ( ! ((Y - X >= NoTrials)));
1065         OneUlp = X * U2;
1066         I = 1;
1067         while (I <= NoTrials) {
1068                 X = X + One;
1069                 SqXMinX (Defect);
1070                 if (J > Zero) break;
1071                 I = I + 1;
1072                 }
1073         printf("Test for sqrt monotonicity.\n");
1074         I = - 1;
1075         X = BMinusU2;
1076         Y = Radix;
1077         Z = Radix + Radix * U2;
1078         NotMonot = False;
1079         Monot = False;
1080         while ( ! (NotMonot || Monot)) {
1081                 I = I + 1;
1082                 X = SQRT(X);
1083                 Q = SQRT(Y);
1084                 Z = SQRT(Z);
1085                 if ((X > Q) || (Q > Z)) NotMonot = True;
1086                 else {
1087                         Q = FLOOR(Q + Half);
1088                         if (!(I > 0 || Radix == Q * Q)) Monot = True;
1089                         else if (I > 0) {
1090                         if (I > 1) Monot = True;
1091                         else {
1092                                 Y = Y * BInvrse;
1093                                 X = Y - U1;
1094                                 Z = Y + U1;
1095                                 }
1096                         }
1097                         else {
1098                                 Y = Q;
1099                                 X = Y - U2;
1100                                 Z = Y + U2;
1101                                 }
1102                         }
1103                 }
1104         if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
1105         else {
1106                 BadCond(Defect, "");
1107                 printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
1108                 }
1109         /*=============================================*/
1110         /*SPLIT
1111         }
1112 #include "paranoia.h"
1113 void part5(VOID){
1114 */
1115         Milestone = 80;
1116         /*=============================================*/
1117         MinSqEr = MinSqEr + Half;
1118         MaxSqEr = MaxSqEr - Half;
1119         Y = (SQRT(One + U2) - One) / U2;
1120         SqEr = (Y - One) + U2 / Eight;
1121         if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1122         SqEr = Y + U2 / Eight;
1123         if (SqEr < MinSqEr) MinSqEr = SqEr;
1124         Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
1125         SqEr = Y + U1 / Eight;
1126         if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1127         SqEr = (Y + One) + U1 / Eight;
1128         if (SqEr < MinSqEr) MinSqEr = SqEr;
1129         OneUlp = U2;
1130         X = OneUlp;
1131         for( Indx = 1; Indx <= 3; ++Indx) {
1132                 Y = SQRT((X + U1 + X) + F9);
1133                 Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
1134                 Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
1135                 SqEr = (Y + Half) + Z;
1136                 if (SqEr < MinSqEr) MinSqEr = SqEr;
1137                 SqEr = (Y - Half) + Z;
1138                 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1139                 if (((Indx == 1) || (Indx == 3)))
1140                         X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
1141                 else {
1142                         OneUlp = U1;
1143                         X = - OneUlp;
1144                         }
1145                 }
1146         /*=============================================*/
1147         Milestone = 85;
1148         /*=============================================*/
1149         SqRWrng = False;
1150         Anomaly = False;
1151         RSqrt = Other; /* ~dgh */
1152         if (Radix != One) {
1153                 printf("Testing whether sqrt is rounded or chopped.\n");
1154                 D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
1155         /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
1156                 X = D / Radix;
1157                 Y = D / A1;
1158                 if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
1159                         Anomaly = True;
1160                         }
1161                 else {
1162                         X = Zero;
1163                         Z2 = X;
1164                         Y = One;
1165                         Y2 = Y;
1166                         Z1 = Radix - One;
1167                         FourD = Four * D;
1168                         do  {
1169                                 if (Y2 > Z2) {
1170                                         Q = Radix;
1171                                         Y1 = Y;
1172                                         do  {
1173                                                 X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
1174                                                 Q = Y1;
1175                                                 Y1 = X1;
1176                                                 } while ( ! (X1 <= Zero));
1177                                         if (Q <= One) {
1178                                                 Z2 = Y2;
1179                                                 Z = Y;
1180                                                 }
1181                                         }
1182                                 Y = Y + Two;
1183                                 X = X + Eight;
1184                                 Y2 = Y2 + X;
1185                                 if (Y2 >= FourD) Y2 = Y2 - FourD;
1186                                 } while ( ! (Y >= D));
1187                         X8 = FourD - Z2;
1188                         Q = (X8 + Z * Z) / FourD;
1189                         X8 = X8 / Eight;
1190                         if (Q != FLOOR(Q)) Anomaly = True;
1191                         else {
1192                                 Break = False;
1193                                 do  {
1194                                         X = Z1 * Z;
1195                                         X = X - FLOOR(X / Radix) * Radix;
1196                                         if (X == One)
1197                                                 Break = True;
1198                                         else
1199                                                 Z1 = Z1 - One;
1200                                         } while ( ! (Break || (Z1 <= Zero)));
1201                                 if ((Z1 <= Zero) && (! Break)) Anomaly = True;
1202                                 else {
1203                                         if (Z1 > RadixD2) Z1 = Z1 - Radix;
1204                                         do  {
1205                                                 NewD();
1206                                                 } while ( ! (U2 * D >= F9));
1207                                         if (D * Radix - D != W - D) Anomaly = True;
1208                                         else {
1209                                                 Z2 = D;
1210                                                 I = 0;
1211                                                 Y = D + (One + Z) * Half;
1212                                                 X = D + Z + Q;
1213                                                 SR3750();
1214                                                 Y = D + (One - Z) * Half + D;
1215                                                 X = D - Z + D;
1216                                                 X = X + Q + X;
1217                                                 SR3750();
1218                                                 NewD();
1219                                                 if (D - Z2 != W - Z2) Anomaly = True;
1220                                                 else {
1221                                                         Y = (D - Z2) + (Z2 + (One - Z) * Half);
1222                                                         X = (D - Z2) + (Z2 - Z + Q);
1223                                                         SR3750();
1224                                                         Y = (One + Z) * Half;
1225                                                         X = Q;
1226                                                         SR3750();
1227                                                         if (I == 0) Anomaly = True;
1228                                                         }
1229                                                 }
1230                                         }
1231                                 }
1232                         }
1233                 if ((I == 0) || Anomaly) {
1234                         BadCond(Failure, "Anomalous arithmetic with Integer < ");
1235                         printf("Radix^Precision = %.7e\n", W);
1236                         printf(" fails test whether sqrt rounds or chops.\n");
1237                         SqRWrng = True;
1238                         }
1239                 }
1240         if (! Anomaly) {
1241                 if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
1242                         RSqrt = Rounded;
1243                         printf("Square root appears to be correctly rounded.\n");
1244                         }
1245                 else  {
1246                         if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
1247                                 || (MinSqEr + Radix < Half)) SqRWrng = True;
1248                         else {
1249                                 RSqrt = Chopped;
1250                                 printf("Square root appears to be chopped.\n");
1251                                 }
1252                         }
1253                 }
1254         if (SqRWrng) {
1255                 printf("Square root is neither chopped nor correctly rounded.\n");
1256                 printf("Observed errors run from %.7e ", MinSqEr - Half);
1257                 printf("to %.7e ulps.\n", Half + MaxSqEr);
1258                 TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
1259                         "sqrt gets too many last digits wrong");
1260                 }
1261         /*=============================================*/
1262         Milestone = 90;
1263         /*=============================================*/
1264         Pause();
1265         printf("Testing powers Z^i for small Integers Z and i.\n");
1266         N = 0;
1267         /* ... test powers of zero. */
1268         I = 0;
1269         Z = -Zero;
1270         M = 3;
1271         Break = False;
1272         do  {
1273                 X = One;
1274                 SR3980();
1275                 if (I <= 10) {
1276                         I = 1023;
1277                         SR3980();
1278                         }
1279                 if (Z == MinusOne) Break = True;
1280                 else {
1281                         Z = MinusOne;
1282                         /* .. if(-1)^N is invalid, replace MinusOne by One. */
1283                         I = - 4;
1284                         }
1285                 } while ( ! Break);
1286         PrintIfNPositive();
1287         N1 = N;
1288         N = 0;
1289         Z = A1;
1290         M = (int)FLOOR(Two * LOG(W) / LOG(A1));
1291         Break = False;
1292         do  {
1293                 X = Z;
1294                 I = 1;
1295                 SR3980();
1296                 if (Z == AInvrse) Break = True;
1297                 else Z = AInvrse;
1298                 } while ( ! (Break));
1299         /*=============================================*/
1300                 Milestone = 100;
1301         /*=============================================*/
1302         /*  Powers of Radix have been tested, */
1303         /*         next try a few primes     */
1304         M = NoTrials;
1305         Z = Three;
1306         do  {
1307                 X = Z;
1308                 I = 1;
1309                 SR3980();
1310                 do  {
1311                         Z = Z + Two;
1312                         } while ( Three * FLOOR(Z / Three) == Z );
1313                 } while ( Z < Eight * Three );
1314         if (N > 0) {
1315                 printf("Errors like this may invalidate financial calculations\n");
1316                 printf("\tinvolving interest rates.\n");
1317                 }
1318         PrintIfNPositive();
1319         N += N1;
1320         if (N == 0) printf("... no discrepancies found.\n");
1321         if (N > 0) Pause();
1322         else printf("\n");
1323         /*=============================================*/
1324         /*SPLIT
1325         }
1326 #include "paranoia.h"
1327 void part6(VOID){
1328 */
1329         Milestone = 110;
1330         /*=============================================*/
1331         printf("Seeking Underflow thresholds UfThold and E0.\n");
1332         D = U1;
1333         if (Precision != FLOOR(Precision)) {
1334                 D = BInvrse;
1335                 X = Precision;
1336                 do  {
1337                         D = D * BInvrse;
1338                         X = X - One;
1339                         } while ( X > Zero);
1340                 }
1341         Y = One;
1342         Z = D;
1343         /* ... D is power of 1/Radix < 1. */
1344         do  {
1345                 C = Y;
1346                 Y = Z;
1347                 Z = Y * Y;
1348                 } while ((Y > Z) && (Z + Z > Z));
1349         Y = C;
1350         Z = Y * D;
1351         do  {
1352                 C = Y;
1353                 Y = Z;
1354                 Z = Y * D;
1355                 } while ((Y > Z) && (Z + Z > Z));
1356         if (Radix < Two) HInvrse = Two;
1357         else HInvrse = Radix;
1358         H = One / HInvrse;
1359         /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
1360         CInvrse = One / C;
1361         E0 = C;
1362         Z = E0 * H;
1363         /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
1364         do  {
1365                 Y = E0;
1366                 E0 = Z;
1367                 Z = E0 * H;
1368                 } while ((E0 > Z) && (Z + Z > Z));
1369         UfThold = E0;
1370         E1 = Zero;
1371         Q = Zero;
1372         E9 = U2;
1373         S = One + E9;
1374         D = C * S;
1375         if (D <= C) {
1376                 E9 = Radix * U2;
1377                 S = One + E9;
1378                 D = C * S;
1379                 if (D <= C) {
1380                         BadCond(Failure, "multiplication gets too many last digits wrong.\n");
1381                         Underflow = E0;
1382                         Y1 = Zero;
1383                         PseudoZero = Z;
1384                         Pause();
1385                         }
1386                 }
1387         else {
1388                 Underflow = D;
1389                 PseudoZero = Underflow * H;
1390                 UfThold = Zero;
1391                 do  {
1392                         Y1 = Underflow;
1393                         Underflow = PseudoZero;
1394                         if (E1 + E1 <= E1) {
1395                                 Y2 = Underflow * HInvrse;
1396                                 E1 = FABS(Y1 - Y2);
1397                                 Q = Y1;
1398                                 if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
1399                                 }
1400                         PseudoZero = PseudoZero * H;
1401                         } while ((Underflow > PseudoZero)
1402                                 && (PseudoZero + PseudoZero > PseudoZero));
1403                 }
1404         /* Comment line 4530 .. 4560 */
1405         if (PseudoZero != Zero) {
1406                 printf("\n");
1407                 Z = PseudoZero;
1408         /* ... Test PseudoZero for "phoney- zero" violates */
1409         /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
1410                    ... */
1411                 if (PseudoZero <= Zero) {
1412                         BadCond(Failure, "Positive expressions can underflow to an\n");
1413                         printf("allegedly negative value\n");
1414                         printf("PseudoZero that prints out as: %g .\n", PseudoZero);
1415                         X = - PseudoZero;
1416                         if (X <= Zero) {
1417                                 printf("But -PseudoZero, which should be\n");
1418                                 printf("positive, isn't; it prints out as  %g .\n", X);
1419                                 }
1420                         }
1421                 else {
1422                         BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
1423                         printf("value PseudoZero that prints out as %g .\n", PseudoZero);
1424                         }
1425                 TstPtUf();
1426                 }
1427         /*=============================================*/
1428         Milestone = 120;
1429         /*=============================================*/
1430         if (CInvrse * Y > CInvrse * Y1) {
1431                 S = H * S;
1432                 E0 = Underflow;
1433                 }
1434         if (! ((E1 == Zero) || (E1 == E0))) {
1435                 BadCond(Defect, "");
1436                 if (E1 < E0) {
1437                         printf("Products underflow at a higher");
1438                         printf(" threshold than differences.\n");
1439                         if (PseudoZero == Zero)
1440                         E0 = E1;
1441                         }
1442                 else {
1443                         printf("Difference underflows at a higher");
1444                         printf(" threshold than products.\n");
1445                         }
1446                 }
1447         printf("Smallest strictly positive number found is E0 = %g .\n", E0);
1448         Z = E0;
1449         TstPtUf();
1450         Underflow = E0;
1451         if (N == 1) Underflow = Y;
1452         I = 4;
1453         if (E1 == Zero) I = 3;
1454         if (UfThold == Zero) I = I - 2;
1455         UfNGrad = True;
1456         switch (I)  {
1457                 case    1:
1458                 UfThold = Underflow;
1459                 if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
1460                         UfThold = Y;
1461                         BadCond(Failure, "Either accuracy deteriorates as numbers\n");
1462                         printf("approach a threshold = %.17e\n", UfThold);;
1463                         printf(" coming down from %.17e\n", C);
1464                         printf(" or else multiplication gets too many last digits wrong.\n");
1465                         }
1466                 Pause();
1467                 break;
1468
1469                 case    2:
1470                 BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
1471                 printf("Q == Y while denying that |Q - Y| == 0; these values\n");
1472                 printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
1473                 printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
1474                 UfThold = Q;
1475                 break;
1476
1477                 case    3:
1478                 X = X;
1479                 break;
1480
1481                 case    4:
1482                 if ((Q == UfThold) && (E1 == E0)
1483                         && (FABS( UfThold - E1 / E9) <= E1)) {
1484                         UfNGrad = False;
1485                         printf("Underflow is gradual; it incurs Absolute Error =\n");
1486                         printf("(roundoff in UfThold) < E0.\n");
1487                         Y = E0 * CInvrse;
1488                         Y = Y * (OneAndHalf + U2);
1489                         X = CInvrse * (One + U2);
1490                         Y = Y / X;
1491                         IEEE = (Y == E0);
1492                         }
1493                 }
1494         if (UfNGrad) {
1495                 printf("\n");
1496                 sigsave = sigfpe;
1497                 if (setjmp(ovfl_buf)) {
1498                         printf("Underflow / UfThold failed!\n");
1499                         R = H + H;
1500                         }
1501                 else R = SQRT(Underflow / UfThold);
1502                 sigsave = 0;
1503                 if (R <= H) {
1504                         Z = R * UfThold;
1505                         X = Z * (One + R * H * (One + H));
1506                         }
1507                 else {
1508                         Z = UfThold;
1509                         X = Z * (One + H * H * (One + H));
1510                         }
1511                 if (! ((X == Z) || (X - Z != Zero))) {
1512                         BadCond(Flaw, "");
1513                         printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
1514                         Z9 = X - Z;
1515                         printf("yet X - Z yields %.17e .\n", Z9);
1516                         printf("    Should this NOT signal Underflow, ");
1517                         printf("this is a SERIOUS DEFECT\nthat causes ");
1518                         printf("confusion when innocent statements like\n");;
1519                         printf("    if (X == Z)  ...  else");
1520                         printf("  ... (f(X) - f(Z)) / (X - Z) ...\n");
1521                         printf("encounter Division by Zero although actually\n");
1522                         sigsave = sigfpe;
1523                         if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
1524                         else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
1525                         sigsave = 0;
1526                         }
1527                 }
1528         printf("The Underflow threshold is %.17e, %s\n", UfThold,
1529                    " below which");
1530         printf("calculation may suffer larger Relative error than ");
1531         printf("merely roundoff.\n");
1532         Y2 = U1 * U1;
1533         Y = Y2 * Y2;
1534         Y2 = Y * U1;
1535         if (Y2 <= UfThold) {
1536                 if (Y > E0) {
1537                         BadCond(Defect, "");
1538                         I = 5;
1539                         }
1540                 else {
1541                         BadCond(Serious, "");
1542                         I = 4;
1543                         }
1544                 printf("Range is too narrow; U1^%d Underflows.\n", I);
1545                 }
1546         /*=============================================*/
1547         /*SPLIT
1548         }
1549 #include "paranoia.h"
1550 void part7(VOID){
1551 */
1552         Milestone = 130;
1553         /*=============================================*/
1554         Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
1555         Y2 = Y + Y;
1556         printf("Since underflow occurs below the threshold\n");
1557         printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
1558         printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n",
1559                 HInvrse, Y2);
1560         printf("actually calculating yields:");
1561         if (setjmp(ovfl_buf)) {
1562                 sigsave = 0;
1563                 BadCond(Serious, "trap on underflow.\n");
1564                 }
1565         else {
1566                 sigsave = sigfpe;
1567                 V9 = POW(HInvrse, Y2);
1568                 sigsave = 0;
1569                 printf(" %.17e .\n", V9);
1570                 if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
1571                         BadCond(Serious, "this is not between 0 and underflow\n");
1572                 printf("   threshold = %.17e .\n", UfThold);
1573                 }
1574                 else if (! (V9 > UfThold * (One + E9)))
1575                         printf("This computed value is O.K.\n");
1576                 else {
1577                         BadCond(Defect, "this is not between 0 and underflow\n");
1578                         printf("   threshold = %.17e .\n", UfThold);
1579                         }
1580                 }
1581         /*=============================================*/
1582         Milestone = 140;
1583         /*=============================================*/
1584         printf("\n");
1585         /* ...calculate Exp2 == exp(2) == 7.389056099... */
1586         X = Zero;
1587         I = 2;
1588         Y = Two * Three;
1589         Q = Zero;
1590         N = 0;
1591         do  {
1592                 Z = X;
1593                 I = I + 1;
1594                 Y = Y / (I + I);
1595                 R = Y + Q;
1596                 X = Z + R;
1597                 Q = (Z - X) + R;
1598                 } while(X > Z);
1599         Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
1600         X = Z * Z;
1601         Exp2 = X * X;
1602         X = F9;
1603         Y = X - U1;
1604         printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
1605                 Exp2);
1606         for(I = 1;;) {
1607                 Z = X - BInvrse;
1608                 Z = (X + One) / (Z - (One - BInvrse));
1609                 Q = POW(X, Z) - Exp2;
1610                 if (FABS(Q) > TwoForty * U2) {
1611                         N = 1;
1612                         V9 = (X - BInvrse) - (One - BInvrse);
1613                         BadCond(Defect, "Calculated");
1614                         printf(" %.17e for\n", POW(X,Z));
1615                         printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
1616                         printf("\tdiffers from correct value by %.17e .\n", Q);
1617                         printf("\tThis much error may spoil financial\n");
1618                         printf("\tcalculations involving tiny interest rates.\n");
1619                         break;
1620                         }
1621                 else {
1622                         Z = (Y - X) * Two + Y;
1623                         X = Y;
1624                         Y = Z;
1625                         Z = One + (X - F9)*(X - F9);
1626                         if (Z > One && I < NoTrials) I++;
1627                         else  {
1628                                 if (X > One) {
1629                                         if (N == 0)
1630                                            printf("Accuracy seems adequate.\n");
1631                                         break;
1632                                         }
1633                                 else {
1634                                         X = One + U2;
1635                                         Y = U2 + U2;
1636                                         Y += X;
1637                                         I = 1;
1638                                         }
1639                                 }
1640                         }
1641                 }
1642         /*=============================================*/
1643         Milestone = 150;
1644         /*=============================================*/
1645         printf("Testing powers Z^Q at four nearly extreme values.\n");
1646         N = 0;
1647         Z = A1;
1648         Q = FLOOR(Half - LOG(C) / LOG(A1));
1649         Break = False;
1650         do  {
1651                 X = CInvrse;
1652                 Y = POW(Z, Q);
1653                 IsYeqX();
1654                 Q = - Q;
1655                 X = C;
1656                 Y = POW(Z, Q);
1657                 IsYeqX();
1658                 if (Z < One) Break = True;
1659                 else Z = AInvrse;
1660                 } while ( ! (Break));
1661         PrintIfNPositive();
1662         if (N == 0) printf(" ... no discrepancies found.\n");
1663         printf("\n");
1664
1665         /*=============================================*/
1666         Milestone = 160;
1667         /*=============================================*/
1668         Pause();
1669         printf("Searching for Overflow threshold:\n");
1670         printf("This may generate an error.\n");
1671         Y = - CInvrse;
1672         V9 = HInvrse * Y;
1673         sigsave = sigfpe;
1674         if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
1675         do {
1676                 V = Y;
1677                 Y = V9;
1678                 V9 = HInvrse * Y;
1679                 } while(V9 < Y);
1680         I = 1;
1681 overflow:
1682         sigsave = 0;
1683         Z = V9;
1684         printf("Can `Z = -Y' overflow?\n");
1685         printf("Trying it on Y = %.17e .\n", Y);
1686         V9 = - Y;
1687         V0 = V9;
1688         if (V - Y == V + V0) printf("Seems O.K.\n");
1689         else {
1690                 printf("finds a ");
1691                 BadCond(Flaw, "-(-Y) differs from Y.\n");
1692                 }
1693         if (Z != Y) {
1694                 BadCond(Serious, "");
1695                 printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
1696                 }
1697         if (I) {
1698                 Y = V * (HInvrse * U2 - HInvrse);
1699                 Z = Y + ((One - HInvrse) * U2) * V;
1700                 if (Z < V0) Y = Z;
1701                 if (Y < V0) V = Y;
1702                 if (V0 - V < V0) V = V0;
1703                 }
1704         else {
1705                 V = Y * (HInvrse * U2 - HInvrse);
1706                 V = V + ((One - HInvrse) * U2) * Y;
1707                 }
1708         printf("Overflow threshold is V  = %.17e .\n", V);
1709         if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
1710         else printf("There is no saturation value because \
1711 the system traps on overflow.\n");
1712         V9 = V * One;
1713         printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
1714         V9 = V / One;
1715         printf("                           nor for V / 1 = %.17e .\n", V9);
1716         printf("Any overflow signal separating this * from the one\n");
1717         printf("above is a DEFECT.\n");
1718         /*=============================================*/
1719         Milestone = 170;
1720         /*=============================================*/
1721         if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
1722                 BadCond(Failure, "Comparisons involving ");
1723                 printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
1724                         V, V0, UfThold);
1725                 }
1726         /*=============================================*/
1727         Milestone = 175;
1728         /*=============================================*/
1729         printf("\n");
1730         for(Indx = 1; Indx <= 3; ++Indx) {
1731                 switch (Indx)  {
1732                         case 1: Z = UfThold; break;
1733                         case 2: Z = E0; break;
1734                         case 3: Z = PseudoZero; break;
1735                         }
1736                 if (Z != Zero) {
1737                         V9 = SQRT(Z);
1738                         Y = V9 * V9;
1739                         if (Y / (One - Radix * E9) < Z
1740                            || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
1741                                 if (V9 > U1) BadCond(Serious, "");
1742                                 else BadCond(Defect, "");
1743                                 printf("Comparison alleges that what prints as Z = %.17e\n", Z);
1744                                 printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
1745                                 }
1746                         }
1747                 }
1748         /*=============================================*/
1749         Milestone = 180;
1750         /*=============================================*/
1751         for(Indx = 1; Indx <= 2; ++Indx) {
1752                 if (Indx == 1) Z = V;
1753                 else Z = V0;
1754                 V9 = SQRT(Z);
1755                 X = (One - Radix * E9) * V9;
1756                 V9 = V9 * X;
1757                 if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
1758                         Y = V9;
1759                         if (X < W) BadCond(Serious, "");
1760                         else BadCond(Defect, "");
1761                         printf("Comparison alleges that Z = %17e\n", Z);
1762                         printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
1763                         }
1764                 }
1765         /*=============================================*/
1766         /*SPLIT
1767         }
1768 #include "paranoia.h"
1769 int part8(VOID){
1770 */
1771         Milestone = 190;
1772         /*=============================================*/
1773         Pause();
1774         X = UfThold * V;
1775         Y = Radix * Radix;
1776         if (X*Y < One || X > Y) {
1777                 if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
1778                 else BadCond(Flaw, "");
1779
1780                 printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
1781                         X, "is too far from 1.\n");
1782                 }
1783         /*=============================================*/
1784         Milestone = 200;
1785         /*=============================================*/
1786         for (Indx = 1; Indx <= 5; ++Indx)  {
1787                 X = F9;
1788                 switch (Indx)  {
1789                         case 2: X = One + U2; break;
1790                         case 3: X = V; break;
1791                         case 4: X = UfThold; break;
1792                         case 5: X = Radix;
1793                         }
1794                 Y = X;
1795                 sigsave = sigfpe;
1796                 if (setjmp(ovfl_buf))
1797                         printf("  X / X  traps when X = %g\n", X);
1798                 else {
1799                         V9 = (Y / X - Half) - Half;
1800                         if (V9 == Zero) continue;
1801                         if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
1802                         else BadCond(Serious, "");
1803                         printf("  X / X differs from 1 when X = %.17e\n", X);
1804                         printf("  instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
1805                         }
1806                 sigsave = 0;
1807                 }
1808         /*=============================================*/
1809         Milestone = 210;
1810         /*=============================================*/
1811         MyZero = Zero;
1812         printf("\n");
1813         printf("What message and/or values does Division by Zero produce?\n") ;
1814 #ifndef NOPAUSE
1815         printf("This can interupt your program.  You can ");
1816         printf("skip this part if you wish.\n");
1817         printf("Do you wish to compute 1 / 0? ");
1818         fflush(stdout);
1819         read (KEYBOARD, ch, 8);
1820         if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1821 #endif
1822                 sigsave = sigfpe;
1823                 printf("    Trying to compute 1 / 0 produces ...");
1824                 if (!setjmp(ovfl_buf)) printf("  %.7e .\n", One / MyZero);
1825                 sigsave = 0;
1826 #ifndef NOPAUSE
1827                 }
1828         else printf("O.K.\n");
1829         printf("\nDo you wish to compute 0 / 0? ");
1830         fflush(stdout);
1831         read (KEYBOARD, ch, 80);
1832         if ((ch[0] == 'Y') || (ch[0] == 'y')) {
1833 #endif
1834                 sigsave = sigfpe;
1835                 printf("\n    Trying to compute 0 / 0 produces ...");
1836                 if (!setjmp(ovfl_buf)) printf("  %.7e .\n", Zero / MyZero);
1837                 sigsave = 0;
1838 #ifndef NOPAUSE
1839                 }
1840         else printf("O.K.\n");
1841 #endif
1842         /*=============================================*/
1843         Milestone = 220;
1844         /*=============================================*/
1845         Pause();
1846         printf("\n");
1847         {
1848                 static char *msg[] = {
1849                         "FAILUREs  encountered =",
1850                         "SERIOUS DEFECTs  discovered =",
1851                         "DEFECTs  discovered =",
1852                         "FLAWs  discovered =" };
1853                 int i;
1854                 for(i = 0; i < 4; i++) if (ErrCnt[i])
1855                         printf("The number of  %-29s %d.\n",
1856                                 msg[i], ErrCnt[i]);
1857                 }
1858         printf("\n");
1859         if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
1860                         + ErrCnt[Flaw]) > 0) {
1861                 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
1862                         Defect] == 0) && (ErrCnt[Flaw] > 0)) {
1863                         printf("The arithmetic diagnosed seems ");
1864                         printf("Satisfactory though flawed.\n");
1865                         }
1866                 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
1867                         && ( ErrCnt[Defect] > 0)) {
1868                         printf("The arithmetic diagnosed may be Acceptable\n");
1869                         printf("despite inconvenient Defects.\n");
1870                         }
1871                 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
1872                         printf("The arithmetic diagnosed has ");
1873                         printf("unacceptable Serious Defects.\n");
1874                         }
1875                 if (ErrCnt[Failure] > 0) {
1876                         printf("Potentially fatal FAILURE may have spoiled this");
1877                         printf(" program's subsequent diagnoses.\n");
1878                         }
1879                 }
1880         else {
1881                 printf("No failures, defects nor flaws have been discovered.\n");
1882                 if (! ((RMult == Rounded) && (RDiv == Rounded)
1883                         && (RAddSub == Rounded) && (RSqrt == Rounded)))
1884                         printf("The arithmetic diagnosed seems Satisfactory.\n");
1885                 else {
1886                         if (StickyBit >= One &&
1887                                 (Radix - Two) * (Radix - Nine - One) == Zero) {
1888                                 printf("Rounding appears to conform to ");
1889                                 printf("the proposed IEEE standard P");
1890                                 if ((Radix == Two) &&
1891                                          ((Precision - Four * Three * Two) *
1892                                           ( Precision - TwentySeven -
1893                                            TwentySeven + One) == Zero))
1894                                         printf("754");
1895                                 else printf("854");
1896                                 if (IEEE) printf(".\n");
1897                                 else {
1898                                         printf(",\nexcept for possibly Double Rounding");
1899                                         printf(" during Gradual Underflow.\n");
1900                                         }
1901                                 }
1902                         printf("The arithmetic diagnosed appears to be Excellent!\n");
1903                         }
1904                 }
1905         if (fpecount)
1906                 printf("\nA total of %d floating point exceptions were registered.\n",
1907                         fpecount);
1908         printf("END OF TEST.\n");
1909         return 0;
1910         }
1911
1912 /*SPLIT subs.c
1913 #include "paranoia.h"
1914 */
1915
1916  FLOAT
1917 Sign (FP X)
1918 #ifdef KR_headers
1919 FLOAT X;
1920 #endif
1921 { return X >= 0. ? 1.0 : -1.0; }
1922
1923  void
1924 Pause(VOID)
1925 {
1926 #ifndef NOPAUSE
1927         char ch[8];
1928
1929         printf("\nTo continue, press RETURN");
1930         fflush(stdout);
1931         read(KEYBOARD, ch, 8);
1932 #endif
1933         printf("\nDiagnosis resumes after milestone Number %d", Milestone);
1934         printf("          Page: %d\n\n", PageNo);
1935         ++Milestone;
1936         ++PageNo;
1937         }
1938
1939  void
1940 TstCond (INT K, INT Valid, CHARP T)
1941 #ifdef KR_headers
1942 int K, Valid;
1943 char *T;
1944 #endif
1945 { if (! Valid) { BadCond(K,T); printf(".\n"); } }
1946
1947  void
1948 BadCond(INT K, CHARP T)
1949 #ifdef KR_headers
1950 int K;
1951 char *T;
1952 #endif
1953 {
1954         static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
1955
1956         ErrCnt [K] = ErrCnt [K] + 1;
1957         printf("%s:  %s", msg[K], T);
1958         }
1959
1960
1961  FLOAT
1962 Random(VOID)
1963 /*  Random computes
1964      X = (Random1 + Random9)^5
1965      Random1 = X - FLOOR(X) + 0.000005 * X;
1966    and returns the new value of Random1
1967 */
1968 {
1969         FLOAT X, Y;
1970
1971         X = Random1 + Random9;
1972         Y = X * X;
1973         Y = Y * Y;
1974         X = X * Y;
1975         Y = X - FLOOR(X);
1976         Random1 = Y + X * 0.000005;
1977         return(Random1);
1978         }
1979
1980  void
1981 SqXMinX (INT ErrKind)
1982 #ifdef KR_headers
1983 int ErrKind;
1984 #endif
1985 {
1986         FLOAT XA, XB;
1987
1988         XB = X * BInvrse;
1989         XA = X - XB;
1990         SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
1991         if (SqEr != Zero) {
1992                 if (SqEr < MinSqEr) MinSqEr = SqEr;
1993                 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
1994                 J = J + 1.0;
1995                 BadCond(ErrKind, "\n");
1996                 printf("sqrt( %.17e) - %.17e  = %.17e\n", X * X, X, OneUlp * SqEr);
1997                 printf("\tinstead of correct value 0 .\n");
1998                 }
1999         }
2000
2001  void
2002 NewD(VOID)
2003 {
2004         X = Z1 * Q;
2005         X = FLOOR(Half - X / Radix) * Radix + X;
2006         Q = (Q - X * Z) / Radix + X * X * (D / Radix);
2007         Z = Z - Two * X * D;
2008         if (Z <= Zero) {
2009                 Z = - Z;
2010                 Z1 = - Z1;
2011                 }
2012         D = Radix * D;
2013         }
2014
2015  void
2016 SR3750(VOID)
2017 {
2018         if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
2019                 I = I + 1;
2020                 X2 = SQRT(X * D);
2021                 Y2 = (X2 - Z2) - (Y - Z2);
2022                 X2 = X8 / (Y - Half);
2023                 X2 = X2 - Half * X2 * X2;
2024                 SqEr = (Y2 + Half) + (Half - X2);
2025                 if (SqEr < MinSqEr) MinSqEr = SqEr;
2026                 SqEr = Y2 - X2;
2027                 if (SqEr > MaxSqEr) MaxSqEr = SqEr;
2028                 }
2029         }
2030
2031  void
2032 IsYeqX(VOID)
2033 {
2034         if (Y != X) {
2035                 if (N <= 0) {
2036                         if (Z == Zero && Q <= Zero)
2037                                 printf("WARNING:  computing\n");
2038                         else BadCond(Defect, "computing\n");
2039                         printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
2040                         printf("\tyielded %.17e;\n", Y);
2041                         printf("\twhich compared unequal to correct %.17e ;\n",
2042                                 X);
2043                         printf("\t\tthey differ by %.17e .\n", Y - X);
2044                         }
2045                 N = N + 1; /* ... count discrepancies. */
2046                 }
2047         }
2048
2049  void
2050 SR3980(VOID)
2051 {
2052         do {
2053                 Q = (FLOAT) I;
2054                 Y = POW(Z, Q);
2055                 IsYeqX();
2056                 if (++I > M) break;
2057                 X = Z * X;
2058                 } while ( X < W );
2059         }
2060
2061  void
2062 PrintIfNPositive(VOID)
2063 {
2064         if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
2065         }
2066
2067  void
2068 TstPtUf(VOID)
2069 {
2070         N = 0;
2071         if (Z != Zero) {
2072                 printf("Since comparison denies Z = 0, evaluating ");
2073                 printf("(Z + Z) / Z should be safe.\n");
2074                 sigsave = sigfpe;
2075                 if (setjmp(ovfl_buf)) goto very_serious;
2076                 Q9 = (Z + Z) / Z;
2077                 printf("What the machine gets for (Z + Z) / Z is  %.17e .\n",
2078                         Q9);
2079                 if (FABS(Q9 - Two) < Radix * U2) {
2080                         printf("This is O.K., provided Over/Underflow");
2081                         printf(" has NOT just been signaled.\n");
2082                         }
2083                 else {
2084                         if ((Q9 < One) || (Q9 > Two)) {
2085 very_serious:
2086                                 N = 1;
2087                                 ErrCnt [Serious] = ErrCnt [Serious] + 1;
2088                                 printf("This is a VERY SERIOUS DEFECT!\n");
2089                                 }
2090                         else {
2091                                 N = 1;
2092                                 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2093                                 printf("This is a DEFECT!\n");
2094                                 }
2095                         }
2096                 sigsave = 0;
2097                 V9 = Z * One;
2098                 Random1 = V9;
2099                 V9 = One * Z;
2100                 Random2 = V9;
2101                 V9 = Z / One;
2102                 if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
2103                         if (N > 0) Pause();
2104                         }
2105                 else {
2106                         N = 1;
2107                         BadCond(Defect, "What prints as Z = ");
2108                         printf("%.17e\n\tcompares different from  ", Z);
2109                         if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
2110                         if (! ((Z == Random2)
2111                                 || (Random2 == Random1)))
2112                                 printf("1 * Z == %g\n", Random2);
2113                         if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
2114                         if (Random2 != Random1) {
2115                                 ErrCnt [Defect] = ErrCnt [Defect] + 1;
2116                                 BadCond(Defect, "Multiplication does not commute!\n");
2117                                 printf("\tComparison alleges that 1 * Z = %.17e\n",
2118                                         Random2);
2119                                 printf("\tdiffers from Z * 1 = %.17e\n", Random1);
2120                                 }
2121                         Pause();
2122                         }
2123                 }
2124         }
2125
2126  void
2127 notify(CHARP s)
2128 #ifdef KR_headers
2129  char *s;
2130 #endif
2131 {
2132         printf("%s test appears to be inconsistent...\n", s);
2133         printf("   PLEASE NOTIFY KARPINKSI!\n");
2134         }
2135
2136 /*SPLIT msgs.c
2137 #include "paranoia.h"
2138 */
2139
2140  void
2141 msglist(CHARPP s)
2142 #ifdef KR_headers
2143 char **s;
2144 #endif
2145 { while(*s) printf("%s\n", *s++); }
2146
2147  void
2148 Instructions(VOID)
2149 {
2150   static char *instr[] = {
2151         "Lest this program stop prematurely, i.e. before displaying\n",
2152         "    `END OF TEST',\n",
2153         "try to persuade the computer NOT to terminate execution when an",
2154         "error like Over/Underflow or Division by Zero occurs, but rather",
2155         "to persevere with a surrogate value after, perhaps, displaying some",
2156         "warning.  If persuasion avails naught, don't despair but run this",
2157         "program anyway to see how many milestones it passes, and then",
2158         "amend it to make further progress.\n",
2159         "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
2160         0};
2161
2162         msglist(instr);
2163         }
2164
2165  void
2166 Heading(VOID)
2167 {
2168   static char *head[] = {
2169         "Users are invited to help debug and augment this program so it will",
2170         "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
2171         "Please send suggestions and interesting results to",
2172         "\tRichard Karpinski",
2173         "\tComputer Center U-76",
2174         "\tUniversity of California",
2175         "\tSan Francisco, CA 94143-0704, USA\n",
2176         "In doing so, please include the following information:",
2177 #ifdef Single
2178         "\tPrecision:\tsingle;",
2179 #else
2180         "\tPrecision:\tdouble;",
2181 #endif
2182         "\tVersion:\t10 February 1989;",
2183         "\tComputer:\n",
2184         "\tCompiler:\n",
2185         "\tOptimization level:\n",
2186         "\tOther relevant compiler options:",
2187         0};
2188
2189         msglist(head);
2190         }
2191
2192  void
2193 Characteristics(VOID)
2194 {
2195         static char *chars[] = {
2196          "Running this program should reveal these characteristics:",
2197         "     Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
2198         "     Precision = number of significant digits carried.",
2199         "     U2 = Radix/Radix^Precision = One Ulp",
2200         "\t(OneUlpnit in the Last Place) of 1.000xxx .",
2201         "     U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
2202         "     Adequacy of guard digits for Mult., Div. and Subt.",
2203         "     Whether arithmetic is chopped, correctly rounded, or something else",
2204         "\tfor Mult., Div., Add/Subt. and Sqrt.",
2205         "     Whether a Sticky Bit used correctly for rounding.",
2206         "     UnderflowThreshold = an underflow threshold.",
2207         "     E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
2208         "     V = an overflow threshold, roughly.",
2209         "     V0  tells, roughly, whether  Infinity  is represented.",
2210         "     Comparisions are checked for consistency with subtraction",
2211         "\tand for contamination with pseudo-zeros.",
2212         "     Sqrt is tested.  Y^X is not tested.",
2213         "     Extra-precise subexpressions are revealed but NOT YET tested.",
2214         "     Decimal-Binary conversion is NOT YET tested for accuracy.",
2215         0};
2216
2217         msglist(chars);
2218         }
2219
2220  void
2221 History(VOID)
2222 { /* History */
2223  /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
2224         with further massaging by David M. Gay. */
2225
2226   static char *hist[] = {
2227         "The program attempts to discriminate among",
2228         "   FLAWs, like lack of a sticky bit,",
2229         "   Serious DEFECTs, like lack of a guard digit, and",
2230         "   FAILUREs, like 2+2 == 5 .",
2231         "Failures may confound subsequent diagnoses.\n",
2232         "The diagnostic capabilities of this program go beyond an earlier",
2233         "program called `MACHAR', which can be found at the end of the",
2234         "book  `Software Manual for the Elementary Functions' (1980) by",
2235         "W. J. Cody and W. Waite. Although both programs try to discover",
2236         "the Radix, Precision and range (over/underflow thresholds)",
2237         "of the arithmetic, this program tries to cope with a wider variety",
2238         "of pathologies, and to say how well the arithmetic is implemented.",
2239         "\nThe program is based upon a conventional radix representation for",
2240         "floating-point numbers, but also allows logarithmic encoding",
2241         "as used by certain early WANG machines.\n",
2242         "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
2243         "see source comments for more history.",
2244         0};
2245
2246         msglist(hist);
2247         }