cgccfloating-pointtype-conversionpelles-c

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


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).


Solution

  • 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.