summaryrefslogblamecommitdiffstats
path: root/testsuites/samples/paranoia/paranoia.c
blob: da03db660640f46018869f01adfdd0698022159c (plain) (tree)
1
2
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
  
        

























































































































































                                                                            



                     






                                                                             
                  





































                                                         
                         



































































                                                                       
                      








                                
             
































































































































































                                                                           
               
        


                           


                                       
                                                                        


                     
                         
                  
                             





                                     
                   

                                        

                                
         
                                  
     
                                                                                     
                                                           
                    


















































                                                                             
                    


























































































































                                                                                    
                                   




























































































































































































































































































































































































































































































































































                                                                               
                                                  











                                                     
                                                
















                                                                      
                                                      























































                                                                           
                                                  










































































                                                                            

                                                    

                      
                  



                                                         
                      




















                                                                                   
                                      










                                                    
                                           





























                                                                              
                     




































































                                                                                  
                            


                                           
                        
                            
                                                    

                        
                                                       










































































































































































































































































                                                                                 
                         





























































































































































































































































































                                                                               
                                                              











































                                                                    
                     





























































































































































































                                                                                      
/*
 *  $Id$
 *
 *      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 (cassidy@cygnus.com)
 *      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");
    /*=============================================*/
    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);
        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 {
        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);
        V = V + ((One - HInvrse) * U2) * Y;
    }
    printf ("Overflow threshold is V  = %.17e .\n", V);
    if (I)
        printf ("Overflow saturates at V0 = %.17e .\n", V0);
    else
        printf ("There is no saturation value because \
the system traps on overflow.\n");
    V9 = V * One;
    printf ("No Overflow should be signaled for V * 1 = %.17e\n", V9);
    V9 = V / One;
    printf ("                           nor for V / 1 = %.17e .\n", V9);
    printf ("Any overflow signal separating this * from the one\n");
    printf ("above is a DEFECT.\n");
    /*=============================================*/
    Milestone = 170;
    /*=============================================*/
    if (!(-V < V && -V0 < V0 && -UfThold < V && UfThold < V)) {
        BadCond (Failure, "Comparisons involving ");
        printf ("+-%g, +-%g\nand +-%g are confused by Overflow.",
            V, V0, UfThold);
    }
    /*=============================================*/
    Milestone = 175;
    /*=============================================*/
    printf ("\n");
    for (Indx = 1; Indx <= 3; ++Indx) {
        switch (Indx) {
        case 1:
            Z = UfThold;
            break;
        case 2:
            Z = E0;
            break;
        case 3:
            Z = PseudoZero;
            break;
        }
        if (Z != Zero) {
            V9 = SQRT (Z);
            Y = V9 * V9;
            if (Y / (One - Radix * E9) < Z
                || Y > (One + Radix * E9) * Z) {        /* dgh: + E9 --> * E9 */
                if (V9 > U1)
                    BadCond (Serious, "");
                else
                    BadCond (Defect, "");
                printf ("Comparison alleges that what prints as Z = %.17e\n", Z);
                printf (" is too far from sqrt(Z) ^ 2 = %.17e .\n", Y);
            }
        }
    }
    /*=============================================*/
    Milestone = 180;
    /*=============================================*/
    for (Indx = 1; Indx <= 2; ++Indx) {
        if (Indx == 1)
            Z = V;
        else
            Z = V0;
        V9 = SQRT (Z);
        X = (One - Radix * E9) * V9;
        V9 = V9 * X;
        if (((V9 < (One - Two * Radix * E9) * Z) || (V9 > Z))) {
            Y = V9;
            if (X < WVar)
                BadCond (Serious, "");
            else
                BadCond (Defect, "");
            printf ("Comparison alleges that Z = %17e\n", Z);
            printf (" is too far from sqrt(Z) ^ 2 (%.17e) .\n", Y);
        }
    }
    /*=============================================*/
    Milestone = 190;
    /*=============================================*/
    Pause ();
    X = UfThold * V;
    Y = Radix * Radix;
    if (X * Y < One || X > Y) {
        if (X * Y < U1 || X > Y / U1)
            BadCond (Defect, "Badly");
        else
            BadCond (Flaw, "");

        printf (" unbalanced range; UfThold * V = %.17e\n\t%s\n",
            X, "is too far from 1.\n");
    }
    /*=============================================*/
    Milestone = 200;
    /*=============================================*/
    for (Indx = 1; Indx <= 5; ++Indx) {
        X = F9;
        switch (Indx) {
        case 2:
            X = One + U2;
            break;
        case 3:
            X = V;
            break;
        case 4:
            X = UfThold;
            break;
        case 5:
            X = Radix;
        }
        Y = X;
        sigsave = _sigfpe;
        if (setjmp (ovfl_buf))
            printf ("  X / X  traps when X = %g\n", X);
        else {
            V9 = (Y / X - Half) - Half;
            if (V9 == Zero)
                continue;
            if (V9 == -U1 && Indx < 5)
                BadCond (Flaw, "");
            else
                BadCond (Serious, "");
            printf ("  X / X differs from 1 when X = %.17e\n", X);
            printf ("  instead, X / X - 1/2 - 1/2 = %.17e .\n", V9);
        }
        sigsave = 0;
    }
    /*=============================================*/
    Milestone = 210;
    /*=============================================*/
    MyZero = Zero;
    printf ("\n");
    printf ("What message and/or values does Division by Zero produce?\n");
#ifndef BATCHMODE
    printf ("This can interupt your program.  You can ");
    printf ("skip this part if you wish.\n");
    printf ("Do you wish to compute 1 / 0? ");
    fflush (stdout);
    read (KEYBOARD, ch, 8);
    if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif /* !BATCHMODE */
        sigsave = _sigfpe;
        printf ("    Trying to compute 1 / 0 produces ...");
        if (!setjmp (ovfl_buf))
            printf ("  %.7e .\n", One / MyZero);
        sigsave = 0;
#ifndef BATCHMODE
    } else
        printf ("O.K.\n");
    printf ("\nDo you wish to compute 0 / 0? ");
    fflush (stdout);
    read (KEYBOARD, ch, 80);
    if ((ch[0] == 'Y') || (ch[0] == 'y')) {
#endif /* !BATCHMODE */
        sigsave = _sigfpe;
        printf ("\n    Trying to compute 0 / 0 produces ...");
        if (!setjmp (ovfl_buf))
            printf ("  %.7e .\n", Zero / MyZero);
        sigsave = 0;
#ifndef BATCHMODE
    } else
        printf ("O.K.\n");
#endif /* !BATCHMODE */
    /*=============================================*/
    Milestone = 220;
    /*=============================================*/

    Pause ();
    printf ("\n");
    {
        static char *msg[] =
        {
            "FAILUREs  encountered =",
            "SERIOUS DEFECTs  discovered =",
            "DEFECTs  discovered =",
            "FLAWs  discovered ="};
        int     i;
        for (i = 0; i < 4; i++)
            if (ErrCnt[i])
                printf ("The number of  %-29s %d.\n",
                    msg[i], ErrCnt[i]);
    }

    printf ("\n");
    if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[Defect]
            + ErrCnt[Flaw]) > 0) {
        if ((ErrCnt[Failure] + ErrCnt[Serious] + ErrCnt[
                    Defect] == 0) && (ErrCnt[Flaw] > 0)) {
            printf ("The arithmetic diagnosed seems ");
            printf ("Satisfactory though flawed.\n");
        }
        if ((ErrCnt[Failure] + ErrCnt[Serious] == 0)
            && (ErrCnt[Defect] > 0)) {
            printf ("The arithmetic diagnosed may be Acceptable\n");
            printf ("despite inconvenient Defects.\n");
        }
        if ((ErrCnt[Failure] + ErrCnt[Serious]) > 0) {
            printf ("The arithmetic diagnosed has ");
            printf ("unacceptable Serious Defects.\n");
        }
        if (ErrCnt[Failure] > 0) {
            printf ("Potentially fatal FAILURE may have spoiled this");
            printf (" program's subsequent diagnoses.\n");
        }
    } else {

        printf ("No failures, defects nor flaws have been discovered.\n");
        if (!((RMult == Rounded) && (RDiv == Rounded)
                && (RAddSub == Rounded) && (RSqrt == Rounded)))
            printf ("The arithmetic diagnosed seems Satisfactory.\n");
        else {
            if (StickyBit >= One &&
                (Radix - Two) * (Radix - Nine - One) == Zero) {
                printf ("Rounding appears to conform to ");
                printf ("the proposed IEEE standard P");
                if ((Radix == Two) &&
                    ((Precision - Four * Three * Two) *
                        (Precision - TwentySeven -
                            TwentySeven + One) == Zero))
                    printf ("754");
                else
                    printf ("854");
                if (IEEE)
                    printf (".\n");
                else {
                    printf (",\nexcept for possibly Double Rounding");
                    printf (" during Gradual Underflow.\n");
                }
            }
            printf ("The arithmetic diagnosed appears to be Excellent!\n");
        }
    }

    if (fpecount)
        printf ("\nA total of %d floating point exceptions were registered.\n",
            fpecount);
    printf ("END OF TEST.\n");
    return 0;
}

FLOAT
Sign (X)
     FLOAT   X;
{
    return X >= 0. ? 1.0 : -1.0;
}

void
Pause ()
{
#ifndef BATCHMODE
    char    ch[8];

    printf ("\nTo continue, press RETURN");
    fflush (stdout);
    read (KEYBOARD, ch, 8);
#endif /* !BATCHMODE */
#ifndef CYGNUS
    printf ("\nDiagnosis resumes after milestone Number %d", Milestone);
    printf ("          Page: %d\n\n", PageNo);
    ++Milestone;
    ++PageNo;
#endif /* !CYGNUS */
}

void
TstCond (K, Valid, T)
     int     K, Valid;
     char   *T;
{
#ifdef CYGNUS
    printf ("TEST: %s\n", T);
#endif /* CYGNUS */
    if (!Valid) {
        BadCond (K, T);
        printf (".\n");
    }
#ifdef CYGNUS
    printf ("PASS: %s\n", T);
#endif /* CYGNUS */
}

void
BadCond (K, T)
     int     K;
     char   *T;
{
    static char *msg[] =
    {"FAILURE", "SERIOUS DEFECT", "DEFECT", "FLAW"};

    ErrCnt[K] = ErrCnt[K] + 1;
#ifndef CYGNUS
    printf ("%s:  %s", msg[K], T);
#else
    printf ("ERROR: Severity: %s:  %s", msg[K], T);
#endif /* CYGNUS */
}

/*
 * Random computes
 *     X = (Random1 + Random9)^5
 *     Random1 = X - FLOOR(X) + 0.000005 * X;
 *   and returns the new value of Random1
*/
FLOAT
Random ()
{
    FLOAT   X, Y;

    X = Random1 + Random9;
    Y = X * X;
    Y = Y * Y;
    X = X * Y;
    Y = X - FLOOR (X);
    Random1 = Y + X * 0.000005;
    return (Random1);
}

void
SqXMinX (ErrKind)
     int     ErrKind;
{
    FLOAT   XA, XB;

    XB = X * BInvrse;
    XA = X - XB;
    SqEr = ((SQRT (X * X) - XB) - XA) / OneUlp;
    if (SqEr != Zero) {
        if (SqEr < MinSqEr)
            MinSqEr = SqEr;
        if (SqEr > MaxSqEr)
            MaxSqEr = SqEr;
        J = J + 1.0;
        BadCond (ErrKind, "\n");
        printf ("sqrt( %.17e) - %.17e  = %.17e\n", X * X, X, OneUlp * SqEr);
        printf ("\tinstead of correct value 0 .\n");
    }
}

void
NewD ()
{
    X = Z1 * Q;
    X = FLOOR (Half - X / Radix) * Radix + X;
    Q = (Q - X * Z) / Radix + X * X * (D / Radix);
    Z = Z - Two * X * D;
    if (Z <= Zero) {
        Z = -Z;
        Z1 = -Z1;
    }
    D = Radix * D;
}

void
SR3750 ()
{
    if (!((X - Radix < Z2 - Radix) || (X - Z2 > WVar - Z2))) {
        I = I + 1;
        X2 = SQRT (X * D);
        Y2 = (X2 - Z2) - (Y - Z2);
        X2 = X8 / (Y - Half);
        X2 = X2 - Half * X2 * X2;
        SqEr = (Y2 + Half) + (Half - X2);
        if (SqEr < MinSqEr)
            MinSqEr = SqEr;
        SqEr = Y2 - X2;
        if (SqEr > MaxSqEr)
            MaxSqEr = SqEr;
    }
}

void
IsYeqX ()
{
    if (Y != X) {
        if (N <= 0) {
            if (Z == Zero && Q <= Zero)
                printf ("WARNING:  computing\n");
            else
                BadCond (Defect, "computing\n");
            printf ("\t(%.17e) ^ (%.17e)\n", Z, Q);
            printf ("\tyielded %.17e;\n", Y);
            printf ("\twhich compared unequal to correct %.17e ;\n",
                X);
            printf ("\t\tthey differ by %.17e .\n", Y - X);
        }
        N = N + 1;              /* ... count discrepancies. */
    }
}

void
SR3980 ()
{
    do {
        Q = (FLOAT) I;
        Y = POW (Z, Q);
        IsYeqX ();
        if (++I > M)
            break;
        X = Z * X;
    }
    while (X < WVar);
}

void
PrintIfNPositive ()
{
    if (N > 0)
        printf ("Similar discrepancies have occurred %d times.\n", N);
}

void
TstPtUf ()
{
    N = 0;
    if (Z != Zero) {
        printf ("Since comparison denies Z = 0, evaluating ");
        printf ("(Z + Z) / Z should be safe.\n");
        sigsave = _sigfpe;
        if (setjmp (ovfl_buf))
            goto very_serious;
        Q9 = (Z + Z) / Z;
        printf ("What the machine gets for (Z + Z) / Z is  %.17e .\n",
            Q9);
        if (FABS (Q9 - Two) < Radix * U2) {
            printf ("This is O.K., provided Over/Underflow");
            printf (" has NOT just been signaled.\n");
        } else {
            if ((Q9 < One) || (Q9 > Two)) {
              very_serious:
                N = 1;
                ErrCnt[Serious] = ErrCnt[Serious] + 1;
                printf ("This is a VERY SERIOUS DEFECT!\n");
            } else {
                N = 1;
                ErrCnt[Defect] = ErrCnt[Defect] + 1;
                printf ("This is a DEFECT!\n");
            }
        }
        sigsave = 0;
        V9 = Z * One;
        Random1 = V9;
        V9 = One * Z;
        Random2 = V9;
        V9 = Z / One;
        if ((Z == Random1) && (Z == Random2) && (Z == V9)) {
            if (N > 0)
                Pause ();
        } else {
            N = 1;
            BadCond (Defect, "What prints as Z = ");
            printf ("%.17e\n\tcompares different from  ", Z);
            if (Z != Random1)
                printf ("Z * 1 = %.17e ", Random1);
            if (!((Z == Random2)
                    || (Random2 == Random1)))
                printf ("1 * Z == %g\n", Random2);
            if (!(Z == V9))
                printf ("Z / 1 = %.17e\n", V9);
            if (Random2 != Random1) {
                ErrCnt[Defect] = ErrCnt[Defect] + 1;
                BadCond (Defect, "Multiplication does not commute!\n");
                printf ("\tComparison alleges that 1 * Z = %.17e\n",
                    Random2);
                printf ("\tdiffers from Z * 1 = %.17e\n", Random1);
            }
            Pause ();
        }
    }
}

void
notify (s)
     char   *s;
{
    printf ("%s test appears to be inconsistent...\n", s);
    printf ("   PLEASE NOTIFY KARPINKSI!\n");
}

void
msglist (s)
     char  **s;
{
    while (*s)
        printf ("%s\n", *s++);
}

void
Instructions ()
{
    static char *instr[] =
    {
        "Lest this program stop prematurely, i.e. before displaying\n",
        "    `END OF TEST',\n",
        "try to persuade the computer NOT to terminate execution when an",
        "error like Over/Underflow or Division by Zero occurs, but rather",
        "to persevere with a surrogate value after, perhaps, displaying some",
        "warning.  If persuasion avails naught, don't despair but run this",
        "program anyway to see how many milestones it passes, and then",
        "amend it to make further progress.\n",
        "Answer questions with Y, y, N or n (unless otherwise indicated).\n",
        0};

    msglist (instr);
}

void
Heading ()
{
    static char *head[] =
    {
        "Users are invited to help debug and augment this program so it will",
        "cope with unanticipated and newly uncovered arithmetic pathologies.\n",
        "Please send suggestions and interesting results to",
        "\tRichard Karpinski",
        "\tComputer Center U-76",
        "\tUniversity of California",
        "\tSan Francisco, CA 94143-0704, USA\n",
        "In doing so, please include the following information:",
#ifdef SINGLE_PRECISION
        "\tPrecision:\tsingle;",
#else /* !SINGLE_PRECISION */
        "\tPrecision:\tdouble;",
#endif /* SINGLE_PRECISION */
        "\tVersion:\t10 February 1989;",
        "\tComputer:\n",
        "\tCompiler:\n",
        "\tOptimization level:\n",
        "\tOther relevant compiler options:",
        0};

    msglist (head);
}

void
Characteristics ()
{
    static char *chars[] =
    {
        "Running this program should reveal these characteristics:",
        "     Radix = 1, 2, 4, 8, 10, 16, 100, 256 ...",
        "     Precision = number of significant digits carried.",
        "     U2 = Radix/Radix^Precision = One Ulp",
        "\t(OneUlpnit in the Last Place) of 1.000xxx .",
        "     U1 = 1/Radix^Precision = One Ulp of numbers a little less than 1.0 .",
        "     Adequacy of guard digits for Mult., Div. and Subt.",
        "     Whether arithmetic is chopped, correctly rounded, or something else",
        "\tfor Mult., Div., Add/Subt. and Sqrt.",
        "     Whether a Sticky Bit used correctly for rounding.",
        "     UnderflowThreshold = an underflow threshold.",
        "     E0 and PseudoZero tell whether underflow is abrupt, gradual, or fuzzy.",
        "     V = an overflow threshold, roughly.",
        "     V0  tells, roughly, whether  Infinity  is represented.",
        "     Comparisions are checked for consistency with subtraction",
        "\tand for contamination with pseudo-zeros.",
        "     Sqrt is tested.  Y^X is not tested.",
        "     Extra-precise subexpressions are revealed but NOT YET tested.",
        "     Decimal-Binary conversion is NOT YET tested for accuracy.",
        0};

    msglist (chars);
}

void
History ()
{                               /* History */
    /* Converted from Brian Wichmann's Pascal version to C by Thos Sumner,
        with further massaging by David M. Gay. */

    static char *hist[] =
    {
        "The program attempts to discriminate among",
        "   FLAWs, like lack of a sticky bit,",
        "   Serious DEFECTs, like lack of a guard digit, and",
        "   FAILUREs, like 2+2 == 5 .",
        "Failures may confound subsequent diagnoses.\n",
        "The diagnostic capabilities of this program go beyond an earlier",
        "program called `MACHAR', which can be found at the end of the",
        "book  `Software Manual for the Elementary Functions' (1980) by",
        "W. J. Cody and W. Waite. Although both programs try to discover",
        "the Radix, Precision and range (over/underflow thresholds)",
        "of the arithmetic, this program tries to cope with a wider variety",
        "of pathologies, and to say how well the arithmetic is implemented.",
        "\nThe program is based upon a conventional radix representation for",
        "floating-point numbers, but also allows logarithmic encoding",
        "as used by certain early WANG machines.\n",
        "BASIC version of this program (C) 1983 by Prof. W. M. Kahan;",
        "see source comments for more history.",
        0};

    msglist (hist);
}