Internal representation of double variables - 3.4.6 vs 4.1.0

Internal representation of double variables - 3.4.6 vs 4.1.0

[max]

Hi gcc developers and users,

I have discovered that my code gives different results if compiled
with different gcc versions, namely 3.4.6 and 4.1.0.
Since I wanted to understand why, I compile my code again w/o any
optimization (-O0) and with debug symbols (-g).
I have found that differences (very small, 10e-12 on a 32-bit machine)
started to appear in the return value of a routine which performs
vector-vector multiplication, i.d.

double vecdot(double *v)
{
  double sum = 0;
  for(i = 0; i < n; i++)
    sum += v[i] * v[i];
  return sum;
}

even if elements of v[] are the same. Do these versions use different
"internal" representation of doubles?
I agree that the sum above is ill-conditioned, but why do different
gccs give (w/o optimization) different results?

Thanks for your help,
Max

Re: Internal representation of double variables - 3.4.6 vs 4.1.0

[Terry Frankcombe]

On Fri, 2007-03-09 at 18:29 +0100, max wrote:
> Hi gcc developers and users,
> 
> I have discovered that my code gives different results if compiled
> with different gcc versions, namely 3.4.6 and 4.1.0.
> Since I wanted to understand why, I compile my code again w/o any
> optimization (-O0) and with debug symbols (-g).
> I have found that differences (very small, 10e-12 on a 32-bit machine)
> started to appear in the return value of a routine which performs
> vector-vector multiplication, i.d.
> 
> double vecdot(double *v)
> {
>   double sum = 0;
>   for(i = 0; i < n; i++)
>     sum += v[i] * v[i];
>   return sum;
> }
> 
> even if elements of v[] are the same. Do these versions use different
> "internal" representation of doubles?
> I agree that the sum above is ill-conditioned, but why do different
> gccs give (w/o optimization) different results?
> 
> Thanks for your help,
> Max

(Slightly off-topic.  But only slightly!)

After reading this, I went off looking for a gcc option enforcing IEEE
floating point behaviour, assuming gcc was like the Intel compilers and
by default sacrificed some exactness in the floating point model for
speed, even with no optimisation.  I could find none.  So, does gcc use
a well-defined and reproducible floating-point model by default?  If
not, can one turn on strict IEEE arithmetic?

Ciao
Terry

-- 
Dr Terry Frankcombe
Physical Chemistry, Department of Chemistry
Göteborgs Universitet
SE-412 96 Göteborg Sweden
Ph: +46 76 224 0887   Skype: terry.frankcombe
<terry@chem.gu.se>

Re: Internal representation of double variables - 3.4.6 vs 4.1.0

[Brian Dessent]

Terry Frankcombe wrote:

> After reading this, I went off looking for a gcc option enforcing IEEE
> floating point behaviour, assuming gcc was like the Intel compilers and
> by default sacrificed some exactness in the floating point model for
> speed, even with no optimisation.  I could find none.  So, does gcc use
> a well-defined and reproducible floating-point model by default?  If
> not, can one turn on strict IEEE arithmetic?

That is not the problem.  The problem is that the intel x86 processor
uses an internal 80 bit representation for 'double', which normally has
only 64 bits.  This causes excess precision, which can appear as
rounding errors in the ULP, but are not really errors, just
misunderstandings of how IEE 754 works.  You can work around this with
-ffloat-store but this incurs a speed penalty as it means temporary
values can't be kept in registers but repeatedly stored and loaded from
memory.  Or you can use SSE for fp math (-mfpmath=sse) instead of 387,
if your architecture supports it.

See gcc bug 323 <http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323> for
examples of this non-bug being reported over and over again dozens of
times over the years.

Brian

Re: Internal representation of double variables - 3.4.6 vs 4.1.0

[Terry Frankcombe]

On Fri, 2007-03-09 at 10:33 -0800, Brian Dessent wrote:
> Terry Frankcombe wrote:
> 
> > After reading this, I went off looking for a gcc option enforcing IEEE
> > floating point behaviour, assuming gcc was like the Intel compilers and
> > by default sacrificed some exactness in the floating point model for
> > speed, even with no optimisation.  I could find none.  So, does gcc use
> > a well-defined and reproducible floating-point model by default?  If
> > not, can one turn on strict IEEE arithmetic?
> 
> That is not the problem.  The problem is that the intel x86 processor
> uses an internal 80 bit representation for 'double', which normally has
> only 64 bits.  This causes excess precision, which can appear as
> rounding errors in the ULP, but are not really errors, just
> misunderstandings of how IEE 754 works.  You can work around this with
> -ffloat-store but this incurs a speed penalty as it means temporary
> values can't be kept in registers but repeatedly stored and loaded from
> memory.  Or you can use SSE for fp math (-mfpmath=sse) instead of 387,
> if your architecture supports it.
> 
> See gcc bug 323 <http://gcc.gnu.org/bugzilla/show_bug.cgi?id=323> for
> examples of this non-bug being reported over and over again dozens of
> times over the years.
> 
> Brian


So the short answers are "No" and "No", right?  ;-)

My understanding of -ffloat-store is that it only works when there's
actually a variable involved.  Is that right?

(I'm used to working with Compaq and Intel Fortran compilers, where you
can set flags that say, effectively, "Use a standard, well-defined and
portable floating point arithmetic model irrespective of what the CPU's
trying to do.")

Ciao
Terry

Re: Internal representation of double variables - 3.4.6 vs 4.1.0

[Brian Dessent]

Terry Frankcombe wrote:

> My understanding of -ffloat-store is that it only works when there's
> actually a variable involved.  Is that right?

Right, it's not a guaranteed fix.

> (I'm used to working with Compaq and Intel Fortran compilers, where you
> can set flags that say, effectively, "Use a standard, well-defined and
> portable floating point arithmetic model irrespective of what the CPU's
> trying to do.")

I'd say you pretty much get that with -mfpmath=sse.  The only reason
anyone should ever use 387 is if binary compatibility with ancient
pre-SSE pentiums/K6s is required.  Otherwise, it's an ancient and creaky
fp architecture that would be best obsoleted (which they finally did in
x86_64 AFAIK.)

Brian

Re: Internal representation of double variables - 3.4.6 vs 4.1.0

[Terry Frankcombe]

On Fri, 2007-03-09 at 10:54 -0800, Brian Dessent wrote:
> Terry Frankcombe wrote:
> 
> > My understanding of -ffloat-store is that it only works when there's
> > actually a variable involved.  Is that right?
> 
> Right, it's not a guaranteed fix.
> 
> > (I'm used to working with Compaq and Intel Fortran compilers, where you
> > can set flags that say, effectively, "Use a standard, well-defined and
> > portable floating point arithmetic model irrespective of what the CPU's
> > trying to do.")
> 
> I'd say you pretty much get that with -mfpmath=sse.  The only reason
> anyone should ever use 387 is if binary compatibility with ancient
> pre-SSE pentiums/K6s is required.  Otherwise, it's an ancient and creaky
> fp architecture that would be best obsoleted (which they finally did in
> x86_64 AFAIK.)
> 
> Brian

For the world at large:

Is there any sensible reason that these options seem to buried in
processor-specific sections?  Surely the hardware floating point model
varies on more than just the x86 family and friends.  Why isn't there a
simple, global option to get a portable floating point model?
(Particularly now that gfortran is becoming more popular!)  Tying the
floating point model to a particular hardware feature/instruction set
(SSE) seems absolutely absurd to me.  It should be abstracted to a much
higher level.

Terry

RE: Internal representation of double variables - 3.4.6 vs 4.1.0

[John (Eljay) Love-Jensen]

Hi Terry,

> Is there any sensible reason that these options seem to buried in processor-specific sections?

My guess is because they are processor-specific options.

> Surely the hardware floating point model varies on more than just the x86 family and friends.

Yes.  For example, not all 680x0 machines has a FPU.

> Why isn't there a simple, global option to get a portable floating point model?

Because there are different FPU facilities on different platforms that may have a different, non-portable floating-point model.

> Particularly now that gfortran is becoming more popular!

If the Fortran users express a desire for some Fortran level floating-point compliance to a portable floating-point model, perhaps the maintainer of the gfortran front-end tool-chain driver will add a flag to that effect.

I imagine such a flag would have these guarantees:
+ behind the scenes, specifies the processor-specific option to insure floating-point portability
+ may impose severe performance penalties (one-to-two orders of magnitude) on some platforms
+ may cause the compile to fail if the floating-point constraint cannot be fulfilled (which is probably a good thing)

> Tying the floating point model to a particular hardware feature/instruction set (SSE) seems absolutely absurd to me.  It should be abstracted to a much higher level.

It's a trade-off, without a doubt.

Java has a facility that allows strong guarantees on floating-point fidelity.  And Mathematica allows very powerful expression of symbolic and mathematical manipulation and expression.  And there's always Lisp.

Sincerely,
--Eljay