source: rtems/testsuites/samples/paranoia/paranoia.c @ 16b1599f

4.115
Last change on this file since 16b1599f was e313551, checked in by Ralf Corsepius <ralf.corsepius@…>, on 02/22/11 at 10:58:44

Add HAVE_CONFIG_H.

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