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