0

I wrote this code to find the xth greatest prime numbers:

for (int i = 3 /* 2 has already been added to the list */; i < maxNumber; i += 2) {
  for (int tested = 0; ; tested++) {
    if (primes[tested] == 0) {
      break;
    }
    if (i % (int) primes[tested] == 0) {
      goto loop;
    }
  }

  count++;
  if (count == primesSize) {
    primesSize += 2000;
    primes = (double*) realloc(primes, sizeof(double) * primesSize);
  }

  primes[count - 1] = i;
  printf("Prime number #%d: %d\n", count, i);
  printf("Prime size: %d\n", primesSize);

  loop: /* statement that does nothing */ if (1) {}
}

However it returned a "Floating point exception" when using big numbers (> 8,000).

When happens here:

  • The user chooses a number.
  • maxNumber is set to the square root of the chosen number.
  • Firstly, a double pointer of size 1000 * sizeof(double) is allocated. It is stored in the primes variable.
  • If a number is found to be prime, it is added to the array represented by the pointer.
    • When the 1,000th number is added to the array, the primes pointer is reallocated to store 2,000 more numbers.

When I used gdb to find out the cause of the error, I found that this part was the cause of problem:

for (int tested = 0; ; tested++) {
  if (primes[tested] == 0) {
    break;
  }
  if (i % (int) primes[tested] == 0 /* breaks here */) {
    goto loop;
  }
}

Update: I thought the first if statement would catch that issue, because printf("%f", primes[tested]) prints 0. However, it doesn't and the "break" is not executed.

When the code broke, tested was 1001. I convert primes[tested] to an integer because the modulo arithmetic operation I use requires ints to work. However, when I print primes[tested] from code it shows 0. If I print the value from gdb, I get 6.1501785659964211e-319.

What am I missing? Should I modify my call to realloc to avoid this exception?

Maxime Launois
  • 928
  • 6
  • 19
  • Possible duplicate of [Why does a%b produce SIGFPE when b is zero?](https://stackoverflow.com/questions/1081250/why-does-ab-produce-sigfpe-when-b-is-zero) – 1201ProgramAlarm Jun 24 '19 at 18:16
  • @1201ProgramAlarm I know that `primes[tested]` is zero when the program crashed. But my question is different: the `if` statement *should* catch this occurrence to prevent such FPEs, but it doesn't. I edited my question to the best of my ability. – Maxime Launois Jun 24 '19 at 18:22
  • 3
    Primality is a property of integers. Why are you using an array of `double` to store them? – dbush Jun 24 '19 at 18:26
  • And what is the type of primes[] array here in your code?? – UkFLSUI Jun 24 '19 at 18:28
  • Continuing on @dbush: `double` for your purposes even is less precise than `uint64_t`: `double` can only store those values exactly of which all the most significant non-zero bits fit into the mantissa. But there can be huge gaps in between two such values, and all those in between are lost for you anyway, as being rounded. Following IEEE754, you only have 53 bits for mantissa, whereas you can use the full range of 64 bits of the integral value... Even better: you don't need to cast to int (losing even more precision) as `uint64_t` supports `%` natively. – Aconcagua Jun 24 '19 at 18:35
  • Off-topic: You might want to consult your favourite search engine for "sieve of erathostenes" for a more efficient algorithm – at least up to the point at which it consumes too much memory... – Aconcagua Jun 24 '19 at 18:40
  • Why are you using doubles at all? The concept of "prime" only applies to integers. The use of double here only makes the code slower, more error-prone, and using more memory. – Lee Daniel Crocker Jun 24 '19 at 20:14

3 Answers3

3

I thought the first if statement would catch that issue, because printf("%f", primes[tested]) prints 0. However, it doesn't and the "break" is not executed.

You test whether primes[tested] == 0, but your code is only valid if ((int)primes[tested]) == 0. These are not at all the same thing. Moreover, printing the value of primes[tested] with format %f does not reliably tell you differently, because it gives you only 6 digits after the decimal point. Try a "%e" format instead, and test the condition you actually require, not a related, weaker one.

But even better, don't use a floating-point type here. FP has no business being used in a discrete math problem, such as you appear to be trying to solve. If primes[tested] in fact holds prime or possibly-prime numbers then unsigned long long int likely has the same size as double, and almost surely can exactly represent a wider range of primes. Or if it just contains flags, such as in a prime number sieve, then anything wider than char is wasteful.

John Bollinger
  • 160,171
  • 8
  • 81
  • 157
2

Floating point numbers which are really close to zero are still not exactly zero. So your check for equals zero fails.

Note that if you get unlucky enough on the type of machine you compile on even

double f = 1.1;
double x = f;
double y = x;

if( y == f )
    puts("This is not always true!");

Floating point math on computers is tricky and doesn't work as you'd expect from writing math in mathematics where x equals y equals f by definition. No, computer floating point works on bit patterns and they have to be exactly the same.

Anyway, to answer your question. Use the exact same int cast on your if statement as in your modulus and it should work.

And also the new memory returned from realloc will not automatically be set to zero.

And also the third: If you had to cast the return from realloc with (double*) then you're in a C++ compiler and should be using std::vector<double>. It is much better really. Otherwise if you're writing C code then write C code.

Zan Lynx
  • 53,022
  • 10
  • 79
  • 131
  • I don't think your example would ever yield inequality (unless there's NaN involved somewhere): You are just assigning the exact bit patterns again and again. Better examples to be found [here](https://stackoverflow.com/questions/588004/is-floating-point-math-broken)... – Aconcagua Jun 24 '19 at 18:48
  • 1
    @Aconcagua Wouldn't you be surprised then if f and x are assigned to 80-bit registers and y is in a 64-bit memory location. – Zan Lynx Jun 24 '19 at 18:49
  • And then? All values of smaller floating point types can be represented exactly in the larger ones – and you can get them back exactly, too... – Aconcagua Jun 24 '19 at 18:54
  • @Aconcagua The C equality operator does not say 80-bits of 0.1 is equal to 64 bits of 0.1. The 80 bit version has extra bits, which makes it larger. – Zan Lynx Jun 24 '19 at 19:05
  • 1
    @Aconcagua Check out FLT_EVAL_METHOD at https://en.wikipedia.org/wiki/C99#IEEE.C2.A0754_floating_point_support and the note about how GCC defaulted to 2, which is to promote everything to the largest type. – Zan Lynx Jun 24 '19 at 19:06
  • And the 5th bullet point here: https://gcc.gnu.org/onlinedocs/gcc-4.1.2/gcc/Disappointments.html – Zan Lynx Jun 24 '19 at 19:08
  • A double residing in a 80-bit register still is a double, the surplus bits should just be set to 0 or be ignored, shouldn't they? Matter changes quickly in the following example: `float x = 0.1; double y = 0.1; double(x) == y`, then I'm immediately with you... – Aconcagua Jun 24 '19 at 19:22
  • Side note: Any pure C compiler would accept casting the result of malloc or realloc, wouldn't it? Any C compiler not accepting the cast wouldn't be standard compliant, would it? So this can pretty well be C code, applying the cast is just a matter of style (like the surplus parentheses when `return`ing)... – Aconcagua Jun 24 '19 at 19:29
  • @Aconcagua, "*the surplus bits should just be set to 0 or be ignored, shouldn't they?*" Ignored - no, that would not be expected. Set to zero - not guaranteed, and historically not always done. One would *like* to be able to depend on that, but no, it is not safe. This is one of the fine details underscoring why comparing FP values for exact equality is rarely the right thing to do. – John Bollinger Jun 24 '19 at 20:21
  • @JohnBollinger This would allow then an even simpler example: `f = 0.1; if(f == 0.1)` Is it at least guaranteed that the surplus bits are filled the same way every time 64 bits are loaded into 80? Otherwise we cannot even rely on `f == f` (NaN test) any more... – Aconcagua Jun 25 '19 at 05:55
  • @JohnBollinger Standard actually doesn't seem to state much about FP. Closest I find is values of smaller type being a subset of larger types (6.2.4.10). We'd need to *deduce* from that, but that gets pretty shaky quickly... 6.2.6.1.8 is interesting, too, would it apply here? 64 bit double stored in 80-bit register having multiple representations (irrelevant bits for IEEE754 double precision set to arbitrary values). Operators (e. g. multiplication) would need to yield the same result, that might require the FPU being aware of if processing double or long double... – Aconcagua Jun 25 '19 at 06:47
  • No, @Aconcagua, the standard indeed doesn't have much to say about FP details, and that;s half the point. We're not engaged in deducing what must happen, but rather considering what *can* happen, and indeed what has happened in does happen in (some) real-world implementations. – John Bollinger Jun 25 '19 at 11:41
0

I thought the first if statement would catch that issue, because printf("%f", primes[tested]) prints 0.

If this print statement is giving right output according to what you have said, then obviously the Floating point exception makes sense.

You are thinking that if statement should handle the exception, but that really doesn't execute like the way you are thinking. At first it executes the i % (int) primes[tested] part then it compares the result with 0. So, obviously the exception occurs even before the if could function. Hope you understand.

FYI, If b = 0, then a % b causes floating point exception.

And if you have more doubt on the execution steps of the if statement, then run this code yourself:

if (printf("Hello") == 5) {
    printf(" World");
}

Then try to understand which printf() executes first.

UkFLSUI
  • 5,509
  • 6
  • 32
  • 47