While Loop Problem in MinGW and only MinGW

[Caglow]

I've got a problem with this piece of iterating code:

/* Iterate to apparent position */
while (a1 != aT) {
	/* <-- assume a1 == aT here after the 3rd iteration */
	a1 = aT;
	/* Recalculate position */
	Venus_C(JDE - aT);
	/* Calculate */
	ep = sqrt(pow(Venus_x(), 2) + pow(Venus_y(), 2)
			+ pow(Venus_z(), 2));
	aT = 0.0057755183 * ep;
}

It's part of an astronomical library which is used in two of my projects. Assume on line 3 that after the 3rd iteration, the two variables, a1 and aT, both of which are type double, are identical. The only problem is...the loop never exits. I even tried putting a printf in there to show that the two variables are the same and yet, it still continues with the loop.

I had no idea how GDB works but when I had it on, I can say that it was of no help at all. Please do not ask why...it just did not help. I decided to use Visual C++'s debugger. That would've helped...except I found out that when I compiled with VC++, it worked as expected. I tried using Cygwin and, remarkably, that also worked. So with that, I pinpointed the problem down to MinGW. Since Cygwin was GCC 4.3.2 rather than 4.4.0, I assumed that GCC 4.4.0 was the culprit. No such luck. Even after downgrading down as far as GCC 3.3.3, I had the exact same problem...the loop just doesn't want to exit.

I would use Cygwin except I can't because one of my 2 projects is not free software (though the other is) along with the fact that I've had NOTHING but trouble when I try to use it with anything else.

Also, it's not just any while loop that fails. ONLY that one fails actually. The following worked perfectly as expected:

double a = 0, b = 20;
while (a != b) {
	a+=0.01;
}

Using Code::Blocks IDE 8.02 on Windows XP SP3 with MinGW 5.1.6 and GCC 4.4.0.

[Dave Sinkula]

http://c-faq.com/fp/fpequal.html

[Caglow]

Doesn't work...I still have the exact same problem as before. Even when I set epsilon to 0.01, a value far too imprecise for my measurements, it still cannot exit when it should.

[Dave Sinkula]

I even tried putting a printf in there to show that the two variables are the same

To how many digits of precision are you printing? Perhaps also consider examining the bytes that contain the double value -- maybe they are indeed different.

But again, further digging into why the floating point equality check fails won't really give you a solution, so I'm back to...

Doesn't work...I still have the exact same problem as before. Even when I set epsilon to 0.01, a value far too imprecise for my measurements, it still cannot exit when it should.

Well, I'd recommend continuing along this line anyway.

Can you create a standalone code snippet that demonstrates the problem you see? Otherwise it's much more difficult to play along.

[Caglow]

I am looking for a minimum of 9-10 digits of final precision (in other words, after everything's completed calculating, I can round off and the last digits will be accurate). Generally, I try to go for as many as the compiler will allow which works in VC++ and in Cygwin but doesn't appear to work in MinGW.

It's not really possible to do a standalone code snippet here. If it was, then I would post it. I will try to gather up a portion of the code and attach it as an attachment here. However, the entire project (75,000 lines of code and terms) is so entwined that it's nearly impossible to take a chunk out and have it still work. One thing relies on another which relies on another which end up relying on 20% of the whole thing. Actually, for this, I'm going to estimate about 3000 lines but...

I am printing 6 digits beyond the decimal.

Ignore what I was saying about that other loop working...it doesn't. Apparently, it never built correctly or something because now, it does the infinite loop as well. So if you wanted the standalone code snippet, I guess the one in the first post would work.

I have no idea how to examine the bytes of a number...I generally work at high level.

[Dave Sinkula]

It's not really possible to do a standalone code snippet here. If it was, then I would post it. I will try to gather up a portion of the code and attach it as an attachment here. However, the entire project (75,000 lines of code and terms) is so entwined that it's nearly impossible to take a chunk out and have it still work. One thing relies on another which relies on another which end up relying on 20% of the whole thing. Actually, for this, I'm going to estimate about 3000 lines but...

I thought you had isolated a few cases in which the values shown in a debugger might have been able to be plucked and put in a little canned function. If that's not that case, then I guess that option is out.

I am printing 6 digits beyond the decimal. [...] I have no idea how to examine the bytes of a number...I generally work at high level.

I'd venture that you need to go well beyond 6 digits to see.
For illustration only:

nclude <stdio.h>
#include <string.h>
#include <limits.h>

char *bits_block(char *dest, const void *object, size_t size);

int main(void)
{
   char result[80];
   double e, d = 123.456;
   printf("d: %17.15f = %10.6f = \"%s\"\n", d, d, bits_block(result, &d, sizeof d));

   memcpy(&e, &d, sizeof e);
   *(unsigned char*)&e &= ~1;
   printf("e: %17.15f = %10.6f = \"%s\"\n", e, e, bits_block(result, &e, sizeof e));
   printf("d and e %s equal\n", d == e ? "are" : "are not");

   *(unsigned char*)&e |= 1;
   printf("e: %17.15f = %g = \"%s\"\n", e, e, bits_block(result, &e, sizeof e));
   printf("d and e %s equal\n", d == e ? "are" : "are not");
   return 0;
}

char *bits_uchar(char *dest, unsigned char d)
{
   char *start = dest;
   unsigned char bit;
   for ( bit = 1 << (CHAR_BIT - 1); bit > 0; bit >>= 1 )
   {
      *dest++ = d & bit ? '1' : '0';
   }
   *dest = 0;
   return start;
}

char *bits_block(char *dest, const void *object, size_t size)
{
   char *start = dest;
   const unsigned char *byte;
   for ( byte = object; size-- > 0; ++byte )
   {
      bits_uchar(dest, *byte);
      dest += CHAR_BIT;
      *dest++ = size > 0 ? '-' : 0;
   }
   return start;
}

/* my output
d: 123.456000000000000 = 123.456000 = "01110111-10111110-10011111-00011010-00101111-11011101-01011110-01000000"
e: 123.455999999999990 = 123.456000 = "01110110-10111110-10011111-00011010-00101111-11011101-01011110-01000000"
d and e are not equal
e: 123.456000000000000 = 123.456 = "01110111-10111110-10011111-00011010-00101111-11011101-01011110-01000000"
d and e are equal
*/

Again, the basic nugget to take from this is "don't check floating point values for equality".

[Caglow]

I went out to 20 digits (I think double only allowed 19 though) and the two numbers are still identical in every aspect.