Very weird issue with double to int conversion inside the loop

Hi,

I have spotted some piece of code which yields different result on addition of integers and doubles depending if it was performed in the loop or manually (+= vs + bare +).

Here is the (not very clean but short code):

Code:

---

#include <stdio.h>
#include <limits.h>
#include <unistd.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

/* Own power function, when used instead of pow()
the loop returns correct result.
*/
int pow2(int a, int b){


    int base = a;
    if (!b)
        return 1;


    while (b--) {
        a *= base;
    }


    return a/base;
}


int main (void) {
   
    int c,d,tmp,tmp2;
    for (d = 0, tmp = 0, tmp2 = tmp+1; tmp < 3; tmp++, tmp2++) {
        c = 9 * (int)pow(10, tmp) * tmp2;
        printf("TEST loop iter %d\n",c);
        d += c;
    }
    // add manually each iteration
    int test = (9 * (int)pow(10, 0) * 1) + (9 * (int)pow(10, 1) * 2) + (9 * (int)pow(10, 2) * 3);
    printf("Sum added manually: %d, sum in the loop: %d \n", test, d);


    return 0;
}

---

So first of all I have tested this code on both linux and windows machines and the result are different! On windows (mingw) the result is incorrect:
Windows:

Code:

---

gcc -Wall -O2 tmp.c -o tmp && ./tmp.exe
TEST loop iter 9
TEST loop iter 180
TEST loop iter 2673
Sum added manually: 2889, sum in the loop: 2862

---

Linux:

Code:

---

gcc -Wall -O2 tmp.c -o tmp -lm && ./tmp
TEST loop iter 9
TEST loop iter 180
TEST loop iter 2700
Sum added manually: 2889, sum in the loop: 2889

---

So the as you can see the addition performed manually yields different result on mingw than when it is performed inside the for loop. Why is that?

Optimize Options (Using the GNU Compiler Collection (GCC))

Quote:

---

Originally Posted by gcc

-ffloat-store

Do not store floating-point variables in registers, and inhibit other options that might change whether a floating-point value is taken from a register or memory.

This option prevents undesirable excess precision on machines such as the 68000 where the floating registers (of the 68881) keep more precision than a double is supposed to have. Similarly for the x86 architecture. For most programs, the excess precision does only good, but a few programs rely on the precise definition of IEEE floating point. Use -ffloat-store for such programs, after modifying them to store all pertinent intermediate computations into variables.

---

Every time there is some variation in the unseen store/load of your intermediate floating point value, it's an opportunity for creeping differences in your accumulated floating point errors.

See also What Every Computer Scientist Should Know About Floating-Point Arithmetic

@Salem, thanks the -ffloat-store flag fixed the issue. Now I wonder why, I read the documentation and it says that it prevent storing doubles in dedicated floating point registers which have more precision that double is supposed to have - OK but how does it affect the cast to integer then?

Also where can I read on the machincs of casting from double to integer - what exactly is happening?

I have my suspicions: Please compile and run this program and post what is printed (on Windows) -- don't forget to link the libm library.

Code:

---

/* test.c */

#include <stdio.h>
#include <math.h>

struct dbl_s {
  unsigned long long int F:52;
  unsigned int E:11;
  unsigned int S:1;
} __attribute__((packed));

static void showDouble( double x )
{
  struct dbl_s *p = (struct dbl_s *)&x;

  printf( "%.50g (F=%llu, E=%d)\n",
    x, (unsigned long long int)p->F, (int)p->E - 1023 );
}

int main( void )
{
  double r;
  int i;

  r = 0;
  for ( i = 0; i < 3; i++ )
  {
    r = pow(10.0f, i);
    showDouble( r );
  }
}

---

Quote:

---

Originally Posted by MartinR

@Salem, thanks the -ffloat-store flag fixed the issue. Now I wonder why, I read the documentation and it says that it prevent storing doubles in dedicated floating point registers which have more precision that double is supposed to have - OK but how does it affect the cast to integer then?

Also where can I read on the machincs of casting from double to integer - what exactly is happening?

---

Have a read of: The Floating-Point Guide - What Every Programmer Should Know About Floating-Point Arithmetic

The issue is that the library is calculating pow(10,2) as something just slightly under 100.0 and it is being truncated down to 99 when it is converted to an integer.

This is also my suspicion... Since x^y (^ as potency) is calculated as e^(y*ln(x)) and, maybe, MinGW implementation of pow() can be a little bit wrong (there's a final approximation that need to be made), one or more of 10^x can be approximated, not exact...

In the latest glibc (linux) the results are correct, but I don't know if mingw libc (libm) is up to date, hence my question about the F and E fields of the double value calculated by pow(10,i), that should be:

For 10^0: F=0, E=0
For 10^1: F=1125899906842624, E=3
For 10^2: F=2533274790395904, E=6

In the equation (since these values are normalized):

Attachment 16358

Quote:

---

Originally Posted by flp1969

In the latest glibc (linux) the results are correct, but I don't know if mingw libc (libm) is up to date, hence my question about the F and E fields of the double value calculated by pow(10,i), that should be:

---

If the user is using mingw, then is it the mingw libc (libm), or is it from the Mickey$oft C Runtime Library that mingw uses?

I don't do programming on Windoze anymore, so I don't know.

Quote:

---

Originally Posted by rstanley

If the user is using mingw, then is it the mingw libc (libm), or is it from the Mickey$oft C Runtime Library that mingw uses?

I don't do programming on Windoze anymore, so I don't know.

---

Neither do I... It is a valid point, since MinGW uses MSVCRT??.dll and maybe pow() is implemented there.

As hamster_nz said, there is a problem of rounding when using floating point, specially with pow(), since the power is calculated as:
Attachment 16359
The natural logarithm is approximated, and the multiplication by y is approximated and the power of e to that expoent is approximated... So, for instance, if you open bc at the terminal and do:

Code:

---

e(2*l(10))

---

You'll get 99.9999999999999999984, not 100.
If pow() function isn't carefully crafted these rounding errors will accumulate and gives wrong results.

It's almost NEVER a good idea to deal with floating point when doing integer calculations.

PS: using (int)pow(10,i) isn't the correct way to do it... You should use (int)rint(pow(10,i)) to round the value to the nearest integer.

Sorry for late response.

Here is the output of your prosed program

Code:

---

$ gcc -Wall -O0 tmp.c -o tmp && time ./tmp.exe
tmp.c: In function 'showDouble':
tmp.c:16:23: warning: unknown conversion type character 'l' in format [-Wformat=]
  printf( "%.50g (F=%llu, E=%d)\n",
                      ^
tmp.c:16:30: warning: format '%d' expects argument of type 'int', but argument 3 has type 'long long unsigned int' [-Wformat=]
  printf( "%.50g (F=%llu, E=%d)\n",
                            ~^
                            %I64d
    x, (unsigned long long int)p->F, (int)p->E - 1023 );
        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
tmp.c:16:11: warning: too many arguments for format [-Wformat-extra-args]
  printf( "%.50g (F=%llu, E=%d)\n",
          ^~~~~~~~~~~~~~~~~~~~~~~~
1 (F=0, E=-180)
10 (F=1125899906842624, E=-180)
100 (F=2533274790395904, E=-180)

---

Any conslusion from that?

If the code is too complex for you to fix for Windows (I don't do Windows in 14 years), then forget it... The explanation in post #9 should be sufficient.

And, yet... you forget, in both codes, to link libm.