nsz Git - libfirm/blob - ir/be/test/extreme/paranoia.c

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