Context Navigation

source: rtems/testsuites/samples/paranoia/paranoia.c @ 3324383c

4.115

Last change on this file since 3324383c was 3324383c, checked in by Joel Sherrill <joel.sherrill@…>, on 05/05/14 at 14:47:30
testsuites: Remove BSP_SMALL_MEMORY
Property mode set to `100644`
File size: 69.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: