source: rtems/testsuites/samples/paranoia/paranoia.c @ 9b4422a2

4.11
Last change on this file since 9b4422a2 was 9b4422a2, checked in by Joel Sherrill <joel.sherrill@…>, on May 3, 2012 at 3:09:24 PM

Remove All CVS Id Strings Possible Using a Script

Script does what is expected and tries to do it as
smartly as possible.

+ remove occurrences of two blank comment lines

next to each other after Id string line removed.

+ remove entire comment blocks which only exited to

contain CVS Ids

+ If the processing left a blank line at the top of

a file, it was removed.

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