real.c implementation

real.c implementation

[Brad Lucier]

Recently, there has been discussion about the floating-point arithmetic
implementation in real.c.  I plan to submit some patches to fix some small
errors soon (possibly this weekend) and adapt Gaston Gonnet's tests
(previously mentioned on this list) for the situation where the target
arithmetic you're testing is not the same as the supported host arithmetic.

However, with Roger Sayle's submission of the sqrt patch, I'd like to
make some general comments about strategy.

There are basically two approaches to implementing proper IEEE arithmetic:

1.	Implement algorithms where the working precision is the target
	precision plus two bits (guard and sticky) and apply whatever
	tricks there are in the literature (Tuckerman's test for 0.5ulp
	accuracy of sqrt, Guy Steele's and Will Clinger's algorithms for
	exact decimal<->binary conversions) to get the correct answers.

2.	Implement somewhat inaccurate algorithms in an intermediate precision
	so large that conversion to the target precision gives correct
	results.

I believe that real.c adopts the second approach.

One may have several goals here:

A.	Implement correct round-to-nearest target precision results.

B.	Implement correct directed-rounding (up, down, or toward zero)
	target precision results.

Now, strategy 1 works (by definition ;-) for both goals A and B.  However,
generally speaking, strategy 2 works only for goal A (i.e., it does not
work with directed roundings), and it works only when the intermediate
precision is "large enough".  Large enough means that to get correctly
rounded results to k bits accuracy in the target precision, one must
calculate intermediate results to at least 2k+2 bits accuracy, and this
bound is sharp.

Based on these theoretical results, I expect that the current 160-bit
intermediate precision should be fine for goal A for single, double,
and extended-double IEEE arithmetic, but not for either IBM 106-bit
precision arithmetic or quad-length IEEE arithmetic (wth 113 bits mantissa).
(It is because of the inadequate length of the intermediate precision for
IEEE-quad arithmetic that I suggested to Roger Sayle that for his patch

http://gcc.gnu.org/ml/gcc-patches/2002-10/msg01033.html

he include the target mode so that one can use Tuckerman's test to
get correctly-rounded results in the target precision.
In other words, I was confused, and conflated the two strategies 1 and 2
in order to get around the limited precision of the intermediate
calculations.)

Having said this, Richard's adaptation of paranoia to use real.c passes
in *all* these precisions.  I believe, however, that Gonnet's implementation
of the lattice algorithm will find examples (after several hours of
computation, unfortunately) where the current real.c will give incorrectly
rounded results, in either IBM 106-bit precision or IEEE-quad precision.

So, what to do?  It depends on answers to the following questions:

(i)	Is goal B ever a possibility?  If so, we'll have to abandon
	strategy 2, any rely solely on strategy 1.  This will require
	a significant rewrite of real.c

(ii)	Should one believe the theorems and increase the intermediate
	precision so that results before rounding are accurate to at
	lease 2*113+2=228 bits?  One would need 256 bits intermediate
	accuracy to do this.  One could either believe the theorems or
	wait for someone like me to adapt Gonnet's tests to come up with
	real examples of erroneous results of the current implementation
	for, e.g., quad-IEEE arithmetic.

In part (ii), one might have the intermediate precision be machine-dependent,
so only target machines with TFMODE arithmetic (is that right?  I mean 128-bit
floating-point numbers) would actually compute intermediate results to
256 bits accuracy, and target machines with only smaller arithmetic could
get by with 160-bit arithmetic.

At any rate, I wanted to bring up these topics while the discussion was
still "warm"; I intend to try to get some results this weekend.

Brad

Re: real.c implementation

[Richard Henderson]

On Thu, Oct 17, 2002 at 02:01:05PM -0500, Brad Lucier wrote:
> 1.	Implement algorithms where the working precision is the target
> 	precision plus two bits (guard and sticky) and apply whatever
> 	tricks there are in the literature (Tuckerman's test for 0.5ulp
> 	accuracy of sqrt, Guy Steele's and Will Clinger's algorithms for
> 	exact decimal<->binary conversions) to get the correct answers.
> 
> 2.	Implement somewhat inaccurate algorithms in an intermediate precision
> 	so large that conversion to the target precision gives correct
> 	results.
> 
> I believe that real.c adopts the second approach.

It was not *intended* that real.c adopt the second approach.

My understanding is that the guard bit is supposed to be the correct
value of the .5ulp bit, and the sticky bit is supposed to be set if
the result is not exact at the .5ulp position.

Thus one should get correct results for any target precision if the
following rules are observed:

	* The intermediate representation has _at least_ two more
	  bits than the target format.  Obviously.

	* The bit 0 of the intermediate precision is sticky; all
	  other bits contain the proper values.

	* When rounding to the target format, the guard bit is
	  read from the target's lsb-1, and the sticky bit is the
	  sum of all bits between lsb-2 and 0.

Perhaps I'm being naieve.  If this is incorrect, stop me now.
An explanation of why this is wrong would also be appreciated.

If this is correct, do we get better results for decimal<->binary
conversion if all computations are correctly rounded for, say,
a 158-bit intermediate representation?

If so, does this allow us to get correct results for the worst-case
113-bit target format with a 126-bit intermediate format?  I.e. can
we remove the extra word added the other week?  I'm guessing not, 
but I don't actually know.


r~

Re: real.c implementation

[Brad Lucier]

> My understanding is that the guard bit is supposed to be the correct
> value of the .5ulp bit, and the sticky bit is supposed to be set if
> the result is not exact at the .5ulp position.
> 
> Thus one should get correct results for any target precision if the
> following rules are observed:
> 
> 	* The intermediate representation has _at least_ two more
> 	  bits than the target format.  Obviously.
> 
> 	* The bit 0 of the intermediate precision is sticky; all
> 	  other bits contain the proper values.
> 
> 	* When rounding to the target format, the guard bit is
> 	  read from the target's lsb-1, and the sticky bit is the
> 	  sum of all bits between lsb-2 and 0.
> 
> Perhaps I'm being naieve.  If this is incorrect, stop me now.
> An explanation of why this is wrong would also be appreciated.

This may very well be correct; I did not really understand what your
goal was before by reading the code.  I'll think about it more this
weekend.  Perhaps this is why paranoia is passing for +,-,*,/.

> If this is correct, do we get better results for decimal<->binary
> conversion if all computations are correctly rounded for, say,
> a 158-bit intermediate representation?

I don't know.  What I do know is that having correct intermediate results
to 2k+2 bits after doing all the conversion stuff does guarantee the
correct target result to k bits.  Something similar will come up with
sqrt.

> If so, does this allow us to get correct results for the worst-case
> 113-bit target format with a 126-bit intermediate format?  I.e. can
> we remove the extra word added the other week?  I'm guessing not, 
> but I don't actually know.

The problem is that unless you want to do some algorithmic tricks (along
the lines I gave in my e-mail) the easiest and safest way to compute
decimal<->binary conversions and sqrt is to calculate results to > 2k+2 bits
and then round to k bits.  I'll try to come up with some tests to see
if 113-bit decimal<->binary conversion works.

Brad

Re: real.c implementation

[Brad Lucier]

> My understanding is that the guard bit is supposed to be the correct
> value of the .5ulp bit, and the sticky bit is supposed to be set if
> the result is not exact at the .5ulp position.
> 
> Thus one should get correct results for any target precision if the
> following rules are observed:
> 
> 	* The intermediate representation has _at least_ two more
> 	  bits than the target format.  Obviously.
> 
> 	* The bit 0 of the intermediate precision is sticky; all
> 	  other bits contain the proper values.
> 
> 	* When rounding to the target format, the guard bit is
> 	  read from the target's lsb-1, and the sticky bit is the
> 	  sum of all bits between lsb-2 and 0.
> 
> Perhaps I'm being naieve.  If this is incorrect, stop me now.
> An explanation of why this is wrong would also be appreciated.

This seems to be correct.  I also think it's a very clever idea.

It may give more confidence if you set up an artificial arithmetic
that has 158-bit accuracy to be tested by paranoia.

> If this is correct, do we get better results for decimal<->binary
> conversion if all computations are correctly rounded for, say,
> a 158-bit intermediate representation?
> 
> If so, does this allow us to get correct results for the worst-case
> 113-bit target format with a 126-bit intermediate format?  I.e. can
> we remove the extra word added the other week?  I'm guessing not, 
> but I don't actually know.

I don't know the answers to these questions.  Steele and White, Clinger,
and David Gay (his code is at netlib.org) worked out provably correct methods
to do exact binary<->decimal conversion.  All of them rely on extra-precision
integer arithmetic for some cases.  I don't know if Moshier's code has this
property, and I don't know if it benefits from extra intermediate precision.

Is exact conversion a goal for gcc?

Brad

Re: real.c implementation

[Richard Henderson]

On Sun, Oct 20, 2002 at 10:00:41PM -0500, Brad Lucier wrote:
> It may give more confidence if you set up an artificial arithmetic
> that has 158-bit accuracy to be tested by paranoia.

Working on it.  There are a couple of failures.

> I don't know the answers to these questions.  Steele and White, Clinger,
> and David Gay (his code is at netlib.org) worked out provably correct methods
> to do exact binary<->decimal conversion.  All of them rely on extra-precision
> integer arithmetic for some cases.  I don't know if Moshier's code has this
> property, and I don't know if it benefits from extra intermediate precision.

Neither the current code nor the previous use arbitrary precision integer
arithmetic.  Both versions, as far as I can tell, rely on being about to
compute 10**N with lg(N) multiplications with "enough" accuracy.

> Is exact conversion a goal for gcc?

I'd say so.

I'm actually *less* concerned about binary->decimal, since for almost
all targets this is only used for debugging dumps -- fp constants are
emitted in the assembly as integer data.  There are a few hold-outs
though, so I don't want to be gratuitously wrong either.


r~

Re: real.c implementation

[Richard Henderson]

On Sun, Oct 20, 2002 at 10:00:41PM -0500, Brad Lucier wrote:
> It may give more confidence if you set up an artificial arithmetic
> that has 158-bit accuracy to be tested by paranoia.

This uncovered a bug in the addition routines.  We'd get the right
result any time the sticky bit was in a different word from the target
format's lsb, which is not the case with the 158-bit artifical format.

With this change, the real_internal format gets 3 DEFECTS related to
sqrt, and one FLAW:

UfThold = 1.2763521207738921031157871691547449085512057e-80807077
V  = 3.0000000000000000000000000000000000000000000e+80807123
FLAW:   unbalanced range; UfThold * V = 4.56719261665907161938651510223838443642478919680e+46
        is too far from 1.

which is simply due to the range I've given the exponent field.
I could make it smaller, but then we'd have a hole in the bitfields,
which would pessimize gcc initializing the things.

Results for all normal fp target formats are unchanged.


r~


        * real.c (sticky_rshift_significand): Return inexact, don't
        or it in immediately.
        (sub_significands): Accept incomming carry.
        (div_significands, rtd_divmod): Update for sub_significands change.
        (round_for_format): Update for sticky_rshift_significand change.
        (do_add): Don't involve the inexact bit in addition, do give the
        inexact bit as the subtraction carry-in.
        (encode_internal, decode_internal, real_internal_format): New.
        * real.h (real_internal_format): Declare.

Index: real.c
===================================================================
RCS file: /cvs/gcc/gcc/gcc/real.c,v
retrieving revision 1.100
diff -c -p -d -u -r1.100 real.c
--- real.c	19 Oct 2002 23:03:21 -0000	1.100
+++ real.c	21 Oct 2002 23:58:39 -0000
@@ -81,7 +81,7 @@ static void get_zero PARAMS ((REAL_VALUE
 static void get_canonical_qnan PARAMS ((REAL_VALUE_TYPE *, int));
 static void get_canonical_snan PARAMS ((REAL_VALUE_TYPE *, int));
 static void get_inf PARAMS ((REAL_VALUE_TYPE *, int));
-static void sticky_rshift_significand PARAMS ((REAL_VALUE_TYPE *,
+static bool sticky_rshift_significand PARAMS ((REAL_VALUE_TYPE *,
 					       const REAL_VALUE_TYPE *,
 					       unsigned int));
 static void rshift_significand PARAMS ((REAL_VALUE_TYPE *,
@@ -97,7 +97,7 @@ static bool add_significands PARAMS ((RE
 				      const REAL_VALUE_TYPE *));
 static bool sub_significands PARAMS ((REAL_VALUE_TYPE *,
 				      const REAL_VALUE_TYPE *,
-				      const REAL_VALUE_TYPE *));
+				      const REAL_VALUE_TYPE *, int));
 static void neg_significand PARAMS ((REAL_VALUE_TYPE *,
 				     const REAL_VALUE_TYPE *));
 static int cmp_significands PARAMS ((const REAL_VALUE_TYPE *,
@@ -182,10 +182,9 @@ get_inf (r, sign)
 
 \f
 /* Right-shift the significand of A by N bits; put the result in the
-   significand of R.  If any one bits are shifted out, set the least
-   significant bit of R.  */
+   significand of R.  If any one bits are shifted out, return true.  */
 
-static void
+static bool
 sticky_rshift_significand (r, a, n)
      REAL_VALUE_TYPE *r;
      const REAL_VALUE_TYPE *a;
@@ -220,7 +219,7 @@ sticky_rshift_significand (r, a, n)
 	r->sig[i] = 0;
     }
 
-  r->sig[0] |= (sticky != 0);
+  return sticky != 0;
 }
 
 /* Right-shift the significand of A by N bits; put the result in the
@@ -327,15 +326,16 @@ add_significands (r, a, b)
   return carry;
 }
 
-/* Subtract the significands of A and B, placing the result in R.
-   Return true if there was carry out of the most significant word.  */
+/* Subtract the significands of A and B, placing the result in R.  CARRY is
+   true if there's a borrow incoming to the least significant word.
+   Return true if there was borrow out of the most significant word.  */
 
 static inline bool
-sub_significands (r, a, b)
+sub_significands (r, a, b, carry)
      REAL_VALUE_TYPE *r;
      const REAL_VALUE_TYPE *a, *b;
+     int carry;
 {
-  bool carry = false;
   int i;
 
   for (i = 0; i < SIGSZ; ++i)
@@ -500,7 +500,7 @@ div_significands (r, a, b)
     start:
       if (msb || cmp_significands (&u, b) >= 0)
 	{
-	  sub_significands (&u, &u, b);
+	  sub_significands (&u, &u, b, 0);
 	  set_significand_bit (r, bit);
 	}
     }
@@ -570,6 +570,7 @@ do_add (r, a, b, subtract_p)
 {
   int dexp, sign, exp;
   REAL_VALUE_TYPE t;
+  bool inexact = false;
 
   /* Determine if we need to add or subtract.  */
   sign = a->sign;
@@ -648,13 +649,13 @@ do_add (r, a, b, subtract_p)
 	  return;
 	}
 
-      sticky_rshift_significand (&t, b, dexp);
+      inexact |= sticky_rshift_significand (&t, b, dexp);
       b = &t;
     }
 
   if (subtract_p)
     {
-      if (sub_significands (r, a, b))
+      if (sub_significands (r, a, b, inexact))
 	{
 	  /* We got a borrow out of the subtraction.  That means that
 	     A and B had the same exponent, and B had the larger
@@ -671,7 +672,7 @@ do_add (r, a, b, subtract_p)
 	  /* We got carry out of the addition.  This means we need to
 	     shift the significand back down one bit and increase the
 	     exponent.  */
-	  sticky_rshift_significand (r, r, 1);
+	  inexact |= sticky_rshift_significand (r, r, 1);
 	  r->sig[SIGSZ-1] |= SIG_MSB;
 	  if (++exp > MAX_EXP)
 	    {
@@ -692,6 +693,8 @@ do_add (r, a, b, subtract_p)
      is positive.  */
   if (r->class == rvc_zero)
     r->sign = 0;
+  else
+    r->sig[0] |= inexact;
 }
 
 /* Return R = A * B.  */
@@ -1430,7 +1433,7 @@ rtd_divmod (num, den)
     start:
       if (msb || cmp_significands (num, den) >= 0)
 	{
-	  sub_significands (num, num, den);
+	  sub_significands (num, num, den, 0);
 	  q |= 1;
 	}
     }
@@ -2329,7 +2332,7 @@ round_for_format (fmt, r)
       if (shift)
 	{
 	  shift = fmt->log2_b - shift;
-	  sticky_rshift_significand (r, r, shift);
+	  r->sig[0] |= sticky_rshift_significand (r, r, shift);
 	  r->exp += shift;
 	}
     }
@@ -2355,7 +2358,7 @@ round_for_format (fmt, r)
 	    goto underflow;
 
 	  /* De-normalize the significand.  */
-	  sticky_rshift_significand (r, r, diff);
+	  r->sig[0] |= sticky_rshift_significand (r, r, diff);
 	  r->exp += diff;
 	}
     }
@@ -2395,7 +2398,7 @@ round_for_format (fmt, r)
 	      if (shift)
 		{
 		  shift = fmt->log2_b - shift;
-		  sticky_rshift_significand (r, r, shift);
+		  rshift_significand (r, r, shift);
 		  r->exp += shift;
 		  if (r->exp > emax2)
 		    goto overflow;
@@ -4307,6 +4310,51 @@ const struct real_format c4x_extended_fo
     false,
     false,
     false
+  };
+
+\f
+/* A synthetic "format" for internal arithmetic.  It's the size of the
+   internal significand minus the two bits needed for proper rounding.
+   The encode and decode routines exist only to satisfy our paranoia
+   harness.  */
+
+static void encode_internal PARAMS ((const struct real_format *fmt,
+				     long *, const REAL_VALUE_TYPE *));
+static void decode_internal PARAMS ((const struct real_format *,
+				     REAL_VALUE_TYPE *, const long *));
+
+static void
+encode_internal (fmt, buf, r)
+     const struct real_format *fmt ATTRIBUTE_UNUSED;
+     long *buf;
+     const REAL_VALUE_TYPE *r;
+{
+  memcpy (buf, r, sizeof (*r));
+}
+
+static void
+decode_internal (fmt, r, buf)
+     const struct real_format *fmt ATTRIBUTE_UNUSED;
+     REAL_VALUE_TYPE *r;
+     const long *buf;
+{
+  memcpy (r, buf, sizeof (*r));
+}
+
+const struct real_format real_internal_format = 
+  {
+    encode_internal,
+    decode_internal,
+    2,
+    1,
+    SIGNIFICAND_BITS - 2,
+    -MAX_EXP,
+    MAX_EXP,
+    true,
+    true,
+    false,
+    true,
+    true 
   };
 \f
 /* Set up default mode to format mapping for IEEE.  Everyone else has

Re: real.c implementation

[Fergus Henderson]

On 21-Oct-2002, Richard Henderson <rth@redhat.com> wrote:
> On Sun, Oct 20, 2002 at 10:00:41PM -0500, Brad Lucier wrote:
> 
> > Is exact conversion a goal for gcc?
> 
> I'd say so.

Deterministic results are a requirement for the Mercury compiler to
conform to the Mercury language specification.  I think that requires
exact conversion.

We got a bug report about this yesterday, so it is something that some of
our users do care about.  (It's not a high priority, though.  The benefit
to effort ratio is too low.)

-- 
Fergus Henderson <fjh@cs.mu.oz.au>  |  "I have always known that the pursuit
The University of Melbourne         |  of excellence is a lethal habit"
WWW: <http://www.cs.mu.oz.au/~fjh>  |     -- the last words of T. S. Garp.

Re: real.c implementation

[Brigitte Verdonk]

> From what I have just seen, the program generates a lot of warning 
> but still compiles and seems to work with g++-3.2. With 3.3 things 
> are more complicated but it should be possible to compile using 
> -fpermissive as the errors are mostly related to strstream.
> Anyway, fixing this is relatively trivial (see below for a 
> fixed version, done very quickly and without being very careful, but I 
> hope that I have not introduced bugs into the program).
> 
> In any case, here are the results for the various versions of gcc I 
> have access to (gcc-3.1, gcc-3.2.1, and gcc-3.3 20021014), and the
> source code apropriate for g++ 3.3 (I have verified that the code
> still works with gcc-3.2, but not with other compilers).
> 

Thanks a lot for your feedback on IeeeCC754 as well as for the updated
version of the source code and the logfiles for different compiler
versions!

If you agree, we would like to make your temporary fix of the source code
available on the IeeeCC754 website. For the moment, we're upgrading the
source code to include some additional functionality and at the same time
take care of all the compiler warnings.

Until our new version is ready, would it be ok that I make your
version available on the IeeeCC754 website? And can I also include the 
Logfiles you sent?

Kind regards,

Brigitte Verdonk


-----------------------------Brigitte Verdonk---------------------------------
Mathematics & Computer Science           Tel  +32 3 820.24.03
University of Antwerp (UIA)              Fax  +32 3 820.24.21
Universiteitsplein 1                     URL: http://www.uia.ac.be/u/verdonk
B2610 Antwerp (Belgium)                  Email: Brigitte.Verdonk@ua.ac.be
------------------------------------------------------------------------------

"Alice laughed: "There's no use trying," she said; "one can't 
believe impossible things."  

"I daresay you haven't had much practice," said the Queen. 
"When I was younger, I always did it for half an hour a day. 
Why, sometimes I've believed as many as six impossible things 
before breakfast."  

Lewis Carroll

Re: real.c implementation

[Stephen L Moshier]

Here is a recently developed set of IEEE test vectors that includes
the set from the thesis of Coonen, some values from ucbtest, and extensions
to quad precision.  Unfortunately the program was written for C++
version n-1 and I have C++ version n, a hopeless situation.  It might
be useful for a C++ enthusiast to work on the program.
   http://win-www.uia.ac.be/u/cant/ieeecc754.html

Testing the decimal constants extracted from the 80-bit and 128-bit
series of these vectors seemed to generate no errors that were outside the
IEEE spec limits.  I also extracted and checked all of the decimal
constants from the quad precision glibc math library, about 1700 numbers;
they all come out correctly rounded now.

Re: real.c implementation

[Stephen L Moshier]

> It failed completely for denormal numbers for which the target
> format also used a 15-bit exponent, i.e 80-bit double-extended
> and 128-bit quad precisions.
>
> we noticed two things, GCC
> did not always get denormals right

I have tried to investigate these statements as bug reports.  The
information provided is sketchy, but I think you may be referring to a
rounding feature that was left incomplete (and was so documented)
because it was never used.  The decimal-to-binary and binary-to-binary
conversions do know the intended output mode and do include the code
needed to round correctly.  On the other hand GCC's REAL_ARITHMETIC
macro should have a mode argument but does not, so the add, subtract,
multiply, divide programs default to the widest available result mode.
That created, in GCC, problems that are similar to the extra-precise
register mysteries of the Intel x86 floating point unit.

The program is used in other packages besides GCC, and its maintenance
will continue.  If you have a bug to report, please provide enough
information to reproduce the problem.

Re: real.c implementation

[Brad Lucier]

> As a practical matter, conversion according to IEEE is
> already correctly rounded over a wide domain that covers the vast
> majority of real programs, and the largest allowable error is less
> than twice the error of correct rounding; so the distinction is minor.

I suppose it's a question of how much one's programs rely on exact conversion.

One man's "defect" or "issue" is another man's disaster.

Brad

Re: real.c implementation

[Robert Dewar]

> They enable the flush to zero mode deliberately because it is faster
> than IEEE mode, or because they are running a VAX program.  In those
> cases they don't want to use mieee.

No one wants to use -mieee in real life, it is a performance disaster.
It may be necessary for strict compliance in some situations, but from
a pragmatic point of view, it is really important that any compiler
work well in practice without -mieee.

Re: real.c implementation

[Stephen L Moshier]

> I don't know the answers to these questions.  Steele and White,
> Clinger, and David Gay (his code is at netlib.org) worked out
> provably correct methods to do exact binary<->decimal conversion.
>  All of them rely on extra-precision integer arithmetic for some
> cases.

Gay's program is in the GCC java library subdirectory.  It used
to be in the C++ library also; I remember chasing a bug in it.
Dewar has mentioned that the GCC Ada front end does correctly rounded
conversions.  Glibc, the GNU system library, has been through
several decimal converters; currently it has a correctly rounded
method that uses the GNU multiple precision library (gmp).  I think
that gmp can also be used as a floating point emulator.

The C language standard references IEEE.  Some months ago the IEEE 754
revision committee posted some messages to the effect that they were
going to require a correctly rounded method, but then they seemed to
withdraw that.  As a practical matter, conversion according to IEEE is
already correctly rounded over a wide domain that covers the vast
majority of real programs, and the largest allowable error is less
than twice the error of correct rounding; so the distinction is minor.

Re: real.c implementation

[Stephen L Moshier]

>On Tue, Oct 22, 2002 at 02:11:27PM -0400, Robert Dewar wrote:
>> > We do now.  We did have code in the compiler to prevent this
>> > (unconditionally, mind!) until recently; I considered this a
>> > bug and removed it.  If you don't want to trap on denormals,
>> > then don't write denormal constants.  Simple as that.
>>
>> Well not quite as simple, you are asking people to do careful
analysis
>> of legacy code and tailor parts of it in a target dependent manner.
>
>No, we were asking them to do that before, since when they
>wrote 2**-149, we *silently* gave them 0.

If you write 10^-4000 to an IEEE double, the nearest representable
value is zero, the compiler generates zero and there is no warning.
In this other case the target machine's arithmetic is flush-to-zero.
In that system, the closest representable number to your example is zero
so the compiler should generate zero.  If the alpha people wanted a
warning, which they did not, check_float_value could give a warning.

>Anyway, if they leave nothing unchanged, they'll get a SIGFPE,
>which can be fixed immediately by using -mieee.

They enable the flush to zero mode deliberately because it is faster
than IEEE mode, or because they are running a VAX program.  In those
cases they don't want to use mieee.

Re: real.c implementation

[Andreas Schwab]

Richard Henderson <rth@redhat.com> writes:

|> On Tue, Oct 22, 2002 at 08:18:37PM +0200, Andreas Schwab wrote:
|> > GAS should be able to handle them correctly.
|> 
|> How about the output of __builtin_nan("0x12345") ?

Haven't tried that yet.

Andreas.

-- 
Andreas Schwab, SuSE Labs, schwab@suse.de
SuSE Linux AG, Deutschherrnstr. 15-19, D-90429 Nürnberg
Key fingerprint = 58CA 54C7 6D53 942B 1756  01D3 44D5 214B 8276 4ED5
"And now for something completely different."

Re: real.c implementation

[Robert Dewar]

> 
> I think the question was more along the lines of why the existing
> infrastructure wasn't improved as opposed to building your own...

The Ada front end requires arbitrary precision real arithmetic to be
implemented. This is pretty Ada specific, but of course the routines could
be used in other languages if anyone decided this was useful (I can't see
any real use).

The Ada floating-point model is pretty specific to Ada and is for instance
quite different from the Java model.

Note that the use of the arbitrary precision arithmetic in Ada has nothing
to do with the code generator per se (obviously any float constants that make
it to the code generator have restricted precision and restricted range).

Re: real.c implementation

[Eric Christopher]

On Tue, 2002-10-22 at 17:50, Robert Dewar wrote:
> > I don't understand why your argument doesn't apply to having Ada a
> > completely separate codebase to GCC?  After all, it has different
> > requirements to C, for example.
> 
> 
> I am sorry, I don't understand your point here. Clearly Ada and C can
> share large parts of a code generator, but there are requirements in
> Ada that are not satisfied by the existing GCC code generator for sure.

I think the question was more along the lines of why the existing
infrastructure wasn't improved as opposed to building your own...

I think...

-eric

-- 
Yeah, I used to play basketball...

Re: real.c implementation

[Robert Dewar]

> I don't understand why your argument doesn't apply to having Ada a
> completely separate codebase to GCC?  After all, it has different
> requirements to C, for example.


I am sorry, I don't understand your point here. Clearly Ada and C can
share large parts of a code generator, but there are requirements in
Ada that are not satisfied by the existing GCC code generator for sure.

Re: real.c implementation

[Richard Henderson]

On Tue, Oct 22, 2002 at 08:18:37PM +0200, Andreas Schwab wrote:
> GAS should be able to handle them correctly.

How about the output of __builtin_nan("0x12345") ?


r~

Re: real.c implementation

[Neil Booth]

Robert Dewar wrote:-

> > This is the part I don't understand.  On a machine with IEEE format
> > arithmetic all langauges benefit if floating-point literals are
> > generated precisely, even if the language specification does not
> > require it.
> 
> First, "generated precisely" is not well defined (and as I have pointed out
> in a previous message, for example, Ada and Java have different requirements
> for what precisely means).

I don't understand why your argument doesn't apply to having Ada a
completely separate codebase to GCC?  After all, it has different
requirements to C, for example.

Neil.

Re: real.c implementation

[Robert Dewar]

> Anyway, if they leave nothing unchanged, they'll get a SIGFPE,
> which can be fixed immediately by using -mieee.

Yes, but -mieee is a disaster on that platform as I am sure you know (the
Berkeley folks regard the Alpha as a useless piece of junk since its fpt
performance is frightful -- they start of course from the point of view
that anything less than full IEEE compliance is junk :-)

I agree it is not too terrible to have the SIGFPE and have to fix it.
Though priobably there should be an option (in particular, I would
qualify behavior based on whether -mieee is set),.

Re: real.c implementation

[Andreas Schwab]

Richard Henderson <rth@redhat.com> writes:

|> On Tue, Oct 22, 2002 at 06:38:09AM -0400, Robert Dewar wrote:
|> > One thing I am not clear on for all ports is whether any of them trust
|> > the assembler to do decimal to float conversion correctly. I know that
|> > this was the case in GCC 2.8.1, and if you do that, the game gets lost
|> > at the assembler level. It is definitely important to generate all
|> > float constants in proper machine form in hex.
|> 
|> The only place this happens is if the target allows fp constants
|> as arguments to an instruction, e.g. on m68k.  Which is probably
|> going to fail if you try to put numbers like "+Inf" or whatnot
|> there.  One could consider that a port bug that it allows this.

GAS should be able to handle them correctly.

Andreas.

-- 
Andreas Schwab, SuSE Labs, schwab@suse.de
SuSE Linux AG, Deutschherrnstr. 15-19, D-90429 Nürnberg
Key fingerprint = 58CA 54C7 6D53 942B 1756  01D3 44D5 214B 8276 4ED5
"And now for something completely different."

Re: real.c implementation

[Richard Henderson]

On Tue, Oct 22, 2002 at 02:11:27PM -0400, Robert Dewar wrote:
> > We do now.  We did have code in the compiler to prevent this
> > (unconditionally, mind!) until recently; I considered this a
> > bug and removed it.  If you don't want to trap on denormals,
> > then don't write denormal constants.  Simple as that.
> 
> Well not quite as simple, you are asking people to do careful analysis
> of legacy code and tailor parts of it in a target dependent manner.

No, we were asking them to do that before, since when they
wrote 2**-149, we *silently* gave them 0.

Anyway, if they leave nothing unchanged, they'll get a SIGFPE,
which can be fixed immediately by using -mieee.


r~

Re: real.c implementation

[Robert Dewar]

> We do now.  We did have code in the compiler to prevent this
> (unconditionally, mind!) until recently; I considered this a
> bug and removed it.  If you don't want to trap on denormals,
> then don't write denormal constants.  Simple as that.

Well not quite as simple, you are asking people to do careful analysis
of legacy code and tailor parts of it in a target dependent manner.

Re: real.c implementation

[Richard Henderson]

On Tue, Oct 22, 2002 at 06:58:31AM -0400, Robert Dewar wrote:
> Well if you define "The Right Thing" to be what Java needs, fine, but then
> this would be wrong for Ada, since Ada requires that float literals be
> rounded away from infinity (0.5 always rounds up).

Easy enough to handle with an argument to the real_from_string
function.

> What about the handling of denormals?

What about them?  They should work as advertized.

The main reason for the recent real.c rewrite was to fix this
very issue for double-extended and quad ieee formats, so if
you remember problems with denormals in the past, they are in
theory fixed now.  You should definitely re-test.

> Also, the issue of what should be done on non-IEEE machines is not clear.
> For example, on the Alpha, you probably do not want to generate denormal
> literals by default.

We do now.  We did have code in the compiler to prevent this
(unconditionally, mind!) until recently; I considered this a
bug and removed it.  If you don't want to trap on denormals,
then don't write denormal constants.  Simple as that.

As for truely non-ieee machines, well, we do the best we can.



r~

Re: real.c implementation

[Richard Henderson]

On Tue, Oct 22, 2002 at 06:38:09AM -0400, Robert Dewar wrote:
> One thing I am not clear on for all ports is whether any of them trust
> the assembler to do decimal to float conversion correctly. I know that
> this was the case in GCC 2.8.1, and if you do that, the game gets lost
> at the assembler level. It is definitely important to generate all
> float constants in proper machine form in hex.

The only place this happens is if the target allows fp constants
as arguments to an instruction, e.g. on m68k.  Which is probably
going to fail if you try to put numbers like "+Inf" or whatnot
there.  One could consider that a port bug that it allows this.

Otherwise, we do emit fp constants as integers.  See assemble_real
in varasm.c.


r~

Re: real.c implementation

[Zack Weinberg]

On Tue, Oct 22, 2002 at 06:38:09AM -0400, Robert Dewar wrote:
> 
> One thing I am not clear on for all ports is whether any of them trust
> the assembler to do decimal to float conversion correctly. I know that
> this was the case in GCC 2.8.1, and if you do that, the game gets lost
> at the assembler level. It is definitely important to generate all
> float constants in proper machine form in hex.

It is my understanding that float constants are now universally emitted as
hexadecimal literals.  However, there could be some dark corners.  Grep for
real_to_decimal.

Re: real.c implementation

[Andrew Haley]

Andrew Haley writes:
 > Matt Thomas writes:
 >  > At 01:55 AM 10/22/2002, Andrew Haley wrote:
 >  > > gcj, if course, requires exact IEEE float conversion and
 >  > >relies on gcc to do it right,
 >  > 
 >  > And what happens if the platform doesn't have IEEE floating point and
 >  > uses a different floating point format?
 > 
 > Then it can't support Java.

Sorry, that's too crude: gcj cannot, in its present state, conform to
the Java Language Specifciation on such a machine.

Andrew.

Re: real.c implementation

[Andrew Haley]

Matt Thomas writes:
 > At 01:55 AM 10/22/2002, Andrew Haley wrote:
 > > gcj, if course, requires exact IEEE float conversion and
 > >relies on gcc to do it right,
 > 
 > And what happens if the platform doesn't have IEEE floating point and
 > uses a different floating point format?

Then it can't support Java.

Andrew.

Re: real.c implementation

[Robert Dewar]

> >I agree.  gcj, if course, requires exact IEEE float conversion and
> >relies on gcc to do it right,
> 
> And what happens if the platform doesn't have IEEE floating point and
> uses a different floating point format?

Obviously an even vaguely correct Java implementation requires a target
that is at least vaguely IEEE compliant. So if you have nothing like
IEEE floating-point in sight, you have to either write an emulator (and
consider the emulator to be the target of code generation), or forget
about having an even vaguely correct Java implementation.

Of course for Ada, the situation is more complex, Ada has a precise
arithmetic model that is parametrized by the target, and carefully designed
to accomodate all usual computer floating-point arithmetic.
So "correct" for Ada is quite target dependent when it comes to floating-point.

Re: real.c implementation

[Matt Thomas]

At 01:55 AM 10/22/2002, Andrew Haley wrote:
>Richard Henderson writes:
>  > On Mon, Oct 21, 2002 at 11:42:56PM -0400, Robert Dewar wrote:
>  > > The same is true of course for Ada, but what we do in Ada is have the
>  > > front end take care of everything, we do not trust the back end of GCC
>  > > or the assembler to do conversions correctly.
>  > >
>  > > I am surprised that the Mercury compiler does not take the same 
> approach.
>  >
>  > That seems like a complete waste of effort.  Why not do it
>  > correctly, once, in the right place?
>
>I agree.  gcj, if course, requires exact IEEE float conversion and
>relies on gcc to do it right,

And what happens if the platform doesn't have IEEE floating point and
uses a different floating point format?


-- 
Matt Thomas               Internet:   matt@3am-software.com
3am Software Foundry      WWW URL:    http://www.3am-software.com/bio/matt/
Cupertino, CA             Disclaimer: I avow all knowledge of this message

Re: real.c implementation

[Robert Dewar]

> This is the part I don't understand.  On a machine with IEEE format
> arithmetic all langauges benefit if floating-point literals are
> generated precisely, even if the language specification does not
> require it.

First, "generated precisely" is not well defined (and as I have pointed out
in a previous message, for example, Ada and Java have different requirements
for what precisely means).

Second, requiring exact rounding through the entire range has very severe
implications. Be sure that is what you want before you require it. This
for instance goes far beyond the requirements of IEEE 754.

Re: real.c implementation

[Robert Dewar]

> a.  Java has somewhat strict requirements for decimal-float conversion, and
> b.  We need gcc to do The Right Thing.
> 
> Neither of these can reasonably be disputed.

Well if you define "The Right Thing" to be what Java needs, fine, but then
this would be wrong for Ada, since Ada requires that float literals be
rounded away from infinity (0.5 always rounds up).

I don't know if any other languages are clearly defined. 

What about the handling of denormals?

Also note that this means that the Java requirements are different from
those of IEEE if you are stating them correctly, since IEEE specifically
does NOT require last bit accuracy for decimal conversions. To be exactly
last bit accurate throughout the range for decimal conversions is tricky,
it requires arbitrary precision arithmetic, far beyond anything that is
in the back end of GCC today.

(at least that's my best understanding from the algorithms available
a few years ago, perhaps there are new discoveries I am unaware of).

Also, the issue of what should be done on non-IEEE machines is not clear.
For example, on the Alpha, you probably do not want to generate denormal
literals by default.

I often find that in the Java world, people think that muttering the magic
incantation "IEEE 754" answers all questions. It does not, as the original
designers of Java (who created a language not practical to implement on the
x86) learned :-)

Re: real.c implementation

[Andrew Haley]

Robert Dewar writes:

 > Second, with regard to decimal constants, we don't want to generate constants
 > in decimal in any case, because assemblers cannot be trusted. So when we did
 > generate in decimal, we noticed two things, GCC did not always get denormals
 > right (there may have been other problems, but I can't remember), and in any
 > case assemblers did not do the conversions last bit accurate with the 
 > required rounding mode.

This is the part I don't understand.  On a machine with IEEE format
arithmetic all langauges benefit if floating-point literals are
generated precisely, even if the language specification does not
require it.

Andrew.

Re: real.c implementation

[Robert Dewar]

> The reason for our surprise is that you chose not to fix gcc but to
> duplicate functionality instead.

Sorry I don't understand.

The Ada front end is required to have much greater capability than the
GCC back end, because it is required to evaluate precisely accurately
all static floating-point expressions (even if the results are far outside
the range of numbers representable in machine floating-point). That's a
semantic requirement in Ada, for example, the following is legal:

   a : constant := 1.3E+786666123;
   b : constant := 1.7E+91234234234;
   c : constant := a/b;

c must be exactly represented at compile time without any overflow.

Second, with regard to decimal constants, we don't want to generate constants
in decimal in any case, because assemblers cannot be trusted. So when we did
generate in decimal, we noticed two things, GCC did not always get denormals
right (there may have been other problems, but I can't remember), and in any
case assemblers did not do the conversions last bit accurate with the 
required rounding mode.

"fix gcc" here is not well defined, there is no clear semantic spec of what
gcc is expected to do here (most languages don't care about precise float
semantics, certainly C and C++ did not in the past, so this is not so
surprising). Even if there were a clear spec, it would be surprising if
it exactly matched the very detailed rules in Annex G of the Ada standard
in all cases.

Re: real.c implementation

[Andrew Haley]

Robert Dewar writes:
 > > I agree.  gcj, if course, requires exact IEEE float conversion and
 > > relies on gcc to do it right,
 > 
 > Well what does this mean?

For me to describe this precidely would require me to quote the Java
Language Specification at considerable length, and there's no point in
me doing that here.

 > Are you talking about the IEEE requirements for float conversion

Yes.

 > Or are you talking about last bit exact conversion (they are not the same,
 > you need arbitrary precision for last bit exact conversion).

 > And what about the issue of biased vs unbiased rounding. The IEEE standard
 > has nothing to say about this (more accurately, it requires the choice to be
 > delayed till runtime so that it is done with the current dynamic rounding,
 > but that's too silly to take seriously).

Java requires round to nearest or even.

 > Just saying "exact" is not good enough here.
 
Good enough for what?  All I was saying is that

a.  Java has somewhat strict requirements for decimal-float conversion, and
b.  We need gcc to do The Right Thing.

Neither of these can reasonably be disputed.

Andrew.

Re: real.c implementation

[Robert Dewar]

> I agree.  gcj, if course, requires exact IEEE float conversion and
> relies on gcc to do it right,

Well what does this mean?

Are you talking about the IEEE requirements for float conversion
Or are you talking about last bit exact conversion (they are not the same,
you need arbitrary precision for last bit exact conversion).

And what about the issue of biased vs unbiased rounding. The IEEE standard
has nothing to say about this (more accurately, it requires the choice to be
delayed till runtime so that it is done with the current dynamic rounding,
but that's too silly to take seriously).

Just saying "exact" is not good enough here.

One thing I am not clear on for all ports is whether any of them trust
the assembler to do decimal to float conversion correctly. I know that
this was the case in GCC 2.8.1, and if you do that, the game gets lost
at the assembler level. It is definitely important to generate all
float constants in proper machine form in hex.

Re: real.c implementation

[Andrew Haley]

Robert Dewar writes:
 > 
 > The reason for my surprise is that if the Mercury spec requires this to be
 > "done right", I would have thought that you would not simply trust the
 > backend of GCC to do things right (in our experience, at least in GCC 2,
 > it for sure does not do the "right thing" for denormals in all situations).

The reason for our surprise is that you chose not to fix gcc but to
duplicate functionality instead.

Andrew.

Re: real.c implementation

[Andrew Haley]

Richard Henderson writes:
 > On Mon, Oct 21, 2002 at 11:42:56PM -0400, Robert Dewar wrote:
 > > The same is true of course for Ada, but what we do in Ada is have the
 > > front end take care of everything, we do not trust the back end of GCC
 > > or the assembler to do conversions correctly.
 > > 
 > > I am surprised that the Mercury compiler does not take the same approach.
 > 
 > That seems like a complete waste of effort.  Why not do it 
 > correctly, once, in the right place?

I agree.  gcj, if course, requires exact IEEE float conversion and
relies on gcc to do it right,

Andrew.

Re: real.c implementation

[Robert Dewar]

> That seems like a complete waste of effort.  Why not do it 
> correctly, once, in the right place?

Well in Ada, you are required in any case to implement multiple precision
real arithmetic at compile time, since it is required that static expressions
of arbitrary precision must be evaluated exactly. Furthremore, the rules
for the exact evaluation of things like the 'Machine attribute are very
precise and "do it correctly" is quite Ada specific. Finally, under no
circumstances can you rely on the assembler to convert decimal float
values accurately, since miscelleaneous assemblers cannot be trusted
to do that right.

> I am surprised that the Mercury compiler does not take the same approach.

The reason for my surprise is that if the Mercury spec requires this to be
"done right", I would have thought that you would not simply trust the
backend of GCC to do things right (in our experience, at least in GCC 2,
it for sure does not do the "right thing" for denormals in all situations).

Of course I am all in favor of doing things as right as possible in GCC.
But it would not surprise me if different front ends had subtly different
requirements. For example, an interesting issue is whether unbiased rounding
is required for literals. In Ada 95, you are required to do biased rounding
(0.5 ULP always rounding up) for literals, but I can easily see some other
language requiring unbiased rounding, and indeed there is a discussion that
this requirement of Ada should be changed :-)

Re: real.c implementation

[Richard Henderson]

On Mon, Oct 21, 2002 at 11:42:56PM -0400, Robert Dewar wrote:
> The same is true of course for Ada, but what we do in Ada is have the
> front end take care of everything, we do not trust the back end of GCC
> or the assembler to do conversions correctly.
> 
> I am surprised that the Mercury compiler does not take the same approach.

That seems like a complete waste of effort.  Why not do it 
correctly, once, in the right place?


r~

Re: real.c implementation

[Robert Dewar]

> Deterministic results are a requirement for the Mercury compiler to
> conform to the Mercury language specification.  I think that requires
> exact conversion.

The same is true of course for Ada, but what we do in Ada is have the
front end take care of everything, we do not trust the back end of GCC
or the assembler to do conversions correctly.

I am surprised that the Mercury compiler does not take the same approach.

Re: real.c implementation

[Roger Sayle]

Brad Lucier writes:
> Richard Henderson writes:
> > Perhaps I'm being naieve.  If this is incorrect, stop me now.
> > An explanation of why this is wrong would also be appreciated.
>
> I don't know.  What I do know is that having correct intermediate
> results to 2k+2 bits after doing all the conversion stuff does guarantee
> the correct target result to k bits.  Something similar will come up
> with sqrt.

A reference on this topic recommended by Brad that I found interesting
(even if I didn't understand much of it) is available on the net:

"High Precision Division and Square Root", Alan H. Karp and Peter
Markstein, HP Labs Technical Report 93-93-42, June 1993.
http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf

Roger
--

Re: real.c implementation

[Richard Henderson]

On Fri, Oct 18, 2002 at 05:25:26PM -0500, Brad Lucier wrote:
> This may very well be correct; I did not really understand what your
> goal was before by reading the code.

Ok, that means more commentary is needed.  I'll see what I can
come up with.


r~

Re: real.c implementation

[Brad Lucier]

(Sorry, this wasn't cc'ed correctly to the gcc list when I replied
to Richard.)

> My understanding is that the guard bit is supposed to be the correct
> value of the .5ulp bit, and the sticky bit is supposed to be set if
> the result is not exact at the .5ulp position.
> 
> Thus one should get correct results for any target precision if the
> following rules are observed:
> 
> 	* The intermediate representation has _at least_ two more
> 	  bits than the target format.  Obviously.
> 
> 	* The bit 0 of the intermediate precision is sticky; all
> 	  other bits contain the proper values.
> 
> 	* When rounding to the target format, the guard bit is
> 	  read from the target's lsb-1, and the sticky bit is the
> 	  sum of all bits between lsb-2 and 0.
> 
> Perhaps I'm being naieve.  If this is incorrect, stop me now.
> An explanation of why this is wrong would also be appreciated.

This may very well be correct; I did not really understand what your
goal was before by reading the code.  I'll think about it more this
weekend.  Perhaps this is why paranoia is passing for +,-,*,/.

> If this is correct, do we get better results for decimal<->binary
> conversion if all computations are correctly rounded for, say,
> a 158-bit intermediate representation?

I don't know.  What I do know is that having correct intermediate results
to 2k+2 bits after doing all the conversion stuff does guarantee the
correct target result to k bits.  Something similar will come up with
sqrt.

> If so, does this allow us to get correct results for the worst-case
> 113-bit target format with a 126-bit intermediate format?  I.e. can
> we remove the extra word added the other week?  I'm guessing not, 
> but I don't actually know.

The problem is that unless you want to do some algorithmic tricks (along
the lines I gave in my e-mail) the easiest and safest way to compute
decimal<->binary conversions and sqrt is to calculate results to > 2k+2 bits
and then round to k bits.  I'll try to come up with some tests to see
if 113-bit decimal<->binary conversion works.

Brad