2

Printing constant 0.8804418 in Fortran

   program hello
      print *, 0.8804418
   end program hello

gives 0.880441785 while the C version

#include <stdio.h>

int main() {
  printf("%.10g", 0.8804418);
}

prints 0.8804418 on the same system (Linux x86-64 with gfortran and gcc). Why is the output different? Note that increasing precision in printf doesn't change output.

This is not a duplicate of Is floating point math broken? or similar. This question is specifically about difference in Fortran and C representation (or formatting).

Community
  • 1
  • 1
vitaut
  • 49,672
  • 25
  • 199
  • 336
  • 4
    How is this a duplicate? – vitaut Apr 28 '16 at 19:38
  • It's not represented exactly. If it was, it would have the same single and double mantissa. Very few fractions that can be written in decimal can also be written exactly in base 2. – DigitalRoss Apr 28 '16 at 19:55
  • @DigitalRoss Good point. Corrected. – vitaut Apr 28 '16 at 19:59
  • 1
    FYI: [binary64](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) `0.8804418` --> `0.8804418000000000521509946338483132421970367431640625`, exactly. binary32 --> `0.88044178485870361328125` – chux - Reinstate Monica Apr 28 '16 at 20:00
  • @chux The number of bits of the mantissa of floating point numbers is **not** equal to the number of decimal places. Simply calculated one decimal place takes up log2(10) bits, so binary64 floats have a precision of about 16 decimal places. – Jan Christoph Terasa Apr 29 '16 at 04:23
  • @Christoph Terasa [Your point](http://stackoverflow.com/questions/36923963/why-is-the-same-floating-point-constant-printed-differently-in-fortran-and-in-c?noredirect=1#comment61419204_36923963) as it applies to my [comment](http://stackoverflow.com/questions/36923963/why-is-the-same-floating-point-constant-printed-differently-in-fortran-and-in-c?noredirect=1#comment61410026_36923963) is unclear. My comment did not allude to precision. It simply stated the exact value of OP’s number once it is converted from text to `double/float`. If you disagree, please provide your exact value. – chux - Reinstate Monica Apr 29 '16 at 13:51
  • @Christoph Terasa Consider `atof("1.0")` which has an _exact_ value of `1.0`. The next higher `double` [binary64](https://en.wikipedia.org/wiki/Double-precision_floating-point_format) has an _exact_ value of `1.0000000000000002220446049250313080847263336181640625` a 53 digit number. The difference in these two is `2/pow(1,53)`. Perhaps [this](http://stackoverflow.com/questions/16839658/printf-width-specifier-to-maintain-precision-of-floating-point-value/19897395#19897395) will help address differences in precision and exact values. – chux - Reinstate Monica Apr 29 '16 at 13:51
  • @chux Sorry for the stark answer, mea culpa on that issue, thanks for clearing that up. After reading some sources on the term "exact" w.r.t float/double, I can still not wrap my head around the concept of "exactly 53 digits" when the data type itself doesn't have that amount of bits to represent that number. The difference is basically the difference between incrementing the LSB of the mantissa? – Jan Christoph Terasa Apr 29 '16 at 15:13
  • @Christoph Terasa Yes, the LSBit difference was `2/pow(1,53)` going from `1.0` to the next higher number. See [ULP](https://en.wikipedia.org/wiki/Unit_in_the_last_place). This may help: Think of a n-bit _binary_ fraction. The MSBit is the 0.5 place. Next bit is the 0.25 place. Next bit is 0.125 place. Next bit is 0.0625 place. The issue is primarily that `double` typically uses binary coding and people often express numbers in decimal. – chux - Reinstate Monica Apr 29 '16 at 16:19

1 Answers1

5

By default, Fortran's REAL number constants are single-precision; In C, however, floating point literals have double precision.

When you translate 0.8804418 to single precision and then print it as a double in C, you get 0.8804417849 printed (demo)

float x = 0.8804418f;
printf("%.10g\n", x);

Fortran's printout appears to be the same number rounded up.

Fortran's syntax for double-precision REAL numbers uses suffix d:

print *, 0.8804418d+0

This prints 0.88044180000000005 (demo).

Sergey Kalinichenko
  • 714,442
  • 84
  • 1,110
  • 1,523
  • Thanks. I was thinking that something like this was going on, but was confused because the full version of Fortran example had a `DOUBLE PRECISION` variable, but I guess there was a conversion to single precision somewhere. – vitaut Apr 28 '16 at 19:51