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

4.104.115
Last change on this file since e1871b9 was e1871b9, checked in by Ralf Corsepius <ralf.corsepius@…>, on 11/24/08 at 16:39:04

2008-11-24 Ralf Corsépius <ralf.corsepius@…>

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