Pelles C and GCC give different result with this C primality test

Question

I compile this code with this compilers. For number i write 18446744073709551615 (2^64-1). Pelles's executable says "18446744073709551615 is prime" but GCC's executable says "18446744073709551615 isn't prime". Why results are different?

#include <stdio.h>
#include <math.h>
int main(void)
{
    unsigned long long number;
    printf("number: ");
    scanf("%llu",&number);
    unsigned long trsq=truncl(sqrtl(number));
    char s=1;
    for(unsigned long i=2;i<=trsq;i++) {
        if (number%i==0) {
            s=0;
            break;
        }
    }
    if (s==1) {
        printf("%llu is prime\n",number);
    } else {
        printf("%llu isn't prime\n",number);
    }
    return 0;
}

Edit:

I've tested and gcc give 12, pelles c give 8 for sizeof(long double).

For what it's worth, GCC is correct, since the example number is obviously not prime - it is at least divisible by 5 and some other number larger than 1. Confirmation that it is indeed not prime: http://www.wolframalpha.com/input/?i=is+18446744073709551615+prime — CmdrMoozy, Jan 30 '14 at 22:05
Pelles says: *warning #2215: Conversion from 'unsigned long long int' to 'long double'; possible loss of data.* and *warning #2215: Conversion from 'long double' to 'unsigned long int'; possible loss of data.* so there is your problem. and if you hard code the value: *warning #2072: Overflow in constant '18446744073709551615'.* — this, Jan 30 '14 at 22:09
Yeah, that's right. I was thinking of a bug I had in a program a while back where `int` and `long int` are the same size on Windows, which is not the case here. @self. has the correct answer for you. — CmdrMoozy, Jan 30 '14 at 22:16
unrelated, but at a minimum you should be starting i at 3 and using += 2 instead of ++ for a 2x speed boost! (you'll need to check even numbers separately). But you can go further than that with (for example) sieves or statistical methods. — Dave, Jan 30 '14 at 22:17
anyway you use long long for your number and only long for your iterator and limit, so I'd guess that's where your problem lies (the compilers may be using different lengths for long, causing an overflow, causing no checks to be performed at all) — Dave, Jan 30 '14 at 22:20
@self If you hard code the value, it should be 18446744073709551615ULL. Otherwise, yeah you will get an overflow if it cannot be represented as a signed int. — , Jan 30 '14 at 22:22
@user3255144 if you use a *unsigned long long trsq* instead of *unsigned long* you will get a correct result. Probably still a coincidence. — this, Jan 30 '14 at 22:23
@self. Actually you probably nailed it. Pelles warns of a loss of precision. If 2^64-1 is rounded to 2^64 when converted to `long double`, then its square root is 2^32, which overflows a 32-bit `unsigned long`. This overflow may well cause `trsq` to have the value 0. If this happens, the loop runs 0 times and the program reports a prime. — Gilles 'SO- stop being evil', Jan 30 '14 at 22:56

score 3 · Accepted Answer · edited May 23 '17 at 12:21

Many things are left unspecified in the definition of the C language, including the size of numeric types and what happens in case of overflow or loss of precision. So it is common to get different results with different implementations (different compilers, different hardware, different operating systems) when you don't check numeric type ranges or when you use floating point.

Given the information provided by self. in two comments, here is a plausible explanation of what's happening. I don't have Pelles C to check.

Pelles warns:

warning #2215: Conversion from 'unsigned long long int' to 'long double'; possible loss of data. warning #2215: Conversion from 'long double' to 'unsigned long int'; possible loss of data

Conjecture #1: the mathematical value 2^64-1 cannot be represented exactly in a long double (this is likely, as 2^64-1 requires 64 bits of mantissa and few implementations have that much). It is rounded to 2^64, which can be represented exactly.

The value of number is 2^64-1, and it is an unsigned long long. Since the function sqrtl expects a long double argument, the value is converted to that type. Given conjecture #1, sqrtl receives a value of 2^64. The result is therefore 2^32. Since this is an integer, truncl returns the same value.

Conjecture #2: unsigned long is a 32-bit type (this is pretty much the norm on 32-bit machines, and is also the case on 64-bit versions of Windows, at least with the Microsoft compiler).

If unsigned long is a 32-bit type, then the value 2^32 overflows it. What happens in case of overflow in a conversion from a floating-point value to an integer value is not defined by the C standard, compilers can choose to do whatever they want.

Conjecture #3: In Pelles C, when a floating-point value is converted to an integer type, it is wrapped modulo the size of the type, like what happens when converting to a smaller integer type.

Under conjecture #3, trying to assign the value 2^32 to trsq, which is of type unsigned long and 32-bit wide, sets it to 0. Thus trsq has the value 0, the for loop runs 0 times, and the program erroneously reports that the number is prime.

An easy fix is to change trsq to be unsigned long long.

Note that your program may report some numbers as prime if their largest prime factor is very close to their square root (e.g. if the number is the square of a prime), because the conversion of number to a floating-point value may round it down, so trsq may end up being less than the square root, even less than the largest integer that is smaller than the square root.

You can avoid all this trouble by performing an integer square root computation.

Pelles C and GCC give different result with this C primality test

1 Answers1