source: rtems/c/src/tests/samples/paranoia/paranoia.c @ ac7d5ef0

4.104.114.84.95
Last change on this file since ac7d5ef0 was ac7d5ef0, checked in by Joel Sherrill <joel.sherrill@…>, on 05/11/95 at 17:39:37

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