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

Subversion Repositories eco32

[/] [eco32/] [trunk/] [fp/] [tests/] [orig/] [paranoia.c] - Blame information for rev 54

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

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

powered by: WebSVN 2.1.0

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