Source Code 1992 March/Source_Code_CD-ROM_Walnut_Creek_March_1992.iso/usenet/altsrcs/2/2334

Internet Message Format | 1990-12-28 | 54.9 KB

From: [e-mail] (Glenn Geers)
Newsgroups: alt.sources
Subject: Alternative 80386 math lib v2.0 part05/06
Message-ID: <[e-mail]>
Date: 16 Dec 90 21:07:40 GMT

Submitted-by: root@trantor
Archive-name: mathlib2.0/part05

---- Cut Here and feed the following to sh ----
#!/bin/sh
# this is mathlib.05 (part 5 of mathlib2.0)
# do not concatenate these parts, unpack them in order with /bin/sh
# file paranoia.c continued
#
if test ! -r _shar_seq_.tmp; then
    echo 'Please unpack part 1 first!'
    exit 1
fi
(read Scheck
 if test "$Scheck" != 5; then
    echo Please unpack part "$Scheck" next!
    exit 1
 else
    exit 0
 fi
) < _shar_seq_.tmp || exit 1
if test ! -f _shar_wnt_.tmp; then
    echo 'x - still skipping paranoia.c'
else
echo 'x - continuing file paranoia.c'
sed 's/^X//' << 'SHAR_EOF' >> 'paranoia.c' &&
X    X = OneAndHalf - U2;
X    Y = OneAndHalf + U2;
X    Z = (X - U2) * Y2;
X    T = Y * Y1;
X    Z = Z - X;
X    T = T - X;
X    X = X * Y2;
X    Y = (Y + U2) * Y1;
X    X = X - OneAndHalf;
X    Y = Y - OneAndHalf;
X    if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
X        X = (OneAndHalf + U2) * Y2;
X        Y = OneAndHalf - U2 - U2;
X        Z = OneAndHalf + U2 + U2;
X        T = (OneAndHalf - U2) * Y1;
X        X = X - (Z + U2);
X        StickyBit = Y * Y1;
X        S = Z * Y2;
X        T = T - Y;
X        Y = (U2 - Y) + StickyBit;
X        Z = S - (Z + U2 + U2);
X        StickyBit = (Y2 + U2) * Y1;
X        Y1 = Y2 * Y1;
X        StickyBit = StickyBit - Y2;
X        Y1 = Y1 - Half;
X        if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
X            && ( StickyBit == Zero) && (Y1 == Half)) {
X            RMult = Rounded;
X            printf("Multiplication appears to round correctly.\n");
X            }
X        else    if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
X                && (T < Zero) && (StickyBit + U2 == Zero)
X                && (Y1 < Half)) {
X                RMult = Chopped;
X                printf("Multiplication appears to chop.\n");
X                }
X            else printf("* is neither chopped nor correctly rounded.\n");
X        if ((RMult == Rounded) && (GMult == No)) notify("Multiplication");
X        }
X    else printf("* is neither chopped nor correctly rounded.\n");
X    /*=============================================*/
X    Milestone = 45;
X    /*=============================================*/
X    Y2 = One + U2;
X    Y1 = One - U2;
X    Z = OneAndHalf + U2 + U2;
X    X = Z / Y2;
X    T = OneAndHalf - U2 - U2;
X    Y = (T - U2) / Y1;
X    Z = (Z + U2) / Y2;
X    X = X - OneAndHalf;
X    Y = Y - T;
X    T = T / Y1;
X    Z = Z - (OneAndHalf + U2);
X    T = (U2 - OneAndHalf) + T;
X    if (! ((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
X        X = OneAndHalf / Y2;
X        Y = OneAndHalf - U2;
X        Z = OneAndHalf + U2;
X        X = X - Y;
X        T = OneAndHalf / Y1;
X        Y = Y / Y1;
X        T = T - (Z + U2);
X        Y = Y - Z;
X        Z = Z / Y2;
X        Y1 = (Y2 + U2) / Y2;
X        Z = Z - OneAndHalf;
X        Y2 = Y1 - Y2;
X        Y1 = (F9 - U1) / F9;
X        if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
X            && (Y2 == Zero) && (Y2 == Zero)
X            && (Y1 - Half == F9 - Half )) {
X            RDiv = Rounded;
X            printf("Division appears to round correctly.\n");
X            if (GDiv == No) notify("Division");
X            }
X        else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
X            && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
X            RDiv = Chopped;
X            printf("Division appears to chop.\n");
X            }
X        }
X    if (RDiv == Other) printf("/ is neither chopped nor correctly rounded.\n");
X    BInvrse = One / Radix;
X    TstCond (Failure, (BInvrse * Radix - Half == Half),
X           "Radix * ( 1 / Radix ) differs from 1");
X    /*=============================================*/
X    /*SPLIT
X    }
#include "paranoia.h"
part4(){
*/
X    Milestone = 50;
X    /*=============================================*/
X    TstCond (Failure, ((F9 + U1) - Half == Half)
X           && ((BMinusU2 + U2 ) - One == Radix - One),
X           "Incomplete carry-propagation in Addition");
X    X = One - U1 * U1;
X    Y = One + U2 * (One - U2);
X    Z = F9 - Half;
X    X = (X - Half) - Z;
X    Y = Y - One;
X    if ((X == Zero) && (Y == Zero)) {
X        RAddSub = Chopped;
X        printf("Add/Subtract appears to be chopped.\n");
X        }
X    if (GAddSub == Yes) {
X        X = (Half + U2) * U2;
X        Y = (Half - U2) * U2;
X        X = One + X;
X        Y = One + Y;
X        X = (One + U2) - X;
X        Y = One - Y;
X        if ((X == Zero) && (Y == Zero)) {
X            X = (Half + U2) * U1;
X            Y = (Half - U2) * U1;
X            X = One - X;
X            Y = One - Y;
X            X = F9 - X;
X            Y = One - Y;
X            if ((X == Zero) && (Y == Zero)) {
X                RAddSub = Rounded;
X                printf("Addition/Subtraction appears to round correctly.\n");
X                if (GAddSub == No) notify("Add/Subtract");
X                }
X            else printf("Addition/Subtraction neither rounds nor chops.\n");
X            }
X        else printf("Addition/Subtraction neither rounds nor chops.\n");
X        }
X    else printf("Addition/Subtraction neither rounds nor chops.\n");
X    S = One;
X    X = One + Half * (One + Half);
X    Y = (One + U2) * Half;
X    Z = X - Y;
X    T = Y - X;
X    StickyBit = Z + T;
X    if (StickyBit != Zero) {
X        S = Zero;
X        BadCond(Flaw, "(X - Y) + (Y - X) is non zero!\n");
X        }
X    StickyBit = Zero;
X    if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
X        && (RMult == Rounded) && (RDiv == Rounded)
X        && (RAddSub == Rounded) && (FLOOR(RadixD2) == RadixD2)) {
X        printf("Checking for sticky bit.\n");
X        X = (Half + U1) * U2;
X        Y = Half * U2;
X        Z = One + Y;
X        T = One + X;
X        if ((Z - One <= Zero) && (T - One >= U2)) {
X            Z = T + Y;
X            Y = Z - X;
X            if ((Z - T >= U2) && (Y - T == Zero)) {
X                X = (Half + U1) * U1;
X                Y = Half * U1;
X                Z = One - Y;
X                T = One - X;
X                if ((Z - One == Zero) && (T - F9 == Zero)) {
X                    Z = (Half - U1) * U1;
X                    T = F9 - Z;
X                    Q = F9 - Y;
X                    if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
X                        Z = (One + U2) * OneAndHalf;
X                        T = (OneAndHalf + U2) - Z + U2;
X                        X = One + Half / Radix;
X                        Y = One + Radix * U2;
X                        Z = X * Y;
X                        if (T == Zero && X + Radix * U2 - Z == Zero) {
X                            if (Radix != Two) {
X                                X = Two + U2;
X                                Y = X / Two;
X                                if ((Y - One == Zero)) StickyBit = S;
X                                }
X                            else StickyBit = S;
X                            }
X                        }
X                    }
X                }
X            }
X        }
X    if (StickyBit == One) printf("Sticky bit apparently used correctly.\n");
X    else printf("Sticky bit used incorrectly or not at all.\n");
X    TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
X            RMult == Other || RDiv == Other || RAddSub == Other),
X        "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");
X    /*=============================================*/
X    Milestone = 60;
X    /*=============================================*/
X    printf("\n");
X    printf("Does Multiplication commute?  ");
X    printf("Testing on %d random pairs.\n", NoTrials);
X    Random9 = SQRT(3.0);
X    Random1 = Third;
X    I = 1;
X    do  {
X        X = Random();
X        Y = Random();
X        Z9 = Y * X;
X        Z = X * Y;
X        Z9 = Z - Z9;
X        I = I + 1;
X        } while ( ! ((I > NoTrials) || (Z9 != Zero)));
X    if (I == NoTrials) {
X        Random1 = One + Half / Three;
X        Random2 = (U2 + U1) + One;
X        Z = Random1 * Random2;
X        Y = Random2 * Random1;
X        Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
X            Three) * ((U2 + U1) + One);
X        }
X    if (! ((I == NoTrials) || (Z9 == Zero)))
X        BadCond(Defect, "X * Y == Y * X trial fails.\n");
X    else printf("     No failures found in %d integer pairs.\n", NoTrials);
X    /*=============================================*/
X    Milestone = 70;
X    /*=============================================*/
X    printf("\nRunning test of square root(x).\n");
X    TstCond (Failure, (Zero == SQRT(Zero))
X           && (- Zero == SQRT(- Zero))
X           && (One == SQRT(One)), "Square root of 0.0, -0.0 or 1.0 wrong");
X    MinSqEr = Zero;
X    MaxSqEr = Zero;
X    J = Zero;
X    X = Radix;
X    OneUlp = U2;
X    SqXMinX (Serious);
X    X = BInvrse;
X    OneUlp = BInvrse * U1;
X    SqXMinX (Serious);
X    X = U1;
X    OneUlp = U1 * U1;
X    SqXMinX (Serious);
X    if (J != Zero) Pause();
X    printf("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
X    J = Zero;
X    X = Two;
X    Y = Radix;
X    if ((Radix != One)) do  {
X        X = Y;
X        Y = Radix * Y;
X        } while ( ! ((Y - X >= NoTrials)));
X    OneUlp = X * U2;
X    I = 1;
X    while (I <= NoTrials) {
X        X = X + One;
X        SqXMinX (Defect);
X        if (J > Zero) break;
X        I = I + 1;
X        }
X    printf("Test for sqrt monotonicity.\n");
X    I = - 1;
X    X = BMinusU2;
X    Y = Radix;
X    Z = Radix + Radix * U2;
X    NotMonot = False;
X    Monot = False;
X    while ( ! (NotMonot || Monot)) {
X        I = I + 1;
X        X = SQRT(X);
X        Q = SQRT(Y);
X        Z = SQRT(Z);
X        if ((X > Q) || (Q > Z)) NotMonot = True;
X        else {
X            Q = FLOOR(Q + Half);
X            if ((I > 0) || (Radix == Q * Q)) Monot = True;
X            else if (I > 0) {
X            if (I > 1) Monot = True;
X            else {
X                Y = Y * BInvrse;
X                X = Y - U1;
X                Z = Y + U1;
X                }
X            }
X            else {
X                Y = Q;
X                X = Y - U2;
X                Z = Y + U2;
X                }
X            }
X        }
X    if (Monot) printf("sqrt has passed a test for Monotonicity.\n");
X    else {
X        BadCond(Defect, "");
X        printf("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
X        }
X    /*=============================================*/
X    /*SPLIT
X    }
#include "paranoia.h"
part5(){
*/
X    Milestone = 80;
X    /*=============================================*/
X    MinSqEr = MinSqEr + Half;
X    MaxSqEr = MaxSqEr - Half;
X    Y = (SQRT(One + U2) - One) / U2;
X    SqEr = (Y - One) + U2 / Eight;
X    if (SqEr > MaxSqEr) MaxSqEr = SqEr;
X    SqEr = Y + U2 / Eight;
X    if (SqEr < MinSqEr) MinSqEr = SqEr;
X    Y = ((SQRT(F9) - U2) - (One - U2)) / U1;
X    SqEr = Y + U1 / Eight;
X    if (SqEr > MaxSqEr) MaxSqEr = SqEr;
X    SqEr = (Y + One) + U1 / Eight;
X    if (SqEr < MinSqEr) MinSqEr = SqEr;
X    OneUlp = U2;
X    X = OneUlp;
X    for( Indx = 1; Indx <= 3; ++Indx) {
X        Y = SQRT((X + U1 + X) + F9);
X        Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
X        Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
X        SqEr = (Y + Half) + Z;
X        if (SqEr < MinSqEr) MinSqEr = SqEr;
X        SqEr = (Y - Half) + Z;
X        if (SqEr > MaxSqEr) MaxSqEr = SqEr;
X        if (((Indx == 1) || (Indx == 3)))
X            X = OneUlp * Sign (X) * FLOOR(Eight / (Nine * SQRT(OneUlp)));
X        else {
X            OneUlp = U1;
X            X = - OneUlp;
X            }
X        }
X    /*=============================================*/
X    Milestone = 85;
X    /*=============================================*/
X    SqRWrng = False;
X    Anomaly = False;
X    RSqrt = Other; /* ~dgh */
X    if (Radix != One) {
X        printf("Testing whether sqrt is rounded or chopped.\n");
X        D = FLOOR(Half + POW(Radix, One + Precision - FLOOR(Precision)));
X    /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
X        X = D / Radix;
X        Y = D / A1;
X        if ((X != FLOOR(X)) || (Y != FLOOR(Y))) {
X            Anomaly = True;
X            }
X        else {
X            X = Zero;
X            Z2 = X;
X            Y = One;
X            Y2 = Y;
X            Z1 = Radix - One;
X            FourD = Four * D;
X            do  {
X                if (Y2 > Z2) {
X                    Q = Radix;
X                    Y1 = Y;
X                    do  {
X                        X1 = FABS(Q + FLOOR(Half - Q / Y1) * Y1);
X                        Q = Y1;
X                        Y1 = X1;
X                        } while ( ! (X1 <= Zero));
X                    if (Q <= One) {
X                        Z2 = Y2;
X                        Z = Y;
X                        }
X                    }
X                Y = Y + Two;
X                X = X + Eight;
X                Y2 = Y2 + X;
X                if (Y2 >= FourD) Y2 = Y2 - FourD;
X                } while ( ! (Y >= D));
X            X8 = FourD - Z2;
X            Q = (X8 + Z * Z) / FourD;
X            X8 = X8 / Eight;
X            if (Q != FLOOR(Q)) Anomaly = True;
X            else {
X                Break = False;
X                do  {
X                    X = Z1 * Z;
X                    X = X - FLOOR(X / Radix) * Radix;
X                    if (X == One)
X                        Break = True;
X                    else
X                        Z1 = Z1 - One;
X                    } while ( ! (Break || (Z1 <= Zero)));
X                if ((Z1 <= Zero) && (! Break)) Anomaly = True;
X                else {
X                    if (Z1 > RadixD2) Z1 = Z1 - Radix;
X                    do  {
X                        NewD();
X                        } while ( ! (U2 * D >= F9));
X                    if (D * Radix - D != W - D) Anomaly = True;
X                    else {
X                        Z2 = D;
X                        I = 0;
X                        Y = D + (One + Z) * Half;
X                        X = D + Z + Q;
X                        SR3750();
X                        Y = D + (One - Z) * Half + D;
X                        X = D - Z + D;
X                        X = X + Q + X;
X                        SR3750();
X                        NewD();
X                        if (D - Z2 != W - Z2) Anomaly = True;
X                        else {
X                            Y = (D - Z2) + (Z2 + (One - Z) * Half);
X                            X = (D - Z2) + (Z2 - Z + Q);
X                            SR3750();
X                            Y = (One + Z) * Half;
X                            X = Q;
X                            SR3750();
X                            if (I == 0) Anomaly = True;
X                            }
X                        }
X                    }
X                }
X            }
X        if ((I == 0) || Anomaly) {
X            BadCond(Failure, "Anomalous arithmetic with Integer < ");
X            printf("Radix^Precision = %.7e\n", W);
X            printf(" fails test whether sqrt rounds or chops.\n");
X            SqRWrng = True;
X            }
X        }
X    if (! Anomaly) {
X        if (! ((MinSqEr < Zero) || (MaxSqEr > Zero))) {
X            RSqrt = Rounded;
X            printf("Square root appears to be correctly rounded.\n");
X            }
X        else  {
X            if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
X                || (MinSqEr + Radix < Half)) SqRWrng = True;
X            else {
X                RSqrt = Chopped;
X                printf("Square root appears to be chopped.\n");
X                }
X            }
X        }
X    if (SqRWrng) {
X        printf("Square root is neither chopped nor correctly rounded.\n");
X        printf("Observed errors run from %.7e ", MinSqEr - Half);
X        printf("to %.7e ulps.\n", Half + MaxSqEr);
X        TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
X            "sqrt gets too many last digits wrong");
X        }
X    /*=============================================*/
X    Milestone = 90;
X    /*=============================================*/
X    Pause();
X    printf("Testing powers Z^i for small Integers Z and i.\n");
X    N = 0;
X    /* ... test powers of zero. */
X    I = 0;
X    Z = -Zero;
X    M = 3.0;
X    Break = False;
X    do  {
X        X = One;
X        SR3980();
X        if (I <= 10) {
X            I = 1023;
X            SR3980();
X            }
X        if (Z == MinusOne) Break = True;
X        else {
X            Z = MinusOne;
X            PrintIfNPositive();
X            N = 0;
X            /* .. if(-1)^N is invalid, replace MinusOne by One. */
X            I = - 4;
X            }
X        } while ( ! Break);
X    PrintIfNPositive();
X    N1 = N;
X    N = 0;
X    Z = A1;
X    M = FLOOR(Two * LOG(W) / LOG(A1));
X    Break = False;
X    do  {
X        X = Z;
X        I = 1;
X        SR3980();
X        if (Z == AInvrse) Break = True;
X        else Z = AInvrse;
X        } while ( ! (Break));
X    /*=============================================*/
X        Milestone = 100;
X    /*=============================================*/
X    /*  Powers of Radix have been tested, */
X    /*         next try a few primes     */
X    M = NoTrials;
X    Z = Three;
X    do  {
X        X = Z;
X        I = 1;
X        SR3980();
X        do  {
X            Z = Z + Two;
X            } while ( Three * FLOOR(Z / Three) == Z );
X        } while ( Z < Eight * Three );
X    if (N > 0) {
X        printf("Errors like this may invalidate financial calculations\n");
X        printf("\tinvolving interest rates.\n");
X        }
X    PrintIfNPositive();
X    N += N1;
X    if (N == 0) printf("... no discrepancis found.\n");
X    if (N > 0) Pause();
X    else printf("\n");
X    /*=============================================*/
X    /*SPLIT
X    }
#include "paranoia.h"
part6(){
*/
X    Milestone = 110;
X    /*=============================================*/
X    printf("Seeking Underflow thresholds UfThold and E0.\n");
X    D = U1;
X    if (Precision != FLOOR(Precision)) {
X        D = BInvrse;
X        X = Precision;
X        do  {
X            D = D * BInvrse;
X            X = X - One;
X            } while ( X > Zero);
X        }
X    Y = One;
X    Z = D;
X    /* ... D is power of 1/Radix < 1. */
X    do  {
X        C = Y;
X        Y = Z;
X        Z = Y * Y;
X        } while ((Y > Z) && (Z + Z > Z));
X    Y = C;
X    Z = Y * D;
X    do  {
X        C = Y;
X        Y = Z;
X        Z = Y * D;
X        } while ((Y > Z) && (Z + Z > Z));
X    if (Radix < Two) HInvrse = Two;
X    else HInvrse = Radix;
X    H = One / HInvrse;
X    /* ... 1/HInvrse == H == Min(1/Radix, 1/2) */
X    CInvrse = One / C;
X    E0 = C;
X    Z = E0 * H;
X    /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
X    do  {
X        Y = E0;
X        E0 = Z;
X        Z = E0 * H;
X        } while ((E0 > Z) && (Z + Z > Z));
X    UfThold = E0;
X    E1 = Zero;
X    Q = Zero;
X    E9 = U2;
X    S = One + E9;
X    D = C * S;
X    if (D <= C) {
X        E9 = Radix * U2;
X        S = One + E9;
X        D = C * S;
X        if (D <= C) {
X            BadCond(Failure, "multiplication gets too many last digits wrong.\n");
X            Underflow = E0;
X            Y1 = Zero;
X            PseudoZero = Z;
X            Pause();
X            }
X        }
X    else {
X        Underflow = D;
X        PseudoZero = Underflow * H;
X        UfThold = Zero;
X        do  {
X            Y1 = Underflow;
X            Underflow = PseudoZero;
X            if (E1 + E1 <= E1) {
X                Y2 = Underflow * HInvrse;
X                E1 = FABS(Y1 - Y2);
X                Q = Y1;
X                if ((UfThold == Zero) && (Y1 != Y2)) UfThold = Y1;
X                }
X            PseudoZero = PseudoZero * H;
X            } while ((Underflow > PseudoZero)
X                && (PseudoZero + PseudoZero > PseudoZero));
X        }
X    /* Comment line 4530 .. 4560 */
X    if (PseudoZero != Zero) {
X        printf("\n");
X        Z = PseudoZero;
X    /* ... Test PseudoZero for "phoney- zero" violates */
X    /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
X           ... */
X        if (PseudoZero <= Zero) {
X            BadCond(Failure, "Positive expressions can underflow to an\n");
X            printf("allegedly negative value\n");
X            printf("PseudoZero that prints out as: %g .\n", PseudoZero);
X            X = - PseudoZero;
X            if (X <= Zero) {
X                printf("But -PseudoZero, which should be\n");
X                printf("positive, isn't; it prints out as  %g .\n", X);
X                }
X            }
X        else {
X            BadCond(Flaw, "Underflow can stick at an allegedly positive\n");
X            printf("value PseudoZero that prints out as %g .\n", PseudoZero);
X            }
X        TstPtUf();
X        }
X    /*=============================================*/
X    Milestone = 120;
X    /*=============================================*/
X    if (CInvrse * Y > CInvrse * Y1) {
X        S = H * S;
X        E0 = Underflow;
X        }
X    if (! ((E1 == Zero) || (E1 == E0))) {
X        BadCond(Defect, "");
X        if (E1 < E0) {
X            printf("Products underflow at a higher");
X            printf(" threshold than differences.\n");
X            if (PseudoZero == Zero)
X            E0 = E1;
X            }
X        else {
X            printf("Difference underflows at a higher");
X            printf(" threshold than products.\n");
X            }
X        }
X    printf("Smallest strictly positive number found is E0 = %g .\n", E0);
X    Z = E0;
X    TstPtUf();
X    Underflow = E0;
X    if (N == 1) Underflow = Y;
X    I = 4;
X    if (E1 == Zero) I = 3;
X    if (UfThold == Zero) I = I - 2;
X    UfNGrad = True;
X    switch (I)  {
X        case    1:
X        UfThold = Underflow;
X        if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
X            UfThold = Y;
X            BadCond(Failure, "Either accuracy deteriorates as numbers\n");
X            printf("approach a threshold = %.17e\n", UfThold);;
X            printf(" coming down from %.17e\n", C);
X            printf(" or else multiplication gets too many last digits wrong.\n");
X            }
X        Pause();
X        break;
X
X        case    2:
X        BadCond(Failure, "Underflow confuses Comparison, which alleges that\n");
X        printf("Q == Y while denying that |Q - Y| == 0; these values\n");
X        printf("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
X        printf ("|Q - Y| = %.17e .\n" , FABS(Q - Y2));
X        UfThold = Q;
X        break;
X
X        case    3:
X        X = X;
X        break;
X
X        case    4:
X        if ((Q == UfThold) && (E1 == E0)
X            && (FABS( UfThold - E1 / E9) <= E1)) {
X            UfNGrad = False;
X            printf("Underflow is gradual; it incurs Absolute Error =\n");
X            printf("(roundoff in UfThold) < E0.\n");
X            Y = E0 * CInvrse;
X            Y = Y * (OneAndHalf + U2);
X            X = CInvrse * (One + U2);
X            Y = Y / X;
X            IEEE = (Y == E0);
X            }
X        }
X    if (UfNGrad) {
X        printf("\n");
X        sigsave = sigfpe;
X        if (setjmp(ovfl_buf)) {
X            printf("Underflow / UfThold failed!\n");
X            R = H + H;
X            }
X        else R = SQRT(Underflow / UfThold);
X        sigsave = 0;
X        if (R <= H) {
X            Z = R * UfThold;
X            X = Z * (One + R * H * (One + H));
X            }
X        else {
X            Z = UfThold;
X            X = Z * (One + H * H * (One + H));
X            }
X        if (! ((X == Z) || (X - Z != Zero))) {
X            BadCond(Flaw, "");
X            printf("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
X            Z9 = X - Z;
X            printf("yet X - Z yields %.17e .\n", Z9);
X            printf("    Should this NOT signal Underflow, ");
X            printf("this is a SERIOUS DEFECT\nthat causes ");
X            printf("confusion when innocent statements like\n");;
X            printf("    if (X == Z)  ...  else");
X            printf("  ... (f(X) - f(Z)) / (X - Z) ...\n");
X            printf("encounter Division by Zero although actually\n");
X            sigsave = sigfpe;
X            if (setjmp(ovfl_buf)) printf("X / Z fails!\n");
X            else printf("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
X            sigsave = 0;
X            }
X        }
X    printf("The Underflow threshold is %.17e, %s\n", UfThold,
X           " below which");
X    printf("calculation may suffer larger Relative error than ");
X    printf("merely roundoff.\n");
X    Y2 = U1 * U1;
X    Y = Y2 * Y2;
X    Y2 = Y * U1;
X    if (Y2 <= UfThold) {
X        if (Y > E0) {
X            BadCond(Defect, "");
X            I = 5;
X            }
X        else {
X            BadCond(Serious, "");
X            I = 4;
X            }
X        printf("Range is too narrow; U1^%d Underflows.\n", I);
X        }
X    /*=============================================*/
X    /*SPLIT
X    }
#include "paranoia.h"
part7(){
*/
X    Milestone = 130;
X    /*=============================================*/
X    Y = - FLOOR(Half - TwoForty * LOG(UfThold) / LOG(HInvrse)) / TwoForty;
X    Y2 = Y + Y;
X    printf("Since underflow occurs below the threshold\n");
X    printf("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
X    printf("should afflict the expression\n\t(%.17e) ^ (%.17e);\n", HInvrse, Y);
X    V9 = POW(HInvrse, Y2);
X    printf("actually calculating yields: %.17e .\n", V9);
X    if (! ((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
X        BadCond(Serious, "this is not between 0 and underflow\n");
X        printf("   threshold = %.17e .\n", UfThold);
X        }
X    else if (! (V9 > UfThold * (One + E9)))
X        printf("This computed value is O.K.\n");
X    else {
X        BadCond(Defect, "this is not between 0 and underflow\n");
X        printf("   threshold = %.17e .\n", UfThold);
X        }
X    /*=============================================*/
X    Milestone = 140;
X    /*=============================================*/
X    printf("\n");
X    /* ...calculate Exp2 == exp(2) == 7.389056099... */
X    X = Zero;
X    I = 2;
X    Y = Two * Three;
X    Q = Zero;
X    N = 0;
X    do  {
X        Z = X;
X        I = I + 1;
X        Y = Y / (I + I);
X        R = Y + Q;
X        X = Z + R;
X        Q = (Z - X) + R;
X        } while(X > Z);
X    Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
X    X = Z * Z;
X    Exp2 = X * X;
X    X = F9;
X    Y = X - U1;
X    printf("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
X        Exp2);
X    for(I = 1;;) {
X        Z = X - BInvrse;
X        Z = (X + One) / (Z - (One - BInvrse));
X        Q = POW(X, Z) - Exp2;
X        if (FABS(Q) > TwoForty * U2) {
X            N = 1;
X             V9 = (X - BInvrse) - (One - BInvrse);
X            BadCond(Defect, "Calculated");
X            printf(" %.17e for\n", POW(X,Z));
X            printf("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
X            printf("\tdiffers from correct value by %.17e .\n", Q);
X            printf("\tThis much error may spoil financial\n");
X            printf("\tcalculations involving tiny interest rates.\n");
X            break;
X            }
X        else {
X            Z = (Y - X) * Two + Y;
X            X = Y;
X            Y = Z;
X            Z = One + (X - F9)*(X - F9);
X            if (Z > One && I < NoTrials) I++;
X            else  {
X                if (X > One) {
X                    if (N == 0)
X                       printf("Accuracy seems adequate.\n");
X                    break;
X                    }
X                else {
X                    X = One + U2;
X                    Y = U2 + U2;
X                    Y += X;
X                    I = 1;
X                    }
X                }
X            }
X        }
X    /*=============================================*/
X    Milestone = 150;
X    /*=============================================*/
X    printf("Testing powers Z^Q at four nearly extreme values.\n");
X    N = 0;
X    Z = A1;
X    Q = FLOOR(Half - LOG(C) / LOG(A1));
X    Break = False;
X    do  {
X        X = CInvrse;
X        Y = POW(Z, Q);
X        IsYeqX();
X        Q = - Q;
X        X = C;
X        Y = POW(Z, Q);
X        IsYeqX();
X        if (Z < One) Break = True;
X        else Z = AInvrse;
X        } while ( ! (Break));
X    PrintIfNPositive();
X    if (N == 0) printf(" ... no discrepancies found.\n");
X    printf("\n");
X
X    /*=============================================*/
X    Milestone = 160;
X    /*=============================================*/
X    Pause();
X    printf("Searching for Overflow threshold:\n");
X    printf("This may generate an error.\n");
X    Y = - CInvrse;
X    V9 = HInvrse * Y;
X    sigsave = sigfpe;
X    if (setjmp(ovfl_buf)) { I = 0; V9 = Y; goto overflow; }
X    do {
X        V = Y;
X        Y = V9;
X        V9 = HInvrse * Y;
X        } while(V9 < Y);
X    I = 1;
overflow:
X    sigsave = 0;
X    Z = V9;
X    printf("Can `Z = -Y' overflow?\n");
X    printf("Trying it on Y = %.17e .\n", Y);
X    V9 = - Y;
X    V0 = V9;
X    if (V - Y == V + V0) printf("Seems O.K.\n");
X    else {
X        printf("finds a ");
X        BadCond(Flaw, "-(-Y) differs from Y.\n");
X        }
X    if (Z != Y) {
X        BadCond(Serious, "");
X        printf("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
X        }
X    if (I) {
X        Y = V * (HInvrse * U2 - HInvrse);
X        Z = Y + ((One - HInvrse) * U2) * V;
X        if (Z < V0) Y = Z;
X        if (Y < V0) V = Y;
X        if (V0 - V < V0) V = V0;
X        }
X    else {
X        V = Y * (HInvrse * U2 - HInvrse);
X        V = V + ((One - HInvrse) * U2) * Y;
X        }
X    printf("Overflow threshold is V  = %.17e .\n", V);
X    if (I) printf("Overflow saturates at V0 = %.17e .\n", V0);
X    else printf("There is no saturation value because \
the system traps on overflow.\n");
X    V9 = V * One;
X    printf("No Overflow should be signaled for V * 1 = %.17e\n", V9);
X    V9 = V / One;
X    printf("                           nor for V / 1 = %.17e .\n", V9);
X    printf("Any overflow signal separating this * from the one\n");
X    printf("above is a DEFECT.\n");
X    /*=============================================*/
X    Milestone = 170;
X    /*=============================================*/
X    if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
X        BadCond(Failure, "Comparisons involving ");
X        printf("+-%g, +-%g\nand +-%g are confused by Overflow.",
X            V, V0, UfThold);
X        }
X    /*=============================================*/
X    Milestone = 175;
X    /*=============================================*/
X    printf("\n");
X    for(Indx = 1; Indx <= 3; ++Indx) {
X        switch (Indx)  {
X            case 1: Z = UfThold; break;
X            case 2: Z = E0; break;
X            case 3: Z = PseudoZero; break;
X            }
X        if (Z != Zero) {
X            V9 = SQRT(Z);
X            Y = V9 * V9;
X            if (Y / (One - Radix * E9) < Z
X               || Y > (One + Radix * E9) * Z) { /* dgh: + E9 --> * E9 */
X                if (V9 > U1) BadCond(Serious, "");
X                else BadCond(Defect, "");
X                printf("Comparison alleges that what prints as Z = %.17e\n", Z);
X                printf(" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
X                }
X            }
X        }
X    /*=============================================*/
X    Milestone = 180;
X    /*=============================================*/
X    for(Indx = 1; Indx <= 2; ++Indx) {
X        if (Indx == 1) Z = V;
X        else Z = V0;
X        V9 = SQRT(Z);
X        X = (One - Radix * E9) * V9;
X        V9 = V9 * X;
X        if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
X            Y = V9;
X            if (X < W) BadCond(Serious, "");
X            else BadCond(Defect, "");
X            printf("Comparison alleges that Z = %17e\n", Z);
X            printf(" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
X            }
X        }
X    /*=============================================*/
X    /*SPLIT
X    }
#include "paranoia.h"
part8(){
*/
X    Milestone = 190;
X    /*=============================================*/
X    Pause();
X    X = UfThold * V;
X    Y = Radix * Radix;
X    if (X*Y < One || X > Y) {
X        if (X * Y < U1 || X > Y/U1) BadCond(Defect, "Badly");
X        else BadCond(Flaw, "");
X
X        printf(" unbalanced range; UfThold * V = %.17e\n\t%s\n",
X            X, "is too far from 1.\n");
X        }
X    /*=============================================*/
X    Milestone = 200;
X    /*=============================================*/
X    for (Indx = 1; Indx <= 5; ++Indx)  {
X        X = F9;
X        switch (Indx)  {
X            case 2: X = One + U2; break;
X            case 3: X = V; break;
X            case 4: X = UfThold; break;
X            case 5: X = Radix;
X            }
X        Y = X;
X        sigsave = sigfpe;
X        if (setjmp(ovfl_buf))
X            printf("  X / X  traps when X = %g\n", X);
X        else {
X            V9 = (Y / X - Half) - Half;
X            if (V9 == Zero) continue;
X            if (V9 == - U1 && Indx < 5) BadCond(Flaw, "");
X            else BadCond(Serious, "");
X            printf("  X / X differs from 1 when X = %.17e\n", X);
X            printf("  instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
X            }
X        sigsave = 0;
X        }
X    /*=============================================*/
X    Milestone = 210;
X    /*=============================================*/
X    MyZero = Zero;
X    printf("\n");
X    printf("What message and/or values does Division by Zero produce?\n") ;
#ifndef NOPAUSE
X    printf("This can interupt your program.  You can ");
X    printf("skip this part if you wish.\n");
X    printf("Do you wish to compute 1 / 0? ");
X    fflush(stdout);
X    read (KEYBOARD, ch, 8);
X    if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif
X        sigsave = sigfpe;
X        printf("    Trying to compute 1 / 0 produces ...");
X        if (!setjmp(ovfl_buf)) printf("  %.7e .\n", One / MyZero);
X        sigsave = 0;
#ifndef NOPAUSE
X        }
X    else printf("O.K.\n");
X    printf("\nDo you wish to compute 0 / 0? ");
X    fflush(stdout);
X    read (KEYBOARD, ch, 80);
X    if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif
X        sigsave = sigfpe;
X        printf("\n    Trying to compute 0 / 0 produces ...");
X        if (!setjmp(ovfl_buf)) printf("  %.7e .\n", Zero / MyZero);
X        sigsave = 0;
#ifndef NOPAUSE
X        }
X    else printf("O.K.\n");
#endif
X    /*=============================================*/
X    Milestone = 220;
X    /*=============================================*/
X    Pause();
X    printf("\n");
X    {
X        static char *msg[] = {
X            "FAILUREs  encountered =",
X            "SERIOUS DEFECTs  discovered =",
X            "DEFECTs  discovered =",
X            "FLAWs  discovered =" };
X        int i;
X        for(i = 0; i < 4; i++) if (ErrCnt[i])
X            printf("The number of  %-29s %d.\n",
X                msg[i], ErrCnt[i]);
X        }
X    printf("\n");
X    if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
X            + ErrCnt[Flaw]) > 0) {
X        if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
X            Defect] == 0) && (ErrCnt[Flaw] > 0)) {
X            printf("The arithmetic diagnosed seems ");
X            printf("Satisfactory though flawed.\n");
X            }
X        if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
X            && ( ErrCnt[Defect] > 0)) {
X            printf("The arithmetic diagnosed may be Acceptable\n");
X            printf("despite inconvenient Defects.\n");
X            }
X        if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
X            printf("The arithmetic diagnosed has ");
X            printf("unacceptable Serious Defects.\n");
X            }
X        if (ErrCnt[Failure] > 0) {
X            printf("Potentially fatal FAILURE may have spoiled this");
X            printf(" program's subsequent diagnoses.\n");
X            }
X        }
X    else {
X        printf("No failures, defects nor flaws have been discovered.\n");
X        if (! ((RMult == Rounded) && (RDiv == Rounded)
X            && (RAddSub == Rounded) && (RSqrt == Rounded)))
X            printf("The arithmetic diagnosed seems Satisfactory.\n");
X        else {
X            if (StickyBit >= One &&
X                (Radix - Two) * (Radix - Nine - One) == Zero) {
X                printf("Rounding appears to conform to ");
X                printf("the proposed IEEE standard P");
X                if ((Radix == Two) &&
X                     ((Precision - Four * Three * Two) *
X                      ( Precision - TwentySeven -
X                       TwentySeven + One) == Zero))
X                    printf("754");
X                else printf("854");
X                if (IEEE) printf(".\n");
X                else {
X                    printf(",\nexcept for possibly Double Rounding");
X                    printf(" during Gradual Underflow.\n");
X                    }
X                }
X            printf("The arithmetic diagnosed appears to be Excellent!\n");
X            }
X        }
X    if (fpecount)
X        printf("\nA total of %d floating point exceptions were registered.\n",
X            fpecount);
X    printf("END OF TEST.\n");
X
#ifdef TEST
X    ieee_retrospective((FILE *)NULL);
#endif
X    asm("fnclex");
X        return(0);
X    }
X
/*SPLIT subs.c
#include "paranoia.h"
*/
X
/* Sign */
X
FLOAT Sign (X)
FLOAT X;
{ return X >= 0. ? 1.0 : -1.0; }
X
/* Pause */
X
Pause()
{
#ifndef NOPAUSE
X    char ch[8];
X
X    printf("\nTo continue, press RETURN");
X    fflush(stdout);
X    read(KEYBOARD, ch, 8);
#endif
X    printf("\nDiagnosis resumes after milestone Number %d", Milestone);
X    printf("          Page: %d\n\n", PageNo);
X    ++Milestone;
X    ++PageNo;
X    }
X
X /* TstCond */
X
TstCond (K, Valid, T)
int K, Valid;
char *T;
{ if (! Valid) { BadCond(K,T); printf(".\n"); } }
X
BadCond(K, T)
int K;
char *T;
{
X    static char *msg[] = { "FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW" };
X
X    ErrCnt [K] = ErrCnt [K] + 1;
X    printf("%s:  %s", msg[K], T);
X    }
X
/* Random */
/*  Random computes
X     X = (Random1 + Random9)^5
X     Random1 = X - FLOOR(X) + 0.000005 * X;
X   and returns the new value of Random1
*/
X
FLOAT Random()
{
X    FLOAT X, Y;
X
X    X = Random1 + Random9;
X    Y = X * X;
X    Y = Y * Y;
X    X = X * Y;
X    Y = X - FLOOR(X);
X    Random1 = Y + X * 0.000005;
X    return(Random1);
X    }
X
/* SqXMinX */
X
SqXMinX (ErrKind)
int ErrKind;
{
X    FLOAT XA, XB;
X
X    XB = X * BInvrse;
X    XA = X - XB;
X    SqEr = ((SQRT(X * X) - XB) - XA) / OneUlp;
X    if (SqEr != Zero) {
X        if (SqEr < MinSqEr) MinSqEr = SqEr;
X        if (SqEr > MaxSqEr) MaxSqEr = SqEr;
X        J = J + 1.0;
X        BadCond(ErrKind, "\n");
X        printf("sqrt( %.17e) - %.17e  = %.17e\n", X * X, X, OneUlp * SqEr);
X        printf("\tinstead of correct value 0 .\n");
X        }
X    }
X
/* NewD */
X
NewD()
{
X    X = Z1 * Q;
X    X = FLOOR(Half - X / Radix) * Radix + X;
X    Q = (Q - X * Z) / Radix + X * X * (D / Radix);
X    Z = Z - Two * X * D;
X    if (Z <= Zero) {
X        Z = - Z;
X        Z1 = - Z1;
X        }
X    D = Radix * D;
X    }
X
/* SR3750 */
X
SR3750()
{
X    if (! ((X - Radix < Z2 - Radix) || (X - Z2 > W - Z2))) {
X        I = I + 1;
X        X2 = SQRT(X * D);
X        Y2 = (X2 - Z2) - (Y - Z2);
X        X2 = X8 / (Y - Half);
X        X2 = X2 - Half * X2 * X2;
X        SqEr = (Y2 + Half) + (Half - X2);
X        if (SqEr < MinSqEr) MinSqEr = SqEr;
X        SqEr = Y2 - X2;
X        if (SqEr > MaxSqEr) MaxSqEr = SqEr;
X        }
X    }
X
/* IsYeqX */
X
IsYeqX()
{
X    if (Y != X) {
X        if (N <= 0) {
X            if (Z == Zero && Q <= Zero)
X                printf("WARNING:  computing\n");
X            else BadCond(Defect, "computing\n");
X            printf("\t(%.17e) ^ (%.17e)\n", Z, Q);
X            printf("\tyielded %.17e;\n", Y);
X            printf("\twhich compared unequal to correct %.17e ;\n",
X                X);
X            printf("\t\tthey differ by %.17e .\n", Y - X);
X            }
X        N = N + 1; /* ... count discrepancies. */
X        }
X    }
X
/* SR3980 */
X
SR3980()
{
X    do {
X        Q = (FLOAT) I;
X        Y = POW(Z, Q);
X        IsYeqX();
X        if (++I > M) break;
X        X = Z * X;
X        } while ( X < W );
X    }
X
/* PrintIfNPositive */
X
PrintIfNPositive()
{
X    if (N > 0) printf("Similar discrepancies have occurred %d times.\n", N);
X    }
X
/* TstPtUf */
X
TstPtUf()
{
X    N = 0;
X    if (Z != Zero) {
X        printf("Since comparison denies Z = 0, evaluating ");
X        printf("(Z + Z) / Z should be safe.\n");
X        sigsave = sigfpe;
X        if (setjmp(ovfl_buf)) goto very_serious;
X        Q9 = (Z + Z) / Z;
X        printf("What the machine gets for (Z + Z) / Z is  %.17e .\n",
X            Q9);
X        if (FABS(Q9 - Two) < Radix * U2) {
X            printf("This is O.K., provided Over/Underflow");
X            printf(" has NOT just been signaled.\n");
X            }
X        else {
X            if ((Q9 < One) || (Q9 > Two)) {
very_serious:
X                N = 1;
X                ErrCnt [Serious] = ErrCnt [Serious] + 1;
X                printf("This is a VERY SERIOUS DEFECT!\n");
X                }
X            else {
X                N = 1;
X                ErrCnt [Defect] = ErrCnt [Defect] + 1;
X                printf("This is a DEFECT!\n");
X                }
X            }
X        sigsave = 0;
X        V9 = Z * One;
X        Random1 = V9;
X        V9 = One * Z;
X        Random2 = V9;
X        V9 = Z / One;
X        if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
X            if (N > 0) Pause();
X            }
X        else {
X            N = 1;
X            BadCond(Defect, "What prints as Z = ");
X            printf("%.17e\n\tcompares different from  ", Z);
X            if (Z != Random1) printf("Z * 1 = %.17e ", Random1);
X            if (! ((Z == Random2)
X                || (Random2 == Random1)))
X                printf("1 * Z == %g\n", Random2);
X            if (! (Z == V9)) printf("Z / 1 = %.17e\n", V9);
X            if (Random2 != Random1) {
X                ErrCnt [Defect] = ErrCnt [Defect] + 1;
X                BadCond(Defect, "Multiplication does not commute!\n");
X                printf("\tComparison alleges that 1 * Z = %.17e\n",
X                    Random2);
X                printf("\tdiffers from Z * 1 = %.17e\n", Random1);
X                }
X            Pause();
X            }
X        }
X    }
X
notify(s)
char *s;
{
X    printf("%s test appears to be inconsistent...\n", s);
X    printf("   PLEASE NOTIFY KARPINKSI!\n");
X    }
X
/*SPLIT msgs.c */
X
/* Instructions */
X
msglist(s)
char **s;
{ while(*s) printf("%s\n", *s++); }
X
Instructions()
{
X  static char *instr[] = {
X    "Lest this program stop prematurely, i.e. before displaying\n",
X    "    `END OF TEST',\n",
X    "try to persuade the computer NOT to terminate execution when an",
X    "error like Over/Underflow or Division by Zero occurs, but rather",
X    "to persevere with a surrogate value after, perhaps, displaying some",
X    "warning.  If persuasion avails naught, don't despair but run this",
X    "program anyway to see how many milestones it passes, and then",
X    "amend it to make further progress.\n",
X    "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
X    0};
X
X    msglist(instr);
X    }
X
/* Heading */
X
Heading()
{
X  static char *head[] = {
X    "Users are invited to help debug and augment this program so it will",
X    "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
X    "Please send suggestions and interesting results to",
X    "\tRichard Karpinski",
X    "\tComputer Center U-76",
X    "\tUniversity of California",
X    "\tSan Francisco, CA 94143-0704, USA\n",
X    "In doing so, please include the following information:",
#ifdef Single
X    "\tPrecision:\tsingle;",
#else
X    "\tPrecision:\tdouble;",
#endif
X    "\tVersion:\t10 February 1989;",
X    "\tComputer:\n",
X    "\tCompiler:\n",
X    "\tOptimization level:\n",
X    "\tOther relevant compiler options:",
X    0};
X
X    msglist(head);
X    }
X
/* Characteristics */
X
Characteristics()
{
X    static char *chars[] = {
X     "Running this program should reveal these characteristics:",
X    "     Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
X    "     Precision = number of significant digits carried.",
X    "     U2 = Radix/Radix^Precision = One Ulp",
X    "\t(OneUlpnit in the Last Place) of 1.000xxx .",
X    "     U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
X    "     Adequacy of guard digits for Mult., Div. and Subt.",
X    "     Whether arithmetic is chopped, correctly rounded, or something else",
X    "\tfor Mult., Div., Add/Subt. and Sqrt.",
X    "     Whether a Sticky Bit used correctly for rounding.",
X    "     UnderflowThreshold = an underflow threshold.",
X    "     E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
X    "     V = an overflow threshold, roughly.",
X    "     V0  tells, roughly, whether  Infinity  is represented.",
X    "     Comparisions are checked for consistency with subtraction",
X    "\tand for contamination with pseudo-zeros.",
X    "     Sqrt is tested.  Y^X is not tested.",
X    "     Extra-precise subexpressions are revealed but NOT YET tested.",
X    "     Decimal-Binary conversion is NOT YET tested for accuracy.",
X    0};
X
X    msglist(chars);
X    }
X
History()
X
{ /* History */
X /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
X    with further massaging by David M. Gay. */
X
X  static char *hist[] = {
X    "The program attempts to discriminate among",
X    "   FLAWs, like lack of a sticky bit,",
X    "   Serious DEFECTs, like lack of a guard digit, and",
X    "   FAILUREs, like 2+2 == 5 .",
X    "Failures may confound subsequent diagnoses.\n",
X    "The diagnostic capabilities of this program go beyond an earlier",
X    "program called `MACHAR', which can be found at the end of the",
X    "book  `Software Manual for the Elementary Functions' (1980) by",
X    "W. J. Cody and W. Waite. Although both programs try to discover",
X    "the Radix, Precision and range (over/underflow thresholds)",
X    "of the arithmetic, this program tries to cope with a wider variety",
X    "of pathologies, and to say how well the arithmetic is implemented.",
X    "\nThe program is based upon a conventional radix representation for",
X    "floating-point numbers, but also allows logarithmic encoding",
X    "as used by certain early WANG machines.\n",
X    "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
X    "see source comments for more history.",
X    0};
X
X    msglist(hist);
X    }
SHAR_EOF
echo 'File paranoia.c is complete' &&
chmod 0644 paranoia.c ||
echo 'restore of paranoia.c failed'
Wc_c="`wc -c < 'paranoia.c'`"
test 57395 -eq "$Wc_c" ||
    echo 'paranoia.c: original size 57395, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= CHANGELOG ==============
if test -f 'CHANGELOG' -a X"$1" != X"-c"; then
    echo 'x - skipping CHANGELOG (File already exists)'
    rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting CHANGELOG (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'CHANGELOG' &&
Thu Nov 15 07:05:39 EDT 1990
Shar'd beta version for release to testers.
X
Sat Nov 17 17:06:17 EDT 1990
Fixed floor and ceil so that they address their local variables at -ve
offsets relative to the frame pointer.
X
Fixed hypot.s so that unnecessary register loads are avoided. I wasn't aware
of the correct assembler mnemonic to multiply %st(0) by the contents of a memory
location.
X
Sat Nov 17 18:49:31 EDT 1990
Copysign rewritten in assembler.
X
Sun Nov 18 11:21:24 EDT 1990
Sinh and cosh now in assembler.
Tanh now in assembler.
Fixed nasty bug in acos.s.
exp__E no longer needed - source deleted.
X
Tue Nov 20 21:24:26 EDT 1990
Sinh and cosh now multiply by 0.5 instead of dividing by 2. It's faster.
Added setinternal for convenience in some cases. Using this makes the
output of paranoia closer to that obtained on a Sun 4.
Added setcont.
Modified the Makefile so that test and paranoia are synonymous.
Modified atan2.c to avoid a floating point register load.
Fixed bug in atan2 - was giving incorrect result when y/x = +/-inf.
X
Thu Nov 22 06:52:47 EDT 1990
Fixed local variable handling in assembler files. Stack needs to be
allocated for local variables.
All assembler routines are now .align 4.
X
Thu Nov 22 19:24:09 EDT 1990
Added MASK_ALL to fpumath.h
X
Fri Nov 23 21:11:33 EDT 1990
Inverse hyperbolics missing from fpumath.h - fixed.
Some definitions missing from atanh.c - didn't work on some systems.
Coded inverse hyperbolics in assembler.
Coded atan2 in assembler.
X
Sat Nov 24 17:12:34 EDT 1990
Modified hypot.s to avoid use of ffree.
atan2 was not working correctly - fixed
Ready for netwide release.
X
Mon Nov 26 19:25:56 EDT 1990
Put in explicit call to cc for compiling lgamma. This means that the
library can link to code not compiled with gcc (not quite true pow.s has to be
changed).
X
Wed Nov 28 22:57:43 EDT 1990
Coded the is??? functions required for IEEE conformance
Coded ieee_retrospective() which prints a summary of IEEE exceptions that are on
when the function is called.
X
Fri Nov 30 17:35:29 EDT 1990
Removed the last vestige of GNU specifics. You can now use the library with 'cc'
although you need gcc/gas to create it.
X
Sat Dec  1 15:48:19 EDT 1990
Found a bug in pow (pow(-x,x) shoudn't work for x non-integral) fixed.
X
Sat Dec  1 21:00:48 EST 1990
Rewrote copysign in a more efficient way.
X
Sun Dec  2 06:51:59 EDT 1990
Added infinity().
Amended README to show the current state.
Changed fpumath.h so as to define HUGE = infinity() if you compile (user
code) with IEEE defined.
X
Sun Dec  2 08:26:20 EDT 1990
Modified log functions and inverse hyperbolics for greater accuracy.
X
Sun Dec  2 10:08:47 EDT 1990
Added MASK_ALL1 to fpumath.h. This enables 64 bit precision significands.
Coded sqrtp to ensure that sqrt is rounded to 53 bits - used in
paranoia.c.
Put sqrtp in fpumath.h and in README.
Made very minor changes to pow.s, floor.s and ceil.s.
ieee_retrospective was printing its leading info for a LOS exception -
fixed.
Rewrote finite so that it does not use the floating point unit and so does
not raise any exceptions.
Produced d2dcomb.summ.
X
Sun Dec  2 18:40:47 EDT 1990
Fixed exp, exp10, expm1, exp2, sinh, cosh, tanh to gain more speed.
X
Wed Dec  5 17:25:50 EDT 1990
Coded max_..., min_..., etc.
X
Wed Dec  5 19:38:56 EDT 1990
Fixed some local parameter handling errors in a few of the assembler
files.
X
Wed Dec  5 20:54:37 EDT 1990
Define MAX(MIN)DOUBLE in terms of max(min)_(sub)normal in fpumath.h.
X
Sat Dec  8 06:09:01 EDT 1990
Fixed bug in isnan 8(%ebp) not $8(%ebp) - typo!
Fixed bug in isinf - not doing correct test.
X
Sat Dec  8 06:36:07 EDT 1990
Coded nextafter.
X
Sat Dec  8 06:46:44 EDT 1990
MAXDOUBLE is wrong in math.h it is correct in fpumath.h
Released to net. Posted to alt.sources.
X
Sun Dec  9 08:26:23 EDT 1990
Added a couple more special cases to nextafter.
X
Sun Dec  9 12:16:04 EDT 1990
Slight modification to Makefile
X
Sun Dec  9 14:33:13 EDT 1990
Coded float version - still testing. You must have an ANSI compiler (gcc is
ok) to use this version.
X
Mon Dec 10 17:21:16 EDT 1990
Float version done.
Fully prototyped all functions in fpumath.h.
Version 2.0 done.
X
Tue Dec 11 06:22:51 EDT 1990
Added a couple of extra cases to nextafter[f].
Made a couple of minor changes to fpumath.h.
Removed AT&T code (j0.c, j1.c, jn.c, erf.c, lgamma.c). These were
avaialble with BSD 4.3 so I had no way of knowing...
X
Wed Dec 12 05:50:31 EDT 1990
Finally got nextafter[f] working correctly.
X
Thu Dec 13 09:44:37 EDT 1990
sinh[f]() now gives good results for all argument ranges. Compute using a
combination of exp() and a ratio of polynomials.
X
Thu Dec 13 18:48:55 EDT 1990
Compute expm1() as expm1() = exp(x)*(1-exp(-x)). This is much better than
before but still not perfect. I need Hart & Cheney...
X
Thu Dec 13 19:38:31 EDT 1990
cosh() was not producing good results for small x. Fixed.
X
Fri Dec 14 08:19:56 EDT 1990
Berkeley provided j0.c, j1.c, erf.c, gamma.c and lgamma.c (no jn.c yet) now
being used. This is *not* BSD code and is not covered by the BSD licence.
***No Berkeley code used.***
mathimpl.h deleted.
README updated.
Defined _getsw() (in ieee_retrospective.c) as a short.
X
Fri Dec 14 11:43:30 EDT 1990
Added shar to Makefile.
Version 2.0 ready for distribution.
SHAR_EOF
chmod 0644 CHANGELOG ||
echo 'restore of CHANGELOG failed'
Wc_c="`wc -c < 'CHANGELOG'`"
test 5231 -eq "$Wc_c" ||
    echo 'CHANGELOG: original size 5231, current size' "$Wc_c"
rm -f _shar_wnt_.tmp
fi
# ============= COPYING ==============
if test -f 'COPYING' -a X"$1" != X"-c"; then
    echo 'x - skipping COPYING (File already exists)'
    rm -f _shar_wnt_.tmp
else
> _shar_wnt_.tmp
echo 'x - extracting COPYING (Text)'
sed 's/^X//' << 'SHAR_EOF' > 'COPYING' &&
X
X            GNU GENERAL PUBLIC LICENSE
X             Version 1, February 1989
X
X Copyright (C) 1989 Free Software Foundation, Inc.
X                    675 Mass Ave, Cambridge, MA 02139, USA
X Everyone is permitted to copy and distribute verbatim copies
X of this license document, but changing it is not allowed.
X
X                Preamble
X
X  The license agreements of most software companies try to keep users
at the mercy of those companies.  By contrast, our General Public
License is intended to guarantee your freedom to share and change free
software--to make sure the software is free for all its users.  The
General Public License applies to the Free Software Foundation's
software and to any other program whose authors commit to using it.
You can use it for your programs, too.
X
X  When we speak of free software, we are referring to freedom, not
price.  Specifically, the General Public License is designed to make
sure that you have the freedom to give away or sell copies of free
software, that you receive source code or can get it if you want it,
that you can change the software or use pieces of it in new free
programs; and that you know you can do these things.
X
X  To protect your rights, we need to make restrictions that forbid
anyone to deny you these rights or to ask you to surrender the rights.
These restrictions translate to certain responsibilities for you if you
distribute copies of the software, or if you modify it.
X
X  For example, if you distribute copies of a such a program, whether
gratis or for a fee, you must give the recipients all the rights that
you have.  You must make sure that they, too, receive or can get the
source code.  And you must tell them their rights.
X
X  We protect your rights with two steps: (1) copyright the software, and
(2) offer you this license which gives you legal permission to copy,
distribute and/or modify the software.
X
X  Also, for each author's protection and ours, we want to make certain
that everyone understands that there is no warranty for this free
software.  If the software is modified by someone else and passed on, we
want its recipients to know that what they have is not the original, so
that any problems introduced by others will not reflect on the original
authors' reputations.
X
X  The precise terms and conditions for copying, distribution and
modification follow.
X
X            GNU GENERAL PUBLIC LICENSE
X   TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
X
X  0. This License Agreement applies to any program or other work which
contains a notice placed by the copyright holder saying it may be
distributed under the terms of this General Public License.  The
"Program", below, refers to any such program or work, and a "work based
on the Program" means either the Program or any work containing the
Program or a portion of it, either verbatim or with modifications.  Each
licensee is addressed as "you".
X
X  1. You may copy and distribute verbatim copies of the Program's source
code as you receive it, in any medium, provided that you conspicuously and
appropriately publish on each copy an appropriate copyright notice and
disclaimer of warranty; keep intact all the notices that refer to this
General Public License and to the absence of any warranty; and give any
other recipients of the Program a copy of this General Public License
along with the Program.  You may charge a fee for the physical act of
transferring a copy.
X
X  2. You may modify your copy or copies of the Program or any portion of
it, and copy and distribute such modifications under the terms of Paragraph
1 above, provided that you also do the following:
X
X    a) cause the modified files to carry prominent notices stating that
X    you changed the files and the date of any change; and
X
X    b) cause the whole of any work that you distribute or publish, that
X    in whole or in part contains the Program or any part thereof, either
X    with or without modifications, to be licensed at no charge to all
X    third parties under the terms of this General Public License (except
X    that you may choose to grant warranty protection to some or all
X    third parties, at your option).
X
X    c) If the modified program normally reads commands interactively when
X    run, you must cause it, when started running for such interactive use
X    in the simplest and most usual way, to print or display an
X    announcement including an appropriate copyright notice and a notice
X    that there is no warranty (or else, saying that you provide a
X    warranty) and that users may redistribute the program under these
X    conditions, and telling the user how to view a copy of this General
X    Public License.
X
X    d) You may charge a fee for the physical act of transferring a
X    copy, and you may at your option offer warranty protection in
X    exchange for a fee.
X
Mere aggregation of another independent work with the Program (or its
derivative) on a volume of a storage or distribution medium does not bring
the other work under the scope of these terms.
X
X  3. You may copy and distribute the Program (or a portion or derivative of
it, under Paragraph 2) in object code or executable form under the terms of
Paragraphs 1 and 2 above provided that you also do one of the following:
X
X    a) accompany it with the complete corresponding machine-readable
X    source code, which must be distributed under the terms of
X    Paragraphs 1 and 2 above; or,
X
X    b) accompany it with a written offer, valid for at least three
X    years, to give any third party free (except for a nominal charge
X    for the cost of distribution) a complete machine-readable copy of the
X    corresponding source code, to be distributed under the terms of
X    Paragraphs 1 and 2 above; or,
X
X    c) accompany it with the information you received as to where the
X    corresponding source code may be obtained.  (This alternative is
X    allowed only for noncommercial distribution and only if you
X    received the program in object code or executable form alone.)
X
Source code for a work means the preferred form of the work for making
modifications to it.  For an executable file, complete source code means
all the source code for all modules it contains; but, as a special
exception, it need not include source code for modules which are standard
libraries that accompany the operating system on which the executable
file runs, or for standard header files or definitions files that
accompany that operating system.
X
X  4. You may not copy, modify, sublicense, distribute or transfer the
Program except as expressly provided under this General Public License.
Any attempt otherwise to copy, modify, sublicense, distribute or transfer
the Program is void, and will automatically terminate your rights to use
the Program under this License.  However, parties who have received
copies, or rights to use copies, from you under this General Public
License will not have their licenses terminated so long as such parties
remain in full compliance.
X
X  5. By copying, distributing or modifying the Program (or any work based
on the Program) you indicate your acceptance of this license to do so,
and all its terms and conditions.
X
X  6. Each time you redistribute the Program (or any work based on the
Program), the recipient automatically receives a license from the original
licensor to copy, distribute or modify the Program subject to these
terms and conditions.  You may not impose any further restrictions on the
recipients' exercise of the rights granted herein.
X
X  7. The Free Software Foundation may publish revised and/or new versions
of the General Public License from time to time.  Such new versions will
be similar in spirit to the present version, but may differ in detail to
address new problems or concerns.
X
Each version is given a distinguishing version number.  If the Program
specifies a version number of the license which applies to it and "any
later version", you have the option of following the terms and conditions
either of that version or of any later version published by the Free
Software Foundation.  If the Program does not specify a version number of
the license, you may choose any version ever published by the Free Software
Foundation.
X
X  8. If you wish to incorporate parts of the Program into other free
programs whose distribution conditions are different, write to the author
to ask for permission.  For software which is copyrighted by the Free
Software Foundation, write to the Free Software Foundation; we sometimes
make exceptions for this.  Our decision will be guided by the two goals
of preserving the free status of all derivatives of our free software and
of promoting the sharing and reuse of software generally.
X
X                NO WARRANTY
X
X  9. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT WHEN
OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES
PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED
OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.  THE ENTIRE RISK AS
TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU.  SHOULD THE
PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,
REPAIR OR CORRECTION.
X
X  10. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES,
INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING
OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED
TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY
YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER
PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE
POSSIBILITY OF SUCH DAMAGES.
X
X             END OF TERMS AND CONDITIONS
X
X    Appendix: How to Apply These Terms to Your New Programs
X
X  If you develop a new program, and you want it to be of the greatest
possible use to humanity, the best way to achieve this is to make it
free software which everyone can redistribute and change under these
terms.
X
X  To do so, attach the following notices to the program.  It is safest to
attach them to the start of each source file to most effectively convey
the exclusion of warranty; and each file should have at least the
"copyright" line and a pointer to where the full notice is found.
X
X    <one line to give the program's name and a brief idea of what it does.>
X    Copyright (C) 19yy  <name of author>
SHAR_EOF
true || echo 'restore of COPYING failed'
fi
echo 'End of mathlib2.0 part 5'
echo 'File COPYING is continued in part 6'
echo 6 > _shar_seq_.tmp
exit 0
--
Glenn Geers                       | "So when it's over, we're back to people.
Department of Theoretical Physics |  Just to prove that human touch can have
The University of Sydney          |  no equal."
Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'