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

Subversion Repositories eco32

[/] [eco32/] [trunk/] [lcc/] [tst/] [paranoia.c] - Blame information for rev 4

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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