numerical error in C program

gopalpathak2 · 08-04-2009, 02:16 AM

if a=1946.498046875 & b=0.00003814697265625
then z=a-b=1946.498046875 why? i.e. in my C program it do not give correct results

neonsignal · 08-04-2009, 03:19 AM

Floating point numbers can only be represented with a certain amount of precision.

In C, you would normally use the 'double' type for floating point numbers, which has more precision than the 'float' type. This will improve the result in your example. (use "%lf" when formatting the output in printf [ignore this formatting, see below, %lf is only for scanf]).

pixellany · 08-04-2009, 08:15 AM

Moved: This thread is more suitable in <Programming> and has been moved accordingly to help your thread/question get the exposure it deserves.

wje_lq · 08-04-2009, 12:42 PM

Quote:

Originally Posted by neonsignal

Floating point numbers can only be represented with a certain amount of precision.

In C, you would normally use the 'double' type for floating point numbers, which has more precision than the 'float' type.

Yes, this fixes the problem.

Quote:

Originally Posted by neonsignal

(use "%lf" when formatting the output in printf).

The man page does not assign any meaning to %lf. %f is already adequate for variables of type double. To demonstrate, you'll get this output:

Code:

fz is 1946.49804687500000000000000000000000000
dz is 1946.49808502197265625000000000000000000

if you run this bash script:

Code:

#!/bin/sh
cat > 1.c <<EOD; gcc -Wall 1.c -o 1; ./1
#include <stdio.h>
int main(void)
{
  float  fa, fb, fz;
  double da, db, dz;

  fa=1946.498046875; fb=0.00003814697265625; fz=fa+fb;
  da=1946.498046875; db=0.00003814697265625; dz=da+db;

  printf("fz is %40.35f\n",fz);
  printf("dz is %40.35f\n",dz);

  return 0;

} /* main() */
EOD

If you need even greater precision, use variables of type long double. Constants should be followed by the letter "L", either upper or lower case. (I recommend you use upper case "L", because lower case "l" looks too much like the numeral 1.) The format is "Lf" (upper case; lower doesn't work).

johnsfine · 08-04-2009, 01:32 PM

Quote:

Originally Posted by wje_lq

If you need even greater precision, use variables of type long double.

IIUC, most architectures including x86_64 don't have long double support any more precision than double.

If you want more precision than double, use the GNU Multiple Precision library:
http://gmplib.org/

Sergei Steshenko · 08-04-2009, 04:22 PM

Quote:

Originally Posted by gopalpathak2

if a=1946.498046875 & b=0.00003814697265625
then z=a-b=1946.498046875 why? i.e. in my C program it do not give correct results

Start from here: http://en.wikipedia.org/wiki/Floating_point and rethink the whole issue of (non)periodic fractions WRT radix.

neonsignal · 08-04-2009, 10:24 PM

Quote:

most architectures including x86_64 don't have long double support any more precision than double

The 8087 architecture supports an 80 bit float, which is used for 'long double' by the better compilers.

johnsfine · 08-05-2009, 06:54 AM

Quote:

Originally Posted by neonsignal

The 8087 architecture supports an 80 bit float, which is used for 'long double' by the better compilers.

Do you know a command line switch or other method to make gcc for x86_64 use the 8087 for 80 bit long double?

If I understood the documentation right, there is no way to get the Intel compiler (that I normally use) to use the 8087 in x86_64 code.

I also failed (but after not trying very hard) to get long doubles to have more precision than doubles using the Microsoft C++ compiler in x86_64 architecture.

Maybe none of those are the "better" compilers.

IIRC, most 32 bit x86 compilers support true 80 bit long doubles. 32 bit x86 may be the most popular architecture, but it isn't what I meant by most architectures.

neonsignal · 08-05-2009, 08:36 AM

Quote:

Do you know a command line switch or other method to make gcc for x86_64 use the 8087 for 80 bit long double?

Not sure that it is possible (or at least, practical), because some of the 64 bit instructions share state with the 8087 stack.

For the Itanium architecture, there are facilities for doing higher precision multiply and accumulate (which would be useful eg for signal processing apps). But you'd have to resort to assembler for that in most compilers. The HP-UX compiler provides an extended floating point type for this purpose (which uses the 82 bit precision in the hardware). There has been discussion about how this might be done in gcc.

As you suggest, perhaps the GMP (or the APREC) libraries are the option for most applications that need higher precision.

Sergei Steshenko · 08-05-2009, 05:06 PM

Quote:

Originally Posted by neonsignal

Not sure that it is possible (or at least, practical), because some of the 64 bit instructions share state with the 8087 stack.

For the Itanium architecture, there are facilities for doing higher precision multiply and accumulate (which would be useful eg for signal processing apps). But you'd have to resort to assembler for that in most compilers. The HP-UX compiler provides an extended floating point type for this purpose (which uses the 82 bit precision in the hardware). There has been discussion about how this might be done in gcc.

As you suggest, perhaps the GMP (or the APREC) libraries are the option for most applications that need higher precision.

It is not possible - hardware (x86_64) doesn't support it.