paranoia.c

RTEMS (Real-Time Executive for Multiprocessor Systems) is a free open source real-time operating sys

        Y = Y - Z;
        Z = Z / Y2;
        Y1 = (Y2 + U2) / Y2;
        Z = Z - OneAndHalf;
        Y2 = Y1 - Y2;
        Y1 = (F9 - U1) / F9;
        if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
            && (Y2 == Zero) && (Y2 == Zero)
            && (Y1 - Half == F9 - Half)) {
            RDiv = Rounded;
            printf ("Division appears to round correctly.\n");
            if (GDiv == No)
                notify ("Division");
        } else if ((X < Zero) && (Y < Zero) && (Z < Zero) && (T < Zero)
            && (Y2 < Zero) && (Y1 - Half < F9 - Half)) {
            RDiv = Chopped;
            printf ("Division appears to chop.\n");
        }
    }
    if (RDiv == Other)
        printf ("/ is neither chopped nor correctly rounded.\n");
    BInvrse = One / Radix;
    TstCond (Failure, (BInvrse * Radix - Half == Half),
        "Radix * ( 1 / Radix ) differs from 1");
    /*=============================================*/
    Milestone = 50;
    /*=============================================*/
    TstCond (Failure, ((F9 + U1) - Half == Half)
        && ((BMinusU2 + U2) - One == Radix - One),
        "Incomplete carry-propagation in Addition");
    X = One - U1 * U1;
    Y = One + U2 * (One - U2);
    Z = F9 - Half;
    X = (X - Half) - Z;
    Y = Y - One;
    if ((X == Zero) && (Y == Zero)) {
        RAddSub = Chopped;
        printf ("Add/Subtract appears to be chopped.\n");
    }
    if (GAddSub == Yes) {
        X = (Half + U2) * U2;
        Y = (Half - U2) * U2;
        X = One + X;
        Y = One + Y;
        X = (One + U2) - X;
        Y = One - Y;
        if ((X == Zero) && (Y == Zero)) {
            X = (Half + U2) * U1;
            Y = (Half - U2) * U1;
            X = One - X;
            Y = One - Y;
            X = F9 - X;
            Y = One - Y;
            if ((X == Zero) && (Y == Zero)) {
                RAddSub = Rounded;
                printf ("Addition/Subtraction appears to round correctly.\n");
                if (GAddSub == No)
                    notify ("Add/Subtract");
            } else
                printf ("Addition/Subtraction neither rounds nor chops.\n");
        } else
            printf ("Addition/Subtraction neither rounds nor chops.\n");
    } else
        printf ("Addition/Subtraction neither rounds nor chops.\n");
    S = One;
    X = One + Half * (One + Half);
    Y = (One + U2) * Half;
    Z = X - Y;
    T = Y - X;
    StickyBit = Z + T;
    if (StickyBit != Zero) {
        S = Zero;
        BadCond (Flaw, "(X - Y) + (Y - X) is non zero!\n");
    }
    StickyBit = Zero;
    if ((GMult == Yes) && (GDiv == Yes) && (GAddSub == Yes)
        && (RMult == Rounded) && (RDiv == Rounded)
        && (RAddSub == Rounded) && (FLOOR (RadixD2) == RadixD2)) {
        printf ("Checking for sticky bit.\n");
        X = (Half + U1) * U2;
        Y = Half * U2;
        Z = One + Y;
        T = One + X;
        if ((Z - One <= Zero) && (T - One >= U2)) {
            Z = T + Y;
            Y = Z - X;
            if ((Z - T >= U2) && (Y - T == Zero)) {
                X = (Half + U1) * U1;
                Y = Half * U1;
                Z = One - Y;
                T = One - X;
                if ((Z - One == Zero) && (T - F9 == Zero)) {
                    Z = (Half - U1) * U1;
                    T = F9 - Z;
                    Q = F9 - Y;
                    if ((T - F9 == Zero) && (F9 - U1 - Q == Zero)) {
                        Z = (One + U2) * OneAndHalf;
                        T = (OneAndHalf + U2) - Z + U2;
                        X = One + Half / Radix;
                        Y = One + Radix * U2;
                        Z = X * Y;
                        if (T == Zero && X + Radix * U2 - Z == Zero) {
                            if (Radix != Two) {
                                X = Two + U2;
                                Y = X / Two;
                                if ((Y - One == Zero))
                                    StickyBit = S;
                            } else
                                StickyBit = S;
                        }
                    }
                }
            }
        }
    }
    if (StickyBit == One)
        printf ("Sticky bit apparently used correctly.\n");
    else
        printf ("Sticky bit used incorrectly or not at all.\n");
    TstCond (Flaw, !(GMult == No || GDiv == No || GAddSub == No ||
            RMult == Other || RDiv == Other || RAddSub == Other),
        "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");
    /*=============================================*/
    Milestone = 60;
    /*=============================================*/
    printf ("\n");
    printf ("Does Multiplication commute?  ");
    printf ("Testing on %d random pairs.\n", NoTrials);
    Random9 = SQRT (3.0);
    Random1 = Third;
    I = 1;
    do {
        X = Random ();
        Y = Random ();
        Z9 = Y * X;
        Z = X * Y;
        Z9 = Z - Z9;
        I = I + 1;
    }
    while (!((I > NoTrials) || (Z9 != Zero)));
    if (I == NoTrials) {
        Random1 = One + Half / Three;
        Random2 = (U2 + U1) + One;
        Z = Random1 * Random2;
        Y = Random2 * Random1;
        Z9 = (One + Half / Three) * ((U2 + U1) + One) - (One + Half /
            Three) * ((U2 + U1) + One);
    }
    if (!((I == NoTrials) || (Z9 == Zero)))
        BadCond (Defect, "X * Y == Y * X trial fails.\n");
    else
        printf ("     No failures found in %d integer pairs.\n", NoTrials);
    /*=============================================*/
    Milestone = 70;
    /*=============================================*/
    printf ("\nRunning test of square root(x).\n");
    TstCond (Failure, (Zero == SQRT (Zero))
        && (-Zero == SQRT (-Zero))
        && (One == SQRT (One)), "Square root of 0.0, -0.0 or 1.0 wrong");
    MinSqEr = Zero;
    MaxSqEr = Zero;
    J = Zero;
    X = Radix;
    OneUlp = U2;
    SqXMinX (Serious);
    X = BInvrse;
    OneUlp = BInvrse * U1;
    SqXMinX (Serious);
    X = U1;
    OneUlp = U1 * U1;
    SqXMinX (Serious);
    if (J != Zero)
        Pause ();
    printf ("Testing if sqrt(X * X) == X for %d Integers X.\n", NoTrials);
    J = Zero;
    X = Two;
    Y = Radix;
    if ((Radix != One))
        do {
            X = Y;
            Y = Radix * Y;
        }
        while (!((Y - X >= NoTrials)));
    OneUlp = X * U2;
    I = 1;
    while (I <= NoTrials) {
        X = X + One;
        SqXMinX (Defect);
        if (J > Zero)
            break;
        I = I + 1;
    }
    printf ("Test for sqrt monotonicity.\n");
    I = -1;
    X = BMinusU2;
    Y = Radix;
    Z = Radix + Radix * U2;
    NotMonot = False;
    Monot = False;
    while (!(NotMonot || Monot)) {
        I = I + 1;
        X = SQRT (X);
        Q = SQRT (Y);
        Z = SQRT (Z);
        if ((X > Q) || (Q > Z))
            NotMonot = True;
        else {
            Q = FLOOR (Q + Half);
            if ((I > 0) || (Radix == Q * Q))
                Monot = True;
            else if (I > 0) {
                if (I > 1)
                    Monot = True;
                else {
                    Y = Y * BInvrse;
                    X = Y - U1;
                    Z = Y + U1;
                }
            } else {
                Y = Q;
                X = Y - U2;
                Z = Y + U2;
            }
        }
    }
    if (Monot)
        printf ("sqrt has passed a test for Monotonicity.\n");
    else {
        BadCond (Defect, "");
        printf ("sqrt(X) is non-monotonic for X near %.7e .\n", Y);
    }
    /*=============================================*/
    Milestone = 80;
    /*=============================================*/
    MinSqEr = MinSqEr + Half;
    MaxSqEr = MaxSqEr - Half;
    Y = (SQRT (One + U2) - One) / U2;
    SqEr = (Y - One) + U2 / Eight;
    if (SqEr > MaxSqEr)
        MaxSqEr = SqEr;
    SqEr = Y + U2 / Eight;
    if (SqEr < MinSqEr)
        MinSqEr = SqEr;
    Y = ((SQRT (F9) - U2) - (One - U2)) / U1;
    SqEr = Y + U1 / Eight;
    if (SqEr > MaxSqEr)
        MaxSqEr = SqEr;
    SqEr = (Y + One) + U1 / Eight;
    if (SqEr < MinSqEr)
        MinSqEr = SqEr;
    OneUlp = U2;
    X = OneUlp;
    for (Indx = 1; Indx <= 3; ++Indx) {
        Y = SQRT ((X + U1 + X) + F9);
        Y = ((Y - U2) - ((One - U2) + X)) / OneUlp;
        Z = ((U1 - X) + F9) * Half * X * X / OneUlp;
        SqEr = (Y + Half) + Z;
        if (SqEr < MinSqEr)
            MinSqEr = SqEr;
        SqEr = (Y - Half) + Z;
        if (SqEr > MaxSqEr)
            MaxSqEr = SqEr;
        if (((Indx == 1) || (Indx == 3)))
            X = OneUlp * Sign (X) * FLOOR (Eight / (Nine * SQRT (OneUlp)));
        else {
            OneUlp = U1;
            X = -OneUlp;
        }
    }
    /*=============================================*/
    Milestone = 85;
    /*=============================================*/
    SqRWrng = False;
    Anomaly = False;
    RSqrt = Other;              /* ~dgh */
    if (Radix != One) {
        printf ("Testing whether sqrt is rounded or chopped.\n");
        D = FLOOR (Half + POW (Radix, One + Precision - FLOOR (Precision)));
        /* ... == Radix^(1 + fract) if (Precision == Integer + fract. */
        X = D / Radix;
        Y = D / A1;
        if ((X != FLOOR (X)) || (Y != FLOOR (Y))) {
            Anomaly = True;
        } else {
            X = Zero;
            Z2 = X;
            Y = One;
            Y2 = Y;
            Z1 = Radix - One;
            FourD = Four * D;
            do {
                if (Y2 > Z2) {
                    Q = Radix;
                    Y1 = Y;
                    do {
                        X1 = FABS (Q + FLOOR (Half - Q / Y1) * Y1);
                        Q = Y1;
                        Y1 = X1;
                    }
                    while (!(X1 <= Zero));
                    if (Q <= One) {
                        Z2 = Y2;
                        Z = Y;
                    }
                }
                Y = Y + Two;
                X = X + Eight;
                Y2 = Y2 + X;
                if (Y2 >= FourD)
                    Y2 = Y2 - FourD;
            }
            while (!(Y >= D));
            X8 = FourD - Z2;
            Q = (X8 + Z * Z) / FourD;
            X8 = X8 / Eight;
            if (Q != FLOOR (Q))
                Anomaly = True;
            else {
                Break = False;
                do {
                    X = Z1 * Z;
                    X = X - FLOOR (X / Radix) * Radix;
                    if (X == One)
                        Break = True;
                    else
                        Z1 = Z1 - One;
                }
                while (!(Break || (Z1 <= Zero)));
                if ((Z1 <= Zero) && (!Break))
                    Anomaly = True;
                else {
                    if (Z1 > RadixD2)
                        Z1 = Z1 - Radix;
                    do {
                        NewD ();
                    }
                    while (!(U2 * D >= F9));
                    if (D * Radix - D != WVar - D)
                        Anomaly = True;
                    else {
                        Z2 = D;
                        I = 0;
                        Y = D + (One + Z) * Half;
                        X = D + Z + Q;
                        SR3750 ();
                        Y = D + (One - Z) * Half + D;
                        X = D - Z + D;
                        X = X + Q + X;
                        SR3750 ();
                        NewD ();
                        if (D - Z2 != WVar - Z2)
                            Anomaly = True;
                        else {
                            Y = (D - Z2) + (Z2 + (One - Z) * Half);
                            X = (D - Z2) + (Z2 - Z + Q);
                            SR3750 ();
                            Y = (One + Z) * Half;
                            X = Q;
                            SR3750 ();
                            if (I == 0)
                                Anomaly = True;
                        }
                    }
                }
            }
        }
        if ((I == 0) || Anomaly) {
            BadCond (Failure, "Anomalous arithmetic with Integer < ");
            printf ("Radix^Precision = %.7e\n", WVar);
            printf (" fails test whether sqrt rounds or chops.\n");
            SqRWrng = True;
        }
    }
    if (!Anomaly) {
        if (!((MinSqEr < Zero) || (MaxSqEr > Zero))) {
            RSqrt = Rounded;
            printf ("Square root appears to be correctly rounded.\n");
        } else {
            if ((MaxSqEr + U2 > U2 - Half) || (MinSqEr > Half)
                || (MinSqEr + Radix < Half))
                SqRWrng = True;
            else {
                RSqrt = Chopped;
                printf ("Square root appears to be chopped.\n");
            }
        }
    }
    if (SqRWrng) {
        printf ("Square root is neither chopped nor correctly rounded.\n");
        printf ("Observed errors run from %.7e ", MinSqEr - Half);
        printf ("to %.7e ulps.\n", Half + MaxSqEr);
        TstCond (Serious, MaxSqEr - MinSqEr < Radix * Radix,
            "sqrt gets too many last digits wrong");
    }
    /*=============================================*/
    Milestone = 90;
    /*=============================================*/
    Pause ();
    printf ("Testing powers Z^i for small Integers Z and i.\n");
    N = 0;
    /* ... test powers of zero. */
    I = 0;
    Z = -Zero;
    M = 3;
    Break = False;
    do {
        X = One;
        SR3980 ();
        if (I <= 10) {
            I = 1023;
            SR3980 ();
        }
        if (Z == MinusOne)
            Break = True;
        else {
            Z = MinusOne;
            /* .. if(-1)^N is invalid, replace MinusOne by One. */
            I = -4;
        }
    }
    while (!Break);
    PrintIfNPositive ();
    N1 = N;
    N = 0;
    Z = A1;
    M = (int) FLOOR (Two * LOG (WVar) / LOG (A1));
    Break = False;
    do {