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