Re: Rounding errors using doubles?

Re: Rounding errors using doubles?

[Ross Smith]

From: Sam Lantinga <slouken@devolution.com>
>
>I have to thank everyone here on this list for your responsiveness.
>The problem was caused by unexpected values in code using doubles
>that was ported from Windows.  Apparently, by default, Windows uses 
>53 bits of precision for it's double operations.  Linux uses 64 bits.
>
>Since the code relies on the exact behavior of doubles quite extensively,
>we are setting the precision of the FPU to 53 bits using the fldcw
>instruction. :)
>
>BTW, does anyone know what precision is used on the PPC and Alpha
>architectures?

It's not really a difference between operating systems or (to a first
approximation) CPUs.

There's a standard called IEC 559 (formerly IEEE 754), which specifies
three standard floating-point arithmetic formats. Two of them are
32-bit and 64-bit formats (with 24-bit and 53-bit precision,
respectively); modern C/C++ compilers almost universally equate these
to float and double. The third is an 80-bit (64-bit precision) format
used for intermediate results in multi-step calculations.

The Intel X86 processors have can perform arithmetic in all three
modes, although it always works internally in 80-bit mode and
inserts automatic conversions for the other two. I don't know
anything about PPCs or Alphas, but the IEC standard is pretty much
universal now, and I'd be mildly amazed if either of them differed
in any important ways.

Exactly how floating-point arithmetic is done is a function of the
compiler, not the operating system. You didn't say which compiler you
were using on Windows, but Microsoft Visual C++ is the most likely.
Both MSVC and EGCS perform intermediate calculations in 80-bit internal
registers wherever possible, but allow values to spill into 64-bit
memory when the compiler can't manage to fit the entire calculation
into registers (a common problem on the register-poor X86
architecture). (The details of exactly how EGCS should handle this
were the subject of heated debate on this list not so long ago.)

It looks like what happened was that, at some critical point in your
program, MSVC allowed a value to spill to 64 bits (thus truncating it
to 53-bit precision) while EGCS was able to keep it in an 80-bit
register (retaining 64-bit precision). (Incidentally, because people
tend to attach far too much importance to this sort of observation, I
hasten to add that this implies nothing about the relative performance
of the two compilers in general; perhaps they were using different
optimisation settings, or perhaps they simple made different decisions
about which intermediate results should be kept in registers and which
sacrificed.)

Both compilers have options to force pure 64-bit arithmetic (53-bit
precision) throughout (at some cost in speed): -ffloat-store on EGCS,
/Op on MSVC.

--
Ross Smith ................................... mailto:ross.s@ihug.co.nz
.............. The Internet Group, Auckland, New Zealand ..............
         "The award for the Most Effective Promotion of Linux
         goes to Microsoft."             -- Nicholas Petreley

Re: Rounding errors using doubles?

[craig]

>Both compilers have options to force pure 64-bit arithmetic (53-bit
>precision) throughout (at some cost in speed): -ffloat-store on EGCS,
>/Op on MSVC.

Take -ffloat-store off that list: it doesn't do the above.  It affects
only explicit *stores* to *variables*, not intermediate coputations.  egcs
provides no way to force 64-bit FP arithmetic.  In fact, I'd be surprised
if -ffloat-store "improved" the behavior of the (buggy) program originally
submitted in this thread.

Other than that, thanks for the decent explanation -- if this isn't
in the FAQ, it should be, and if it is, we should refrain from
repeating it, and point people there instead!

        tq vm, (burley)

Re: Rounding errors using doubles?

[craig]

BTW, <slouken@devolution.com> wasn't on the recipients list for Ross's
or my previous email.  I took the liberty of forwarding both of them to
him, and adding his email address to *this* email's list of recipients

So, future responders to this thread should ensure his address is
on the list of recipients, e.g. by starting the reply to this email,
or (presumably) one prior to Ross's, but not Ross's or my previous one.

(<slouken@devolution.com> is the originator of this thread.)

        tq vm, (burley)

Re: Rounding errors using doubles?

[craig]

BTW, <slouken@devolution.com> wasn't on the recipients list for Ross's
or my previous email.  I took the liberty of forwarding both of them to
him, and adding his email address to *this* email's list of recipients

So, future responders to this thread should ensure his address is
on the list of recipients, e.g. by starting the reply to this email,
or (presumably) one prior to Ross's, but not Ross's or my previous one.

(<slouken@devolution.com> is the originator of this thread.)

        tq vm, (burley)

Re: Rounding errors using doubles?

[craig]

>Both compilers have options to force pure 64-bit arithmetic (53-bit
>precision) throughout (at some cost in speed): -ffloat-store on EGCS,
>/Op on MSVC.

Take -ffloat-store off that list: it doesn't do the above.  It affects
only explicit *stores* to *variables*, not intermediate coputations.  egcs
provides no way to force 64-bit FP arithmetic.  In fact, I'd be surprised
if -ffloat-store "improved" the behavior of the (buggy) program originally
submitted in this thread.

Other than that, thanks for the decent explanation -- if this isn't
in the FAQ, it should be, and if it is, we should refrain from
repeating it, and point people there instead!

        tq vm, (burley)

Re: Rounding errors using doubles?

[Ross Smith]

From: Sam Lantinga <slouken@devolution.com>
>
>I have to thank everyone here on this list for your responsiveness.
>The problem was caused by unexpected values in code using doubles
>that was ported from Windows.  Apparently, by default, Windows uses 
>53 bits of precision for it's double operations.  Linux uses 64 bits.
>
>Since the code relies on the exact behavior of doubles quite extensively,
>we are setting the precision of the FPU to 53 bits using the fldcw
>instruction. :)
>
>BTW, does anyone know what precision is used on the PPC and Alpha
>architectures?

It's not really a difference between operating systems or (to a first
approximation) CPUs.

There's a standard called IEC 559 (formerly IEEE 754), which specifies
three standard floating-point arithmetic formats. Two of them are
32-bit and 64-bit formats (with 24-bit and 53-bit precision,
respectively); modern C/C++ compilers almost universally equate these
to float and double. The third is an 80-bit (64-bit precision) format
used for intermediate results in multi-step calculations.

The Intel X86 processors have can perform arithmetic in all three
modes, although it always works internally in 80-bit mode and
inserts automatic conversions for the other two. I don't know
anything about PPCs or Alphas, but the IEC standard is pretty much
universal now, and I'd be mildly amazed if either of them differed
in any important ways.

Exactly how floating-point arithmetic is done is a function of the
compiler, not the operating system. You didn't say which compiler you
were using on Windows, but Microsoft Visual C++ is the most likely.
Both MSVC and EGCS perform intermediate calculations in 80-bit internal
registers wherever possible, but allow values to spill into 64-bit
memory when the compiler can't manage to fit the entire calculation
into registers (a common problem on the register-poor X86
architecture). (The details of exactly how EGCS should handle this
were the subject of heated debate on this list not so long ago.)

It looks like what happened was that, at some critical point in your
program, MSVC allowed a value to spill to 64 bits (thus truncating it
to 53-bit precision) while EGCS was able to keep it in an 80-bit
register (retaining 64-bit precision). (Incidentally, because people
tend to attach far too much importance to this sort of observation, I
hasten to add that this implies nothing about the relative performance
of the two compilers in general; perhaps they were using different
optimisation settings, or perhaps they simple made different decisions
about which intermediate results should be kept in registers and which
sacrificed.)

Both compilers have options to force pure 64-bit arithmetic (53-bit
precision) throughout (at some cost in speed): -ffloat-store on EGCS,
/Op on MSVC.

--
Ross Smith ................................... mailto:ross.s@ihug.co.nz
.............. The Internet Group, Auckland, New Zealand ..............
         "The award for the Most Effective Promotion of Linux
         goes to Microsoft."             -- Nicholas Petreley

Re: Rounding errors using doubles?

[Sam Lantinga]

> >Both compilers have options to force pure 64-bit arithmetic (53-bit
> >precision) throughout (at some cost in speed): -ffloat-store on EGCS,
> >/Op on MSVC.

> Take -ffloat-store off that list: it doesn't do the above.  It affects
> only explicit *stores* to *variables*, not intermediate coputations.  egcs
> provides no way to force 64-bit FP arithmetic.  In fact, I'd be surprised
> if -ffloat-store "improved" the behavior of the (buggy) program originally
> submitted in this thread.

Interestingly enough, the original project file (MS VC++ 5.0) did have the
/Op option, and -ffloat-store didn't force 64-bit FP arithmetic.  The
following code does force 64-bit arithmetic on Pentium style CPUs:

// Initialize the FPU at 53-bit precision to mimic Windows /Op behavior
// This code doesn't parse with the G++ parser because of the "::", but
// it works fine as a C file.
void Setup_FPU(void)
{
#ifdef i386
    unsigned short flags;

    asm volatile ("fnstcw %0":"=m" (flags));
    flags &= ~0x0300;       /* Clear precision flags */
    flags |=  0x0200;       /* Set 53-bit precision */
    asm volatile ("fldcw %0"::"m" (flags));
#endif
}

If this doesn't work on other x86 CPUs, I'd like to know.
Linux saves and restores the FPU flags on context switch, so this is
safe to call in an individual program.

BTW, I am on the list. :)

	-Sam Lantinga				(slouken@devolution.com)

Lead Programmer, Loki Entertainment Software
--
Author of Simple DirectMedia Layer -
	http://www.devolution.com/~slouken/SDL/
--

Re: Rounding errors using doubles?

[craig]

>If this doesn't work on other x86 CPUs, I'd like to know.
>Linux saves and restores the FPU flags on context switch, so this is
>safe to call in an individual program.

I've seen that kind of solution recommended before, and most/all
people who've used it seem to be happy with it.

The reason I don't like it as a general solution is that we don't
really know whether's there's any low-level numerical code lurking
down there that *assumes* the FPU modes are set to 80 bits, and
won't work correctly otherwise.  It's therefore not a robust solution,
even though it might happen to be 100% effective.  Whereas, an
option that's like -ffloat-store but affects *all* intermediate
computations *would* be robust, even if not 100% effective (in that
code not compiled with that option would still use 80-bit computations).

>BTW, I am on the list. :)

That's good to know: people can ignore my earlier email then.

        tq vm, (burley)

Re: Rounding errors using doubles?

[craig]

>If this doesn't work on other x86 CPUs, I'd like to know.
>Linux saves and restores the FPU flags on context switch, so this is
>safe to call in an individual program.

I've seen that kind of solution recommended before, and most/all
people who've used it seem to be happy with it.

The reason I don't like it as a general solution is that we don't
really know whether's there's any low-level numerical code lurking
down there that *assumes* the FPU modes are set to 80 bits, and
won't work correctly otherwise.  It's therefore not a robust solution,
even though it might happen to be 100% effective.  Whereas, an
option that's like -ffloat-store but affects *all* intermediate
computations *would* be robust, even if not 100% effective (in that
code not compiled with that option would still use 80-bit computations).

>BTW, I am on the list. :)

That's good to know: people can ignore my earlier email then.

        tq vm, (burley)

Re: Rounding errors using doubles?

[Sam Lantinga]

> >Both compilers have options to force pure 64-bit arithmetic (53-bit
> >precision) throughout (at some cost in speed): -ffloat-store on EGCS,
> >/Op on MSVC.

> Take -ffloat-store off that list: it doesn't do the above.  It affects
> only explicit *stores* to *variables*, not intermediate coputations.  egcs
> provides no way to force 64-bit FP arithmetic.  In fact, I'd be surprised
> if -ffloat-store "improved" the behavior of the (buggy) program originally
> submitted in this thread.

Interestingly enough, the original project file (MS VC++ 5.0) did have the
/Op option, and -ffloat-store didn't force 64-bit FP arithmetic.  The
following code does force 64-bit arithmetic on Pentium style CPUs:

// Initialize the FPU at 53-bit precision to mimic Windows /Op behavior
// This code doesn't parse with the G++ parser because of the "::", but
// it works fine as a C file.
void Setup_FPU(void)
{
#ifdef i386
    unsigned short flags;

    asm volatile ("fnstcw %0":"=m" (flags));
    flags &= ~0x0300;       /* Clear precision flags */
    flags |=  0x0200;       /* Set 53-bit precision */
    asm volatile ("fldcw %0"::"m" (flags));
#endif
}

If this doesn't work on other x86 CPUs, I'd like to know.
Linux saves and restores the FPU flags on context switch, so this is
safe to call in an individual program.

BTW, I am on the list. :)

	-Sam Lantinga				(slouken@devolution.com)

Lead Programmer, Loki Entertainment Software
--
Author of Simple DirectMedia Layer -
	http://www.devolution.com/~slouken/SDL/
--

Re: Rounding errors using doubles?

[Sam Lantinga]

> On Thu, Mar 18, 1999 at 03:41:32PM -0800, Sam Lantinga wrote:
> > On a Pentium II, the following code gives an incorrect result when using
> > doubles, but a correct result when using floats:

> I think this is normal.  You should probably learn why at
> http://www.linuxsupportline.com/~billm/ .

I have to thank everyone here on this list for your responsiveness.
The problem was caused by unexpected values in code using doubles
that was ported from Windows.  Apparently, by default, Windows uses 
53 bits of precision for it's double operations.  Linux uses 64 bits.

Since the code relies on the exact behavior of doubles quite extensively,
we are setting the precision of the FPU to 53 bits using the fldcw
instruction. :)

BTW, does anyone know what precision is used on the PPC and Alpha
architectures?

Thanks for your help!

	-Sam Lantinga				(slouken@devolution.com)

Lead Programmer, Loki Entertainment Software
--
Author of Simple DirectMedia Layer -
	http://www.devolution.com/~slouken/SDL/
--

Re: Rounding errors using doubles?

[Dima Volodin]

Sam Lantinga wrote:

> I have to thank everyone here on this list for your responsiveness.
> The problem was caused by unexpected values in code using doubles
> that was ported from Windows.  Apparently, by default, Windows uses
> 53 bits of precision for it's double operations.  Linux uses 64 bits.
>
> Since the code relies on the exact behavior of doubles quite extensively,
> we are setting the precision of the FPU to 53 bits using the fldcw
> instruction. :)
>
> BTW, does anyone know what precision is used on the PPC and Alpha
> architectures?

Man, I'm telling you - your code is _broken_. It's just sheer luck you see it
working at all.

>         -Sam Lantinga                           (slouken@devolution.com)

Dima

Re: Rounding errors using doubles?

[Dima Volodin]

Sam Lantinga wrote:

> I have to thank everyone here on this list for your responsiveness.
> The problem was caused by unexpected values in code using doubles
> that was ported from Windows.  Apparently, by default, Windows uses
> 53 bits of precision for it's double operations.  Linux uses 64 bits.
>
> Since the code relies on the exact behavior of doubles quite extensively,
> we are setting the precision of the FPU to 53 bits using the fldcw
> instruction. :)
>
> BTW, does anyone know what precision is used on the PPC and Alpha
> architectures?

Man, I'm telling you - your code is _broken_. It's just sheer luck you see it
working at all.

>         -Sam Lantinga                           (slouken@devolution.com)

Dima

Re: Rounding errors using doubles?

[David Edelsohn]

>>>>> Sam Lantinga writes:

Sam> BTW, does anyone know what precision is used on the PPC and Alpha
Sam> architectures?

	Are you referring to wordsize or significand bits?  I think that
the following table will be useful:

/* The number of useful digits depends upon the size of the     */
/* hardware floating point word.  Several word sizes            */
/* and the corresponding digits are:                            */
/*     Word size   Significand bits   Decimal digits            */
/*        32            24                 9                    */
/*        64            53                17                    */
/*        80            64                21                    */
/*       128           106                33                    */
/*       128           113                36                    */

	PowerPC supports 32-bit wordsize floats and 64-bit wordsize
doubles.  I assume that your reference to 53-bits describes significand
bits in which case PowerPC supports that with its 64-bit wordsize.

David

Re: Rounding errors using doubles?

[David Edelsohn]

>>>>> Sam Lantinga writes:

Sam> BTW, does anyone know what precision is used on the PPC and Alpha
Sam> architectures?

	Are you referring to wordsize or significand bits?  I think that
the following table will be useful:

/* The number of useful digits depends upon the size of the     */
/* hardware floating point word.  Several word sizes            */
/* and the corresponding digits are:                            */
/*     Word size   Significand bits   Decimal digits            */
/*        32            24                 9                    */
/*        64            53                17                    */
/*        80            64                21                    */
/*       128           106                33                    */
/*       128           113                36                    */

	PowerPC supports 32-bit wordsize floats and 64-bit wordsize
doubles.  I assume that your reference to 53-bits describes significand
bits in which case PowerPC supports that with its 64-bit wordsize.

David

Re: Rounding errors using doubles?

[Sam Lantinga]

> On Thu, Mar 18, 1999 at 03:41:32PM -0800, Sam Lantinga wrote:
> > On a Pentium II, the following code gives an incorrect result when using
> > doubles, but a correct result when using floats:

> I think this is normal.  You should probably learn why at
> http://www.linuxsupportline.com/~billm/ .

I have to thank everyone here on this list for your responsiveness.
The problem was caused by unexpected values in code using doubles
that was ported from Windows.  Apparently, by default, Windows uses 
53 bits of precision for it's double operations.  Linux uses 64 bits.

Since the code relies on the exact behavior of doubles quite extensively,
we are setting the precision of the FPU to 53 bits using the fldcw
instruction. :)

BTW, does anyone know what precision is used on the PPC and Alpha
architectures?

Thanks for your help!

	-Sam Lantinga				(slouken@devolution.com)

Lead Programmer, Loki Entertainment Software
--
Author of Simple DirectMedia Layer -
	http://www.devolution.com/~slouken/SDL/
--

Rounding errors using doubles?

[Sam Lantinga]

On a Pentium II, the following code gives an incorrect result when using
doubles, but a correct result when using floats:

#include <stdio.h>

#ifdef USE_FLOAT
/* This works */
typedef float percentage_type;
#else
/* This doesn't */
typedef double percentage_type;
#endif

main()
{
        int value = 86;
        percentage_type percentage;

        printf("Starting Percentage is: %d\n", value);
        percentage = ((percentage_type)value)/100.0;
        printf("Ending Percentage is: %d\n", (int)(percentage * 100.0));
}

Incorrect output:
Starting Percentage is: 86
Ending Percentage is: 85

This is VERY bad for our software.  An analysis of the offending assembly
output shows only one difference between the broken and the working versions:

--- t-bad.s     Thu Mar 18 13:32:32 1999
+++ t.s Thu Mar 18 13:33:40 1999
@@ -28,8 +28,8 @@
        fstp %st(0)
        fldl .LC1
        fdivrp %st,%st(1)
-       fstpl -12(%ebp)
-       fldl -12(%ebp)
+       fstps -8(%ebp)
+       flds -8(%ebp)
        fldl .LC1
        fstp %st(0)
        fldl .LC1

This is using egcs 1.1.1.  Any ideas?

Thanks,
	-Sam Lantinga				(slouken@devolution.com)

Lead Programmer, Loki Entertainment Software
--
Author of Simple DirectMedia Layer -
	http://www.devolution.com/~slouken/SDL/
--

Re: Rounding errors using doubles?

[Dima Volodin]

Sam Lantinga wrote:

> On a Pentium II, the following code gives an incorrect result when using
> doubles, but a correct result when using floats:
>
> #include <stdio.h>
>
> #ifdef USE_FLOAT
> /* This works */
> typedef float percentage_type;
> #else
> /* This doesn't */
> typedef double percentage_type;
> #endif
>
> main()
> {
>         int value = 86;
>         percentage_type percentage;
>
>         printf("Starting Percentage is: %d\n", value);
>         percentage = ((percentage_type)value)/100.0;
>         printf("Ending Percentage is: %d\n", (int)(percentage * 100.0));

Your code is broken here. What you need is

printf("Ending Percentage is: %d\n", (int) floor(percentage * 100.0 + 0.5));

And, BTW, converting to int is _not_ rounding.

>         -Sam Lantinga

Dima

Re: Rounding errors using doubles?

[Dima Volodin]

Sam Lantinga wrote:

> On a Pentium II, the following code gives an incorrect result when using
> doubles, but a correct result when using floats:
>
> #include <stdio.h>
>
> #ifdef USE_FLOAT
> /* This works */
> typedef float percentage_type;
> #else
> /* This doesn't */
> typedef double percentage_type;
> #endif
>
> main()
> {
>         int value = 86;
>         percentage_type percentage;
>
>         printf("Starting Percentage is: %d\n", value);
>         percentage = ((percentage_type)value)/100.0;
>         printf("Ending Percentage is: %d\n", (int)(percentage * 100.0));

Your code is broken here. What you need is

printf("Ending Percentage is: %d\n", (int) floor(percentage * 100.0 + 0.5));

And, BTW, converting to int is _not_ rounding.

>         -Sam Lantinga

Dima

Re: Rounding errors using doubles?

[Sylvain Pion]

On Thu, Mar 18, 1999 at 03:41:32PM -0800, Sam Lantinga wrote:
> On a Pentium II, the following code gives an incorrect result when using
> doubles, but a correct result when using floats:

I think this is normal.  You should probably learn why at
http://www.linuxsupportline.com/~billm/ .

-- 
Sylvain

Re: Rounding errors using doubles?

[Sylvain Pion]

On Thu, Mar 18, 1999 at 03:41:32PM -0800, Sam Lantinga wrote:
> On a Pentium II, the following code gives an incorrect result when using
> doubles, but a correct result when using floats:

I think this is normal.  You should probably learn why at
http://www.linuxsupportline.com/~billm/ .

-- 
Sylvain

Re: Rounding errors using doubles?

[Joe Buck]

> On a Pentium II, the following code gives an incorrect result when using
> doubles, but a correct result when using floats:

Take it to a beginner's programming group.  (Hint: .86 can't be
represented exactly in floating point, the nearest IEEE float is
a hair above .86, the nearest IEEE double is a hair below .86;
in C, conversion from float or double to integer truncates, it does not
round).

> This is VERY bad for our software.

Any compiler will be very bad for your software, until you learn how
to use floating point.

Re: Rounding errors using doubles?

[Joe Buck]

I wrote:

> Take it to a beginner's programming group....
> Any compiler will be very bad for your software, until you learn how
> to use floating point.

I apologize, my language was unnecessarily harsh.  (It does bug me when
the compiler gets blamed when the user is at fault, but providers of
compilers should be used to that).

Re: Rounding errors using doubles?

[Joe Buck]

I wrote:

> Take it to a beginner's programming group....
> Any compiler will be very bad for your software, until you learn how
> to use floating point.

I apologize, my language was unnecessarily harsh.  (It does bug me when
the compiler gets blamed when the user is at fault, but providers of
compilers should be used to that).

Re: Rounding errors using doubles?

[Joe Buck]

> On a Pentium II, the following code gives an incorrect result when using
> doubles, but a correct result when using floats:

Take it to a beginner's programming group.  (Hint: .86 can't be
represented exactly in floating point, the nearest IEEE float is
a hair above .86, the nearest IEEE double is a hair below .86;
in C, conversion from float or double to integer truncates, it does not
round).

> This is VERY bad for our software.

Any compiler will be very bad for your software, until you learn how
to use floating point.

Rounding errors using doubles?

[Sam Lantinga]

On a Pentium II, the following code gives an incorrect result when using
doubles, but a correct result when using floats:

#include <stdio.h>

#ifdef USE_FLOAT
/* This works */
typedef float percentage_type;
#else
/* This doesn't */
typedef double percentage_type;
#endif

main()
{
        int value = 86;
        percentage_type percentage;

        printf("Starting Percentage is: %d\n", value);
        percentage = ((percentage_type)value)/100.0;
        printf("Ending Percentage is: %d\n", (int)(percentage * 100.0));
}

Incorrect output:
Starting Percentage is: 86
Ending Percentage is: 85

This is VERY bad for our software.  An analysis of the offending assembly
output shows only one difference between the broken and the working versions:

--- t-bad.s     Thu Mar 18 13:32:32 1999
+++ t.s Thu Mar 18 13:33:40 1999
@@ -28,8 +28,8 @@
        fstp %st(0)
        fldl .LC1
        fdivrp %st,%st(1)
-       fstpl -12(%ebp)
-       fldl -12(%ebp)
+       fstps -8(%ebp)
+       flds -8(%ebp)
        fldl .LC1
        fstp %st(0)
        fldl .LC1

This is using egcs 1.1.1.  Any ideas?

Thanks,
	-Sam Lantinga				(slouken@devolution.com)

Lead Programmer, Loki Entertainment Software
--
Author of Simple DirectMedia Layer -
	http://www.devolution.com/~slouken/SDL/
--