real.c fails floating point tests

[Stephen L Moshier <steve@moshier.net>]

Here are some initial results of testing the new version of real.c.
Exhaustive testing of floating point arithmetic operations is
impractical except in the case of some unary single precision
operations.  We rely on test programs that probe for implementation
errors that have been found and debugged in the past, and on programs
that use random number or deterministic techniques to generate input values
that we think are likely to uncover problems.

GCC uses a software floating point arithmetic for binary <-> decimal
conversion and constant folding.  Previously it passed various test
programs for floating point arithmetic, thus giving some degree of
assurance that the arithmetic was correct.  Virtually all of the
floating point arithmetic test programs that have been developed since
1984 are for IEEE 754 arithmetic, and the programs do not work with
non-IEEE arithmetic.  The previous, working, GCC arithmetic has been
replaced by another arithmetic that does not pass such tests, and
cannot even run them without first going through a considerable
programming and debugging effort to modify the test code.

Decimal <-> binary conversion is the most important floating point
operation for the compiler, since there is no portable alternative but
to write decimal numbers (the ugly C99 hexadecimal input will become
generally available only in the distant future) and hope that they
will be converted with sufficient accuracy into binary machine
numbers.  I have reported some errors for this in the current GCC, and
have also noticed that the decimal conversion method suffers from
double-rounding error (test case to follow later).

Full implementation of C99 floating point pragmas will require that
the GCC arithmetic maintain the IEEE exception flags.  That means the
arithmetic has to be strictly IEEE compliant in each supported
precision.  For the flags to be correct, each operation must be
rounded off directly to target mode without double-rounding.  GCC
previously had a provision for that, but it seems to have been
removed.  Currently GCC calculates an extended precision result and
then double-rounds to the target precision.

It is hard to see how users are going to develop a degree of assurance
that the current GCC arithmetic is correct.  It is perplexing, why it
was so important to remove a working arithmetic and substitute a
non-working arithmetic.  There was a configuration problem for Cray
computers, which do not have a 16-bit integer type.  To solve Cray's
problem, the radix 65536 arithmetic code that depended on 16-bit
integers might have been replaced by the radix 2 arithmetic that was
used in gcc-2.4, and which was kept available as a fallback through
gcc-2.8.


Here is the output of the PARANOIA test program for the current
real.c arithmetic.

Notes:

 This arithmetic has severe underflow and overflow bugs; consequently
the test generates unexpected results while probing for the overflow
and underflow thresholds.

 There is no guard digit, consequently programs that depend on small
differences can fail.  In particular the test function pow2test had to
be omitted because that subtest induces a GCC abort while GCC is attempting
to convert a large number to a decimal string.  The GCC abort is a
consequence of the guard bit defect, which then triggers the overflow bug.

 This test printout is for i386-linux.  I am still trying to get the program
to run on a sparcstation.

----
This version of paranoia omits test for extra precise subexpressions
and includes a few additional tests.

0 = 0.0
00000000 00000000 00000000 00000000 00000000 00000000
1 = 1.00000000000000000000000000000000000000000000000e+0
00000009 00000000 00000000 00000000 00000000 80000000
1+1 = 2.00000000000000000000000000000000000000000000000e+0
00000011 00000000 00000000 00000000 00000000 80000000
Program is now RUNNING tests on small integers:
-1, 0, 1/2, 1, 2, 3, 4, 5, 9, 27, 32 & 240 are O.K.

Searching for Radix and Precision.
Radix = 2.00000000000000000000000000000000000000000000000e+0
00000011 00000000 00000000 00000000 00000000 80000000
Closest relative separation found is U 1 = 6.84227765783602085411977335590779360976690401307e-49
fffffb09 00000000 00000000 00000000 00000000 80000000
Recalculating radix and precision.confirms closest relative separation U1 .
Radix confirmed.
The number of significant digits of the Radix is 1.60000000000000000000000000000000000000000000000e+2
00000041 00000000 00000000 00000000 00000000 a0000000
Subtraction appears to be normalized, as it should be.
Checking for guard digit in *, /, and -.
Serious: - lacks Guard Digit, so cancellation is obscured
Checking rounding on multiply, divide and add/subtract.
* is neither chopped nor correctly rounded.
/ is neither chopped nor correctly rounded.
Addition/Subtraction neither rounds nor chops.
Sticky bit used incorrectly or not at all.
lack(s) of guard digits or failure(s) to correctly round or chop
(noted above) count as one flaw in the final tally below

Does Multiplication commute?  Testing on 20 random pairs.
X * Y == Y * X trial fails.
X=3.63050464352284616725107325384666317732717329962e-1
fffffff9 61a89aaf bf4e4c2a 81e04512 61a05152 b9e1c01e
Y=3.67074300256508518327699029234143153463545945784e-1
fffffff9 738367dd 3ef51b2c 262de0d2 9b8a2a7a bbf129a5

Running test of square root(x).
Testing if sqrt(X * X) == X for 20 Integers X.
Test for sqrt monotonicity.
sqrt has passed a test for Monotonicity.
Testing whether sqrt is rounded or chopped.
Square root is neither chopped nor correctly rounded.
Observed errors run from -6.04938271604938271604938717428160477599954110124e-1
00000005 3acd47eb 232dcd13 74f0bc8b a4587e6b 9add3c0c
	to 5.00000000000000000000000000000000000000000000000e-1
00000001 00000000 00000000 00000000 00000000 80000000
ulps.
Testing powers Z^i for small Integers Z and i.
... no discrepancies found.

Seeking Underflow thresholds UfThold and E0.
Smallest strictly positive number found is E0 = 6.98652447960225958099589122020050053280205923037e-80807125
80000009 00000000 00000000 00000000 00000000 80000000
Z = 6.98652447960225958099589122020050053280205923037e-80807125
80000009 00000000 00000000 00000000 00000000 80000000
Since comparison denies Z = 0, evaluating (Z + Z) / Z should be safe.
What the machine gets for (Z + Z) / Z is 2.00000000000000000000000000000000000000000000000e+0
00000011 00000000 00000000 00000000 00000000 80000000
This is O.K., provided Over/Underflow has NOT just been signaled.
Underflow is gradual; it incurs Absolute Error =
(roundoff in UfThold) < E0.
((CInvrse E0) (1.5+U2)) / (CInvrse (1+U2)) != E0
CInvrse = 1.23611798402889138811011757425565984235889332296e+80807047
7ffff809 00000000 00000000 00000000 00000000 80000000
E0 = 6.98652447960225958099589122020050053280205923037e-80807125
80000009 00000000 00000000 00000000 00000000 80000000
U2 = 1.36845553156720417082395467118155872195338080261e-48
fffffb11 00000000 00000000 00000000 00000000 80000000
X = 1.23611798402889138811011757425565984235889332296e+80807047
7ffff809 00000001 00000000 00000000 00000000 80000000
Y = 1.29542528326416669380795277942005993566740005467e-77
fffff809 00000001 00000000 00000000 00000000 c0000000
Y/X = 1.04795735075967608845681472584251343288499149546e-80807124
80000009 ffffffff ffffffff ffffffff ffffffff bfffffff
The Underflow threshold is 5.10540848309556841246314867661897963420482280488e-80807077
80000501 00000001 00000000 00000000 00000000 80000000
below which calculation may suffer larger Relative error than merely roundoff.
Since underflow occurs below the threshold
UfThold = 2.00000000000000000000000000000000000000000000000e+0
00000011 00000000 00000000 00000000 00000000 80000000
	to the power  -2.68435297000000000000000000000000000000000000000e+8
000000e5 00000000 00000000 00000000 00000000 fffff610
only underflow should afflict the expression 2.00000000000000000000000000000000000000000000000e+0
00000011 00000000 00000000 00000000 00000000 80000000
	to the power  -2.68435298000000000000000000000000000000000000000e+8
000000e5 00000000 00000000 00000000 00000000 fffff620
Actually calculating yields: 2.55270424154778420623157433830948981710241140244e-80807077
800004f9 00000000 00000000 00000000 00000000 80000000
This computed value is O.K.
Searching for Overflow threshold:
This may generate an error.
2.00000000000000000000000000000000000000000000000e+0
00000011 00000000 00000000 00000000 00000000 80000000
	times -Inf
00000006 00000000 00000000 00000000 00000000 00000000
	equals -Inf
00000006 00000000 00000000 00000000 00000000 00000000
Can `Z = -Y' overflow?
Trying it on Y = -Inf
00000006 00000000 00000000 00000000 00000000 00000000
Seems O.K.
Overflow threshold is V  = 3.57831709786311968119428155838269724514906882400e+80807123
7ffffff1 ffffffff ffffffff ffffffff ffffffff ffffffff
Overflow saturates at V0 = +Inf
00000002 00000000 00000000 00000000 00000000 00000000
No Overflow should be signaled for V * 1 = 3.57831709786311968119428155838269724514906882400e+80807123
7ffffff1 ffffffff ffffffff ffffffff ffffffff ffffffff
                           nor for V / 1 = 3.57831709786311968119428155838269724514906882400e+80807123
7ffffff1 ffffffff ffffffff ffffffff ffffffff ffffffff
Any overflow signal separating this * from the one
above is a DEFECT.

 unbalanced range; UfThold * V = 1.82687704666362864775460604089535377456991567872e+47
000004f1 00000001 00000000 00000000 00000000 80000000
	is too far from 1.

The number of  SERIOUS DEFECTs  discovered = 1.
The number of  DEFECTs  discovered =         1.
The number of  FLAWs  discovered =           2.

The arithmetic diagnosed has unacceptable serious defects.
END OF TEST.