paranoia.c

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

/*
 *  $Id: paranoia.c,v 1.7 2000/10/31 16:40:03 joel Exp $
 *
 *      A C version of Kahan's Floating Point Test "Paranoia"
 *
 *                      Thos Sumner, UCSF, Feb. 1985
 *                      David Gay, BTL, Jan. 1986
 *
 *      This is a rewrite from the Pascal version by
 *
 *                      B. A. Wichmann, 18 Jan. 1985
 *
 *      (and does NOT exhibit good C programming style).
 *
 *      Sun May 16 18:21:51 MDT 1993 Jeffrey Wheat ([e-mail])
 *      Removed KR_headers defines, removed ANSI prototyping
 *      Cleaned up and reformated code. Added special CYGNUS
 *      "verbose" mode type messages (-DCYGNUS).
 *      Note: This code is VERY NASTY.
 *
 *      Adjusted to use Standard C headers 19 Jan. 1992 (dmg);
 *      compile with -DKR_headers or insert
 * #define KR_headers
 *      at the beginning if you have an old-style C compiler.
 *
 * (C) Apr 19 1983 in BASIC version by:
 *      Professor W. M. Kahan,
 *      567 Evans Hall
 *      Electrical Engineering & Computer Science Dept.
 *      University of California
 *      Berkeley, California 94720
 *      USA
 *
 * converted to Pascal by:
 *      B. A. Wichmann
 *      National Physical Laboratory
 *      Teddington Middx
 *      TW11 OLW
 *      UK
 *
 * converted to C by:
 *
 *      David M. Gay            and     Thos Sumner
 *      AT&T Bell Labs                  Computer Center, Rm. U-76
 *      600 Mountain Avenue             University of California
 *      Murray Hill, NJ 07974           San Francisco, CA 94143
 *      USA                             USA
 *
 * with simultaneous corrections to the Pascal source (reflected
 * in the Pascal source available over netlib).
 * [A couple of bug fixes from dgh = sun!dhough incorporated 31 July 1986.]
 *
 * Reports of results on various systems from all the versions
 * of Paranoia are being collected by Richard Karpinski at the
 * same address as Thos Sumner.  This includes sample outputs,
 * bug reports, and criticisms.
 *
 * You may copy this program freely if you acknowledge its source.
 * Comments on the Pascal version to NPL, please.
 *
 *
 * The C version catches signals from floating-point exceptions.
 * If signal(SIGFPE,...) is unavailable in your environment, you may
 * #define NOSIGNAL to comment out the invocations of signal.
 *
 * This source file is too big for some C compilers, but may be split
 * into pieces.  Comments containing "SPLIT" suggest convenient places
 * for this splitting.  At the end of these comments is an "ed script"
 * (for the UNIX(tm) editor ed) that will do this splitting.
 *
 * By #defining SINGLE_PRECISION when you compile this source, you may
 * obtain a single-precision C version of Paranoia.
 *
 * The following is from the introductory commentary from Wichmann's work:
 *
 * The BASIC program of Kahan is written in Microsoft BASIC using many
 * facilities which have no exact analogy in Pascal.  The Pascal
 * version below cannot therefore be exactly the same.  Rather than be
 * a minimal transcription of the BASIC program, the Pascal coding
 * follows the conventional style of block-structured languages.  Hence
 * the Pascal version could be useful in producing versions in other
 * structured languages.
 *
 * Rather than use identifiers of minimal length (which therefore have
 * little mnemonic significance), the Pascal version uses meaningful
 * identifiers as follows [Note: A few changes have been made for C]:
 *
 *
 * BASIC   C               BASIC   C               BASIC   C
 *
 *   A                       J                       S    StickyBit
 *   A1   AInverse           J0   NoErrors           T
 *   B    Radix                    [Failure]         T0   Underflow
 *   B1   BInverse           J1   NoErrors           T2   ThirtyTwo
 *   B2   RadixD2                  [SeriousDefect]   T5   OneAndHalf
 *   B9   BMinusU2           J2   NoErrors           T7   TwentySeven
 *   C                             [Defect]          T8   TwoForty
 *   C1   CInverse           J3   NoErrors           U    OneUlp
 *   D                             [Flaw]            U0   UnderflowThreshold
 *   D4   FourD              K    PageNo             U1
 *   E0                      L    Milestone          U2
 *   E1                      M                       V
 *   E2   Exp2               N                       V0
 *   E3                      N1                      V8
 *   E5   MinSqEr            O    Zero               V9
 *   E6   SqEr               O1   One                W
 *   E7   MaxSqEr            O2   Two                X
 *   E8                      O3   Three              X1
 *   E9                      O4   Four               X8
 *   F1   MinusOne           O5   Five               X9   Random1
 *   F2   Half               O8   Eight              Y
 *   F3   Third              O9   Nine               Y1
 *   F6                      P    Precision          Y2
 *   F9                      Q                       Y9   Random2
 *   G1   GMult              Q8                      Z
 *   G2   GDiv               Q9                      Z0   PseudoZero
 *   G3   GAddSub            R                       Z1
 *   H                       R1   RMult              Z2
 *   H1   HInverse           R2   RDiv               Z9
 *   I                       R3   RAddSub
 *   IO   NoTrials           R4   RSqrt
 *   I3   IEEE               R9   Random9
 *
 * SqRWrng
 *
 * All the variables in BASIC are true variables and in consequence,
 * the program is more difficult to follow since the "constants" must
 * be determined (the glossary is very helpful).  The Pascal version
 * uses Real constants, but checks are added to ensure that the values
 * are correctly converted by the compiler.
 *
 * The major textual change to the Pascal version apart from the
 * identifiersis that named procedures are used, inserting parameters
 * wherehelpful.  New procedures are also introduced.  The
 * correspondence is as follows:
 *
 *
 * BASIC       Pascal
 * lines
 *
 *   90- 140   Pause
 *  170- 250   Instructions
 *  380- 460   Heading
 *  480- 670   Characteristics
 *  690- 870   History
 * 2940-2950   Random
 * 3710-3740   NewD
 * 4040-4080   DoesYequalX
 * 4090-4110   PrintIfNPositive
 * 4640-4850   TestPartialUnderflow
 *
*/

#include <stdio.h>
#include <string.h>

#if defined(solaris2)
#include <math.h>
#endif

/*
 * To compile this on host using only libm from newlib (and using host libc)
 * do:
 *       gcc -g -DNEED_REENT -DCYGNUS paranoia.c .../newlib-1.6/newlib/libm.a
 */

#ifdef NEED_REENT
#include <reent.h>
struct _reent libm_reent = _REENT_INIT(libm_reent);
struct _reent *_impure_ptr = &libm_reent;
#endif

#ifndef NOSIGNAL
#include <signal.h>
#include <setjmp.h>
#else   /* NOSIGNAL */
#define longjmp(e,v)
#define setjmp(e)       0
#define jmp_buf  int
#endif /* NOSIGNAL */

#ifdef SINGLE_PRECISION
#define FLOAT float
#define FABS(x) (float)fabs((double)(x))
#define FLOOR(x) (float)floor((double)(x))
#define LOG(x) (float)log((double)(x))
#define POW(x,y) (float)pow((double)(x),(double)(y))
#define SQRT(x) (float)sqrt((double)(x))
#else /* !SINGLE_PRECISION */
#define FLOAT double
#define FABS(x) fabs(x)
#define FLOOR(x) floor(x)
#define LOG(x) log(x)
#define POW(x,y) pow(x,y)
#define SQRT(x) sqrt(x)
#endif /* SINGLE_PRECISION */

jmp_buf ovfl_buf;
extern double fabs (), floor (), log (), pow (), sqrt ();
extern void exit ();
typedef void (*Sig_type) ();
FLOAT   Sign (), Random ();
extern void BadCond ();
extern void SqXMinX ();
extern void TstCond ();
extern void notify ();
/* extern int read (); */
extern void Characteristics ();
extern void Heading ();
extern void History ();
extern void Instructions ();
extern void IsYeqX ();
extern void NewD ();
extern void Pause ();
extern void PrintIfNPositive ();
extern void SR3750 ();
extern void SR3980 ();
extern void TstPtUf ();

Sig_type sigsave;

#define KEYBOARD 0

FLOAT   Radix, BInvrse, RadixD2, BMinusU2;

/*Small floating point constants.*/
FLOAT   Zero = 0.0;
FLOAT   Half = 0.5;
FLOAT   One = 1.0;
FLOAT   Two = 2.0;
FLOAT   Three = 3.0;
FLOAT   Four = 4.0;
FLOAT   Five = 5.0;
FLOAT   Eight = 8.0;
FLOAT   Nine = 9.0;
FLOAT   TwentySeven = 27.0;
FLOAT   ThirtyTwo = 32.0;
FLOAT   TwoForty = 240.0;
FLOAT   MinusOne = -1.0;
FLOAT   OneAndHalf = 1.5;

/*Integer constants*/
int     NoTrials = 20;          /*Number of tests for commutativity. */
#define False 0
#define True 1

/*
 * Definitions for declared types
 *      Guard == (Yes, No);
 *      Rounding == (Chopped, Rounded, Other);
 *      Message == packed array [1..40] of char;
 *      Class == (Flaw, Defect, Serious, Failure);
 */
#define Yes 1
#define No  0
#define Chopped 2
#define Rounded 1
#define Other   0
#define Flaw    3
#define Defect  2
#define Serious 1
#define Failure 0
typedef int Guard, Rounding, Class;
typedef char Message;

/* Declarations of Variables */
int     Indx;
char    ch[8];
FLOAT   AInvrse, A1;
FLOAT   C, CInvrse;
FLOAT   D, FourD;
FLOAT   E0, E1, Exp2, E3, MinSqEr;
FLOAT   SqEr, MaxSqEr, E9;
FLOAT   Third;
FLOAT   F6, F9;
FLOAT   HVar, HInvrse;
int     I;
FLOAT   StickyBit, J;
FLOAT   MyZero;
FLOAT   Precision;
FLOAT   Q, Q9;
FLOAT   R, Random9;
FLOAT   T, Underflow, S;
FLOAT   OneUlp, UfThold, U1, U2;
FLOAT   V, V0, V9;
FLOAT   WVar;
FLOAT   X, X1, X2, X8, Random1;
FLOAT   Y, Y1, Y2, Random2;
FLOAT   Z, PseudoZero, Z1, Z2, Z9;
int     ErrCnt[4];
int     fpecount;
int     Milestone;
int     PageNo;
int     M, N, N1;
Guard   GMult, GDiv, GAddSub;
Rounding RMult, RDiv, RAddSub, RSqrt;
int     Break, Done, NotMonot, Monot, Anomaly, IEEE, SqRWrng, UfNGrad;

/* Computed constants.
 *   U1 gap below 1.0, i.e, 1.0 - U1 is next number below 1.0
 *   U2 gap above 1.0, i.e, 1.0 + U2 is next number above 1.0
 */

int     batchmode;              /* global batchmode test */

/* program name and version variables and macro */
char   *temp;
char   *program_name;
char   *program_vers;

#ifndef VERSION
#define VERSION "1.1 [cygnus]"
#endif /* VERSION */

#define basename(cp)    ((temp=(char *)strrchr((cp), '/')) ? temp+1 : (cp))

#ifndef BATCHMODE
# ifdef CYGNUS
#  define BATCHMODE
# endif
#endif

/* floating point exception receiver */
void
_sigfpe (x)
int x;
{
    fpecount++;
    printf ("\n* * * FLOATING-POINT ERROR %d * * *\n", x);
    fflush (stdout);
    if (sigsave) {
#ifndef NOSIGNAL
        signal (SIGFPE, sigsave);
#endif /* NOSIGNAL */
        sigsave = 0;
        longjmp (ovfl_buf, 1);
    }
    exit (1);
}

#ifdef NOMAIN
#define main paranoia
#endif

int
main (argc, argv)
int argc;
char **argv;
{
    /* First two assignments use integer right-hand sides. */
    Zero = 0;
    One = 1;
    Two = One + One;
    Three = Two + One;
    Four = Three + One;
    Five = Four + One;
    Eight = Four + Four;
    Nine = Three * Three;
    TwentySeven = Nine * Three;
    ThirtyTwo = Four * Eight;
    TwoForty = Four * Five * Three * Four;
    MinusOne = -One;
    Half = One / Two;
    OneAndHalf = One + Half;
    ErrCnt[Failure] = 0;
    ErrCnt[Serious] = 0;
    ErrCnt[Defect] = 0;
    ErrCnt[Flaw] = 0;
    PageNo = 1;

#ifdef BATCHMODE
    batchmode = 1;              /* run test in batchmode? */
#else /* !BATCHMODE */
    batchmode = 0;              /* run test interactively */
#endif /* BATCHMODE */

    program_name = basename (argv[0]);
    program_vers = VERSION;

    printf ("%s version %s\n", program_name, program_vers);

    /*=============================================*/
    Milestone = 0;
    /*=============================================*/
#ifndef NOSIGNAL
    signal (SIGFPE, _sigfpe);
#endif

    if (!batchmode) {
        Instructions ();
        Pause ();
        Heading ();
        Instructions ();
        Pause ();
        Heading ();
        Pause ();
        Characteristics ();
        Pause ();
        History ();
        Pause ();
    }

    /*=============================================*/
    Milestone = 7;
    /*=============================================*/
    printf ("Program is now RUNNING tests on small integers:\n");
    TstCond (Failure, (Zero + Zero == Zero) && (One - One == Zero)
        && (One > Zero) && (One + One == Two),
        "0+0 != 0, 1-1 != 0, 1 <= 0, or 1+1 != 2");
    Z = -Zero;
    if (Z != 0.0) {
        ErrCnt[Failure] = ErrCnt[Failure] + 1;
        printf ("Comparison alleges that -0.0 is Non-zero!\n");
        U1 = 0.001;
        Radix = 1;
        TstPtUf ();
    }
    TstCond (Failure, (Three == Two + One) && (Four == Three + One)
        && (Four + Two * (-Two) == Zero)
        && (Four - Three - One == Zero),
        "3 != 2+1, 4 != 3+1, 4+2*(-2) != 0, or 4-3-1 != 0");
    TstCond (Failure, (MinusOne == (0 - One))
        && (MinusOne + One == Zero) && (One + MinusOne == Zero)
        && (MinusOne + FABS (One) == Zero)
        && (MinusOne + MinusOne * MinusOne == Zero),
        "-1+1 != 0, (-1)+abs(1) != 0, or -1+(-1)*(-1) != 0");
    TstCond (Failure, Half + MinusOne + Half == Zero,
        "1/2 + (-1) + 1/2 != 0");
    /*=============================================*/