paranoia.c

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

        X = Z;
        I = 1;
        SR3980 ();
        if (Z == AInvrse)
            Break = True;
        else
            Z = AInvrse;
    }
    while (!(Break));
    /*=============================================*/
    Milestone = 100;
    /*=============================================*/
    /*  Powers of Radix have been tested, */
    /*         next try a few primes     */
    M = NoTrials;
    Z = Three;
    do {
        X = Z;
        I = 1;
        SR3980 ();
        do {
            Z = Z + Two;
        }
        while (Three * FLOOR (Z / Three) == Z);
    }
    while (Z < Eight * Three);
    if (N > 0) {
        printf ("Errors like this may invalidate financial calculations\n");
        printf ("\tinvolving interest rates.\n");
    }
    PrintIfNPositive ();
    N += N1;
    if (N == 0)
        printf ("... no discrepancies found.\n");
    if (N > 0)
        Pause ();
    else
        printf ("\n");
    /*=============================================*/
    Milestone = 110;
    /*=============================================*/
    printf ("Seeking Underflow thresholds UfThold and E0.\n");
    D = U1;
    if (Precision != FLOOR (Precision)) {
        D = BInvrse;
        X = Precision;
        do {
            D = D * BInvrse;
            X = X - One;
        }
        while (X > Zero);
    }
    Y = One;
    Z = D;
    /* ... D is power of 1/Radix < 1. */
    do {
        C = Y;
        Y = Z;
        Z = Y * Y;
    }
    while ((Y > Z) && (Z + Z > Z));
    Y = C;
    Z = Y * D;
    do {
        C = Y;
        Y = Z;
        Z = Y * D;
    }
    while ((Y > Z) && (Z + Z > Z));
    if (Radix < Two)
        HInvrse = Two;
    else
        HInvrse = Radix;
    HVar = One / HInvrse;
    /* ... 1/HInvrse == HVar == Min(1/Radix, 1/2) */
    CInvrse = One / C;
    E0 = C;
    Z = E0 * HVar;
    /* ...1/Radix^(BIG Integer) << 1 << CInvrse == 1/C */
    do {
        Y = E0;
        E0 = Z;
        Z = E0 * HVar;
    }
    while ((E0 > Z) && (Z + Z > Z));
    UfThold = E0;
    E1 = Zero;
    Q = Zero;
    E9 = U2;
    S = One + E9;
    D = C * S;
    if (D <= C) {
        E9 = Radix * U2;
        S = One + E9;
        D = C * S;
        if (D <= C) {
            BadCond (Failure, "multiplication gets too many last digits wrong.\n");
            Underflow = E0;
            Y1 = Zero;
            PseudoZero = Z;
            Pause ();
        }
    } else {
        Underflow = D;
        PseudoZero = Underflow * HVar;
        UfThold = Zero;
        do {
            Y1 = Underflow;
            Underflow = PseudoZero;
            if (E1 + E1 <= E1) {
                Y2 = Underflow * HInvrse;
                E1 = FABS (Y1 - Y2);
                Q = Y1;
                if ((UfThold == Zero) && (Y1 != Y2))
                    UfThold = Y1;
            }
            PseudoZero = PseudoZero * HVar;
        }
        while ((Underflow > PseudoZero)
            && (PseudoZero + PseudoZero > PseudoZero));
    }
    /* Comment line 4530 .. 4560 */
    if (PseudoZero != Zero) {
        printf ("\n");
        Z = PseudoZero;
        /* ... Test PseudoZero for "phoney- zero" violates */
        /* ... PseudoZero < Underflow or PseudoZero < PseudoZero + PseudoZero
                   ... */
        if (PseudoZero <= Zero) {
            BadCond (Failure, "Positive expressions can underflow to an\n");
            printf ("allegedly negative value\n");
            printf ("PseudoZero that prints out as: %g .\n", PseudoZero);
            X = -PseudoZero;
            if (X <= Zero) {
                printf ("But -PseudoZero, which should be\n");
                printf ("positive, isn't; it prints out as  %g .\n", X);
            }
        } else {
            BadCond (Flaw, "Underflow can stick at an allegedly positive\n");
            printf ("value PseudoZero that prints out as %g .\n", PseudoZero);
        }
        TstPtUf ();
    }
    /*=============================================*/
    Milestone = 120;
    /*=============================================*/
    if (CInvrse * Y > CInvrse * Y1) {
        S = HVar * S;
        E0 = Underflow;
    }
    if (!((E1 == Zero) || (E1 == E0))) {
        BadCond (Defect, "");
        if (E1 < E0) {
            printf ("Products underflow at a higher");
            printf (" threshold than differences.\n");
            if (PseudoZero == Zero)
                E0 = E1;
        } else {
            printf ("Difference underflows at a higher");
            printf (" threshold than products.\n");
        }
    }
    printf ("Smallest strictly positive number found is E0 = %g .\n", E0);
    Z = E0;
    TstPtUf ();
    Underflow = E0;
    if (N == 1)
        Underflow = Y;
    I = 4;
    if (E1 == Zero)
        I = 3;
    if (UfThold == Zero)
        I = I - 2;
    UfNGrad = True;
    switch (I) {
    case 1:
        UfThold = Underflow;
        if ((CInvrse * Q) != ((CInvrse * Y) * S)) {
            UfThold = Y;
            BadCond (Failure, "Either accuracy deteriorates as numbers\n");
            printf ("approach a threshold = %.17e\n", UfThold);;
            printf (" coming down from %.17e\n", C);
            printf (" or else multiplication gets too many last digits wrong.\n");
        }
        Pause ();
        break;

    case 2:
        BadCond (Failure, "Underflow confuses Comparison, which alleges that\n");
        printf ("Q == Y while denying that |Q - Y| == 0; these values\n");
        printf ("print out as Q = %.17e, Y = %.17e .\n", Q, Y2);
        printf ("|Q - Y| = %.17e .\n", FABS (Q - Y2));
        UfThold = Q;
        break;

    case 3:
        X = X;
        break;

    case 4:
        if ((Q == UfThold) && (E1 == E0)
            && (FABS (UfThold - E1 / E9) <= E1)) {
            UfNGrad = False;
            printf ("Underflow is gradual; it incurs Absolute Error =\n");
            printf ("(roundoff in UfThold) < E0.\n");
            Y = E0 * CInvrse;
            Y = Y * (OneAndHalf + U2);
            X = CInvrse * (One + U2);
            Y = Y / X;
            IEEE = (Y == E0);
        }
    }
    if (UfNGrad) {
        printf ("\n");
        sigsave = _sigfpe;
        if (setjmp (ovfl_buf)) {
            printf ("Underflow / UfThold failed!\n");
            R = HVar + HVar;
        } else
            R = SQRT (Underflow / UfThold);
        sigsave = 0;
        if (R <= HVar) {
            Z = R * UfThold;
            X = Z * (One + R * HVar * (One + HVar));
        } else {
            Z = UfThold;
            X = Z * (One + HVar * HVar * (One + HVar));
        }
        if (!((X == Z) || (X - Z != Zero))) {
            BadCond (Flaw, "");
            printf ("X = %.17e\n\tis not equal to Z = %.17e .\n", X, Z);
            Z9 = X - Z;
            printf ("yet X - Z yields %.17e .\n", Z9);
            printf ("    Should this NOT signal Underflow, ");
            printf ("this is a SERIOUS DEFECT\nthat causes ");
            printf ("confusion when innocent statements like\n");;
            printf ("    if (X == Z)  ...  else");
            printf ("  ... (f(X) - f(Z)) / (X - Z) ...\n");
            printf ("encounter Division by Zero although actually\n");
            sigsave = _sigfpe;
            if (setjmp (ovfl_buf))
                printf ("X / Z fails!\n");
            else
                printf ("X / Z = 1 + %g .\n", (X / Z - Half) - Half);
            sigsave = 0;
        }
    }
    printf ("The Underflow threshold is %.17e, %s\n", UfThold,
        " below which");
    printf ("calculation may suffer larger Relative error than ");
    printf ("merely roundoff.\n");
    Y2 = U1 * U1;
    Y = Y2 * Y2;
    Y2 = Y * U1;
    if (Y2 <= UfThold) {
        if (Y > E0) {
            BadCond (Defect, "");
            I = 5;
        } else {
            BadCond (Serious, "");
            I = 4;
        }
        printf ("Range is too narrow; U1^%d Underflows.\n", I);
    }
    /*=============================================*/
    Milestone = 130;
    /*=============================================*/
    Y = -FLOOR (Half - TwoForty * LOG (UfThold) / LOG (HInvrse)) / TwoForty;
    Y2 = Y + Y;
    printf ("Since underflow occurs below the threshold\n");
    printf ("UfThold = (%.17e) ^ (%.17e)\nonly underflow ", HInvrse, Y);
    printf ("should afflict the expression\n\t(%.17e) ^ (%.17e);\n",
        HInvrse, Y2);
    printf ("actually calculating yields:");
    if (setjmp (ovfl_buf)) {
        sigsave = 0;
        BadCond (Serious, "trap on underflow.\n");
    } else {
        sigsave = _sigfpe;
        V9 = POW (HInvrse, Y2);
        sigsave = 0;
        printf (" %.17e .\n", V9);
        if (!((V9 >= Zero) && (V9 <= (Radix + Radix + E9) * UfThold))) {
            BadCond (Serious, "this is not between 0 and underflow\n");
            printf ("   threshold = %.17e .\n", UfThold);
        } else if (!(V9 > UfThold * (One + E9)))
            printf ("This computed value is O.K.\n");
        else {
            BadCond (Defect, "this is not between 0 and underflow\n");
            printf ("   threshold = %.17e .\n", UfThold);
        }
    }
    /*=============================================*/
    Milestone = 140;
    /*=============================================*/
    printf ("\n");
    /* ...calculate Exp2 == exp(2) == 7.389056099... */
    X = Zero;
    I = 2;
    Y = Two * Three;
    Q = Zero;
    N = 0;
    do {
        Z = X;
        I = I + 1;
        Y = Y / (I + I);
        R = Y + Q;
        X = Z + R;
        Q = (Z - X) + R;
    }
    while (X > Z);
    Z = (OneAndHalf + One / Eight) + X / (OneAndHalf * ThirtyTwo);
    X = Z * Z;
    Exp2 = X * X;
    X = F9;
    Y = X - U1;
    printf ("Testing X^((X + 1) / (X - 1)) vs. exp(2) = %.17e as X -> 1.\n",
        Exp2);
    for (I = 1;;) {
        Z = X - BInvrse;
        Z = (X + One) / (Z - (One - BInvrse));
        Q = POW (X, Z) - Exp2;
        if (FABS (Q) > TwoForty * U2) {
            N = 1;
            V9 = (X - BInvrse) - (One - BInvrse);
            BadCond (Defect, "Calculated");
            printf (" %.17e for\n", POW (X, Z));
            printf ("\t(1 + (%.17e) ^ (%.17e);\n", V9, Z);
            printf ("\tdiffers from correct value by %.17e .\n", Q);
            printf ("\tThis much error may spoil financial\n");
            printf ("\tcalculations involving tiny interest rates.\n");
            break;
        } else {
            Z = (Y - X) * Two + Y;
            X = Y;
            Y = Z;
            Z = One + (X - F9) * (X - F9);
            if (Z > One && I < NoTrials)
                I++;
            else {
                if (X > One) {
                    if (N == 0)
                        printf ("Accuracy seems adequate.\n");
                    break;
                } else {
                    X = One + U2;
                    Y = U2 + U2;
                    Y += X;
                    I = 1;
                }
            }
        }
    }
    /*=============================================*/
    Milestone = 150;
    /*=============================================*/
    printf ("Testing powers Z^Q at four nearly extreme values.\n");
    N = 0;
    Z = A1;
    Q = FLOOR (Half - LOG (C) / LOG (A1));
    Break = False;
    do {
        X = CInvrse;
        Y = POW (Z, Q);
        IsYeqX ();
        Q = -Q;
        X = C;
        Y = POW (Z, Q);
        IsYeqX ();
        if (Z < One)
            Break = True;
        else
            Z = AInvrse;
    }
    while (!(Break));
    PrintIfNPositive ();
    if (N == 0)
        printf (" ... no discrepancies found.\n");
    printf ("\n");

    /*=============================================*/
    Milestone = 160;
    /*=============================================*/
    Pause ();
    printf ("Searching for Overflow threshold:\n");
    printf ("This may generate an error.\n");
    Y = -CInvrse;
    V9 = HInvrse * Y;
    sigsave = _sigfpe;
    if (setjmp (ovfl_buf)) {
        I = 0;
        V9 = Y;
        goto overflow;
    }
    do {
        V = Y;
        Y = V9;
        V9 = HInvrse * Y;
    }
    while (V9 < Y);
    I = 1;
  overflow:
    sigsave = 0;
    Z = V9;
    printf ("Can `Z = -Y' overflow?\n");
    printf ("Trying it on Y = %.17e .\n", Y);
    V9 = -Y;
    V0 = V9;
    if (V - Y == V + V0)
        printf ("Seems O.K.\n");
    else {
        printf ("finds a ");
        BadCond (Flaw, "-(-Y) differs from Y.\n");
    }
    if (Z != Y) {
        BadCond (Serious, "");
        printf ("overflow past %.17e\n\tshrinks to %.17e .\n", Y, Z);
    }
    if (I) {
        Y = V * (HInvrse * U2 - HInvrse);
        Z = Y + ((One - HInvrse) * U2) * V;
        if (Z < V0)
            Y = Z;
        if (Y < V0)
            V = Y;
        if (V0 - V < V0)
            V = V0;
    } else {
        V = Y * (HInvrse * U2 - HInvrse);