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