How can I test for primality?

Question

I am writing a little library with some prime number related methods. As I've done the groundwork (aka working methods) and now I'm looking for some optimization. Ofcourse the internet is an excellent place to do so. I've, however, stumbled upon a rounding problem and I was wondering how to solve this.

In the loop I use to test a number for it's primality it's more efficient to search until sqrt(n) instead of n/2 or even n - 1. But due to rounding problems some number get skipped and thus some primes are skipped! For example, the 10000th prime should be: 104729, but the 'optimized' version ends up with: 103811.

Some code (it's open for more optimization, I know, but I can handle only one thing at a time):

/// <summary>
/// Method for testing the primality of a number e.g.: return IsPrime(29);
/// History:
/// 1. Initial version, most basic form of testing: m smaller then n -1
/// 2. Implemented m smaller then sqrt(n), optimization due to prime factoring
/// </summary>
/// <param name="test">Number to be tested on primality</param>
/// <returns>True if the number is prime, false otherwise</returns>
public static bool IsPrime(int test)
{
    // 0 and 1 are not prime numbers
    if (test == 0 || test == 1) return false;

    // 2 and 3 are prime numbers
    if (test == 2) return true;

    // all even numbers, save 2, are not prime
    if (test % 2 == 0) return false;

    double squared = Math.Sqrt(test);
    int flooredAndSquared = Convert.ToInt32(Math.Floor(squared));

    // start with 5, make increments of 2, even numbers do not need to be tested
    for (int idx = 3; idx < flooredAndSquared; idx++)
    {
        if (test % idx == 0)
        {
            return false;
        }
    }
    return true;
}

I know the squared part fails me (or I fail), tried Math.Ceiling as well, with about the same results.

Answer 1

Try the sieve of eratosthenes -- that should take out getting the root and floating point issues.

As for the floor, you will be better served by taking the ceiling.

Answer 2

Here is a halfway decent function I wrote to solve one of the Euler:

private static long IsPrime(long input)
{
    if ((input % 2) == 0)
    {
        return 2;
    }
    else if ((input == 1))
    {
        return 1;
    }
    else
    {
        long threshold = (Convert.ToInt64(Math.Sqrt(input)));
        long tryDivide = 3;
        while (tryDivide < threshold)
        {
            if ((input % tryDivide) == 0)
            {
                Console.WriteLine("Found a factor: " + tryDivide);
                return tryDivide;
            }
            tryDivide += 2;
        }
        Console.WriteLine("Found a factor: " + input);
        return -1;
    }
}

Answer 3

I guess this is your problem:

for (int idx = 3; idx < flooredAndSquared; idx++)

This should be

for (int idx = 3; idx <= flooredAndSquared; idx++)

so you don't get square numbers as primes. Also, you can use "idx += 2" instead of "idx++" because you only have to test odd numbers (as you wrote in the comment directly above...).

Answer 4

Try this...

if (testVal == 2) return true;
if (testVal % 2 == 0) return false;

for (int i = 3; i <= Math.Ceiling(Math.Sqrt(testVal)); i += 2)
{
   if (testVal % i == 0)
       return false;
}

return true;

Ive used this quite a few times.. Not as fast as a sieve.. but it works.

Answer 5

I'll make a sieve as well then. And schnaader, my comments weren't 100% in sync as you noticed :)

Thanks for the replies :)

Answer 6

I posted a class that uses the sieve or Eratosthenes to calculate prime numbers here:

http://stackoverflow.com/questions/573692/is-the-size-of-an-array-constrained-by-the-upper-limit-of-int-2147483647/573770#573770

Answer 7

I don't know if this is quite what you are looking for but if you are really concerned about speed then you should look into probablistic methods for testing primality rather than using a sieve. Rabin-Miller is a probabilistic primality test used by Mathematica.

Answer 8

Your for loop should look like this:

for (int idx = 3; idx * idx <= test; idx++) { ... }

That way, you avoid floating-point computation. Should run faster and it'll be more accurate. This is why the for statement conditional is simply a boolean expression: it makes things like this possible.

Answer 9

You might want to look into Fermat's little theorem.

Here is the pseudo code from the book Algorithms.

Pick a positive integer a < N at random
if a^n-1 is equivalent to 1 (mod n)
   return yes
else
   return no

Implementing Fermat's theorem should be faster then the sieve solution. However, there are Carmichael numbers that pass Fermat's test and are NOT prime. There are workarounds for this. I recommend consulting Section 1.3 in the fore mentioned book. It is all about primality testing and might be helpful for you.

Answer 10

private static bool IsPrime(int number) {
    if (number <= 3)
        return true;
    if ((number & 1) == 0)
        return false;
    int x = (int)Math.Sqrt(number) + 1;
    for (int i = 3; i < x; i += 2) {
        if ((number % i) == 0)
            return false;
    }
    return true;
}

I can't get it any faster...

Answer 11

As Mark said, the Miller-Rabin test is actually a very good way to go. An additional reference (with pseudo-code) is the Wikipedia article about it.

It should be noted that while it is probabilistic, by testing just a very small number of cases, you can determine whether a number is prime for numbers in the int (and nearly long) range. See this part of that Wikipedia article, or the Primality Proving reference for more details.

I would also recommend reading this article on modular exponentiation, as otherwise you're going to be dealing with very very large numbers when trying to do the Miller-Rabin test...

Answer 12

Sadly, I haven't tried the algorithmic approaches before. But if you want to implement your approach efficiently, I'd suggest doing some caching. Create an array to store all prime numbers less than a defined threshold, fill this array, and search within/using it.

In the following example, finding whether a number is prime is O(1) in the best case (namely, when the number is less than or equal to maxPrime, which is 821,461 for a 64K buffer), and is somewhat optimized for other cases (by checking mod over only 64K numbers out of the first 820,000 -- about 8%).

(Note: Don't take this answer as the "optimal" approach. It's more of an example on how to optimize your implementation.)

public static class PrimeChecker
{
    private const int BufferSize = 64 * 1024; // 64K * sizeof(int) == 256 KB

    private static int[] primes;
    public static int MaxPrime { get; private set; }

    public static bool IsPrime(int value)
    {
        if (value <= MaxPrime)
        {
            return Array.BinarySearch(primes, value) >= 0;
        }
        else
        {
            return IsPrime(value, primes.Length) && IsLargerPrime(value);
        }
    }

    static PrimeChecker()
    {
        primes = new int[BufferSize];
        primes[0] = 2;
        for (int i = 1, x = 3; i < primes.Length; x += 2)
        {
            if (IsPrime(x, i))
                primes[i++] = x;
        }
        MaxPrime = primes[primes.Length - 1];
    }

    private static bool IsPrime(int value, int primesLength)
    {
        for (int i = 0; i < primesLength; ++i)
        {
            if (value % primes[i] == 0)
                return false;
        }
        return true;
    }

    private static bool IsLargerPrime(int value)
    {
        int max = (int)Math.Sqrt(value);
        for (int i = MaxPrime + 2; i <= max; i += 2)
        {
            if (value % i == 0)
                return false;
        }
        return true;
    }
}

Answer 13

I thought Prime numbers and primality testing was useful and the AKS algorithm sounds interesting even if it isn't particularly practical compared to a probabily based tests.