git.sur5r.net Git - cc65/blob - test/ref/paranoia.c

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