OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [rtos/] [rtems/] [c/] [src/] [tests/] [samples/] [paranoia/] [paranoia.c] - Blame information for rev 593

Go to most recent revision | Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.