paranoia.c

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

    Milestone = 10;
    /*=============================================*/
    TstCond (Failure, (Nine == Three * Three)
        && (TwentySeven == Nine * Three) && (Eight == Four + Four)
        && (ThirtyTwo == Eight * Four)
        && (ThirtyTwo - TwentySeven - Four - One == Zero),
        "9 != 3*3, 27 != 9*3, 32 != 8*4, or 32-27-4-1 != 0");
    TstCond (Failure, (Five == Four + One) &&
        (TwoForty == Four * Five * Three * Four)
        && (TwoForty / Three - Four * Four * Five == Zero)
        && (TwoForty / Four - Five * Three * Four == Zero)
        && (TwoForty / Five - Four * Three * Four == Zero),
        "5 != 4+1, 240/3 != 80, 240/4 != 60, or 240/5 != 48");
    if (ErrCnt[Failure] == 0) {
        printf ("-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.\n");
        printf ("\n");
    }
    printf ("Searching for Radix and Precision.\n");
    WVar = One;
    do {
        WVar = WVar + WVar;
        Y = WVar + One;
        Z = Y - WVar;
        Y = Z - One;
    }
    while (MinusOne + FABS (Y) < Zero);
    /*.. now WVar is just big enough that |((WVar+1)-WVar)-1| >= 1 ...*/
    Precision = Zero;
    Y = One;
    do {
        Radix = WVar + Y;
        Y = Y + Y;
        Radix = Radix - WVar;
    }
    while (Radix == Zero);
    if (Radix < Two)
        Radix = One;
    printf ("Radix = %f .\n", Radix);
    if (Radix != 1) {
        WVar = One;
        do {
            Precision = Precision + One;
            WVar = WVar * Radix;
            Y = WVar + One;
        }
        while ((Y - WVar) == One);
    }
    /*... now WVar == Radix^Precision is barely too big to satisfy (WVar+1)-WVar == 1
                                                      ...*/
    U1 = One / WVar;
    U2 = Radix * U1;
    printf ("Closest relative separation found is U1 = %.7e .\n\n", U1);
    printf ("Recalculating radix and precision\n ");

    /*save old values*/
    E0 = Radix;
    E1 = U1;
    E9 = U2;
    E3 = Precision;

    X = Four / Three;
    Third = X - One;
    F6 = Half - Third;
    X = F6 + F6;
    X = FABS (X - Third);
    if (X < U2)
        X = U2;

    /*... now X = (unknown no.) ulps of 1+...*/
    do {
        U2 = X;
        Y = Half * U2 + ThirtyTwo * U2 * U2;
        Y = One + Y;
        X = Y - One;
    }
    while (!((U2 <= X) || (X <= Zero)));

    /*... now U2 == 1 ulp of 1 + ... */
    X = Two / Three;
    F6 = X - Half;
    Third = F6 + F6;
    X = Third - Half;
    X = FABS (X + F6);
    if (X < U1)
        X = U1;

    /*... now  X == (unknown no.) ulps of 1 -... */
    do {
        U1 = X;
        Y = Half * U1 + ThirtyTwo * U1 * U1;
        Y = Half - Y;
        X = Half + Y;
        Y = Half - X;
        X = Half + Y;
    }
    while (!((U1 <= X) || (X <= Zero)));
    /*... now U1 == 1 ulp of 1 - ... */
    if (U1 == E1)
        printf ("confirms closest relative separation U1 .\n");
    else
        printf ("gets better closest relative separation U1 = %.7e .\n", U1);
    WVar = One / U1;
    F9 = (Half - U1) + Half;
    Radix = FLOOR (0.01 + U2 / U1);
    if (Radix == E0)
        printf ("Radix confirmed.\n");
    else
        printf ("MYSTERY: recalculated Radix = %.7e .\n", Radix);
    TstCond (Defect, Radix <= Eight + Eight,
        "Radix is too big: roundoff problems");
    TstCond (Flaw, (Radix == Two) || (Radix == 10)
        || (Radix == One), "Radix is not as good as 2 or 10");
    /*=============================================*/
    Milestone = 20;
    /*=============================================*/
    TstCond (Failure, F9 - Half < Half,
        "(1-U1)-1/2 < 1/2 is FALSE, prog. fails?");
    X = F9;
    I = 1;
    Y = X - Half;
    Z = Y - Half;
    TstCond (Failure, (X != One)
        || (Z == Zero), "Comparison is fuzzy,X=1 but X-1/2-1/2 != 0");
    X = One + U2;
    I = 0;
    /*=============================================*/
    Milestone = 25;
    /*=============================================*/
    /*... BMinusU2 = nextafter(Radix, 0) */
    BMinusU2 = Radix - One;
    BMinusU2 = (BMinusU2 - U2) + One;
    /* Purify Integers */
    if (Radix != One) {
        X = -TwoForty * LOG (U1) / LOG (Radix);
        Y = FLOOR (Half + X);
        if (FABS (X - Y) * Four < One)
            X = Y;
        Precision = X / TwoForty;
        Y = FLOOR (Half + Precision);
        if (FABS (Precision - Y) * TwoForty < Half)
            Precision = Y;
    }
    if ((Precision != FLOOR (Precision)) || (Radix == One)) {
        printf ("Precision cannot be characterized by an Integer number\n");
        printf ("of significant digits but, by itself, this is a minor flaw.\n");
    }
    if (Radix == One)
        printf ("logarithmic encoding has precision characterized solely by U1.\n");
    else
        printf ("The number of significant digits of the Radix is %f .\n",
            Precision);
    TstCond (Serious, U2 * Nine * Nine * TwoForty < One,
        "Precision worse than 5 decimal figures  ");
    /*=============================================*/
    Milestone = 30;
    /*=============================================*/
    /* Test for extra-precise subepressions */
    X = FABS (((Four / Three - One) - One / Four) * Three - One / Four);
    do {
        Z2 = X;
        X = (One + (Half * Z2 + ThirtyTwo * Z2 * Z2)) - One;
    }
    while (!((Z2 <= X) || (X <= Zero)));
    X = Y = Z = FABS ((Three / Four - Two / Three) * Three - One / Four);
    do {
        Z1 = Z;
        Z = (One / Two - ((One / Two - (Half * Z1 + ThirtyTwo * Z1 * Z1))
                + One / Two)) + One / Two;
    }
    while (!((Z1 <= Z) || (Z <= Zero)));
    do {
        do {
            Y1 = Y;
            Y = (Half - ((Half - (Half * Y1 + ThirtyTwo * Y1 * Y1)) + Half
                )) + Half;
        }
        while (!((Y1 <= Y) || (Y <= Zero)));
        X1 = X;
        X = ((Half * X1 + ThirtyTwo * X1 * X1) - F9) + F9;
    }
    while (!((X1 <= X) || (X <= Zero)));
    if ((X1 != Y1) || (X1 != Z1)) {
        BadCond (Serious, "Disagreements among the values X1, Y1, Z1,\n");
        printf ("respectively  %.7e,  %.7e,  %.7e,\n", X1, Y1, Z1);
        printf ("are symptoms of inconsistencies introduced\n");
        printf ("by extra-precise evaluation of arithmetic subexpressions.\n");
        notify ("Possibly some part of this");
        if ((X1 == U1) || (Y1 == U1) || (Z1 == U1))
            printf (
                "That feature is not tested further by this program.\n");
    } else {
        if ((Z1 != U1) || (Z2 != U2)) {
            if ((Z1 >= U1) || (Z2 >= U2)) {
                BadCond (Failure, "");
                notify ("Precision");
                printf ("\tU1 = %.7e, Z1 - U1 = %.7e\n", U1, Z1 - U1);
                printf ("\tU2 = %.7e, Z2 - U2 = %.7e\n", U2, Z2 - U2);
            } else {
                if ((Z1 <= Zero) || (Z2 <= Zero)) {
                    printf ("Because of unusual Radix = %f", Radix);
                    printf (", or exact rational arithmetic a result\n");
                    printf ("Z1 = %.7e, or Z2 = %.7e ", Z1, Z2);
                    notify ("of an\nextra-precision");
                }
                if (Z1 != Z2 || Z1 > Zero) {
                    X = Z1 / U1;
                    Y = Z2 / U2;
                    if (Y > X)
                        X = Y;
                    Q = -LOG (X);
                    printf ("Some subexpressions appear to be calculated extra\n");
                    printf ("precisely with about %g extra B-digits, i.e.\n",
                        (Q / LOG (Radix)));
                    printf ("roughly %g extra significant decimals.\n",
                        Q / LOG (10.));
                }
                printf ("That feature is not tested further by this program.\n");
            }
        }
    }
    Pause ();
    /*=============================================*/
    Milestone = 35;
    /*=============================================*/
    if (Radix >= Two) {
        X = WVar / (Radix * Radix);
        Y = X + One;
        Z = Y - X;
        T = Z + U2;
        X = T - Z;
        TstCond (Failure, X == U2,
            "Subtraction is not normalized X=Y,X+Z != Y+Z!");
        if (X == U2)
            printf (
                "Subtraction appears to be normalized, as it should be.");
    }
    printf ("\nChecking for guard digit in *, /, and -.\n");
    Y = F9 * One;
    Z = One * F9;
    X = F9 - Half;
    Y = (Y - Half) - X;
    Z = (Z - Half) - X;
    X = One + U2;
    T = X * Radix;
    R = Radix * X;
    X = T - Radix;
    X = X - Radix * U2;
    T = R - Radix;
    T = T - Radix * U2;
    X = X * (Radix - One);
    T = T * (Radix - One);
    if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero))
        GMult = Yes;
    else {
        GMult = No;
        TstCond (Serious, False,
            "* lacks a Guard Digit, so 1*X != X");
    }
    Z = Radix * U2;
    X = One + Z;
    Y = FABS ((X + Z) - X * X) - U2;
    X = One - U2;
    Z = FABS ((X - U2) - X * X) - U1;
    TstCond (Failure, (Y <= Zero)
        && (Z <= Zero), "* gets too many final digits wrong.\n");
    Y = One - U2;
    X = One + U2;
    Z = One / Y;
    Y = Z - X;
    X = One / Three;
    Z = Three / Nine;
    X = X - Z;
    T = Nine / TwentySeven;
    Z = Z - T;
    TstCond (Defect, X == Zero && Y == Zero && Z == Zero,
        "Division lacks a Guard Digit, so error can exceed 1 ulp\n\
or  1/3  and  3/9  and  9/27 may disagree");
    Y = F9 / One;
    X = F9 - Half;
    Y = (Y - Half) - X;
    X = One + U2;
    T = X / One;
    X = T - X;
    if ((X == Zero) && (Y == Zero) && (Z == Zero))
        GDiv = Yes;
    else {
        GDiv = No;
        TstCond (Serious, False,
            "Division lacks a Guard Digit, so X/1 != X");
    }
    X = One / (One + U2);
    Y = X - Half - Half;
    TstCond (Serious, Y < Zero,
        "Computed value of 1/1.000..1 >= 1");
    X = One - U2;
    Y = One + Radix * U2;
    Z = X * Radix;
    T = Y * Radix;
    R = Z / Radix;
    StickyBit = T / Radix;
    X = R - X;
    Y = StickyBit - Y;
    TstCond (Failure, X == Zero && Y == Zero,
        "* and/or / gets too many last digits wrong");
    Y = One - U1;
    X = One - F9;
    Y = One - Y;
    T = Radix - U2;
    Z = Radix - BMinusU2;
    T = Radix - T;
    if ((X == U1) && (Y == U1) && (Z == U2) && (T == U2))
        GAddSub = Yes;
    else {
        GAddSub = No;
        TstCond (Serious, False,
            "- lacks Guard Digit, so cancellation is obscured");
    }
    if (F9 != One && F9 - One >= Zero) {
        BadCond (Serious, "comparison alleges  (1-U1) < 1  although\n");
        printf ("  subtraction yields  (1-U1) - 1 = 0 , thereby vitiating\n");
        printf ("  such precautions against division by zero as\n");
        printf ("  ...  if (X == 1.0) {.....} else {.../(X-1.0)...}\n");
    }
    if (GMult == Yes && GDiv == Yes && GAddSub == Yes)
        printf (
            "     *, /, and - appear to have guard digits, as they should.\n");
    /*=============================================*/
    Milestone = 40;
    /*=============================================*/
    Pause ();
    printf ("Checking rounding on multiply, divide and add/subtract.\n");
    RMult = Other;
    RDiv = Other;
    RAddSub = Other;
    RadixD2 = Radix / Two;
    A1 = Two;
    Done = False;
    do {
        AInvrse = Radix;
        do {
            X = AInvrse;
            AInvrse = AInvrse / A1;
        }
        while (!(FLOOR (AInvrse) != AInvrse));
        Done = (X == One) || (A1 > Three);
        if (!Done)
            A1 = Nine + One;
    }
    while (!(Done));
    if (X == One)
        A1 = Radix;
    AInvrse = One / A1;
    X = A1;
    Y = AInvrse;
    Done = False;
    do {
        Z = X * Y - Half;
        TstCond (Failure, Z == Half,
            "X * (1/X) differs from 1");
        Done = X == Radix;
        X = Radix;
        Y = One / X;
    }
    while (!(Done));
    Y2 = One + U2;
    Y1 = One - U2;
    X = OneAndHalf - U2;
    Y = OneAndHalf + U2;
    Z = (X - U2) * Y2;
    T = Y * Y1;
    Z = Z - X;
    T = T - X;
    X = X * Y2;
    Y = (Y + U2) * Y1;
    X = X - OneAndHalf;
    Y = Y - OneAndHalf;
    if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T <= Zero)) {
        X = (OneAndHalf + U2) * Y2;
        Y = OneAndHalf - U2 - U2;
        Z = OneAndHalf + U2 + U2;
        T = (OneAndHalf - U2) * Y1;
        X = X - (Z + U2);
        StickyBit = Y * Y1;
        S = Z * Y2;
        T = T - Y;
        Y = (U2 - Y) + StickyBit;
        Z = S - (Z + U2 + U2);
        StickyBit = (Y2 + U2) * Y1;
        Y1 = Y2 * Y1;
        StickyBit = StickyBit - Y2;
        Y1 = Y1 - Half;
        if ((X == Zero) && (Y == Zero) && (Z == Zero) && (T == Zero)
            && (StickyBit == Zero) && (Y1 == Half)) {
            RMult = Rounded;
            printf ("Multiplication appears to round correctly.\n");
        } else if ((X + U2 == Zero) && (Y < Zero) && (Z + U2 == Zero)
                && (T < Zero) && (StickyBit + U2 == Zero)
            && (Y1 < Half)) {
            RMult = Chopped;
            printf ("Multiplication appears to chop.\n");
        } else
            printf ("* is neither chopped nor correctly rounded.\n");
        if ((RMult == Rounded) && (GMult == No))
            notify ("Multiplication");
    } else
        printf ("* is neither chopped nor correctly rounded.\n");
    /*=============================================*/
    Milestone = 45;
    /*=============================================*/
    Y2 = One + U2;
    Y1 = One - U2;
    Z = OneAndHalf + U2 + U2;
    X = Z / Y2;
    T = OneAndHalf - U2 - U2;
    Y = (T - U2) / Y1;
    Z = (Z + U2) / Y2;
    X = X - OneAndHalf;
    Y = Y - T;
    T = T / Y1;
    Z = Z - (OneAndHalf + U2);
    T = (U2 - OneAndHalf) + T;
    if (!((X > Zero) || (Y > Zero) || (Z > Zero) || (T > Zero))) {
        X = OneAndHalf / Y2;
        Y = OneAndHalf - U2;
        Z = OneAndHalf + U2;
        X = X - Y;
        T = OneAndHalf / Y1;
        Y = Y / Y1;
        T = T - (Z + U2);