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

4.11
Last change on this file since 3324383c was 3324383c, checked in by Joel Sherrill <joel.sherrill@…>, on May 5, 2014 at 2:47:30 PM

testsuites: Remove BSP_SMALL_MEMORY

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