Subversion Repositories eco32

/*      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).

        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 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

=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=*=

Below is an "ed script" that splits para.c into 10 files

of the form part[1-8].c, subs.c, and msgs.c, plus a header

file, paranoia.h, that these files require.

r paranoia.c

$

?SPLIT

 .d

+d

-,$w msgs.c

-,$d

?SPLIT

 .d

+d

-,$w subs.c

-,$d

?part8

+d

?include

 .,$w part8.c

 .,$d

-d

?part7

+d

?include

 .,$w part7.c

 .,$d

-d

?part6

+d

?include

 .,$w part6.c

 .,$d

-d

?part5

+d

?include

 .,$w part5.c

 .,$d

-d

?part4

+d

?include

 .,$w part4.c

 .,$d

-d

?part3

+d

?include

 .,$w part3.c

 .,$d

-d

?part2

+d

?include

 .,$w part2.c

 .,$d

?SPLIT

 .d

1,/^#include/-1d

1,$w part1.c

/Computed constants/,$d

1,$s/^int/extern &/

1,$s/^FLOAT/extern &/

1,$s/^char/extern &/

1,$s! = .*!;!

/^Guard/,/^Round/s/^/extern /

/^jmp_buf/s/^/extern /

/^Sig_type/s/^/extern /

s/$/\

extern void sigfpe(INT);/

w paranoia.h

q

*/

#include <stdio.h>

#ifndef NOSIGNAL

#include <signal.h>

#endif

#include <setjmp.h>

#ifdef Single

#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

#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

jmp_buf ovfl_buf;

#ifdef KR_headers

#define VOID /* void */

#define INT /* int */

#define FP /* FLOAT */

#define CHARP /* char * */

#define CHARPP /* char ** */

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();

#else

#define VOID void

#define INT int

#define FP FLOAT

#define CHARP char *

#define CHARPP char **

#ifdef __STDC__

#include <stdlib.h>

#include <math.h>

#else

#ifdef __cplusplus

extern "C" {

#endif

extern double fabs(double), floor(double), log(double);

extern double pow(double,double), sqrt(double);

extern void exit(INT);

#ifdef __cplusplus

        }

#endif

#endif

typedef void (*Sig_type)(int);

FLOAT Sign(FLOAT), Random(void);

extern void BadCond(int, char*);

extern void SqXMinX(int);

extern void TstCond(int, int, char*);

extern void notify(char*);

extern int read(int, char*, int);

#endif

#undef V9

extern void Characteristics(VOID);

extern void Heading(VOID);

extern void History(VOID);

extern void Instructions(VOID);

extern void IsYeqX(VOID);

extern void NewD(VOID);

extern void Pause(VOID);

extern void PrintIfNPositive(VOID);

extern void SR3750(VOID);

extern void SR3980(VOID);

extern void TstPtUf(VOID);

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 H, 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 W;

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 */

/* floating point exception receiver */

 void

sigfpe(INT x)

{

        fpecount++;

        printf("\n* * * FLOATING-POINT ERROR %d * * *\n", x);

        fflush(stdout);

        if (sigsave) {

#ifndef NOSIGNAL

                signal(SIGFPE, sigsave);

#endif

                sigsave = 0;

                longjmp(ovfl_buf, 1);

                }

        exit(1);

}

main(VOID)

{

        /* 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;

        /*=============================================*/

        Milestone = 0;

        /*=============================================*/

#ifndef NOSIGNAL

        signal(SIGFPE, sigfpe);

#endif

        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");

                U2 = 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");

        /*=============================================*/

        /*SPLIT

        {

                extern void part2(VOID), part3(VOID), part4(VOID),

                        part5(VOID), part6(VOID), part7(VOID);

                int part8(VOID);

                part2();

                part3();

                part4();

                part5();

                part6();

                part7();

                return part8();

                }

        }

#include "paranoia.h"

void part2(VOID){

*/

        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");

        W = One;

        do  {

                W = W + W;

                Y = W + One;

                Z = Y - W;

                Y = Z - One;

                } while (MinusOne + FABS(Y) < Zero);

        /*.. now W is just big enough that |((W+1)-W)-1| >= 1 ...*/

        Precision = Zero;

        Y = One;

        do  {

                Radix = W + Y;

                Y = Y + Y;

                Radix = Radix - W;

                } while ( Radix == Zero);

        if (Radix < Two) Radix = One;

        printf("Radix = %f .\n", Radix);

        if (Radix != 1) {

                W = One;

                do  {

                        Precision = Precision + One;

                        W = W * Radix;

                        Y = W + One;

                        } while ((Y - W) == One);

                }

        /*... now W == Radix^Precision is barely too big to satisfy (W+1)-W == 1

                                                      ...*/

        U1 = One / W;

        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);

        W = 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();

        /*=============================================*/

        /*SPLIT