Testing primality of long doubles

[psiblade]

I've been tearing my hair out over this for awhile, and I have a feeling I'm going to kick myself when someone shows me what I'm doing wrong. I've written a very basic boolian function to test input for primality (full function in the attached header). It works perfectly for all integer values I can throw at it, but when I start using floating-point numbers (which I've tried to design it to handle), I start running into weird problems.

For example, the number 9.9999999978e10, which is clearly divisible by 2, should be caught here:

  if (static_cast<int>( fmodl(n, 2) ) == 0)             // tests for divisibility by 2
   {return false;}
   else

The variable n in this case being a long double, the input to the prime() function.

As far as I know, typecasting the result of the modulus will drop anything after the decimal point, (and indeed, a simple test program running static_cast<int>(fmodl(9.9999999978e10,2)) will return 0), eliminating the weirdness that I've heard comes with comparisons with floating-point numbers.

Oddly, prime() is always returning true with non-integer numbers, and if the problem is with this comparison or somewhere else in the function, for the life of me I can't figure out why.

Any insight is much appreciated.

~Jap

[TheDominis]

Last time I checked int ranges from -2147483648 (-2^31) to 2147483647 (2^31 - 1).

If D is 9.9999999978e10 a static cast to an int would yield -2147483648

#include <iostream>
#include <cmath>

bool prime(long double n)
{
	if(n < 0.0)
	{
		n = fabs(n);
	}

	if(n == 0.0 || n == 1.0)
	{
		return false;
	}

	if(n < 4.0)
	{
		return true;
	}

	if(fmodl(n, 2.0) == 0.0)
	{
		return false;
	}

	if(fmodl(n, 3.0) == 0.0)
	{
		return false;
	}

	if(n < 9.0)
	{
		return true;
	}

	for (long double i = 5.0; i <= sqrt(n); i += 6.0)
	{
		if(fmodl(n, i) == 0.0)	
		{
			return false;
		}

		if(fmodl(n, i + 2.0))
		{
			return false;
		}
	}
}


int main(int argc, char* argv[])
{
	if(prime(9.9999999978e10))
	{
		std::cout << "prime\n";
	}
	else
	{
		std::cout << "not prime\n";
	}

	std::cin.get();

	return 0;
}

[psiblade]

*Facepalms* fabsl()... Of course static casts to big numbers are going to show up as negative. Christ, I wouldn't have even looked that direction.

Part of the problem is this is a header from my high-school days, and for some reason weird axioms like "Don't use == to compare floating point values" stuck with me, but for the life of me I can't remember the reasoning behind it. I guess I was lead to believe that for some reason floating point numbers could be a tiny bit inaccurate in their final digits, making 1.0 into 1.00000...000001 for instance, thus producing weirdness when you tried to compare a float to an int. Now that I think about it I don't know why this would happen, and I don't think it was ever really explained to me in the first place. The purpose of all the static casts was to drop everything after the decimal point so that this wouldn't happen.

At any rate, thanks, it works properly now that all the casting garbage has been done away with

[monarch_dodra]

If all you are doing in your "isPrime" function is cast your double to int, why even write a function that accepts double?

If your function accepts ints (or even better, unsigned ints), then writing something like

bool isPrime(unsigned int i);

int main()
{
    double n = some_complex_value;
    return isPrime(n); //implicit cast
}

This will work just fine, and will make things less confusing for your users.

Of course, if a double is too big for the int, you will have a problem, but there is no way of solving it anyways. Quite frankly (IMO), testing primality on a double is just plain wrong to begin with, especially, if it can't be accurately represented.

If you want to check big numbers, consider using unsigned long long int, which ranges to about 4e18. If you need any bigger, find yourself a bigInt/BigNum object.

Last time I checked int ranges from -2147483648 (-2^31) to 2147483647 (2^31 - 1).

If D is 9.9999999978e10 a static cast to an int would yield -2147483648

I'm pretty sure the actal value will be between -2147483648 and 2147483647, but guessing the actual value is probably close to impossible.

[Paul McKenzie]

I've been tearing my hair out over this for awhile, and I have a feeling I'm going to kick myself when someone shows me what I'm doing wrong.

What you're doing wrong is assuming that doubles can be exactly represented in binary.

Please read the following:

http://www.parashift.com/c++-faq-lit...html#faq-29.16
http://www.parashift.com/c++-faq-lit...html#faq-29.17

Regards,

Paul McKenzie

[TheDominis]

If all you are doing in your "isPrime" function is cast your double to int, why even write a function that accepts double?

If your function accepts ints (or even better, unsigned ints), then writing something like

bool isPrime(unsigned int i);

int main()
{
    double n = some_complex_value;
    return isPrime(n); //implicit cast
}

This will work just fine, and will make things less confusing for your users.

Of course, if a double is too big for the int, you will have a problem, but there is no way of solving it anyways. Quite frankly (IMO), testing primality on a double is just plain wrong to begin with, especially, if it can't be accurately represented.

If you want to check big numbers, consider using unsigned long long int, which ranges to about 4e18. If you need any bigger, find yourself a bigInt/BigNum object.


I'm pretty sure the actal value will be between -2147483648 and 2147483647, but guessing the actual value is probably close to impossible.

Cast from double to int on my computer yielded -2147483648.

[monarch_dodra]

Cast from double to int on my computer yielded -2147483648.

int main()
{
    double n = 9.9999999978e10;
    signed int a = static_cast<int>(9.9999999978e10);
    unsigned int b = static_cast<int>(9.9999999978e10);
    //29 D:\*\main.cpp [Warning] integer overflow in expression
    
    std::cout << "double n = " << n << std::endl;
    std::cout << "signed int a = " << a << std::endl;
    std::cout << "unsigned int b = " << b << std::endl;

    system("PAUSE");
    return EXIT_SUCCESS;
}

double n = 1e+011
signed int a = 2147483647
unsigned int b = 2147483647
Press any key to continue . . .

Mine yielded 2147483647, signed or unsigned.
However, I'm pretty sure in my case the cast is done during compilation. I can try to do a run-time cast, to see what happens, but one thing for sure is that we have no guarantees.

[TheDominis]

int main()
{
    double n = 9.9999999978e10;
    signed int a = static_cast<int>(9.9999999978e10);
    unsigned int b = static_cast<int>(9.9999999978e10);
    //29 D:\*\main.cpp [Warning] integer overflow in expression
    
    std::cout << "double n = " << n << std::endl;
    std::cout << "signed int a = " << a << std::endl;
    std::cout << "unsigned int b = " << b << std::endl;

    system("PAUSE");
    return EXIT_SUCCESS;
}

double n = 1e+011
signed int a = 2147483647
unsigned int b = 2147483647
Press any key to continue . . .

Mine yielded 2147483647, signed or unsigned.
However, I'm pretty sure in my case the cast is done during compilation. I can try to do a run-time cast, to see what happens, but one thing for sure is that we have no guarantees.

I didn't want to say with any certainty that a cast from 9.9999999978e10 to an int would always yield -2147483648.

[jefranki]

It's not generally a good idea to test floating point numbers for primality because of round off error. To illustrate:

In a 32-bit floating point number, the maximum value you can exactly convert to an (32-bit) int is 2^24 - 1 ==16777215. Why? Using the IEEE 754 Standard (Wikipedia article) for storing floating points, there are 23 bits allocated for the significand (significant bits), with an implied extra leading 1, so 24 significant bits total. Converting floating point numbers higher than 1.6777215e7 to ints will result in round off error. To contrast, the max value of a signed 32-bit int is 2147483647 (~2.1e9)

The same principle applies to 64-bit and 128-bit values, but the limits are higher.

[Lindley]

You can use boost::numeric_cast if you need to way to easily catch unexpected behavior resulting from overflow.

[psiblade]

What you're doing wrong is assuming that doubles can be exactly represented in binary.

Please read the following:

http://www.parashift.com/c++-faq-lit...html#faq-29.16
http://www.parashift.com/c++-faq-lit...html#faq-29.17

Thanks, this is exactly what I needed to read. To clarify my reasoning a bit, though, I was only casting the entire number in the case of testing negativity (which I see is a mistake, and fabsl() is a much simpler choice), and in the case of testing for single-digit numbers. Sure, a double may cast to 2147483647, but when testing to see if the input is 3, for instance, it obviously isn't 3.

In the case of the loop, where accuracy is important, I was not trying to cast n itself as an int, but only the result of the modulus, the idea being that if the mod came out as 0.0.....001, anything after the decimal point would simply be dropped before being compared to an int value. Now, I realize that results of mods are not always integers, which is why I've added a few lines to my latest revision to return false for any decimal input before we ever get to this point.

I'll tell you how it goes once I can put it through some decent testing.

~Jap

EDIT: I've also taken out the "Test negative numbers as positive" bit entirely, because for the uses this function will receive all negatives should be considered non-prime anyway.

[monarch_dodra]

In the case of the loop, where accuracy is important, I was not trying to cast n itself as an int, but only the result of the modulus, the idea being that if the mod came out as 0.0.....001, anything after the decimal point would simply be dropped before being compared to an int value. Now, I realize that results of mods are not always integers, which is why I've added a few lines to my latest revision to return false for any decimal input before we ever get to this point.

It doesn't matter how you handle your floating point number. If at any one point, your double is not an exact representation of a whole number, the whole concept of primality is flawed.

If you want something safe and/or simple, you should do these tests:

bool isPrimeDouble(double i)
{
    static const double maxDouble = pow(2,std::numeric_limits<double>::digits);
    static const unsigned int maxInt = std::numeric_limits<unsigned int>::max();
    
    if ( i < 0 )
    {
        //i is negative
        std::cout << "negative" << std::endl;
        return false;
    }
    else if ( i > maxDouble )
    {
        //i too big for exact representation
        std::cout << "inexact" << std::endl;
        return false;
    }
    else if ( i > maxInt )
    {
        //i too big for integer cast
        std::cout << "too big" << std::endl;
        return false;
    }
    else if ( (i - static_cast<int>(i)) != 0 )
    {
        //i is not whole
        std::cout << "not whole" << std::endl;
        return false;
    }
    else
    {
        return isPrime(static_cast<int>(i));
    }
}

If the double does not respect any of the above conditions, then there is no point doing a test anyways.

These tests are basically just checking if your double can be an exact cast. But I'm pretty sure boost::numeric_cast can do better.