source: rtems/testsuites/samples/paranoia/paranoia.c @ f59df6b2

4.104.114.84.95
Last change on this file since f59df6b2 was 2936b425, checked in by Joel Sherrill <joel.sherrill@…>, on Jan 23, 1998 at 5:45:05 PM

Solaris port updates from Chris Johns

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