Return to paranoia.c CVS log

Up to [Research Unix] / researchv10no / cmd / lcc / tst

Annotation of researchv10no/cmd/lcc/tst/paranoia.c, revision 1.1

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

This archive runs on limited infrastructure. Preserving old code on modern bandwidth. Automated agents are requested to crawl responsibly.