test/math/eparanoi.c

   1 /* paranoia.c arithmetic tester\r
   2  *\r
   3  * This is an implementation of the PARANOIA program.  It substitutes\r
   4  * subroutine calls for ALL floating point arithmetic operations.\r
   5  * This permits you to substitute your own experimental versions of\r
   6  * arithmetic routines.  It also defeats compiler optimizations,\r
   7  * so for native arithmetic you can be pretty sure you are testing\r
   8  * the arithmetic and not the compiler.\r
   9  *\r
  10  * This version of PARANOIA omits the display of division by zero.\r
  11  * It also omits the test for extra precise subexpressions, since\r
  12  * they cannot occur in this context.  Otherwise it includes all the\r
  13  * tests of the 27 Jan 86 distribution, plus a few additional tests.\r
  14  * Commentary has been reduced to a minimum in order to make the program\r
  15  * smaller.\r
  16  *\r
  17  * The original PARANOIA program, written by W. Kahan, C version\r
  18  * by Thos Sumner and David Gay, can be downloaded free from the\r
  19  * Internet NETLIB.  An MSDOS disk can be obtained for $15 from:\r
  20  *   Richard Karpinski\r
  21  *   6521 Raymond Street\r
  22  *   Oakland, CA 94609\r
  23  *\r
  24  * Steve Moshier, 28 Oct 88\r
  25  * last rev: 23 May 92\r
  26  */\r
  27 \r
  28 #define DEBUG 0\r
  29 \r
  30 /* To use the native arithmetic of the computer, define NATIVE\r
  31  * to be 1.  To use your own supplied arithmetic routines, NATIVE is 0.\r
  32  */\r
  33 #define NATIVE 0\r
  34 \r
  35 /* gcc real.c interface */\r
  36 #define L128DOUBLE 0\r
  37 \r
  38 #include <stdio.h>\r
  39 \r
  40 \r
  41 \r
  42 \r
  43 /* Data structure of floating point number.  If NATIVE was\r
  44  * selected above, you can define LDOUBLE 1 to test 80-bit long double\r
  45  * precision or define it 0 to test 64-bit double precision.\r
  46 */\r
  47 #define LDOUBLE 0\r
  48 #if NATIVE\r
  49 \r
  50 #define NE 1\r
  51 #if LDOUBLE\r
  52 #define FSIZE long double\r
  53 #define FLOAT(x) FSIZE x[NE]\r
  54 static FSIZE eone[NE] = {1.0L}; /* The constant 1.0 */\r
  55 #define ZSQRT sqrtl\r
  56 #define ZLOG logl\r
  57 #define ZFLOOR floorl\r
  58 #define ZPOW powl\r
  59 long double sqrtl(), logl(), floorl(), powl();\r
  60 #define FSETUP einit\r
  61 #else /* not LDOUBLE */\r
  62 #define FSIZE double\r
  63 #define FLOAT(x) FSIZE x[NE]\r
  64 static FSIZE eone[NE] = {1.0};  /* The constant 1.0 */\r
  65 #define ZSQRT sqrt\r
  66 #define ZLOG log\r
  67 #define ZFLOOR floor\r
  68 #define ZPOW pow\r
  69 double sqrt(), log(), floor(), pow();\r
  70 /* Coprocessor initialization,\r
  71  * defeat underflow trap or what have you.\r
  72  * This is required mainly on i386 and 68K processors.\r
  73  */\r
  74 #define FSETUP dprec\r
  75 #endif /* double, not LDOUBLE */\r
  76 \r
  77 #else /* not NATIVE */\r
  78 \r
  79 /* Setup for extended double type.\r
  80  * Put NE = 10 for real.c operating with TFmode support (16-byte reals)\r
  81  * Put NE = 6 for real.c operating with XFmode support (10- or 12-byte reals)\r
  82  * The value of NE must agree with that in ehead.h, if ieee.c is used.\r
  83  */\r
  84 #define NE 6\r
  85 #define FSIZE unsigned short\r
  86 #define FLOAT(x) unsigned short x[NE]\r
  87 extern unsigned short eone[];\r
  88 #define FSETUP einit\r
  89 \r
  90 /* default for FSETUP */\r
  91 /*\r
  92 einit()\r
  93 {}\r
  94 */\r
  95 \r
  96 error(s)\r
  97 char *s;\r
  98 {\r
  99 printf( "error: %s\n", s );\r
 100 }\r
 101 \r
 102 #endif  /* not NATIVE */\r
 103 \r
 104 \r
 105 \r
 106 #if L128DOUBLE\r
 107 /* real.c interface */\r
 108 \r
 109 #undef FSETUP\r
 110 #define FSETUP efsetup\r
 111 \r
 112 FLOAT(enone);\r
 113 \r
 114 #define ONE enone\r
 115 \r
 116 /* Use emov to convert from widest type to widest type, ... */\r
 117 /*\r
 118 #define ENTOE emov\r
 119 #define ETOEN emov\r
 120 */\r
 121 \r
 122 /*                 ... else choose e24toe, e53toe, etc. */\r
 123 #define ENTOE e64toe\r
 124 #define ETOEN etoe64\r
 125 #define NNBITS 64\r
 126 \r
 127 #define NIBITS ((NE-1)*16)\r
 128 extern int rndprc;\r
 129 \r
 130 efsetup()\r
 131 {\r
 132 rndprc = NNBITS;\r
 133 ETOEN(eone, enone);\r
 134 }\r
 135 \r
 136 add(a,b,c)\r
 137 FLOAT(a);\r
 138 FLOAT(b);\r
 139 FLOAT(c);\r
 140 {\r
 141 unsigned short aa[10], bb[10], cc[10];\r
 142 \r
 143 ENTOE(a,aa);\r
 144 ENTOE(b,bb);\r
 145 eadd(aa,bb,cc);\r
 146 ETOEN(cc,c);\r
 147 }\r
 148 \r
 149 sub(a,b,c)\r
 150 FLOAT(a);\r
 151 FLOAT(b);\r
 152 FLOAT(c);\r
 153 {\r
 154 unsigned short aa[10], bb[10], cc[10];\r
 155 \r
 156 ENTOE(a,aa);\r
 157 ENTOE(b,bb);\r
 158 esub(aa,bb,cc);\r
 159 ETOEN(cc,c);\r
 160 }\r
 161 \r
 162 mul(a,b,c)\r
 163 FLOAT(a);\r
 164 FLOAT(b);\r
 165 FLOAT(c);\r
 166 {\r
 167 unsigned short aa[10], bb[10], cc[10];\r
 168 \r
 169 ENTOE(a,aa);\r
 170 ENTOE(b,bb);\r
 171 emul(aa,bb,cc);\r
 172 ETOEN(cc,c);\r
 173 }\r
 174 \r
 175 div(a,b,c)\r
 176 FLOAT(a);\r
 177 FLOAT(b);\r
 178 FLOAT(c);\r
 179 {\r
 180 unsigned short aa[10], bb[10], cc[10];\r
 181 \r
 182 ENTOE(a,aa);\r
 183 ENTOE(b,bb);\r
 184 ediv(aa,bb,cc);\r
 185 ETOEN(cc,c);\r
 186 }\r
 187 \r
 188 int cmp(a,b)\r
 189 FLOAT(a);\r
 190 FLOAT(b);\r
 191 {\r
 192 unsigned short aa[10], bb[10];\r
 193 int c;\r
 194 int ecmp();\r
 195 \r
 196 ENTOE(a,aa);\r
 197 ENTOE(b,bb);\r
 198 c = ecmp(aa,bb);\r
 199 return(c);\r
 200 }\r
 201 \r
 202 mov(a,b)\r
 203 FLOAT(a);\r
 204 FLOAT(b);\r
 205 {\r
 206 int i;\r
 207 \r
 208 for( i=0; i<NE; i++ )\r
 209         b[i] = a[i];\r
 210 }\r
 211 \r
 212 \r
 213 neg(a)\r
 214 FLOAT(a);\r
 215 {\r
 216 unsigned short aa[10];\r
 217 \r
 218 ENTOE(a,aa);\r
 219 eneg(aa);\r
 220 ETOEN(aa,a);\r
 221 }\r
 222 \r
 223 clear(a)\r
 224 FLOAT(a);\r
 225 {\r
 226 int i;\r
 227 \r
 228 for( i=0; i<NE; i++ )\r
 229         a[i] = 0;\r
 230 }\r
 231 \r
 232 FABS(a)\r
 233 FLOAT(a);\r
 234 {\r
 235 unsigned short aa[10];\r
 236 \r
 237 ENTOE(a,aa);\r
 238 eabs(aa);\r
 239 ETOEN(aa,a);\r
 240 }\r
 241 \r
 242 FLOOR(a,b)\r
 243 FLOAT(a);\r
 244 FLOAT(b);\r
 245 {\r
 246 unsigned short aa[10], bb[10];\r
 247 \r
 248 ENTOE(a,aa);\r
 249 efloor(aa,bb);\r
 250 ETOEN(bb,b);\r
 251 }\r
 252 \r
 253 LOG(a,b)\r
 254 FLOAT(a);\r
 255 FLOAT(b);\r
 256 {\r
 257 unsigned short aa[10], bb[10];\r
 258 int rndsav;\r
 259 \r
 260 ENTOE(a,aa);\r
 261 rndsav = rndprc;\r
 262 rndprc = NIBITS;\r
 263 elog(aa,bb);\r
 264 rndprc = rndsav;\r
 265 ETOEN(bb,b);\r
 266 }\r
 267 \r
 268 POW(a,b,c)\r
 269 FLOAT(a);\r
 270 FLOAT(b);\r
 271 FLOAT(c);\r
 272 {\r
 273 unsigned short aa[10], bb[10], cc[10];\r
 274 int rndsav;\r
 275 \r
 276 ENTOE(a,aa);\r
 277 ENTOE(b,bb);\r
 278 rndsav = rndprc;\r
 279 rndprc = NIBITS;\r
 280 epow(aa,bb,cc);\r
 281 rndprc = rndsav;\r
 282 ETOEN(cc,c);\r
 283 }\r
 284 \r
 285 SQRT(a,b)\r
 286 FLOAT(a);\r
 287 FLOAT(b);\r
 288 {\r
 289 unsigned short aa[10], bb[10];\r
 290 \r
 291 ENTOE(a,aa);\r
 292 esqrt(aa,bb);\r
 293 ETOEN(bb,b);\r
 294 }\r
 295 \r
 296 FTOL(x,ip,f)\r
 297 FLOAT(x);\r
 298 long *ip;\r
 299 FLOAT(f);\r
 300 {\r
 301 unsigned short xx[10], ff[10];\r
 302 \r
 303 ENTOE(x,xx);\r
 304 eifrac(xx,ip,ff);\r
 305 ETOEN(ff,f);\r
 306 }\r
 307 \r
 308 LTOF(ip,x)\r
 309 long *ip;\r
 310 FLOAT(x);\r
 311 {\r
 312 unsigned short xx[10];\r
 313 ltoe(ip,xx);\r
 314 ETOEN(xx,x);\r
 315 }\r
 316 \r
 317 TOASC(a,b,c)\r
 318 FLOAT(a);\r
 319 int b;\r
 320 char *c;\r
 321 {\r
 322 unsigned short xx[10];\r
 323 \r
 324 ENTOE(a,xx);\r
 325 etoasc(xx,b,c);\r
 326 }\r
 327 \r
 328 #else /* not L128DOUBLE */\r
 329 \r
 330 #define ONE eone\r
 331 \r
 332 /* Note all arguments of operation subroutines are pointers. */\r
 333 /* c = b + a */\r
 334 #define add(a,b,c) eadd(a,b,c)\r
 335 /* c = b - a */\r
 336 #define sub(a,b,c) esub(a,b,c)\r
 337 /* c = b * a */\r
 338 #define mul(a,b,c) emul(a,b,c)\r
 339 /* c = b / a */\r
 340 #define div(a,b,c) ediv(a,b,c)\r
 341 /* 1 if a>b, 0 if a==b, -1 if a<b */\r
 342 #define cmp(a,b) ecmp(a,b)\r
 343 /* b = a */\r
 344 #define mov(a,b) emov(a,b)\r
 345 /* a = -a */\r
 346 #define neg(a) eneg(a)\r
 347 /* a = 0 */\r
 348 #define clear(a) eclear(a)\r
 349 \r
 350 #define FABS(x) eabs(x)\r
 351 #define FLOOR(x,y) efloor(x,y)\r
 352 #define LOG(x,y) elog(x,y)\r
 353 #define POW(x,y,z) epow(x,y,z)\r
 354 #define SQRT(x,y) esqrt(x,y)\r
 355 \r
 356 /* x = &FLOAT input, i = &long integer part, f = &FLOAT fractional part */\r
 357 #define FTOL(x,i,f) eifrac(x,i,f)\r
 358 \r
 359 /* i = &long integer input, x = &FLOAT output */\r
 360 #define LTOF(i,x) ltoe(i,x)\r
 361 \r
 362 /* Convert FLOAT a to decimal ASCII string with b digits */\r
 363 #define TOASC(a,b,c) etoasc(a,b,c)\r
 364 #endif /* not L128DOUBLE */\r
 365 \r
 366 \r
 367 \r
 368 /* The following subroutines are implementations of the above\r
 369  * named functions, using the native or default arithmetic.\r
 370  */\r
 371 #if NATIVE\r
 372 eadd(a,b,c)\r
 373 FSIZE *a, *b, *c;\r
 374 {\r
 375 *c = *b + *a;\r
 376 }\r
 377 \r
 378 esub(a,b,c)\r
 379 FSIZE *a, *b, *c;\r
 380 {\r
 381 *c = *b - *a;\r
 382 }\r
 383 \r
 384 emul(a,b,c)\r
 385 FSIZE *a, *b, *c;\r
 386 {\r
 387 *c = (*b) * (*a);\r
 388 }\r
 389 \r
 390 ediv(a,b,c)\r
 391 FSIZE *a, *b, *c;\r
 392 {\r
 393 *c = (*b) / (*a);\r
 394 }\r
 395 \r
 396 \r
 397 /* Important note: comparison can be done by subracting\r
 398  * or by a compare instruction that may or may not be\r
 399  * equivalent to subtracting.\r
 400  */\r
 401 ecmp(a,b)\r
 402 FSIZE *a, *b;\r
 403 {\r
 404 if( (*a) > (*b) )\r
 405         return( 1 );\r
 406 if( (*a) < (*b) )\r
 407         return( -1 );\r
 408 if( (*a) != (*b) )\r
 409         goto cmpf;\r
 410 if( (*a) == (*b) )\r
 411         return( 0 );\r
 412 cmpf:\r
 413 printf( "Compare fails\n" );\r
 414 return(0);\r
 415 }\r
 416 \r
 417 \r
 418 emov( a, b )\r
 419 FSIZE *a, *b;\r
 420 {\r
 421 *b = *a;\r
 422 }\r
 423 \r
 424 eneg( a )\r
 425 FSIZE *a;\r
 426 {\r
 427 *a = -(*a);\r
 428 }\r
 429 \r
 430 eclear(a)\r
 431 FSIZE *a;\r
 432 {\r
 433 *a = 0.0;\r
 434 }\r
 435 \r
 436 eabs(x)\r
 437 FSIZE *x;\r
 438 {\r
 439 if( (*x) < 0.0 )\r
 440         *x = -(*x);\r
 441 }\r
 442 \r
 443 efloor(x,y)\r
 444 FSIZE *x, *y;\r
 445 {\r
 446 \r
 447 *y = (FSIZE )ZFLOOR( *x );\r
 448 }\r
 449 \r
 450 elog(x,y)\r
 451 FSIZE *x, *y;\r
 452 {\r
 453 \r
 454 *y = (FSIZE )ZLOG( *x );\r
 455 }\r
 456 \r
 457 epow(x,y,z)\r
 458 FSIZE *x, *y, *z;\r
 459 {\r
 460 \r
 461 *z = (FSIZE )ZPOW( *x, *y );\r
 462 }\r
 463 \r
 464 esqrt(x,y)\r
 465 FSIZE *x, *y;\r
 466 {\r
 467 \r
 468 *y = (FSIZE )ZSQRT( *x );\r
 469 }\r
 470 \r
 471 \r
 472 eifrac(x,i,f)\r
 473 FSIZE *x;\r
 474 long *i;\r
 475 FSIZE *f;\r
 476 {\r
 477 FSIZE y;\r
 478 \r
 479 y = (FSIZE )ZFLOOR( *x );\r
 480 if( y < 0.0 )\r
 481         {\r
 482         *f = y - *x;\r
 483         *i = -y;\r
 484         }\r
 485 else\r
 486         {\r
 487         *f = *x - y;\r
 488         *i = y;\r
 489         }\r
 490 }\r
 491 \r
 492 \r
 493 ltoe(i,x)\r
 494 long *i;\r
 495 FSIZE *x;\r
 496 {\r
 497 *x = *i;\r
 498 }\r
 499 \r
 500 \r
 501 etoasc(a,str,n)\r
 502 FSIZE *a;\r
 503 char *str;\r
 504 int n;\r
 505 {\r
 506 double x;\r
 507 \r
 508 x = (double )(*a);\r
 509 sprintf( str, " %.17e ", x );\r
 510 }\r
 511 \r
 512 /* default for FSETUP */\r
 513 einit()\r
 514 {}\r
 515 \r
 516 #endif  /* NATIVE */\r
 517 \r
 518 \r
 519 \r
 520 \r
 521 FLOAT(Radix);\r
 522 FLOAT(BInvrse);\r
 523 FLOAT(RadixD2);\r
 524 FLOAT(BMinusU2);\r
 525 /*Small floating point constants.*/\r
 526 FLOAT(Zero);\r
 527 FLOAT(Half);\r
 528 FLOAT(One);\r
 529 FLOAT(Two);\r
 530 FLOAT(Three);\r
 531 FLOAT(Four);\r
 532 FLOAT(Five);\r
 533 FLOAT(Six);\r
 534 FLOAT(Eight);\r
 535 FLOAT(Nine);\r
 536 FLOAT(Ten);\r
 537 FLOAT(TwentySeven);\r
 538 FLOAT(ThirtyTwo);\r
 539 FLOAT(TwoForty);\r
 540 FLOAT(MinusOne );\r
 541 FLOAT(OneAndHalf);\r
 542 \r
 543 /*Integer constants*/\r
 544 int NoTrials = 20; /*Number of tests for commutativity. */\r
 545 #define False 0\r
 546 #define True 1\r
 547 \r
 548 /* Definitions for declared types \r
 549         Guard == (Yes, No);\r
 550         Rounding == (Chopped, Rounded, Other);\r
 551         Message == packed array [1..40] of char;\r
 552         Class == (Flaw, Defect, Serious, Failure);\r
 553           */\r
 554 #define Yes 1\r
 555 #define No  0\r
 556 #define Chopped 2\r
 557 #define Rounded 1\r
 558 #define Other   0\r
 559 #define Flaw    3\r
 560 #define Defect  2\r
 561 #define Serious 1\r
 562 #define Failure 0\r
 563 \r
 564 typedef int Guard, Rounding, Class;\r
 565 typedef char Message;\r
 566 \r
 567 /* Declarations of Variables */\r
 568 FLOAT(AInvrse);\r
 569 FLOAT(A1);\r
 570 FLOAT(C);\r
 571 FLOAT(CInvrse);\r
 572 FLOAT(D);\r
 573 FLOAT(FourD);\r
 574 FLOAT(E0);\r
 575 FLOAT(E1);\r
 576 FLOAT(Exp2);\r
 577 FLOAT(E3);\r
 578 FLOAT(MinSqEr);\r
 579 FLOAT(SqEr);\r
 580 FLOAT(MaxSqEr);\r
 581 FLOAT(E9);\r
 582 FLOAT(Third);\r
 583 FLOAT(F6);\r
 584 FLOAT(F9);\r
 585 FLOAT(H);\r
 586 FLOAT(HInvrse);\r
 587 FLOAT(StickyBit);\r
 588 FLOAT(J);\r
 589 FLOAT(MyZero);\r
 590 FLOAT(Precision);\r
 591 FLOAT(Q);\r
 592 FLOAT(Q9);\r
 593 FLOAT(R);\r
 594 FLOAT(Random9);\r
 595 FLOAT(T);\r
 596 FLOAT(Underflow);\r
 597 FLOAT(S);\r
 598 FLOAT(OneUlp);\r
 599 FLOAT(UfThold);\r
 600 FLOAT(U1);\r
 601 FLOAT(U2);\r
 602 FLOAT(V);\r
 603 FLOAT(V0);\r
 604 FLOAT(V9);\r
 605 FLOAT(W);\r
 606 FLOAT(X);\r
 607 FLOAT(X1);\r
 608 FLOAT(X2);\r
 609 FLOAT(X8);\r
 610 FLOAT(Random1);\r
 611 FLOAT(Y);\r
 612 FLOAT(YY1);\r
 613 FLOAT(Y2);\r
 614 FLOAT(Random2);\r
 615 FLOAT(Z);\r
 616 FLOAT(PseudoZero);\r
 617 FLOAT(Z1);\r
 618 FLOAT(Z2);\r
 619 FLOAT(Z9);\r
 620 static FLOAT(t);\r
 621 FLOAT(t2);\r
 622 FLOAT(Sqarg);\r
 623 int ErrCnt[4];\r
 624 int fpecount;\r
 625 int Milestone;\r
 626 int PageNo;\r
 627 int I, M, N, N1, stkflg;\r
 628 Guard GMult, GDiv, GAddSub;\r
 629 Rounding RMult, RDiv, RAddSub, RSqrt;\r
 630 int Break, Done, NotMonot, Monot, Anomaly, IEEE;\r
 631 int SqRWrng, UfNGrad;\r
 632 int k, k2;\r
 633 int Indx;\r
 634 char ch[8];\r
 635 \r
 636 long lngint, lng2; /* intermediate for conversion between int and FLOAT */\r
 637 \r
 638 /* Computed constants. */\r
 639 /*U1  gap below 1.0, i.e, 1.0-U1 is next number below 1.0 */\r
 640 /*U2  gap above 1.0, i.e, 1.0+U2 is next number above 1.0 */\r
 641 \r
 642 \r
 643 show( x )\r
 644 short x[];\r
 645 {\r
 646 int i;\r
 647 char s[80];\r
 648 \r
 649 /* Number of 16-bit groups to display */\r
 650 #if NATIVE\r
 651 #if LDOUBLE\r
 652 #define NPRT (sizeof( long double )/2)\r
 653 #else\r
 654 #define NPRT (sizeof( double )/2)\r
 655 #endif\r
 656 #else\r
 657 #define NPRT NE\r
 658 #endif\r
 659 \r
 660 TOASC( x, s, 70 );\r
 661 printf( "%s\n", s );\r
 662 for( i=0; i<NPRT; i++ )\r
 663         printf( "%04x ", x[i] & 0xffff );\r
 664 printf( "\n" );\r
 665 }\r
 666 \r
 667 /* define NOSIGNAL */\r
 668 #ifndef NOSIGNAL\r
 669 #include <signal.h>\r
 670 #endif\r
 671 #include <setjmp.h>\r
 672 jmp_buf ovfl_buf;\r
 673 /*typedef int (*Sig_type)();*/\r
 674 typedef void (*Sig_type)();\r
 675 Sig_type sigsave;\r
 676 \r
 677 /* Floating point exception receiver */\r
 678 void sigfpe()\r
 679 {\r
 680 fpecount++;\r
 681 printf( "\n* * * FLOATING-POINT ERROR * * *\n" );\r
 682 /* reinitialize the floating point unit */\r
 683 FSETUP();\r
 684 fflush(stdout);\r
 685 if( sigsave )\r
 686         {\r
 687 #ifndef NOSIGNAL\r
 688         signal( SIGFPE, sigsave );\r
 689 #endif\r
 690         sigsave = 0;\r
 691         longjmp( ovfl_buf, 1 );\r
 692         }\r
 693 abort();\r
 694 }\r
 695 \r
 696 \r
 697 main()\r
 698 {\r
 699 \r
 700 /* Do coprocessor or other initializations */\r
 701 FSETUP();\r
 702 \r
 703 printf(\r
 704  "This version of paranoia omits test for extra precise subexpressions\n" );\r
 705 printf( "and includes a few additional tests.\n" );\r
 706 \r
 707 clear(Zero);\r
 708 printf( "0 = " );\r
 709 show( Zero );\r
 710 mov( ONE, One);\r
 711 printf( "1 = " );\r
 712 show( One );\r
 713 add( One, One, Two );\r
 714 printf( "1+1 = " );\r
 715 show( Two );\r
 716 add( Two, One, Three );\r
 717 add( Three, One, Four );\r
 718 add( Four, One, Five );\r
 719 add( Five, One, Six );\r
 720 add( Four, Four, Eight );\r
 721 mul( Three, Three, Nine );\r
 722 add( Nine, One, Ten );\r
 723 mul( Nine, Three, TwentySeven );\r
 724 mul( Four, Eight, ThirtyTwo );\r
 725 mul( Four, Five, t );\r
 726 mul( t, Three, t );\r
 727 mul( t, Four, TwoForty );\r
 728 mov( One, MinusOne );\r
 729 neg( MinusOne );\r
 730 div( Two, One, Half );\r
 731 add( One, Half, OneAndHalf );\r
 732 ErrCnt[Failure] = 0;\r
 733 ErrCnt[Serious] = 0;\r
 734 ErrCnt[Defect] = 0;\r
 735 ErrCnt[Flaw] = 0;\r
 736 PageNo = 1;\r
 737 #ifndef NOSIGNAL\r
 738 signal( SIGFPE, sigfpe );\r
 739 #endif\r
 740 printf("Program is now RUNNING tests on small integers:\n");\r
 741 \r
 742 add( Zero, Zero, t );\r
 743 if( cmp( t, Zero ) != 0)\r
 744         {\r
 745         printf( "0+0 != 0\n" );\r
 746         ErrCnt[Failure] += 1;\r
 747         }\r
 748 sub( One, One, t );\r
 749 if( cmp( t, Zero ) != 0 )\r
 750         {\r
 751         printf( "1-1 != 0\n" );\r
 752         ErrCnt[Failure] += 1;\r
 753         }\r
 754 if( cmp( One, Zero ) <= 0 )\r
 755         {\r
 756         printf( "1 <= 0\n" );\r
 757         ErrCnt[Failure] += 1;\r
 758         }\r
 759 add( One, One, t );\r
 760 if( cmp( t, Two ) != 0 )\r
 761         {\r
 762         printf( "1+1 != 2\n" );\r
 763         ErrCnt[Failure] += 1;\r
 764         }\r
 765 mov( Zero, Z );\r
 766 neg( Z );\r
 767 FLOOR( Z, t );\r
 768 if( cmp(t,Zero) != 0 )\r
 769         {\r
 770         ErrCnt[Serious] += 1;\r
 771         printf( "FLOOR(-0) should equal 0, is = " );\r
 772         show( t );\r
 773         }\r
 774 if( cmp(Z, Zero) != 0)\r
 775         {\r
 776         ErrCnt[Failure] += 1;\r
 777         printf("Comparison alleges that -0.0 is Non-zero!\n");\r
 778         }\r
 779 else\r
 780         {\r
 781         div( TwoForty, One, U1 ); /* U1 = 0.001 */\r
 782         mov( One, Radix );\r
 783         TstPtUf();\r
 784         }\r
 785 add( Two, One, t );\r
 786 if( cmp( t, Three ) != 0 )\r
 787         {\r
 788         printf( "2+1 != 3\n" );\r
 789         ErrCnt[Failure] += 1;\r
 790         }\r
 791 add( Three, One, t );\r
 792 if( cmp( t, Four ) != 0 )\r
 793         {\r
 794         printf( "3+1 != 4\n" );\r
 795         ErrCnt[Failure] += 1;\r
 796         }\r
 797 mov( Two, t );\r
 798 neg( t );\r
 799 mul( Two, t, t );\r
 800 add( Four, t, t );\r
 801 if( cmp( t, Zero ) != 0 )\r
 802         {\r
 803         printf( "4+2*(-2) != 0\n" );\r
 804         ErrCnt[Failure] += 1;\r
 805         }\r
 806 sub( Three, Four, t );\r
 807 sub( One, t, t );\r
 808 if( cmp( t, Zero ) != 0 )\r
 809         {\r
 810         printf( "4-3-1 != 0\n" );\r
 811         ErrCnt[Failure] += 1;\r
 812         }\r
 813         sub( One, Zero, t );\r
 814 if( cmp( t, MinusOne ) != 0 )\r
 815         {\r
 816         printf( "-1 != 0-1\n" );\r
 817         ErrCnt[Failure] += 1;\r
 818         }\r
 819 add( One, MinusOne, t );\r
 820 if( cmp( t, Zero ) != 0 )\r
 821         {\r
 822         printf( "1+(-1) != 0\n" );\r
 823         ErrCnt[Failure] += 1;\r
 824         }\r
 825 mov( One, t );\r
 826 FABS( t );\r
 827 add( MinusOne, t, t );\r
 828 if( cmp( t, Zero ) != 0 )\r
 829         {\r
 830         printf( "-1+abs(1) != 0\n" );\r
 831         ErrCnt[Failure] += 1;\r
 832         }\r
 833 mul( MinusOne, MinusOne, t );\r
 834 add( MinusOne, t, t );\r
 835 if( cmp( t, Zero ) != 0 )\r
 836         {\r
 837         printf( "-1+(-1)*(-1) != 0\n" );\r
 838         ErrCnt[Failure] += 1;\r
 839         }\r
 840 add( Half, MinusOne, t );\r
 841 add( Half, t, t );\r
 842 if( cmp( t, Zero ) != 0 )\r
 843         {\r
 844         printf( "1/2 + (-1) + 1/2 != 0\n" );\r
 845         ErrCnt[Failure] += 1;\r
 846         }\r
 847 Milestone = 10;\r
 848 mul( Three, Three, t );\r
 849 if( cmp( t, Nine ) != 0 )\r
 850         {\r
 851         printf( "3*3 != 9\n" );\r
 852         ErrCnt[Failure] += 1;\r
 853         }\r
 854 mul( Nine, Three, t );\r
 855 if( cmp( t, TwentySeven ) != 0 )\r
 856         {\r
 857         printf( "3*9 != 27\n" );\r
 858         ErrCnt[Failure] += 1;\r
 859         }\r
 860 add( Four, Four, t );\r
 861 if( cmp( t, Eight ) != 0 )\r
 862         {\r
 863         printf( "4+4 != 8\n" );\r
 864         ErrCnt[Failure] += 1;\r
 865         }\r
 866 mul( Eight, Four, t );\r
 867 if( cmp( t, ThirtyTwo ) != 0 )\r
 868         {\r
 869         printf( "8*4 != 32\n" );\r
 870         ErrCnt[Failure] += 1;\r
 871         }\r
 872 sub( TwentySeven, ThirtyTwo, t );\r
 873 sub( Four, t, t );\r
 874 sub( One, t, t );\r
 875 if( cmp( t, Zero ) != 0 )\r
 876         {\r
 877         printf( "32-27-4-1 != 0\n" );\r
 878         ErrCnt[Failure] += 1;\r
 879         }\r
 880 add( Four, One, t );\r
 881 if( cmp( t, Five ) != 0 )\r
 882         {\r
 883         printf( "4+1 != 5\n" );\r
 884         ErrCnt[Failure] += 1;\r
 885         }\r
 886 mul( Four, Five, t );\r
 887 mul( Three, t, t );\r
 888 mul( Four, t, t );\r
 889 if( cmp( t, TwoForty ) != 0 )\r
 890         {\r
 891         printf( "4*5*3*4 != 240\n" );\r
 892         ErrCnt[Failure] += 1;\r
 893         }\r
 894 div( Three, TwoForty, t );\r
 895 mul( Four, Four, t2 );\r
 896 mul( Five, t2, t2 );\r
 897 sub( t2, t2, t );\r
 898 if( cmp( t, Zero ) != 0 )\r
 899         {\r
 900         printf( "240/3 - 4*4*5 != 0\n" );\r
 901         ErrCnt[Failure] += 1;\r
 902         }\r
 903 div( Four, TwoForty, t );\r
 904 mul( Five, Three, t2 );\r
 905 mul( Four, t2, t2 );\r
 906 sub( t2, t, t );\r
 907 if( cmp( t, Zero ) != 0 )\r
 908         {\r
 909         printf( "240/4 - 5*3*4 != 0\n" );\r
 910         ErrCnt[Failure] += 1;\r
 911         }\r
 912 div( Five, TwoForty, t );\r
 913 mul( Four, Three, t2 );\r
 914 mul( Four, t2, t2 );\r
 915 sub( t2, t, t );\r
 916 if( cmp( t, Zero ) != 0 )\r
 917         {\r
 918         printf( "240/5 - 4*3*4 != 0\n" );\r
 919         ErrCnt[Failure] += 1;\r
 920         }\r
 921 if(ErrCnt[Failure] == 0)\r
 922         {\r
 923 printf("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n\n");\r
 924         }\r
 925 printf("Searching for Radix and Precision.\n");\r
 926 mov( One, W );\r
 927 do\r
 928         {\r
 929         add( W, W, W );\r
 930         add( W, One, Y );\r
 931         sub( W, Y, Z );\r
 932         sub( One, Z, Y );\r
 933         mov( Y, t );\r
 934         FABS(t);\r
 935         add( MinusOne, t, t );\r
 936         k = cmp( t, Zero );\r
 937         }\r
 938 while( k < 0 );\r
 939 /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/\r
 940 mov( Zero, Precision );\r
 941 mov( One, Y );\r
 942 do\r
 943         {\r
 944         add( W, Y, Radix );\r
 945         add( Y, Y, Y );\r
 946         sub( W, Radix, Radix );\r
 947         k = cmp( Radix, Zero );\r
 948         }\r
 949 while( k == 0);\r
 950 \r
 951 if( cmp(Radix, Two) < 0 )\r
 952         mov( One, Radix );\r
 953 printf("Radix = " );\r
 954 show( Radix );\r
 955 if( cmp(Radix, One) != 0)\r
 956         {\r
 957         mov( One, W );\r
 958         do\r
 959                 {\r
 960                 add( One, Precision, Precision );\r
 961                 mul( W, Radix, W );\r
 962                 add( W, One, Y );\r
 963                 sub( W, Y, t );\r
 964                 k = cmp( t, One );\r
 965                 }\r
 966         while( k == 0 );\r
 967         }\r
 968 /* now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1 */\r
 969 div( W, One, U1 );\r
 970 mul( Radix, U1, U2 );\r
 971 printf( "Closest relative separation found is U 1 = " );\r
 972 show( U1 );\r
 973 printf( "Recalculating radix and precision." );\r
 974         \r
 975 /*save old values*/\r
 976 mov( Radix, E0 );\r
 977 mov( U1, E1 );\r
 978 mov( U2, E9 );\r
 979 mov( Precision, E3 );\r
 980         \r
 981 div( Three, Four, X );\r
 982 sub( One, X, Third );\r
 983 sub( Third, Half, F6 );\r
 984 add( F6, F6, X );\r
 985 sub( Third, X, X );\r
 986 FABS( X );\r
 987 if( cmp(X, U2) < 0 )\r
 988         mov( U2, X );\r
 989         \r
 990 /*... now X = (unknown no.) ulps of 1+...*/\r
 991 do\r
 992         {\r
 993         mov( X, U2 );\r
 994 /* Y = Half * U2 + ThirtyTwo * U2 * U2; */\r
 995         mul( ThirtyTwo, U2, t );\r
 996         mul( t, U2, t );\r
 997         mul( Half, U2, Y );\r
 998         add( t, Y, Y );\r
 999         add( One, Y, Y );\r
1000         sub( One, Y, X );\r
1001         k = cmp( U2, X );\r
1002         k2 = cmp( X, Zero );\r
1003         }\r
1004 while ( ! ((k <= 0) || (k2 <= 0)));\r
1005         \r
1006 /*... now U2 == 1 ulp of 1 + ... */\r
1007 div( Three, Two, X );\r
1008 sub( Half, X, F6 );\r
1009 add( F6, F6, Third );\r
1010 sub( Half, Third, X );\r
1011 add( F6, X, X );\r
1012 FABS( X );\r
1013 if( cmp(X, U1) < 0 )\r
1014         mov( U1, X );\r
1015         \r
1016 /*... now  X == (unknown no.) ulps of 1 -... */\r
1017 do\r
1018         {\r
1019         mov( X, U1 );\r
1020  /* Y = Half * U1 + ThirtyTwo * U1 * U1;*/\r
1021         mul( ThirtyTwo, U1, t );\r
1022         mul( U1, t, t );\r
1023         mul( Half, U1, Y );\r
1024         add( t, Y, Y );\r
1025         sub( Y, Half, Y );\r
1026         add( Half, Y, X );\r
1027         sub( X, Half, Y );\r
1028         add( Half, Y, X );\r
1029         k = cmp( U1, X );\r
1030         k2 = cmp( X, Zero );\r
1031         } while ( ! ((k <= 0) || (k2 <= 0)));\r
1032 /*... now U1 == 1 ulp of 1 - ... */\r
1033 if( cmp( U1, E1 ) == 0 )\r
1034         printf("confirms closest relative separation U1 .\n");\r
1035 else\r
1036         {\r
1037         printf("gets better closest relative separation U1 = " );\r
1038         show( U1 );\r
1039         }\r
1040 div( U1, One, W );\r
1041 sub( U1, Half, F9 );\r
1042 add( F9, Half, F9 );\r
1043 div( U1, U2, t );\r
1044 div( TwoForty, One, t2 );\r
1045 add( t2, t, t );\r
1046 FLOOR( t, Radix );\r
1047 if( cmp(Radix, E0) == 0 )\r
1048         printf("Radix confirmed.\n");\r
1049 else\r
1050         {\r
1051         printf("MYSTERY: recalculated Radix = " );\r
1052         show( Radix );\r
1053         mov( E0, Radix );\r
1054         }\r
1055 add( Eight, Eight, t );\r
1056 if( cmp( Radix, t ) > 0 )\r
1057         {\r
1058         printf( "Radix is too big: roundoff problems\n" );\r
1059         ErrCnt[Defect] += 1;\r
1060         }\r
1061 k = 1;\r
1062 if( cmp( Radix, Two ) == 0 )\r
1063         k = 0;\r
1064 if( cmp( Radix, Ten ) == 0 )\r
1065         k = 0;\r
1066 if( cmp( Radix, One ) == 0 )\r
1067         k = 0;\r
1068 if( k != 0 )\r
1069         {\r
1070         printf( "Radix is not as good as 2 or 10\n" );\r
1071         ErrCnt[Flaw] += 1;\r
1072         }\r
1073 /*=============================================*/\r
1074 Milestone = 20;\r
1075 /*=============================================*/\r
1076 sub( Half, F9, t );\r
1077 if( cmp( t, Half ) >= 0 )\r
1078         {\r
1079         printf( "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?\n" );\r
1080         ErrCnt[Failure] += 1;\r
1081         }\r
1082 mov( F9, X );\r
1083 I = 1;\r
1084 sub( Half, X, Y );\r
1085 sub( Half, Y, Z );\r
1086 if( (cmp( X, One ) == 0) && (cmp( Z, Zero) != 0) )\r
1087         {\r
1088         printf( "Comparison is fuzzy ,X=1 but X-1/2-1/2 != 0\n" );\r
1089         ErrCnt[Failure] += 1;\r
1090         }\r
1091 add( One, U2, X );\r
1092 I = 0;\r
1093 /*=============================================*/\r
1094 Milestone = 25;\r
1095 /*=============================================*/\r
1096 /*... BMinusU2 = nextafter(Radix, 0) */\r
1097 \r
1098 sub( One, Radix, BMinusU2 );\r
1099 sub( U2, BMinusU2, t );\r
1100 add( One, t, BMinusU2 );\r
1101 /* Purify Integers */\r
1102 if( cmp(Radix,One) != 0 )\r
1103         {\r
1104 /*X = - TwoForty * LOG(U1) / LOG(Radix);*/\r
1105         LOG( U1, X );\r
1106         LOG( Radix, t );\r
1107         div( t, X, X );\r
1108         mul( TwoForty, X, X );\r
1109         neg( X );       \r
1110 \r
1111         add( Half, X, Y );\r
1112         FLOOR( Y, Y );\r
1113         sub( Y, X, t );\r
1114         FABS( t );\r
1115         mul( Four, t, t );\r
1116         if( cmp( t, One ) < 0 )\r
1117                 mov( Y, X );\r
1118         div( TwoForty, X, Precision );\r
1119         add( Half, Precision, Y );\r
1120         FLOOR( Y, Y );\r
1121         sub( Y, Precision, t );\r
1122         FABS( t );\r
1123         mul( TwoForty, t, t );\r
1124         if( cmp( t, Half ) < 0 )\r
1125                 mov( Y, Precision );\r
1126         }\r
1127 FLOOR( Precision, t );\r
1128 if( (cmp( Precision, t ) != 0) || (cmp( Radix, One ) == 0) )\r
1129         {\r
1130         printf("Precision cannot be characterized by an Integer number\n");\r
1131         printf("of significant digits but, by itself, this is a minor flaw.\n");\r
1132         }\r
1133 if( cmp(Radix, One) == 0 ) \r
1134         printf("logarithmic encoding has precision characterized solely by U1.\n");\r
1135 else\r
1136         {\r
1137         printf("The number of significant digits of the Radix is " );\r
1138         show( Precision );\r
1139         }\r
1140 mul( U2, Nine, t );\r
1141 mul( Nine, t, t );\r
1142 mul( TwoForty, t, t );\r
1143 if( cmp( t, One ) >= 0 )\r
1144         {\r
1145         printf( "Precision worse than 5 decimal figures\n" );\r
1146         ErrCnt[Serious] += 1;\r
1147         }\r
1148 /*=============================================*/\r
1149 Milestone = 30;\r
1150 /*=============================================*/\r
1151 /* Test for extra-precise subepressions has been deleted. */\r
1152 Milestone = 35;\r
1153 /*=============================================*/\r
1154 if( cmp(Radix,Two) >= 0 )\r
1155         {\r
1156         mul( Radix, Radix, t );\r
1157         div( t, W, X );\r
1158         add( X, One, Y );\r
1159         sub( X, Y, Z );\r
1160         add( Z, U2, T );\r
1161         sub( Z, T, X );\r
1162         if( cmp( X, U2 ) != 0 )\r
1163                 {\r
1164                 printf( "Subtraction is not normalized X=Y,X+Z != Y+Z!\n" );\r
1165                 ErrCnt[Failure] += 1;\r
1166                 }\r
1167         if( cmp(X,U2) == 0 )\r
1168          printf("Subtraction appears to be normalized, as it should be.");\r
1169         }\r
1170 \r
1171 printf("\nChecking for guard digit in *, /, and -.\n");\r
1172 mul( F9, One, Y );\r
1173 mul( One, F9, Z );\r
1174 sub( Half, F9, X );\r
1175 sub( Half, Y, Y );\r
1176 sub( X, Y, Y );\r
1177 sub( Half, Z, Z );\r
1178 sub( X, Z, Z );\r
1179 add( One, U2, X );\r
1180 mul( X, Radix, T );\r
1181 mul( Radix, X, R );\r
1182 sub( Radix, T, X );\r
1183 mul( Radix, U2, t );\r
1184 sub( t, X, X );\r
1185 sub( Radix, R, T );\r
1186 mul( Radix, U2, t );\r
1187 sub( t, T, T );\r
1188 sub( One, Radix, t );\r
1189 mul( t, X, X );\r
1190 sub( One, Radix, t );\r
1191 mul( t, T, T );\r
1192 \r
1193 k = cmp(X,Zero);\r
1194 k |= cmp(Y,Zero);\r
1195 k |= cmp(Z,Zero);\r
1196 k |= cmp(T,Zero);\r
1197 if( k == 0 )\r
1198         GMult = Yes;\r
1199 else\r
1200         {\r
1201         GMult = No;\r
1202         ErrCnt[Serious] += 1;\r
1203         printf( "* lacks a Guard Digit, so 1*X != X\n" );\r
1204         }\r
1205 mul( Radix, U2, Z );\r
1206 add( One, Z, X );\r
1207 add( X, Z, Y );\r
1208 mul( X, X, t );\r
1209 sub( t, Y, Y );\r
1210 FABS( Y );\r
1211 sub( U2, Y, Y );\r
1212 sub( U2, One, X );\r
1213 sub( U2, X, Z );\r
1214 mul( X, X, t );\r
1215 sub( t, Z, Z );\r
1216 FABS( Z );\r
1217 sub( U1, Z, Z );\r
1218 if( (cmp(Y,Zero) > 0) || (cmp(Z,Zero) > 0) )\r
1219         {\r
1220         ErrCnt[Failure] += 1;\r
1221         printf( "* gets too many final digits wrong.\n" );\r
1222         }\r
1223 sub( U2, One, Y );\r
1224 add( One, U2, X );\r
1225 div( Y, One, Z );\r
1226 sub( X, Z, Y );\r
1227 div( Three, One, X );\r
1228 div( Nine, Three, Z );\r
1229 sub( Z, X, X );\r
1230 div( TwentySeven, Nine, T );\r
1231 sub( T, Z, Z );\r
1232 k = cmp( X, Zero );\r
1233 k |= cmp( Y, Zero );\r
1234 k |= cmp( Z, Zero );\r
1235 if( k )\r
1236         {\r
1237         ErrCnt[Defect] += 1;\r
1238 printf( "Division lacks a Guard Digit, so error can exceed 1 ulp\n" );\r
1239 printf( "or  1/3  and  3/9  and  9/27 may disagree\n" );\r
1240         }\r
1241 div( One, F9, Y );\r
1242 sub( Half, F9, X );\r
1243 sub( Half, Y, Y );\r
1244 sub( X, Y, Y );\r
1245 add( One, U2, X );\r
1246 div( One, X, T );\r
1247 sub( X, T, X );\r
1248 k = cmp( X, Zero );\r
1249 k |= cmp( Y, Zero );\r
1250 k |= cmp( Z, Zero );\r
1251 if( k == 0 )\r
1252         GDiv = Yes;\r
1253 else\r
1254         {\r
1255         GDiv = No;\r
1256         ErrCnt[Serious] += 1;\r
1257         printf( "Division lacks a Guard Digit, so X/1 != X\n" );\r
1258         }\r
1259 add( One, U2, X );\r
1260 div( X, One, X );\r
1261 sub( Half, X, Y );\r
1262 sub( Half, Y, Y );\r
1263 if( cmp(Y,Zero) >= 0 )\r
1264         {\r
1265         ErrCnt[Serious] += 1;\r
1266         printf( "Computed value of 1/1.000..1 >= 1\n" );\r
1267         }\r
1268 sub( U2, One, X );\r
1269 mul( Radix, U2, Y );\r
1270 add( One, Y, Y );\r
1271 mul( X, Radix, Z );\r
1272 mul( Y, Radix, T );\r
1273 div( Radix, Z, R );\r
1274 div( Radix, T, StickyBit );\r
1275 sub( X, R, X );\r
1276 sub( Y, StickyBit, Y );\r
1277 k = cmp( X, Zero );\r
1278 k |= cmp( Y, Zero );\r
1279 if( k )\r
1280         {\r
1281         ErrCnt[Failure] += 1;\r
1282         printf( "* and/or / gets too many last digits wrong\n" );\r
1283         }\r
1284 sub( U1, One, Y );\r
1285 sub( F9, One, X );\r
1286 sub( Y, One, Y );\r
1287 sub( U2, Radix, T );\r
1288 sub( BMinusU2, Radix, Z );\r
1289 sub( T, Radix, T );\r
1290 k = cmp( X, U1 );\r
1291 k |= cmp( Y, U1 );\r
1292 k |= cmp( Z, U2 );\r
1293 k |= cmp( T, U2 );\r
1294 if( k == 0 )\r
1295         GAddSub = Yes;\r
1296 else\r
1297         {\r
1298         GAddSub = No;\r
1299         ErrCnt[Serious] += 1;\r
1300         printf( "- lacks Guard Digit, so cancellation is obscured\n" );\r
1301         }\r
1302 sub( One, F9, t );\r
1303 if( (cmp(F9,One) != 0) && (cmp(t,Zero) >= 0) )\r
1304         {\r
1305         ErrCnt[Serious] += 1;\r
1306         printf("comparison alleges  (1-U1) < 1  although\n");\r
1307         printf("  subtration yields  (1-U1) - 1 = 0 , thereby vitiating\n");\r
1308         printf("  such precautions against division by zero as\n");\r
1309         printf("  ...  if (X == 1.0) {.....} else {.../(X-1.0)...}\n");\r
1310         }\r
1311 if (GMult == Yes && GDiv == Yes && GAddSub == Yes)\r
1312         printf(" *, /, and - appear to have guard digits, as they should.\n");\r
1313 /*=============================================*/\r
1314 Milestone = 40;\r
1315 /*=============================================*/\r
1316 printf("Checking rounding on multiply, divide and add/subtract.\n");\r
1317 RMult = Other;\r
1318 RDiv = Other;\r
1319 RAddSub = Other;\r
1320 div( Two, Radix, RadixD2 );\r
1321 mov( Two, A1 );\r
1322 Done = False;\r
1323 do\r
1324         {\r
1325         mov( Radix, AInvrse );\r
1326         do\r
1327                 {\r
1328                 mov( AInvrse, X );\r
1329                 div( A1, AInvrse, AInvrse );\r
1330                 FLOOR( AInvrse, t );\r
1331                 k = cmp( t, AInvrse );\r
1332                 }\r
1333         while( ! (k != 0 ) );\r
1334         k = cmp( X, One );\r
1335         k2 = cmp( A1, Three );\r
1336         Done = (k == 0) || (k2 > 0);\r
1337         if(! Done)\r
1338                 add( Nine, One, A1 );\r
1339         }\r
1340 while( ! (Done));\r
1341 if( cmp(X, One) == 0 )\r
1342         mov( Radix, A1 );\r
1343 div( A1, One, AInvrse );\r
1344 mov( A1, X );\r
1345 mov( AInvrse, Y );\r
1346 Done = False;\r
1347 do\r
1348         {\r
1349         mul( X, Y, Z );\r
1350         sub( Half, Z, Z );\r
1351         if( cmp( Z, Half ) != 0 )\r
1352                 {\r
1353                 ErrCnt[Failure] += 1;\r
1354                 printf( "X * (1/X) differs from 1\n" );\r
1355                 }\r
1356         k = cmp( X, Radix );\r
1357         Done = (k == 0);\r
1358         mov( Radix, X );\r
1359         div( X, One, Y );\r
1360         }\r
1361 while( ! (Done));\r
1362 \r
1363 add( One, U2, Y2 );\r
1364 sub( U2, One, YY1 );\r
1365 sub( U2, OneAndHalf, X );\r
1366 add( OneAndHalf, U2, Y );\r
1367 sub( U2, X, Z );\r
1368 mul( Z, Y2, Z );\r
1369 mul( Y, YY1, T );\r
1370 sub( X, Z, Z );\r
1371 sub( X, T, T );\r
1372 mul( X, Y2, X );\r
1373 add( Y, U2, Y );\r
1374 mul( Y, YY1, Y );\r
1375 sub( OneAndHalf, X, X );\r
1376 sub( OneAndHalf, Y, Y );\r
1377 k = cmp( X, Zero );\r
1378 k |= cmp( Y, Zero );\r
1379 k |= cmp( Z, Zero );\r
1380 if( cmp( T, Zero ) > 0 )\r
1381         k = 1;\r
1382 if( k == 0 )\r
1383         {\r
1384         add( OneAndHalf, U2, X );\r
1385         mul( X, Y2, X );\r
1386         sub( U2, OneAndHalf, Y );\r
1387         sub( U2, Y, Y );\r
1388         add( OneAndHalf, U2, Z );\r
1389         add( U2, Z, Z );\r
1390         sub( U2, OneAndHalf, T );\r
1391         mul( T, YY1, T );\r
1392         add( Z, U2, t );\r
1393         sub( t, X, X );\r
1394         mul( Y, YY1, StickyBit );\r
1395         mul( Z, Y2, S );\r
1396         sub( Y, T, T );\r
1397         sub( Y, U2, Y );\r
1398         add( StickyBit, Y, Y );\r
1399 /* Z = S - (Z + U2 + U2); */\r
1400         add( Z, U2, t );\r
1401         add( t, U2, t );\r
1402         sub( t, S, Z );\r
1403         add( Y2, U2, t );\r
1404         mul( t, YY1, StickyBit );\r
1405         mul( Y2, YY1, YY1 );\r
1406         sub( Y2, StickyBit, StickyBit );\r
1407         sub( Half, YY1, YY1 );\r
1408         k = cmp( X, Zero );\r
1409         k |= cmp( Y, Zero );\r
1410         k |= cmp( Z, Zero );\r
1411         k |= cmp( T, Zero );\r
1412         k |= cmp( StickyBit, Zero );\r
1413         k |= cmp( YY1, Half );\r
1414         if( k == 0 )\r
1415                 {\r
1416                 RMult = Rounded;\r
1417                 printf("Multiplication appears to round correctly.\n");\r
1418                 }\r
1419         else\r
1420                 {\r
1421                 add( X, U2, t );\r
1422                 k = cmp( t, Zero );\r
1423                 if( cmp( Y, Zero ) >= 0 )\r
1424                         k |= 1;\r
1425                 add( Z, U2, t );\r
1426                 k |= cmp( t, Zero );\r
1427                 if( cmp( T, Zero ) >= 0 )\r
1428                         k |= 1;\r
1429                 add( StickyBit, U2, t );\r
1430                 k |= cmp( t, Zero );\r
1431                 if( cmp(YY1, Half) >= 0 )\r
1432                         k |= 1;\r
1433                 if( k == 0 )\r
1434                         {\r
1435                         printf("Multiplication appears to chop.\n");\r
1436                         }\r
1437                 else\r
1438                         {\r
1439                 printf("* is neither chopped nor correctly rounded.\n");\r
1440                         }\r
1441                 if( (RMult == Rounded) && (GMult == No) )\r
1442                         printf("Multiplication has inconsistent result");\r
1443                 }\r
1444         }\r
1445 else\r
1446         printf("* is neither chopped nor correctly rounded.\n");\r
1447 \r
1448 /*=============================================*/\r
1449 Milestone = 45;\r
1450 /*=============================================*/\r
1451 add( One, U2, Y2 );\r
1452 sub( U2, One, YY1 );\r
1453 add( OneAndHalf, U2, Z );\r
1454 add( Z, U2, Z );\r
1455 div( Y2, Z, X );\r
1456 sub( U2, OneAndHalf, T );\r
1457 sub( U2, T, T );\r
1458 sub( U2, T, Y );\r
1459 div( YY1, Y, Y );\r
1460 add( Z, U2, Z );\r
1461 div( Y2, Z, Z );\r
1462 sub( OneAndHalf, X, X );\r
1463 sub( T, Y, Y );\r
1464 div( YY1, T, T );\r
1465 add( OneAndHalf, U2, t );\r
1466 sub( t, Z, Z );\r
1467 sub( OneAndHalf, U2, t );\r
1468 add( t, T, T );\r
1469 k = 0;\r
1470 if( cmp( X, Zero ) > 0 )\r
1471         k = 1;\r
1472 if( cmp( Y, Zero ) > 0 )\r
1473         k = 1;\r
1474 if( cmp( Z, Zero ) > 0 )\r
1475         k = 1;\r
1476 if( cmp( T, Zero ) > 0 )\r
1477         k = 1;\r
1478 if( k == 0 )\r
1479         {\r
1480         div( Y2, OneAndHalf, X );\r
1481         sub( U2, OneAndHalf, Y );\r
1482         add( U2, OneAndHalf, Z );\r
1483         sub( Y, X, X );\r
1484         div( YY1, OneAndHalf, T );\r
1485         div( YY1, Y, Y );\r
1486         add( Z, U2, t );\r
1487         sub( t, T, T );\r
1488         sub( Z, Y, Y );\r
1489         div( Y2, Z, Z );\r
1490         add( Y2, U2, YY1 );\r
1491         div( Y2, YY1, YY1 );\r
1492         sub( OneAndHalf, Z, Z );\r
1493         sub( Y2, YY1, Y2 );\r
1494         sub( U1, F9, YY1 );\r
1495         div( F9, YY1, YY1 );\r
1496         k = cmp( X, Zero );\r
1497         k |= cmp( Y, Zero );\r
1498         k |= cmp( Z, Zero );\r
1499         k |= cmp( T, Zero );\r
1500         k |= cmp( Y2, Zero );\r
1501         sub( Half, YY1, t );\r
1502         sub( Half, F9, t2 );\r
1503         k |= cmp( t, t2 );\r
1504         if( k == 0 )\r
1505                 {\r
1506                 RDiv = Rounded;\r
1507                 printf("Division appears to round correctly.\n");\r
1508                 if(GDiv == No)\r
1509                         printf("Division test inconsistent\n");\r
1510                 }\r
1511         else\r
1512                 {\r
1513                 k = 0;\r
1514                 if( cmp( X, Zero ) >= 0 )\r
1515                         k = 1;\r
1516                 if( cmp( Y, Zero ) >= 0 )\r
1517                         k = 1;\r
1518                 if( cmp( Z, Zero ) >= 0 )\r
1519                         k = 1;\r
1520                 if( cmp( T, Zero ) >= 0 )\r
1521                         k = 1;\r
1522                 if( cmp( Y, Zero ) >= 0 )\r
1523                         k = 1;\r
1524                 sub( Half, YY1, t );\r
1525                 sub( Half, F9, t2 );\r
1526                 if( cmp( t, t2 ) >= 0 )\r
1527                         k = 1;\r
1528                 if( k == 0 )\r
1529                         {\r
1530                         RDiv = Chopped;\r
1531                         printf("Division appears to chop.\n");\r
1532                         }\r
1533                 }\r
1534         }\r
1535 if(RDiv == Other)\r
1536         printf("/ is neither chopped nor correctly rounded.\n");\r
1537 div( Radix, One, BInvrse );\r
1538 mul( BInvrse, Radix, t );\r
1539 sub( Half, t, t );\r
1540 if( cmp( t, Half ) != 0 )\r
1541         {\r
1542         ErrCnt[Failure] += 1;\r
1543         printf( "Radix * ( 1 / Radix ) differs from 1\n" );\r
1544         }\r
1545 \r
1546 Milestone = 50;\r
1547 /*=============================================*/\r
1548 add( F9, U1, t );\r
1549 sub( Half, t, t );\r
1550 k = cmp( t, Half );\r
1551 add( BMinusU2, U2, t );\r
1552 sub( One, t, t );\r
1553 sub( One, Radix, t2 );\r
1554 k |= cmp( t, t2 );\r
1555 if( k != 0 )\r
1556         {\r
1557         ErrCnt[Failure] += 1;\r
1558         printf( "Incomplete carry-propagation in Addition\n" );\r
1559         }\r
1560 mul( U1, U1, X );\r
1561 sub( X, One, X );\r
1562 sub( U2, One, Y );\r
1563 mul( U2, Y, Y );\r
1564 add( One, Y, Y );\r
1565 sub( Half, F9, Z );\r
1566 sub( Half, X, X );\r
1567 sub( Z, X, X );\r
1568 sub( One, Y, Y );\r
1569 if( (cmp(X,Zero) == 0) && (cmp(Y,Zero) == 0) )\r
1570         {\r
1571         RAddSub = Chopped;\r
1572         printf("Add/Subtract appears to be chopped.\n");\r
1573         }\r
1574 if(GAddSub == Yes)\r
1575         {\r
1576         add( Half, U2, X );\r
1577         mul( X, U2, X );\r
1578         sub( U2, Half, Y );\r
1579         mul( Y, U2, Y );\r
1580         add( One, X, X );\r
1581         add( One, Y, Y );\r
1582         add( One, U2, t );\r
1583         sub( X, t, X );\r
1584         sub( Y, One, Y );\r
1585         k = cmp(X,Zero);\r
1586         if( k )\r
1587                 printf( "1+U2-[u2(1/2+U2)+1] != 0\n" );\r
1588         k2 = cmp(Y,Zero);\r
1589         if( k2 )\r
1590                 printf( "1-[U2(1/2-U2)+1] != 0\n" );\r
1591         k |= k2;\r
1592         if( k == 0 )\r
1593                 {\r
1594                 add( Half, U2, X );\r
1595                 mul( X, U1, X );\r
1596                 sub( U2, Half, Y );\r
1597                 mul( Y, U1, Y );\r
1598                 sub( X, One, X );\r
1599                 sub( Y, One, Y );\r
1600                 sub( X, F9, X );\r
1601                 sub( Y, One, Y );\r
1602                 k = cmp(X,Zero);\r
1603                 if( k )\r
1604                         printf( "F9-[1-U1(1/2+U2)] != 0\n" );\r
1605                 k2 = cmp(Y,Zero);\r
1606                 if( k2 )\r
1607                         printf( "1-[1-U1(1/2-U2)] != 0\n" );\r
1608                 k |= k2;\r
1609                 if( k == 0 )\r
1610                         {\r
1611                         RAddSub = Rounded;\r
1612                 printf("Addition/Subtraction appears to round correctly.\n");\r
1613                         if(GAddSub == No)\r
1614                                 printf( "Add/Subtract test inconsistent\n");\r
1615                         }\r
1616                 else\r
1617                         {\r
1618                  printf("Addition/Subtraction neither rounds nor chops.\n");\r
1619                         }\r
1620                 }\r
1621         else\r
1622                 printf("Addition/Subtraction neither rounds nor chops.\n");\r
1623         }\r
1624 else\r
1625         printf("Addition/Subtraction neither rounds nor chops.\n");\r
1626 \r
1627 mov( One, S );\r
1628 add( One, Half, X );\r
1629 mul( Half, X, X );\r
1630 add( One, X, X );\r
1631 add( One, U2, Y );\r
1632 mul( Y, Half, Y );\r
1633 sub( Y, X, Z );\r
1634 sub( X, Y, T );\r
1635 add( Z, T, StickyBit );\r
1636 if( cmp(StickyBit, Zero) != 0 )\r
1637         {\r
1638         mov( Zero, S );\r
1639         ErrCnt[Flaw] += 1;\r
1640         printf( "(X - Y) + (Y - X) is non zero!\n" );\r
1641         }\r
1642 mov( Zero, StickyBit );\r
1643 FLOOR( RadixD2, t );\r
1644 k2 = cmp( t, RadixD2 );\r
1645 k = 1;\r
1646 if( (GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)\r
1647         && (RMult == Rounded) && (RDiv == Rounded)\r
1648         && (RAddSub == Rounded) && (k2 == 0) )\r
1649         {\r
1650         printf("Checking for sticky bit.\n");\r
1651         k = 0;\r
1652         add( Half, U1, X );\r
1653         mul( X, U2, X );\r
1654         mul( Half, U2, Y );\r
1655         add( One, Y, Z );\r
1656         add( One, X, T );\r
1657         sub( One, Z, t );\r
1658         sub( One, T, t2 );\r
1659         if( cmp(t,Zero) > 0 )\r
1660                 {\r
1661                 k = 1;\r
1662                 printf( "[1+(1/2)U2]-1 > 0\n" );\r
1663                 }\r
1664         if( cmp(t2,U2) < 0 )\r
1665                 {\r
1666                 k = 1;\r
1667                 printf( "[1+U2(1/2+U1)]-1 < U2\n" );\r
1668                 }\r
1669         add( T, Y, Z );\r
1670         sub( X, Z, Y );\r
1671         sub( T, Z, t );\r
1672         sub( T, Y, t2 );\r
1673         if( cmp(t,U2) < 0 )\r
1674                 {\r
1675                 k = 1;\r
1676                 printf( "[[1+U2(1/2+U1)]+(1/2)U2]-[1+U2(1/2+U1)] < U2\n" );\r
1677                 }\r
1678         if( cmp(t2,Zero) != 0 )\r
1679                 {\r
1680                 k = 1;\r
1681                 printf( "(1/2)U2-[1+U2(1/2+U1)] != 0\n" );\r
1682                 }\r
1683         add( Half, U1, X );\r
1684         mul( X, U1, X );\r
1685         mul( Half, U1, Y );\r
1686         sub( Y, One, Z );\r
1687         sub( X, One, T );\r
1688         sub( One, Z, t );\r
1689         sub( F9, T, t2 );\r
1690         if( cmp(t,Zero) != 0 )\r
1691                 {\r
1692                 k = 1;\r
1693                 printf( "(1-(1/2)U1)-1 != 0\n" );\r
1694                 }\r
1695         if( cmp(t2,Zero) != 0 )\r
1696                 {\r
1697                 k = 1;\r
1698                 printf( "[1-U1(1/2+U1)]-F9 != 0\n" );\r
1699                 }\r
1700         sub( U1, Half, Z );\r
1701         mul( Z, U1, Z );\r
1702         sub( Z, F9, T );\r
1703         sub( Y, F9, Q );\r
1704         sub( F9, T, t );\r
1705         if( cmp( t, Zero ) != 0 )\r
1706                 {\r
1707                 k = 1;\r
1708                 printf( "[F9-U1(1/2-U1)]-F9 != 0\n" );\r
1709                 }\r
1710         sub( U1, F9, t );\r
1711         sub( Q, t, t );\r
1712         if( cmp( t, Zero ) != 0 )\r
1713                 {\r
1714                 k = 1;\r
1715                 printf( "(F9-U1)-(F9-(1/2)U1) != 0\n" );\r
1716                 }\r
1717         add( One, U2, Z );\r
1718         mul( Z, OneAndHalf, Z );\r
1719         add( OneAndHalf, U2, T );\r
1720         sub( Z, T, T );\r
1721         add( U2, T, T );\r
1722         div( Radix, Half, X );\r
1723         add( One, X, X );\r
1724         mul( Radix, U2, Y );\r
1725         add( One, Y, Y );\r
1726         mul( X, Y, Z );\r
1727         if( cmp( T, Zero ) != 0 )\r
1728                 {\r
1729                 k = 1;\r
1730                 printf( "(3/2+U2)-3/2(1+U2)+U2 != 0\n" );\r
1731                 }\r
1732         mul( Radix, U2, t );\r
1733         add( X, t, t );\r
1734         sub( Z, t, t );\r
1735         if( cmp( t, Zero ) != 0 )\r
1736                 {\r
1737                 k = 1;\r
1738         printf( "(1+1/2Radix)+Radix*U2-[1+1/(2Radix)][1+Radix*U2] != 0\n" );\r
1739                 }\r
1740         if( cmp(Radix, Two) != 0 )\r
1741                 {\r
1742                 add( Two, U2, X );\r
1743                 div( Two, X, Y );\r
1744                 sub( One, Y, t );\r
1745                 if( cmp( t, Zero) != 0 )\r
1746                         k = 1;\r
1747                 }\r
1748         }\r
1749 if( k == 0 )\r
1750         {\r
1751         printf("Sticky bit apparently used correctly.\n");\r
1752         mov( One, StickyBit );\r
1753         }\r
1754 else\r
1755         {\r
1756         printf("Sticky bit used incorrectly or not at all.\n");\r
1757         }\r
1758 \r
1759 if( GMult == No || GDiv == No || GAddSub == No ||\r
1760                 RMult == Other || RDiv == Other || RAddSub == Other)\r
1761         {\r
1762         ErrCnt[Flaw] += 1;\r
1763  printf("lack(s) of guard digits or failure(s) to correctly round or chop\n");\r
1764 printf( "(noted above) count as one flaw in the final tally below\n" );\r
1765         }\r
1766 /*=============================================*/\r
1767 Milestone = 60;\r
1768 /*=============================================*/\r
1769 printf("\n");\r
1770 printf("Does Multiplication commute?  ");\r
1771 printf("Testing on %d random pairs.\n", NoTrials);\r
1772 SQRT( Three, Random9 );\r
1773 mov( Third, Random1 );\r
1774 I = 1;\r
1775 do\r
1776         {\r
1777         Random();\r
1778         mov( Random1, X );\r
1779         Random();\r
1780         mov( Random1, Y );\r
1781         mul( Y, X, Z9 );\r
1782         mul( X, Y, Z );\r
1783         sub( Z9, Z, Z9 );\r
1784         I = I + 1;\r
1785         }\r
1786 while ( ! ((I > NoTrials) || (cmp(Z9,Zero) != 0)));\r
1787 if(I == NoTrials)\r
1788         {\r
1789         div( Three, Half, t );\r
1790         add( One, t, Random1 );\r
1791         add( U2, U1, t );\r
1792         add( t, One, Random2 );\r
1793         mul( Random1, Random2, Z );\r
1794         mul( Random2, Random1, Y );\r
1795 /* Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /\r
1796  *                      Three) * ((U2 + U1) + One);\r
1797  */\r
1798         div( Three, Half, t2 );\r
1799         add( One, t2, t2 );\r
1800         add( U2, U1, t );\r
1801         add( t, One, t );\r
1802         mul( t2, t, Z9 );\r
1803         mul( t2, t, t );\r
1804         sub( t, Z9, Z9 );\r
1805         }\r
1806 if(! ((I == NoTrials) || (cmp(Z9,Zero) == 0)))\r
1807         {\r
1808         ErrCnt[Defect] += 1;\r
1809         printf( "X * Y == Y * X trial fails.\n");\r
1810         }\r
1811 else\r
1812         {\r
1813         printf("     No failures found in %d integer pairs.\n", NoTrials);\r
1814         }\r
1815 /*=============================================*/\r
1816 Milestone = 70;\r
1817 /*=============================================*/\r
1818 sqtest();\r
1819 Milestone = 90;\r
1820 pow1test();\r
1821 \r
1822 Milestone = 110;\r
1823 \r
1824 printf("Seeking Underflow thresholds UfThold and E0.\n");\r
1825 mov( U1, D );\r
1826 FLOOR( Precision, t );\r
1827 if( cmp(Precision, t) != 0 )\r
1828         {\r
1829         mov( BInvrse, D );\r
1830         mov( Precision, X );\r
1831         do\r
1832                 {\r
1833                 mul( D, BInvrse, D );\r
1834                 sub( One, X, X );\r
1835                 }\r
1836         while( cmp(X, Zero) > 0 );\r
1837         }\r
1838 mov( One, Y );\r
1839 mov( D, Z );\r
1840 /* ... D is power of 1/Radix < 1. */\r
1841 sigsave = sigfpe;\r
1842 if( setjmp(ovfl_buf) )\r
1843         goto under0;\r
1844 do\r
1845         {\r
1846         mov( Y, C );\r
1847         mov( Z, Y );\r
1848         mul( Y, Y, Z );\r
1849         add( Z, Z, t );\r
1850         }\r
1851 while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) );\r
1852 \r
1853 under0:\r
1854 sigsave = 0;\r
1855 \r
1856 mov( C, Y );\r
1857 mul( Y, D, Z );\r
1858 sigsave = sigfpe;\r
1859 if( setjmp(ovfl_buf) )\r
1860         goto under1;\r
1861 do\r
1862         {\r
1863         mov( Y, C );\r
1864         mov( Z, Y );\r
1865         mul( Y, D, Z );\r
1866         add( Z, Z, t );\r
1867         }\r
1868 while( (cmp(Y,Z) > 0) && (cmp(t,Z) > 0) );\r
1869 \r
1870 under1:\r
1871 sigsave = 0;\r
1872 \r
1873 if( cmp(Radix,Two) < 0 )\r
1874         mov( Two, HInvrse );\r
1875 else\r
1876         mov( Radix, HInvrse );\r
1877 div( HInvrse, One, H );\r
1878 /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */\r
1879 div( C, One, CInvrse );\r
1880 mov( C, E0 );\r
1881 mul( E0, H, Z );\r
1882 /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */\r
1883 sigsave = sigfpe;\r
1884 if( setjmp(ovfl_buf) )\r
1885         goto under2;\r
1886 do\r
1887         {\r
1888         mov( E0, Y );\r
1889         mov( Z, E0 );\r
1890         mul( E0, H, Z );\r
1891         add( Z, Z, t );\r
1892         }\r
1893 while( (cmp(E0,Z) > 0) && (cmp(t,Z) > 0) );\r
1894 \r
1895 under2:\r
1896 sigsave = 0;\r
1897 \r
1898 mov( E0, UfThold );\r
1899 mov( Zero, E1 );\r
1900 mov( Zero, Q );\r
1901 mov( U2, E9 );\r
1902 add( One, E9, S );\r
1903 mul( C, S, D );\r
1904 if( cmp(D,C) <= 0 )\r
1905         {\r
1906         mul( Radix, U2, E9 );\r
1907         add( One, E9, S );\r
1908         mul( C, S, D );\r
1909         if( cmp(D, C) <= 0 )\r
1910                 {\r
1911                 ErrCnt[Failure] += 1;\r
1912                 printf( "multiplication gets too many last digits wrong.\n" );\r
1913                 mov( E0, Underflow );\r
1914                 mov( Zero, YY1 );\r
1915                 mov( Z, PseudoZero );\r
1916                 }\r
1917         }\r
1918 else\r
1919         {\r
1920         mov( D, Underflow );\r
1921         mul( Underflow, H, PseudoZero );\r
1922         mov( Zero, UfThold );\r
1923         do\r
1924                 {\r
1925                 mov( Underflow, YY1 );\r
1926                 mov( PseudoZero, Underflow );\r
1927                 add( E1, E1, t );\r
1928                 if( cmp(t, E1) <= 0)\r
1929                         {\r
1930                         mul( Underflow, HInvrse, Y2 );\r
1931                         sub( Y2, YY1, E1 );\r
1932                         FABS( E1 );\r
1933                         mov( YY1, Q );\r
1934                         if( (cmp( UfThold, Zero ) == 0)\r
1935                                 && (cmp(YY1, Y2) != 0) )\r
1936                                 mov( YY1, UfThold );\r
1937                         }\r
1938                 mul( PseudoZero, H, PseudoZero );\r
1939                 add( PseudoZero, PseudoZero, t );\r
1940                 }\r
1941         while( (cmp(Underflow, PseudoZero) > 0)\r
1942                 && (cmp(t, PseudoZero) > 0) );\r
1943         }\r
1944 /* Comment line 4530 .. 4560 */\r
1945 if( cmp(PseudoZero, Zero) != 0 )\r
1946         {\r
1947         printf("\n");\r
1948         mov(PseudoZero, Z );\r
1949 /* ... Test PseudoZero for "phoney- zero" violates */\r
1950 /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero\r
1951                    ... */\r
1952         if( cmp(PseudoZero, Zero) <= 0 )\r
1953                 {\r
1954                 ErrCnt[Failure] += 1;\r
1955                 printf("Positive expressions can underflow to an\n");\r
1956                 printf("allegedly negative value\n");\r
1957                 printf("PseudoZero that prints out as: " );\r
1958                 show( PseudoZero );\r
1959                 mov( PseudoZero, X );\r
1960                 neg( X );\r
1961                 if( cmp(X, Zero) <= 0 )\r
1962                         {\r
1963                         printf("But -PseudoZero, which should be\n");\r
1964                         printf("positive, isn't; it prints out as " );\r
1965                         show( X );\r
1966                         }\r
1967                 }\r
1968         else\r
1969                 {\r
1970                 ErrCnt[Flaw] += 1;\r
1971                 printf( "Underflow can stick at an allegedly positive\n");\r
1972                 printf("value PseudoZero that prints out as " );\r
1973                 show( PseudoZero );\r
1974                 }\r
1975 /*      TstPtUf();*/\r
1976         }\r
1977 \r
1978 /*=============================================*/\r
1979 Milestone = 120;\r
1980 /*=============================================*/\r
1981 mul( CInvrse, Y, t );\r
1982 mul( CInvrse, YY1, t2 );\r
1983 if( cmp(t,t2) > 0 )\r
1984         {\r
1985         mul( H, S, S );\r
1986         mov( Underflow, E0 );\r
1987         }\r
1988 if(! ((cmp(E1,Zero) == 0) || (cmp(E1,E0) == 0)) )\r
1989         {\r
1990         ErrCnt[Defect] += 1;\r
1991         if( cmp(E1,E0) < 0 )\r
1992                 {\r
1993                 printf("Products underflow at a higher");\r
1994                 printf(" threshold than differences.\n");\r
1995                 if( cmp(PseudoZero,Zero) == 0 ) \r
1996                         mov( E1, E0 );\r
1997                 }\r
1998         else\r
1999                 {\r
2000                 printf("Difference underflows at a higher");\r
2001                 printf(" threshold than products.\n");\r
2002                 }\r
2003         }\r
2004 printf("Smallest strictly positive number found is E0 = " );\r
2005 show( E0 );\r
2006 mov( E0, Z );\r
2007 TstPtUf();\r
2008 mov( E0, Underflow );\r
2009 if(N == 1)\r
2010         mov( Y, Underflow );\r
2011 I = 4;\r
2012 if( cmp(E1,Zero) == 0 )\r
2013         I = 3;\r
2014 if( cmp( UfThold,Zero) == 0 )\r
2015         I = I - 2;\r
2016 UfNGrad = True;\r
2017 switch(I)\r
2018         {\r
2019         case 1:\r
2020         mov( Underflow, UfThold );\r
2021         mul( CInvrse, Q, t );\r
2022         mul( CInvrse, Y, t2 );\r
2023         mul( t2, S, t2 );\r
2024         if( cmp( t, t2 ) != 0 )\r
2025                 {\r
2026                 mov( Y, UfThold );\r
2027                 ErrCnt[Failure] += 1;\r
2028                 printf( "Either accuracy deteriorates as numbers\n");\r
2029                 printf("approach a threshold = " );\r
2030                 show( UfThold );\r
2031                 printf(" coming down from " );\r
2032                 show( C );\r
2033         printf(" or else multiplication gets too many last digits wrong.\n");\r
2034                 }\r
2035         break;\r
2036         \r
2037         case    2:\r
2038         ErrCnt[Failure] += 1;\r
2039         printf( "Underflow confuses Comparison which alleges that\n");\r
2040         printf("Q == Y while denying that |Q - Y| == 0; these values\n");\r
2041         printf("print out as Q = " );\r
2042         show( Q );\r
2043         printf( ", Y = " );\r
2044         show( Y );\r
2045         sub( Y2, Q, t );\r
2046         FABS(t);\r
2047         printf ("|Q - Y| = " );\r
2048         show( t );\r
2049         mov( Q, UfThold );\r
2050         break;\r
2051         \r
2052         case 3:\r
2053         mov( X, X );\r
2054         break;\r
2055         \r
2056         case 4:\r
2057         div( E9, E1, t );\r
2058         sub( t, UfThold, t );\r
2059         FABS(t);\r
2060         if( (cmp(Q,UfThold) == 0) && (cmp(E1,E0) == 0)\r
2061                 && (cmp(t,E1) <= 0) )\r
2062                 {\r
2063                 UfNGrad = False;\r
2064                 printf("Underflow is gradual; it incurs Absolute Error =\n");\r
2065                 printf("(roundoff in UfThold) < E0.\n");\r
2066                 mul( E0, CInvrse, Y );\r
2067                 add( OneAndHalf, U2, t );\r
2068                 mul( Y, t, Y );\r
2069                 add( One, U2, X );\r
2070                 mul( CInvrse, X, X );\r
2071                 div( X, Y, t );\r
2072                 IEEE = (cmp(t,E0) == 0);\r
2073                 if( IEEE == 0 )\r
2074                         {\r
2075                 printf( "((CInvrse E0) (1.5+U2)) / (CInvrse (1+U2)) != E0\n" );\r
2076                         printf( "CInvrse = " );\r
2077                         show( CInvrse );\r
2078                         printf( "E0 = " );\r
2079                         show( E0 );\r
2080                         printf( "U2 = " );\r
2081                         show( U2 );\r
2082                         printf( "X = " );\r
2083                         show(X);\r
2084                         printf( "Y = " );\r
2085                         show(Y);\r
2086                         printf( "Y/X = " );\r
2087                         show(t);\r
2088                         }\r
2089                 }\r
2090         }\r
2091 if(UfNGrad)\r
2092         {\r
2093         printf("\n");\r
2094         div( UfThold, Underflow, R );\r
2095         SQRT( R, R );\r
2096         if( cmp(R,H) <= 0)\r
2097                 {\r
2098                 mul( R, UfThold, Z );\r
2099 /* X = Z * (One + R * H * (One + H));*/\r
2100                 add( One, H, X );\r
2101                 mul( H, X, X );\r
2102                 mul( R, X, X );\r
2103                 add( One, X, X );\r
2104                 mul( Z, X, X );\r
2105                 }\r
2106         else\r
2107                 {\r
2108                 mov( UfThold, Z );\r
2109 /*X = Z * (One + H * H * (One + H));*/\r
2110                 add( One, H, X );\r
2111                 mul( H, X, X );\r
2112                 mul( H, X, X );\r
2113                 add( One, X, X );\r
2114                 mul( Z, X, X );\r
2115                 }\r
2116         sub( Z, X, t );\r
2117 /*      if(! ((cmp(X,Z) == 0) || (cmp(t,Zero) != 0)) )*/\r
2118         if( (cmp(X,Z) != 0) && (cmp(t,Zero) == 0) )\r
2119                 {\r
2120 /*              ErrCnt[Flaw] += 1;*/\r
2121                 ErrCnt[Serious] += 1;\r
2122                 printf("X = " );\r
2123                 show( X );\r
2124                 printf( "\tis not equal to Z = " );\r
2125                 show( Z );\r
2126 /*              sub( Z, X, Z9 );*/\r
2127                 printf("yet X - Z yields " );\r
2128                 show( t );\r
2129                 printf("which compares equal to " );\r
2130                 show( Zero );\r
2131                 printf("    Should this NOT signal Underflow, ");\r
2132                 printf("this is a SERIOUS DEFECT\nthat causes ");\r
2133                 printf("confusion when innocent statements like\n");;\r
2134                 printf("    if (X == Z)  ...  else");\r
2135                 printf("  ... (f(X) - f(Z)) / (X - Z) ...\n");\r
2136                 printf("encounter Division by Zero although actually\n");\r
2137                 printf("X / Z = 1 + " );\r
2138                 div( Z, X, t );\r
2139                 sub( Half, t, t );\r
2140                 sub( Half, t, t );\r
2141                 show(t);\r
2142                 }\r
2143         }\r
2144 printf("The Underflow threshold is " );\r
2145 show( UfThold );\r
2146 printf( "below which calculation may suffer larger Relative error than" );\r
2147 printf( " merely roundoff.\n");\r
2148 mul( U1, U1, Y2 );\r
2149 mul( Y2, Y2, Y );\r
2150 mul( Y, U1, Y2 );\r
2151 if( cmp( Y2,UfThold) <= 0 )\r
2152         {\r
2153         if( cmp(Y,E0) > 0 )\r
2154                 {\r
2155                 ErrCnt[Defect] += 1;\r
2156                 I = 5;\r
2157                 }\r
2158         else\r
2159                 {\r
2160                 ErrCnt[Serious] += 1;\r
2161                 I = 4;\r
2162                 }\r
2163         printf("Range is too narrow; U1^%d Underflows.\n", I);\r
2164         }\r
2165 Milestone = 130;\r
2166 \r
2167 /*Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;*/\r
2168 LOG( UfThold, Y );\r
2169 LOG( HInvrse, t );\r
2170 div( t, Y, Y );\r
2171 mul( TwoForty, Y, Y );\r
2172 sub( Y, Half, Y );\r
2173 FLOOR( Y, Y );\r
2174 div( TwoForty, Y, Y );\r
2175 neg(Y);\r
2176 sub( One, Y, Y2 ); /* ***** changed from Y2 = Y + Y */\r
2177 printf("Since underflow occurs below the threshold\n");\r
2178 printf("UfThold = " ); \r
2179 show( HInvrse );\r
2180 printf( "\tto the power  " );\r
2181 show( Y );\r
2182 printf( "only underflow should afflict the expression " );\r
2183 show( HInvrse );\r
2184 printf( "\tto the power  " );\r
2185 show( Y2 );\r
2186 POW( HInvrse, Y2, V9 );\r
2187 printf("Actually calculating yields: " );\r
2188 show( V9 );\r
2189 add( Radix, Radix, t );\r
2190 add( t, E9, t );\r
2191 mul( t, UfThold, t );\r
2192 if( (cmp(V9,Zero) < 0) || (cmp(V9,t) > 0) )\r
2193         {\r
2194         ErrCnt[Serious] += 1;\r
2195         printf( "this is not between 0 and underflow\n");\r
2196         printf("   threshold = " );\r
2197         show( UfThold );\r
2198         }\r
2199 else\r
2200         {\r
2201         add( One, E9, t );\r
2202         mul( UfThold, t, t );\r
2203         if( cmp(V9,t) <= 0 )\r
2204                 printf("This computed value is O.K.\n");\r
2205         else\r
2206                 {\r
2207                 ErrCnt[Defect] += 1;\r
2208                 printf( "this is not between 0 and underflow\n");\r
2209                 printf("   threshold = " );\r
2210                 show( UfThold );\r
2211                 }\r
2212         }\r
2213 \r
2214 Milestone = 140;\r
2215 \r
2216 pow2test();\r
2217         \r
2218 /*=============================================*/\r
2219 Milestone = 160;\r
2220 /*=============================================*/\r
2221 Pause();\r
2222 printf("Searching for Overflow threshold:\n");\r
2223 printf("This may generate an error.\n");\r
2224 sigsave = sigfpe;\r
2225 I = 0;\r
2226 mov( CInvrse, Y ); /* a large power of 2 */\r
2227 neg(Y);\r
2228 mul( HInvrse, Y, V9 ); /* HInvrse = 2 */\r
2229 if (setjmp(ovfl_buf))\r
2230         goto overflow;\r
2231 do\r
2232         {\r
2233         mov( Y, V );\r
2234         mov( V9, Y );\r
2235         mul( HInvrse, Y, V9 );\r
2236         }\r
2237 while( cmp(V9,Y) < 0 ); /* V9 = 2 * Y */\r
2238 I = 1;\r
2239 \r
2240 overflow:\r
2241 \r
2242 show( HInvrse );\r
2243 printf( "\ttimes " );\r
2244 show( Y );\r
2245 printf( "\tequals " );\r
2246 show( V9 );\r
2247 \r
2248 mov( V9, Z );\r
2249 printf("Can `Z = -Y' overflow?\n");\r
2250 printf("Trying it on Y = " );\r
2251 show(Y);\r
2252 mov( Y, V9 );\r
2253 neg( V9 );\r
2254 mov( V9, V0 );\r
2255 sub( Y, V, t );\r
2256 add( V, V0, t2 );\r
2257 if( cmp(t,t2) == 0 )\r
2258         printf("Seems O.K.\n");\r
2259 else\r
2260         {\r
2261         printf("finds a Flaw, -(-Y) differs from Y.\n");\r
2262         printf( "V-Y=t:" );\r
2263         show(V);\r
2264         show(Y);\r
2265         show(t);\r
2266         printf( "V+V0=t2:" );\r
2267         show(V);\r
2268         show(V0);\r
2269         show(t2);\r
2270         ErrCnt[Flaw] += 1;\r
2271         }\r
2272 if( (cmp(Z, Y) != 0) && (I != 0) )\r
2273         {\r
2274         ErrCnt[Serious] += 1;\r
2275         printf("overflow past " );\r
2276         show( Y );\r
2277         printf( "\tshrinks to " );\r
2278         show( Z );\r
2279         printf( "= Y * " );\r
2280         show( HInvrse );\r
2281         }\r
2282 /*Y = V * (HInvrse * U2 - HInvrse);*/\r
2283 mul( HInvrse, U2, Y );\r
2284 sub( HInvrse, Y, Y );\r
2285 mul( V, Y, Y );\r
2286 /*Z = Y + ((One - HInvrse) * U2) * V;*/\r
2287 sub( HInvrse, One, Z );\r
2288 mul( Z, U2, Z );\r
2289 mul( Z, V, Z );\r
2290 add( Y, Z, Z );\r
2291 if( cmp(Z,V0) < 0 )\r
2292         mov( Z, Y );\r
2293 if( cmp(Y,V0) < 0)\r
2294         mov( Y, V );\r
2295 sub( V, V0, t );\r
2296 if( cmp(t,V0) < 0 )\r
2297         mov( V0, V );\r
2298 printf("Overflow threshold is V  = " );\r
2299 show( V );\r
2300 if(I)\r
2301         {\r
2302         printf("Overflow saturates at V0 = " );\r
2303         show( V0 );\r
2304         }\r
2305 else\r
2306 printf("There is no saturation value because the system traps on overflow.\n");\r
2307 \r
2308 mul( V, One, V9 );\r
2309 printf("No Overflow should be signaled for V * 1 = " );\r
2310 show( V9 );\r
2311 div( One, V, V9 );\r
2312         printf("                           nor for V / 1 = " );\r
2313         show( V9 );\r
2314         printf("Any overflow signal separating this * from the one\n");\r
2315         printf("above is a DEFECT.\n");\r
2316 /*=============================================*/\r
2317 Milestone = 170;\r
2318 /*=============================================*/\r
2319 mov( V, t );\r
2320 neg( t );\r
2321 k = 0;\r
2322 if( cmp(t,V) >= 0 )\r
2323         k = 1;\r
2324 mov( V0, t );\r
2325 neg( t );\r
2326 if( cmp(t,V0) >= 0 )\r
2327         k = 1;\r
2328 mov( UfThold, t );\r
2329 neg(t);\r
2330 if( cmp(t,V) >= 0 )\r
2331         k = 1;\r
2332 if( cmp(UfThold,V) >= 0 )\r
2333         k = 1;\r
2334 if( k != 0 )\r
2335         {\r
2336         ErrCnt[Failure] += 1;\r
2337         printf( "Comparisons involving +-");\r
2338         show( V );\r
2339         show( V0 );\r
2340         show( UfThold );\r
2341         printf("are confused by Overflow." );\r
2342         }\r
2343 /*=============================================*/\r
2344 Milestone = 175;\r
2345 /*=============================================*/\r
2346 printf("\n");\r
2347 for(Indx = 1; Indx <= 3; ++Indx) {\r
2348         switch(Indx)\r
2349                 {\r
2350                 case 1: mov(UfThold, Z); break;\r
2351                 case 2: mov( E0, Z); break;\r
2352                 case 3: mov(PseudoZero, Z); break;\r
2353                 }\r
2354 if( cmp(Z, Zero) != 0 )\r
2355         {\r
2356         SQRT( Z, V9 );\r
2357         mul( V9, V9, Y );\r
2358         mul( Radix, E9, t );\r
2359         sub( t, One, t );\r
2360         div( t, Y, t );\r
2361         add( One, Radix, t2 );\r
2362         add( t2, E9, t2 );\r
2363         mul( t2, Z, t2 );\r
2364         if( (cmp(t,Z) < 0) || (cmp(Y,t2) > 0) )\r
2365                 {\r
2366                 if( cmp(V9,U1) > 0 )\r
2367                         ErrCnt[Serious] += 1;\r
2368                 else\r
2369                         ErrCnt[Defect] += 1;\r
2370                 printf("Comparison alleges that what prints as Z = " );\r
2371                 show( Z );\r
2372                 printf(" is too far from sqrt(Z) ^ 2 = " );\r
2373                 show( Y );\r
2374                 }\r
2375         }\r
2376 }\r
2377 \r
2378 Milestone = 180;\r
2379 \r
2380 for(Indx = 1; Indx <= 2; ++Indx)\r
2381         {\r
2382         if(Indx == 1)\r
2383                 mov( V, Z );\r
2384         else\r
2385                 mov( V0, Z );\r
2386         SQRT( Z, V9 );\r
2387         mul( Radix, E9, X );\r
2388         sub( X, One, X );\r
2389         mul( X, V9, X );\r
2390         mul( V9, X, V9 );\r
2391         mul( Two, Radix, t );\r
2392         mul( t, E9, t );\r
2393         sub( t, One, t );\r
2394         mul( t, Z, t );\r
2395         if( (cmp(V9,t) < 0) || (cmp(V9,Z) > 0) )\r
2396                 {\r
2397                 mov( V9, Y );\r
2398                 if( cmp(X,W) <  0 )\r
2399                         ErrCnt[Serious] += 1;\r
2400                 else\r
2401                         ErrCnt[Defect] += 1;\r
2402                 printf("Comparison alleges that Z = " );\r
2403                 show( Z );\r
2404                 printf(" is too far from sqrt(Z) ^ 2 :" );\r
2405                 show( Y );\r
2406                 }\r
2407         }\r
2408 \r
2409 Milestone = 190;\r
2410 \r
2411 Pause();\r
2412 mul( UfThold, V, X ); \r
2413 mul( Radix, Radix, Y );\r
2414 mul( X, Y, t );\r
2415 if( (cmp(t,One) < 0) || (cmp(X,Y) > 0) )\r
2416         {\r
2417         mul( X, Y, t );\r
2418         div( U1, Y, t2 );\r
2419         if( (cmp(t,U1) < 0) || (cmp(X,t2) > 0) )\r
2420                 {\r
2421                 ErrCnt[Defect] += 1;\r
2422                 printf( "Badly " );\r
2423                 }\r
2424         else\r
2425                 {\r
2426                 ErrCnt[Flaw] += 1;\r
2427                 }\r
2428         printf(" unbalanced range; UfThold * V = " );\r
2429         show( X );\r
2430         printf( "\tis too far from 1.\n");\r
2431         }\r
2432 Milestone = 200;\r
2433 \r
2434 for(Indx = 1; Indx <= 5; ++Indx)\r
2435         {\r
2436         mov( F9, X );\r
2437         switch(Indx)\r
2438                 {\r
2439                 case 2: add( One, U2, X ); break;\r
2440                 case 3: mov( V, X ); break;\r
2441                 case 4: mov(UfThold,X); break;\r
2442                 case 5: mov(Radix,X);\r
2443                 }\r
2444         mov( X, Y );\r
2445 \r
2446         sigsave = sigfpe;\r
2447         if (setjmp(ovfl_buf))\r
2448                 {\r
2449                 printf("  X / X  traps when X = " );\r
2450                 show( X );\r
2451                 }\r
2452         else\r
2453                 {\r
2454 /*V9 = (Y / X - Half) - Half;*/\r
2455                 div( X, Y, t );\r
2456                 sub( Half, t, t );\r
2457                 sub( Half, t, V9 );\r
2458                 if( cmp(V9,Zero) == 0 )\r
2459                         continue;\r
2460                 mov( U1, t );\r
2461                 neg(t);\r
2462                 if( (cmp(V9,t) == 0) && (Indx < 5) )\r
2463                         {\r
2464                         ErrCnt[Flaw] += 1;\r
2465                         }\r
2466                 else\r
2467                         {\r
2468                         ErrCnt[Serious] += 1;\r
2469                         }\r
2470                 printf("  X / X differs from 1 when X = " );\r
2471                 show( X );\r
2472                 printf("  instead, X / X - 1/2 - 1/2 = " );\r
2473                 show( V9 );\r
2474                 }\r
2475         }\r
2476 \r
2477         Pause();\r
2478         printf("\n");\r
2479         {\r
2480                 static char *msg[] = {\r
2481                         "FAILUREs  encountered =",\r
2482                         "SERIOUS DEFECTs  discovered =",\r
2483                         "DEFECTs  discovered =",\r
2484                         "FLAWs  discovered =" };\r
2485                 int i;\r
2486                 for(i = 0; i < 4; i++) if (ErrCnt[i])\r
2487                         printf("The number of  %-29s %d.\n",\r
2488                                 msg[i], ErrCnt[i]);\r
2489                 }\r
2490         printf("\n");\r
2491         if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]\r
2492                         + ErrCnt[Flaw]) > 0) {\r
2493                 if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[\r
2494                         Defect] == 0) && (ErrCnt[Flaw] > 0)) {\r
2495                         printf("The arithmetic diagnosed seems ");\r
2496                         printf("satisfactory though flawed.\n");\r
2497                         }\r
2498                 if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)\r
2499                         && ( ErrCnt[Defect] > 0)) {\r
2500                         printf("The arithmetic diagnosed may be acceptable\n");\r
2501                         printf("despite inconvenient Defects.\n");\r
2502                         }\r
2503                 if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {\r
2504                         printf("The arithmetic diagnosed has ");\r
2505                         printf("unacceptable serious defects.\n");\r
2506                         }\r
2507                 if (ErrCnt[Failure] > 0) {\r
2508                         printf("Fatal FAILURE may have spoiled this");\r
2509                         printf(" program's subsequent diagnoses.\n");\r
2510                         }\r
2511                 }\r
2512         else {\r
2513                 printf("No failures, defects nor flaws have been discovered.\n");\r
2514                 if (! ((RMult == Rounded) && (RDiv == Rounded)\r
2515                         && (RAddSub == Rounded) && (RSqrt == Rounded))) \r
2516                         printf("The arithmetic diagnosed seems satisfactory.\n");\r
2517                 else {\r
2518                         k = 0;\r
2519                         if( cmp( Radix, Two ) == 0 )\r
2520                                 k = 1;\r
2521                         if( cmp( Radix, Ten ) == 0 )\r
2522                                 k = 1;\r
2523                         if( (cmp(StickyBit,One) >= 0) && (k == 1) )\r
2524                                 {\r
2525                                 printf("Rounding appears to conform to ");\r
2526                                 printf("the proposed IEEE standard P");\r
2527                                 k = 0;\r
2528                                 k |= cmp( Radix, Two );\r
2529                                 mul( Four, Three, t );\r
2530                                 mul( t, Two, t );\r
2531                                 sub( t, Precision, t );\r
2532                                 sub( TwentySeven, Precision, t2 );\r
2533                                 sub( TwentySeven, t2, t2 );\r
2534                                 add( t2, One, t2 );\r
2535                                 mul( t2, t, t );\r
2536                                 if( (cmp(Radix,Two) == 0)\r
2537                                         && (cmp(t,Zero) == 0) )\r
2538                                         printf("754");\r
2539                                 else\r
2540                                         printf("854");\r
2541                                 if(IEEE)\r
2542                                         printf(".\n");\r
2543                                 else\r
2544                                         {\r
2545                         printf(",\nexcept for possibly Double Rounding");\r
2546                         printf(" during Gradual Underflow.\n");\r
2547                                         }\r
2548                                 }\r
2549                 printf("The arithmetic diagnosed appears to be excellent!\n");\r
2550                         }\r
2551                 }\r
2552         if (fpecount)\r
2553                 printf("\nA total of %d floating point exceptions were registered.\n",\r
2554                         fpecount);\r
2555         printf("END OF TEST.\n");\r
2556         }\r
2557 \r
2558 \r
2559 /* Random */\r
2560 /*  Random computes\r
2561      X = (Random1 + Random9)^5\r
2562      Random1 = X - FLOOR(X) + 0.000005 * X;\r
2563    and returns the new value of Random1\r
2564 */\r
2565 \r
2566 \r
2567 static int randflg = 0;\r
2568 FLOAT(C5em6);\r
2569 \r
2570 Random()\r
2571 {\r
2572 \r
2573 if( randflg == 0 )\r
2574         {\r
2575         mov( Six, t );\r
2576         neg(t);\r
2577         POW( Ten, t, t );\r
2578         mul( Five, t, C5em6 );\r
2579         randflg = 1;\r
2580         }\r
2581 add( Random1, Random9, t );\r
2582 mul( t, t, t2 );\r
2583 mul( t2, t2, t2 );\r
2584 mul( t, t2, t );\r
2585 FLOOR(t, t2 );\r
2586 sub( t2, t, t2 );\r
2587 mul( t, C5em6, t );\r
2588 add( t, t2, Random1 );\r
2589 /*return(Random1);*/\r
2590 }\r
2591 \r
2592 /* SqXMinX */\r
2593 \r
2594 SqXMinX( ErrKind )\r
2595 int ErrKind;\r
2596 {\r
2597 mul( X, BInvrse, t2 );\r
2598 sub( t2, X, t );\r
2599 /*SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;*/\r
2600 mul( X, X, Sqarg );\r
2601 SQRT( Sqarg, SqEr );\r
2602 sub( t2, SqEr, SqEr );\r
2603 sub( t, SqEr, SqEr );\r
2604 div( OneUlp, SqEr, SqEr );\r
2605 if( cmp(SqEr,Zero) != 0)\r
2606         {\r
2607         Showsq( 0 );\r
2608         add( J, One, J );\r
2609         ErrCnt[ErrKind] += 1;\r
2610         printf("sqrt of " );\r
2611         mul( X, X, t );\r
2612         show( t );\r
2613         printf( "minus " );\r
2614         show( X );\r
2615         printf( "equals " );\r
2616         mul( OneUlp, SqEr, t );\r
2617         show( t );\r
2618         printf("\tinstead of correct value 0 .\n");\r
2619         }\r
2620 }\r
2621 \r
2622 /* NewD */\r
2623 \r
2624 NewD()\r
2625 {\r
2626 mul( Z1, Q, X );\r
2627 /*X = FLOOR(Half - X / Radix) * Radix + X;*/\r
2628 div( Radix, X, t );\r
2629 sub( t, Half, t );\r
2630 FLOOR( t, t );\r
2631 mul( t, Radix, t );\r
2632 add( t, X, X );\r
2633 /*Q = (Q - X * Z) / Radix + X * X * (D / Radix);*/\r
2634 mul( X, Z, t );\r
2635 sub( t, Q, t );\r
2636 div( Radix, t, t );\r
2637 div( Radix, D, t2 );\r
2638 mul( X, t2, t2 );\r
2639 mul( X, t2, t2 );\r
2640 add( t, t2, Q );\r
2641 /*Z = Z - Two * X * D;*/\r
2642 mul( Two, X, t );\r
2643 mul( t, D, t );\r
2644 sub( t, Z, Z );\r
2645 \r
2646 if( cmp(Z,Zero) <= 0)\r
2647         {\r
2648         neg(Z);\r
2649         neg(Z1);\r
2650         }\r
2651 mul( Radix, D, D );\r
2652 }\r
2653 \r
2654 /* SR3750 */\r
2655 \r
2656 SR3750()\r
2657 {\r
2658 sub( Radix, X, t );\r
2659 sub( Radix, Z2, t2 );\r
2660 k = 0;\r
2661 if( cmp(t,t2) < 0 )\r
2662         k = 1;\r
2663 sub( Z2, X, t );\r
2664 sub( Z2, W, t2 );\r
2665 if( cmp(t,t2) > 0 )\r
2666         k = 1;\r
2667 /*if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {*/\r
2668 if( k == 0 )\r
2669         {\r
2670         I = I + 1;\r
2671         mul( X, D, X2 );\r
2672         mov( X2, Sqarg );\r
2673         SQRT( X2, X2 );\r
2674 /*Y2 = (X2 - Z2) - (Y - Z2);*/\r
2675         sub( Z2, X2, Y2 );\r
2676         sub( Z2, Y, t );\r
2677         sub( t, Y2, Y2 );\r
2678         sub( Half, Y, X2 );\r
2679         div( X2, X8, X2 );\r
2680         mul( Half, X2, t );\r
2681         mul( t, X2, t );\r
2682         sub( t, X2, X2 );\r
2683 /*SqEr = (Y2 + Half) + (Half - X2);*/\r
2684         add( Y2, Half, SqEr );\r
2685         sub( X2, Half, t );\r
2686         add( t, SqEr, SqEr );\r
2687         Showsq( -1 );\r
2688         sub( X2, Y2, SqEr );\r
2689         Showsq( 1 );\r
2690         }\r
2691 }\r
2692 \r
2693 /* IsYeqX */\r
2694 \r
2695 IsYeqX()\r
2696 {\r
2697 if( cmp(Y,X) != 0 )\r
2698         {\r
2699         if (N <= 0)\r
2700                 {\r
2701                 if( (cmp(Z,Zero) == 0) && (cmp(Q,Zero) <= 0) )\r
2702                         printf("WARNING:  computing\n");\r
2703                 else\r
2704                         {\r
2705                         ErrCnt[Defect] += 1;\r
2706                         printf( "computing\n");\r
2707                         }\r
2708                 show( Z );\r
2709                 printf( "\tto the power " );\r
2710                 show( Q );\r
2711                 printf("\tyielded " );\r
2712                 show( Y );\r
2713                 printf("\twhich compared unequal to correct " );\r
2714                 show( X );\r
2715                 sub( X, Y, t );\r
2716                 printf("\t\tthey differ by " );\r
2717                 show( t );\r
2718                 }\r
2719         N = N + 1; /* ... count discrepancies. */\r
2720         }\r
2721 }\r
2722 \r
2723 /* SR3980 */\r
2724 \r
2725 SR3980()\r
2726 {\r
2727 long li;\r
2728 \r
2729 do\r
2730         {\r
2731 /*Q = (FLOAT) I;*/\r
2732         li = I;\r
2733         LTOF( &li, Q );\r
2734         POW( Z, Q, Y );\r
2735         IsYeqX();\r
2736         if(++I > M)\r
2737                 break;\r
2738         mul( Z, X, X );\r
2739         }\r
2740 while( cmp(X,W) < 0 );\r
2741 }\r
2742 \r
2743 /* PrintIfNPositive */\r
2744 \r
2745 PrintIfNPositive()\r
2746 {\r
2747 if(N > 0)\r
2748         printf("Similar discrepancies have occurred %d times.\n", N);\r
2749 }\r
2750 \r
2751 \r
2752 /* TstPtUf */\r
2753 \r
2754 TstPtUf()\r
2755 {\r
2756 N = 0;\r
2757 if( cmp(Z,Zero) != 0)\r
2758         {\r
2759         printf( "Z = " );\r
2760         show(Z);\r
2761         printf("Since comparison denies Z = 0, evaluating ");\r
2762         printf("(Z + Z) / Z should be safe.\n");\r
2763         sigsave = sigfpe;\r
2764         if (setjmp(ovfl_buf))\r
2765                 goto very_serious;\r
2766         add( Z, Z, Q9 );\r
2767         div( Z, Q9, Q9 );\r
2768         printf("What the machine gets for (Z + Z) / Z is " );\r
2769         show( Q9 );\r
2770         sub( Two, Q9, t );\r
2771         FABS(t);\r
2772         mul( Radix, U2, t2 );\r
2773         if( cmp(t,t2) < 0 )\r
2774                 {\r
2775                 printf("This is O.K., provided Over/Underflow");\r
2776                 printf(" has NOT just been signaled.\n");\r
2777                 }\r
2778         else\r
2779                 {\r
2780                 if( (cmp(Q9,One) < 0) || (cmp(Q9,Two) > 0) )\r
2781                         {\r
2782 very_serious:\r
2783                         N = 1;\r
2784                         ErrCnt [Serious] = ErrCnt [Serious] + 1;\r
2785                         printf("This is a VERY SERIOUS DEFECT!\n");\r
2786                         }\r
2787                 else\r
2788                         {\r
2789                         N = 1;\r
2790                         ErrCnt[Defect] += 1;\r
2791                         printf("This is a DEFECT!\n");\r
2792                         }\r
2793                 }\r
2794         mul( Z, One, V9 );\r
2795         mov( V9, Random1 );\r
2796         mul( One, Z, V9 );\r
2797         mov( V9, Random2 );\r
2798         div( One, Z, V9 );\r
2799         if( (cmp(Z,Random1) == 0) && (cmp(Z,Random2) == 0)\r
2800                 && (cmp(Z,V9) == 0) )\r
2801                 {\r
2802                 if (N > 0)\r
2803                         Pause();\r
2804                 }\r
2805         else\r
2806                 {\r
2807                 N = 1;\r
2808                 ErrCnt[Defect] += 1;\r
2809                 printf( "What prints as Z = ");\r
2810                 show( Z );\r
2811                 printf( "\tcompares different from " );\r
2812                 if( cmp(Z,Random1) != 0)\r
2813                         {\r
2814                         printf("Z * 1 = " );\r
2815                         show( Random1 );\r
2816                         }\r
2817                 if( (cmp(Z,Random2) != 0)\r
2818                         || (cmp(Random2,Random1) != 0) )\r
2819                         {\r
2820                         printf("1 * Z == " );\r
2821                         show( Random2 );\r
2822                         }\r
2823                 if( cmp(Z,V9) != 0 )\r
2824                         {\r
2825                         printf("Z / 1 = " );\r
2826                         show( V9 );\r
2827                         }\r
2828                 if( cmp(Random2,Random1) != 0 )\r
2829                         {\r
2830                         ErrCnt[Defect] += 1;\r
2831                         printf( "Multiplication does not commute!\n");\r
2832                         printf("\tComparison alleges that 1 * Z = " );\r
2833                         show(Random2);\r
2834                         printf("\tdiffers from Z * 1 = " );\r
2835                         show(Random1);\r
2836                         }\r
2837                 Pause();\r
2838                 }\r
2839         }\r
2840 }\r
2841 \r
2842 Pause()\r
2843 {\r
2844 }\r
2845 \r
2846 Sign( x, y )\r
2847 FSIZE *x, *y;\r
2848 {\r
2849 \r
2850 if( cmp( x, Zero ) < 0 )\r
2851         {\r
2852         mov( One, y );\r
2853         neg( y );\r
2854         }\r
2855 else\r
2856         {\r
2857         mov( One, y );\r
2858         }\r
2859 }\r
2860 \r
2861 sqtest()\r
2862 {\r
2863 printf("\nRunning test of square root(x).\n");\r
2864 \r
2865 RSqrt = Other;\r
2866 k = 0;\r
2867 SQRT( Zero, t );\r
2868 k |= cmp( Zero, t );\r
2869 mov( Zero, t );\r
2870 neg(t);\r
2871 SQRT( t, t2 );\r
2872 k |= cmp( t, t2 );\r
2873 SQRT( One, t );\r
2874 k |= cmp( One, t );\r
2875 if( k != 0 )\r
2876         {\r
2877         ErrCnt[Failure] += 1;\r
2878         printf( "Square root of 0.0, -0.0 or 1.0 wrong\n");\r
2879         }\r
2880 mov( Zero, MinSqEr );\r
2881 mov( Zero, MaxSqEr );\r
2882 mov( Zero, J );\r
2883 mov( Radix, X );\r
2884 mov( U2, OneUlp );\r
2885 SqXMinX( Serious );\r
2886 mov( BInvrse, X );\r
2887 mul( BInvrse, U1, OneUlp );\r
2888 SqXMinX( Serious );\r
2889 mov( U1, X );\r
2890 mul( U1, U1, OneUlp );\r
2891 SqXMinX( Serious );\r
2892 if( cmp(J,Zero) != 0)\r
2893         Pause();\r
2894 printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);\r
2895 mov( Zero, J );\r
2896 mov( Two, X );\r
2897 mov( Radix, Y );\r
2898 if( cmp(Radix,One) != 0 )\r
2899         {\r
2900         lngint = NoTrials;\r
2901         LTOF( &lngint, t );\r
2902         FTOL( t, &lng2, X );\r
2903         if( lngint != lng2 )\r
2904                 {\r
2905                 printf( "Integer conversion error\n" );\r
2906                 exit(1);\r
2907                 }\r
2908         do\r
2909                 {\r
2910                 mov( Y, X );\r
2911                 mul( Radix, Y, Y );\r
2912                 sub( X, Y, t2 );\r
2913                 }\r
2914         while( ! (cmp(t2,t) >= 0) );\r
2915         }\r
2916 mul( X, U2, OneUlp );\r
2917 I = 1;\r
2918 while(I < 10)\r
2919         {\r
2920         add( X, One, X );\r
2921         SqXMinX( Defect );\r
2922         if( cmp(J,Zero) > 0 )\r
2923                 break;\r
2924         I = I + 1;\r
2925         }\r
2926 printf("Test for sqrt monotonicity.\n");\r
2927 I = - 1;\r
2928 mov( BMinusU2, X );\r
2929 mov( Radix, Y );\r
2930 mul( Radix, U2, Z );\r
2931 add( Radix, Z, Z );\r
2932 NotMonot = False;\r
2933 Monot = False;\r
2934 while( ! (NotMonot || Monot))\r
2935         {\r
2936         I = I + 1;\r
2937         SQRT(X, X);\r
2938         SQRT(Y,Q);\r
2939         SQRT(Z,Z);\r
2940         if( (cmp(X,Q) > 0) || (cmp(Q,Z) > 0) )\r
2941                 NotMonot = True;\r
2942         else\r
2943                 {\r
2944                 add( Q, Half, Q );\r
2945                 FLOOR( Q, Q );\r
2946                 mul( Q, Q, t );\r
2947                 if( (I > 0) || (cmp(Radix,t) == 0) )\r
2948                         Monot = True;\r
2949                 else if (I > 0)\r
2950                         {\r
2951                         if(I > 1)\r
2952                                 Monot = True;\r
2953                         else\r
2954                                 {\r
2955                                 mul( Y, BInvrse, Y );\r
2956                                 sub( U1, Y, X );\r
2957                                 add( Y, U1, Z );\r
2958                                 }\r
2959                         }\r
2960                 else\r
2961                         {\r
2962                         mov( Q, Y );\r
2963                         sub( U2, Y, X );\r
2964                         add( Y, U2, Z );\r
2965                         }\r
2966                 }\r
2967         }\r
2968 if( Monot )\r
2969         printf("sqrt has passed a test for Monotonicity.\n");\r
2970 else\r
2971         {\r
2972         ErrCnt[Defect] += 1;\r
2973         printf("sqrt(X) is non-monotonic for X near " );\r
2974         show(Y);\r
2975         }\r
2976 /*=============================================*/\r
2977 Milestone = 80;\r
2978 /*=============================================*/\r
2979 add( MinSqEr, Half, MinSqEr );\r
2980 sub( Half, MaxSqEr, MaxSqEr);\r
2981 /*Y = (SQRT(One + U2) - One) / U2;*/\r
2982 add( One, U2, Sqarg );\r
2983 SQRT( Sqarg, Y );\r
2984 sub( One, Y, Y );\r
2985 div( U2, Y, Y );\r
2986 /*SqEr = (Y - One) + U2 / Eight;*/\r
2987 sub( One, Y, t );\r
2988 div( Eight, U2, SqEr );\r
2989 add( t, SqEr, SqEr );\r
2990 Showsq( 1 );\r
2991 div( Eight, U2, SqEr );\r
2992 add( Y, SqEr, SqEr );\r
2993 Showsq( -1 );\r
2994 /*Y = ((SQRT(F9) - U2) - (One - U2)) / U1;*/\r
2995 mov( F9, Sqarg );\r
2996 SQRT( Sqarg, Y );\r
2997 sub( U2, Y, Y );\r
2998 sub( U2, One, t );\r
2999 sub( t, Y, Y );\r
3000 div( U1, Y, Y );\r
3001 div( Eight, U1, SqEr );\r
3002 add( Y, SqEr, SqEr );\r
3003 Showsq( 1 );\r
3004 /*SqEr = (Y + One) + U1 / Eight;*/\r
3005 div( Eight, U1, t );\r
3006 add( Y, One, SqEr );\r
3007 add( SqEr, t, SqEr );\r
3008 Showsq( -1 );\r
3009 mov( U2, OneUlp );\r
3010 mov( OneUlp, X );\r
3011 for( Indx = 1; Indx <= 3; ++Indx)\r
3012         {\r
3013 /*Y = SQRT((X + U1 + X) + F9);*/\r
3014         add( X, U1, Y );\r
3015         add( Y, X, Y );\r
3016         add( Y, F9, Y );\r
3017         mov( Y, Sqarg );\r
3018         SQRT( Sqarg, Y );\r
3019 /*Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;*/\r
3020         sub( U2, One, t );\r
3021         add( t, X, t );\r
3022         sub( U2, Y, Y );\r
3023         sub( t, Y, Y );\r
3024         div( OneUlp, Y, Y );\r
3025 /*Z = ((U1 - X) + F9) * Half * X * X / OneUlp;*/\r
3026         sub( X, U1, t );\r
3027         add( t, F9, t );\r
3028         mul( t, Half, t );\r
3029         mul( t, X, t );\r
3030         mul( t, X, t );\r
3031         div( OneUlp, t, Z );\r
3032         add( Y, Half, SqEr );\r
3033         add( SqEr, Z, SqEr );\r
3034         Showsq( -1 );\r
3035         sub( Half, Y, SqEr );\r
3036         add( SqEr, Z, SqEr );\r
3037         Showsq( 1 );\r
3038         if(((Indx == 1) || (Indx == 3))) \r
3039                 {\r
3040 /*X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));*/\r
3041                 mov( OneUlp, Sqarg );\r
3042                 SQRT( Sqarg, t );\r
3043                 mul( Nine, t, t );\r
3044                 div( t, Eight, t );\r
3045                 FLOOR( t, t );\r
3046                 Sign( X, t2 );\r
3047                 mul( t2, t, t );\r
3048                 mul( OneUlp, t, X );\r
3049                 }\r
3050         else\r
3051                 {\r
3052                 mov( U1, OneUlp );\r
3053                 mov( OneUlp, X );\r
3054                 neg( X );\r
3055                 }\r
3056         }\r
3057 /*=============================================*/\r
3058 Milestone = 85;\r
3059 /*=============================================*/\r
3060 SqRWrng = False;\r
3061 Anomaly = False;\r
3062 if( cmp(Radix,One) != 0 )\r
3063         {\r
3064         printf("Testing whether sqrt is rounded or chopped.\n");\r
3065 /*D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));*/\r
3066         FLOOR( Precision, t2 );\r
3067         add( One, Precision, t );\r
3068         sub( t2, t, t );\r
3069         POW( Radix, t, D );\r
3070         add( Half, D, D );\r
3071         FLOOR( D, D );\r
3072 /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */\r
3073         div( Radix, D, X );\r
3074         div( A1, D, Y );\r
3075         FLOOR( X, t );\r
3076         FLOOR( Y, t2 );\r
3077         if( (cmp(X,t) != 0) || (cmp(Y,t2) != 0) )\r
3078                 {\r
3079                 Anomaly = True;\r
3080                 printf( "Anomaly 1\n" );\r
3081                 }\r
3082         else\r
3083                 {\r
3084                 mov( Zero, X );\r
3085                 mov( X, Z2 );\r
3086                 mov( One, Y );\r
3087                 mov( Y, Y2 );\r
3088                 sub( One, Radix, Z1 );\r
3089                 mul( Four, D, FourD );\r
3090                 do\r
3091                         {\r
3092                         if( cmp(Y2,Z2) >0 )\r
3093                                 {\r
3094                                 mov( Radix, Q );\r
3095                                 mov( Y, YY1 );\r
3096                                 do\r
3097                                         {\r
3098 /*X1 = FABS(Q + FLOOR(Half - Q / YY1) * YY1);*/\r
3099                                         div( YY1, Q, t );\r
3100                                         sub( t, Half, t );\r
3101                                         FLOOR( t, t );\r
3102                                         mul( t, YY1, t );\r
3103                                         add( Q, t, X1 );\r
3104                                         FABS( X1 );\r
3105                                         mov( YY1, Q );\r
3106                                         mov( X1, YY1 );\r
3107                                         }\r
3108                                 while( ! (cmp(X1,Zero) <= 0) );\r
3109                                 if( cmp(Q,One) <= 0 )\r
3110                                         {\r
3111                                         mov( Y2, Z2 );\r
3112                                         mov( Y, Z );\r
3113                                         }\r
3114                                 }\r
3115                         add( Y, Two, Y );\r
3116                         add( X, Eight, X );\r
3117                         add( Y2, X, Y2 );\r
3118                         if( cmp(Y2,FourD) >= 0 )\r
3119                                 sub( FourD, Y2, Y2 );\r
3120                         }\r
3121                 while( ! (cmp(Y,D) >= 0) );\r
3122                 sub( Z2, FourD, X8 );\r
3123                 mul( Z, Z, Q );\r
3124                 add( X8, Q, Q );\r
3125                 div( FourD, Q, Q );\r
3126                 div( Eight, X8, X8 );\r
3127                 FLOOR( Q, t );\r
3128                 if( cmp(Q,t) != 0 )\r
3129                         {\r
3130                         Anomaly = True;\r
3131                         printf( "Anomaly 2\n" );\r
3132                         }\r
3133                 else\r
3134                         {\r
3135                         Break = False;\r
3136                         do\r
3137                                 {\r
3138                                 mul( Z1, Z, X );\r
3139 /*X = X - FLOOR(X / Radix) * Radix;*/\r
3140                                 div( Radix, X, t );\r
3141                                 FLOOR( t, t );\r
3142                                 mul( t, Radix, t );\r
3143                                 sub( t, X, X );\r
3144                                 if( cmp(X,One) == 0 ) \r
3145                                         Break = True;\r
3146                                 else\r
3147                                         sub( One, Z1, Z1 );\r
3148                                 }\r
3149                         while( ! (Break || (cmp(Z1,Zero) <= 0)) );\r
3150                         if( (cmp(Z1,Zero) <= 0) && (! Break))\r
3151                                 {\r
3152                                 printf( "Anomaly 3\n" );\r
3153                                 Anomaly = True;\r
3154                                 }\r
3155                         else\r
3156                                 {\r
3157                                 if( cmp(Z1,RadixD2) > 0)\r
3158                                         sub( Radix, Z1, Z1 );\r
3159                                 do\r
3160                                         {\r
3161                                         NewD();\r
3162                                         mul( U2, D, t );\r
3163                                         }\r
3164                                 while( ! (cmp(t,F9) >= 0) );\r
3165                                 mul( D, Radix, t );\r
3166                                 sub( D, t, t );\r
3167                                 sub( D, W, t2 );\r
3168                                 if (cmp(t,t2) != 0 )\r
3169                                         {\r
3170                                         printf( "Anomaly 4\n" );\r
3171                                         Anomaly = True;\r
3172                                         }\r
3173                                 else\r
3174                                         {\r
3175                                         mov( D, Z2 );\r
3176                                         I = 0;\r
3177                                         add( One, Z, t );\r
3178                                         mul( t, Half, t );\r
3179                                         add( D, t, Y );\r
3180                                         add( D, Z, X );\r
3181                                         add( X, Q, X );\r
3182                                         SR3750();\r
3183                                         sub( Z, One, t );\r
3184                                         mul( t, Half, t );\r
3185                                         add( D, t, Y );\r
3186                                         add( Y, D, Y );\r
3187                                         sub( Z, D, X );\r
3188                                         add( X, D, X );\r
3189                                         add( X, Q, t );\r
3190                                         add( t, X, X );\r
3191                                         SR3750();\r
3192                                         NewD();\r
3193                                         sub( Z2, D, t );\r
3194                                         sub( Z2, W, t2 );\r
3195                                         if(cmp(t,t2) != 0 )\r
3196                                                 {\r
3197                                                 printf( "Anomaly 5\n" );\r
3198                                                 Anomaly = True;\r
3199                                                 }\r
3200                                         else\r
3201                                                 {\r
3202 /*Y = (D - Z2) + (Z2 + (One - Z) * Half);*/\r
3203                                                 sub( Z, One, t );\r
3204                                                 mul( t, Half, t );\r
3205                                                 add( Z2, t, t );\r
3206                                                 sub( Z2, D, Y );\r
3207                                                 add( Y, t, Y );\r
3208 /*X = (D - Z2) + (Z2 - Z + Q);*/\r
3209                                                 sub( Z, Z2, t );\r
3210                                                 add( t, Q, t );\r
3211                                                 sub( Z2, D, X );\r
3212                                                 add( X, t, X );\r
3213                                                 SR3750();\r
3214                                                 add( One, Z, Y );\r
3215                                                 mul( Y, Half, Y );\r
3216                                                 mov( Q, X );\r
3217                                                 SR3750();\r
3218                                                 if(I == 0)\r
3219                                                         {\r
3220                                                         printf( "Anomaly 6\n" );\r
3221                                                         Anomaly = True;\r
3222                                                         }\r
3223                                                 }\r
3224                                         }\r
3225                                 }\r
3226                         }\r
3227                 }\r
3228         if ((I == 0) || Anomaly)\r
3229                 {\r
3230                 ErrCnt[Failure] += 1;\r
3231                 printf( "Anomalous arithmetic with Integer < \n");\r
3232                 printf("Radix^Precision = " );\r
3233                 show( W );\r
3234                 printf(" fails test whether sqrt rounds or chops.\n");\r
3235                 SqRWrng = True;\r
3236                 }\r
3237         }\r
3238 if(! Anomaly)\r
3239         {\r
3240         if(! ((cmp(MinSqEr,Zero) < 0) || (cmp(MaxSqEr,Zero) > 0))) {\r
3241         RSqrt = Rounded;\r
3242         printf("Square root appears to be correctly rounded.\n");\r
3243         }\r
3244         else\r
3245                 {\r
3246                 k = 0;\r
3247                 add( MaxSqEr, U2, t );\r
3248                 sub( Half, U2, t2 );\r
3249                 if( cmp(t,t2) > 0 )\r
3250                         k = 1;\r
3251                 if( cmp( MinSqEr, Half ) > 0 )\r
3252                         k = 1;\r
3253                 add( MinSqEr, Radix, t );\r
3254                 if( cmp( t, Half ) < 0 )\r
3255                         k = 1;\r
3256                 if( k == 1 )\r
3257                         SqRWrng = True;\r
3258                 else\r
3259                         {\r
3260                         RSqrt = Chopped;\r
3261                         printf("Square root appears to be chopped.\n");\r
3262                         }\r
3263                 }\r
3264         }\r
3265 if( SqRWrng )\r
3266         {\r
3267         printf("Square root is neither chopped nor correctly rounded.\n");\r
3268         printf("Observed errors run from " );\r
3269         sub( Half, MinSqEr, t );\r
3270         show( t );\r
3271         printf("\tto " );\r
3272         add( Half, MaxSqEr, t );\r
3273         show( t );\r
3274         printf( "ulps.\n" );\r
3275         sub( MinSqEr, MaxSqEr, t );\r
3276         mul( Radix, Radix, t2 );\r
3277         if( cmp( t, t2 ) >= 0 )\r
3278                 {\r
3279                 ErrCnt[Serious] += 1;\r
3280                 printf( "sqrt gets too many last digits wrong\n");\r
3281                 }\r
3282         }\r
3283 }\r
3284 \r
3285 Showsq( arg )\r
3286 int arg;\r
3287 {\r
3288 \r
3289 k = 0;\r
3290 if( arg <= 0 )\r
3291         {\r
3292         if( cmp(SqEr,MinSqEr) < 0 )\r
3293                 {\r
3294                 k = 1;\r
3295                 mov( SqEr, MinSqEr );\r
3296                 }\r
3297         }\r
3298 if( arg >= 0 )\r
3299         {\r
3300         if( cmp(SqEr,MaxSqEr) > 0 )\r
3301                 {\r
3302                 k = 2;\r
3303                 mov( SqEr, MaxSqEr );\r
3304                 }\r
3305         }\r
3306 #if DEBUG\r
3307 if( k != 0 )\r
3308         {\r
3309         printf( "Square root of " );\r
3310         show( arg );\r
3311         printf( "\tis in error by " );\r
3312         show( SqEr );\r
3313         }\r
3314 #endif\r
3315 }\r
3316 \r
3317 \r
3318 pow1test()\r
3319 {\r
3320 \r
3321 /*=============================================*/\r
3322 Milestone = 90;\r
3323 /*=============================================*/\r
3324 Pause();\r
3325 printf("Testing powers Z^i for small Integers Z and i.\n");\r
3326 N = 0;\r
3327 /* ... test powers of zero. */\r
3328 I = 0;\r
3329 mov( Zero, Z );\r
3330 neg(Z);\r
3331 M = 3;\r
3332 Break = False;\r
3333 do\r
3334         {\r
3335         mov( One, X );\r
3336         SR3980();\r
3337         if(I <= 10)\r
3338                 {\r
3339                 I = 1023;\r
3340                 SR3980();\r
3341                 }\r
3342         if( cmp(Z,MinusOne) == 0 )\r
3343                 Break = True;\r
3344         else\r
3345                 {\r
3346                 mov( MinusOne, Z );\r
3347                 PrintIfNPositive();\r
3348                 N = 0;\r
3349 /* .. if(-1)^N is invalid, replace MinusOne by One. */\r
3350                 I = - 4;\r
3351                 }\r
3352         }\r
3353 while( ! Break );\r
3354 PrintIfNPositive();\r
3355 N1 = N;\r
3356 N = 0;\r
3357 mov( A1, Z );\r
3358 /*M = FLOOR(Two * LOG(W) / LOG(A1));*/\r
3359 LOG( W, t );\r
3360 mul( Two, t, t );\r
3361 FLOOR( t, t );\r
3362 LOG( A1, t2 );\r
3363 div( t2, t, t );\r
3364 FTOL( t, &lngint, t2 );\r
3365 M = lngint;\r
3366 Break = False;\r
3367 do\r
3368         {\r
3369         mov( Z, X );\r
3370         I = 1;\r
3371         SR3980();\r
3372         if( cmp(Z,AInvrse) == 0 )\r
3373                 Break = True;\r
3374         else\r
3375                  mov( AInvrse, Z );\r
3376         }\r
3377 while( ! (Break) );\r
3378 /*=============================================*/\r
3379 Milestone = 100;\r
3380 /*=============================================*/\r
3381 /*  Powers of Radix have been tested, */\r
3382 /*         next try a few primes     */\r
3383 \r
3384 M = NoTrials;\r
3385 \r
3386 mov( Three, Z );\r
3387 do\r
3388         {\r
3389         mov( Z, X );\r
3390         I = 1;\r
3391         SR3980();\r
3392         do\r
3393                 {\r
3394                 add( Z, Two, Z );\r
3395                 div( Three, Z, t );\r
3396                 FLOOR( t, t );\r
3397                 mul( Three, t, t );\r
3398                 }\r
3399         while( cmp(t,Z) == 0 );\r
3400         mul( Eight, Three, t );\r
3401         }\r
3402 while( cmp(Z,t) < 0 );\r
3403 \r
3404 if(N > 0)\r
3405         {\r
3406         printf("Errors like this may invalidate financial calculations\n");\r
3407         printf("\tinvolving interest rates.\n");\r
3408         }\r
3409 PrintIfNPositive();\r
3410 N += N1;\r
3411 if(N == 0)\r
3412         printf("... no discrepancies found.\n");\r
3413 if(N > 0)\r
3414         Pause();\r
3415 else printf("\n");\r
3416 }\r
3417 \r
3418 \r
3419 \r
3420 pow2test()\r
3421 {\r
3422 printf("\n");\r
3423 /* ...calculate Exp2 == exp(2) == 7.38905 60989 30650 22723 04275-... */\r
3424 mov( Zero, X );\r
3425 mov( Two, t2 ); /*I = 2;*/\r
3426 \r
3427 mul( Two, Three, Y );\r
3428 mov( Zero, Q );\r
3429 N = 0;\r
3430 do\r
3431         {\r
3432         mov( X, Z );\r
3433         add( t2, One, t2 ); /*I = I + 1;*/\r
3434         add( t2, t2, t );\r
3435         div( t, Y, Y ); /*Y = Y / (I + I);*/\r
3436         add( Y, Q, R );\r
3437         add( Z, R, X );\r
3438         sub( X, Z, Q );\r
3439         add( Q, R, Q );\r
3440         }\r
3441 while( cmp(X,Z) > 0 );\r
3442 \r
3443 /*Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);*/\r
3444 div( Eight, One, t );\r
3445 add( OneAndHalf, t, Z );\r
3446 mul( OneAndHalf, ThirtyTwo, t );\r
3447 div( t, X, t );\r
3448 add( Z, t, Z );\r
3449 mul( Z, Z, X );\r
3450 mul( X, X, Exp2 );\r
3451 mov( F9, X );\r
3452 sub( U1, X, Y );\r
3453 printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = " );\r
3454 show( Exp2 );\r
3455 printf( "\tas X -> 1.\n" );\r
3456 for(I = 1;;)\r
3457         {\r
3458         sub( BInvrse, X, Z );\r
3459 /*Z = (X + One) / (Z - (One - BInvrse));*/\r
3460         add( X, One, t2 );\r
3461         sub( BInvrse, One, t );\r
3462         sub( t, Z, t );\r
3463         div( t, t2, Z );\r
3464         POW( X, Z, Sqarg );\r
3465         sub( Exp2, Sqarg, Q );\r
3466         mov( Q, t );\r
3467         FABS( t );\r
3468         mul( TwoForty, U2, t2 );\r
3469         if( cmp( t, t2 ) > 0 )\r
3470                 {\r
3471                 N = 1;\r
3472                 sub( BInvrse, X, V9 );\r
3473                 sub( BInvrse, One, t );\r
3474                 sub( t, V9, V9 );\r
3475                 ErrCnt[Defect] += 1;\r
3476                 printf( "Calculated " );\r
3477                 show( Sqarg );\r
3478                 printf(" for \t(1 + " );\r
3479                 show( V9 );\r
3480                 printf( "\tto the power " );\r
3481                 show( Z );\r
3482                 printf("\tdiffers from correct value by " );\r
3483                 show( Q );\r
3484                 printf("\tThis much error may spoil financial\n");\r
3485                 printf("\tcalculations involving tiny interest rates.\n");\r
3486                 break;\r
3487                 }\r
3488         else\r
3489                 {\r
3490                 sub( X, Y, Z );\r
3491                 mul( Z, Two, Z );\r
3492                 add( Z, Y, Z );\r
3493                 mov( Y, X );\r
3494                 mov( Z, Y );\r
3495                 sub( F9, X, Z );\r
3496                 mul( Z, Z, Z );\r
3497                 add( Z, One, Z );\r
3498                 if( (cmp(Z,One) > 0) && (I < NoTrials) )\r
3499                         I++;\r
3500                 else\r
3501                         {\r
3502                         if( cmp(X,One) > 0 )\r
3503                                 {\r
3504                                 if(N == 0)\r
3505                                         printf("Accuracy seems adequate.\n");\r
3506                                 break;\r
3507                                 }\r
3508                         else\r
3509                                 {\r
3510                                 add( One, U2, X );\r
3511                                 add( U2, U2, Y );\r
3512                                 add( X, Y, Y );\r
3513                                 I = 1;\r
3514                                 }\r
3515                         }\r
3516                 }\r
3517         }\r
3518 /*=============================================*/\r
3519 Milestone = 150;\r
3520 /*=============================================*/\r
3521 printf("Testing powers Z^Q at four nearly extreme values.\n");\r
3522 N = 0;\r
3523 mov( A1, Z );\r
3524 /*Q = FLOOR(Half - LOG(C) / LOG(A1));*/\r
3525 LOG( C, t );\r
3526 LOG( A1, t2 );\r
3527 div( t2, t, t );\r
3528 sub( t, Half, t );\r
3529 FLOOR( t, Q );\r
3530 Break = False;\r
3531 do\r
3532         {\r
3533         mov( CInvrse, X );\r
3534         POW( Z, Q, Y );\r
3535         IsYeqX();\r
3536         neg(Q);\r
3537         mov( C, X );\r
3538         POW( Z, Q, Y );\r
3539         IsYeqX();\r
3540         if( cmp(Z,One) < 0 )\r
3541                 Break = True;\r
3542         else\r
3543                 mov( AInvrse, Z );\r
3544         }\r
3545 while( ! (Break));\r
3546 PrintIfNPositive();\r
3547 if(N == 0)\r
3548         printf(" ... no discrepancies found.\n");\r
3549 printf("\n");\r
3550 }\r