14

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.

3
  • 1
    Your for loop seems to start at 3, and increment by 1; your comment states that it starts at 5 and increments by 2.
    – Phil H
    Mar 10, 2009 at 10:19
  • 1
    This question appears to be off-topic because it is about number theory. try math.stackexchange.com.
    – jww
    Feb 17, 2014 at 2:09
  • squared is not the right variable name for the result of a square root. "Squared" means raised to the second power; square root is raising to the 1/2 power. Maybe call it sqrt_test or something. Jul 28, 2014 at 5:09

16 Answers 16

10

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

10

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.

8

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;
    }
}
1
  • 6
    This technique is called memoization, in case anyone wants to search for it.
    – rvighne
    Sep 3, 2014 at 1:31
5

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

Is the size of an array constrained by the upper limit of int (2147483647)?

3
  • The sieve of Eratosthenes is very fast, but you can only use it if you know where the upperbound of the primaries to test will be.
    – Xn0vv3r
    Mar 10, 2009 at 8:48
  • Not at all, take a look at the code. It expands the range in steps, so it's only limited by the capacty of the data type used to store the primes. In this case a long, so the limit is 9223372036854775807.
    – Guffa
    Mar 10, 2009 at 9:24
  • Not at all: you can create a list of lists of sieve approximations, and take "as much as you need". You'll need lazy evaluation of course, like in a functional language. I have written an implementation of this in C# using yield statements which as far as I know worked fine. Don't have laptop with me, so I'd have to get back to this later to post the actual answer. (If you all want.)
    – peSHIr
    Nov 2, 2009 at 10:52
4

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

3

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

Here is the pseudo code from the book Algorithms by S. Dasgupta, C.H. Papadimitriou, and U.V. Vazirani, where n is the number you are testing for primality.

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.

6
  • You want to do that several times to have any real confidence. Mar 9, 2009 at 21:44
  • Yes absolutely, but it is fast enough that you can do that. I edited the answer to mention Carmichael's numbers. Mar 9, 2009 at 21:48
  • 1
    Soloway-Strassen and Miller-Rabin primality testing are both superior to Fermat's Little Theorem in just about every way; both can be trivially extended to deterministic (not just probabilistic) tests, though runtime isn't optimal. Don't bother with FLT.
    – kquinn
    Mar 10, 2009 at 3:21
  • @kquinn, your technically correct, but I brought up FLT because it is the basis of Miller-Rabin. The book I linked goes on to explain the weakness of FLT and then expand it into Miller-Rabin. I also didn't see the Mark posted about Miller-Rabin until after I posted this answer. Mar 10, 2009 at 14:59
  • 1
    What is n (lowercase) and what is N (uppercase)? Jul 28, 2014 at 5:04
1

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.

1
  • 1
    I believe i < Math.Ceiling (or i <= Math.Floor) is enough. You don't need i <= Math.Ceiling.
    – Hosam Aly
    Mar 10, 2009 at 10:52
1

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;
    }
}
3
  • Same error as the OP - this should be "tryDivide <= threshold" or you miss the square numbers.
    – schnaader
    Mar 9, 2009 at 18:34
  • Point taken, I have adjusted back to my original answer. Mar 9, 2009 at 21:19
  • sorry, the <= might still be necessary if the sqrt happens to return a bit "less" then needed (like, n = 9, sqrt(n) == 2.99999999, floor -> 2, algo thinks it's prime... I was kinda mixed up. sorry
    – mitchnull
    Mar 9, 2009 at 21:25
1

First and foremost, primes start from 2. 2 and 3 are primes. Prime should not be dividable by 2 or 3. The rest of the primes are in the form of 6k-1 and 6k+1. Note that you should check the numbers up to SQRT(input). This approach is very efficient. I hope it helps.

public class Prime {

    public static void main(String[] args) {
        System.out.format("%d is prime: %s.\n", 199, isPrime(199)); // Prime
        System.out.format("%d is prime: %s.\n", 198, isPrime(198)); // Not prime
        System.out.format("%d is prime: %s.\n", 104729, isPrime(104729)); // Prime
        System.out.format("%d is prime: %s.\n", 104727, isPrime(982443529)); // Prime
    }

    /**
     * Tells if a number is prime or not.
     *
     * @param input the input
     * @return If the input is prime or not
     */
    private boolean isPrime(long input) {
    if (input <= 1) return false; // Primes start from 2
    if (input <= 3) return true; // 2 and 3 are primes
    if (input % 2 == 0 || input % 3 == 0) return false; // Not prime if dividable by 2 or 3
    // The rest of the primes are in the shape of 6k-1 and 6k+1
    for (long i = 5; i <= Math.sqrt(input); i += 6) if (input % i == 0 || input % (i + 2) == 0) return false;
    return true;
    }

}
1
  • FWIW, I've converted this to VBA for Excel. If anyone is interested, I've posted it below. Dec 26, 2017 at 9:07
1

Repeat mode operations will run very slowly. Use the eratosthenes grid to get the prime list in order.

/*
The Sieve Algorithm
http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes
*/
numbers = new MyBitArray(limit, true);
for (long i = 2; i < limit; i++)
    if (numbers[i])
        for (long j = i * 2; j < limit; j += i)
            numbers[j] = false;
        }

public class MyBitArray: IDisposable
    {
        byte[] bytes;
        public MyBitArray(long limit, bool defaultValue = false)
        {
            long byteCount = (limit & 7) == 0 ? limit >> 3 : (limit >> 3) + 1;
            this.bytes = new byte[byteCount];
            for(long i = 0; i < byteCount; i++)
            {
                bytes[i] = (defaultValue == true ? (byte)0xFF : (byte)0x00);
            }
            this.limit = limit;
        }

        public MyBitArray(long limit, byte[] bytes)
        {
            this.limit = limit;
            this.bytes = bytes;
        }

        public bool this[long index]
        {
            get
            {
                return getValue(index);
            }
            set
            {
                setValue(index, value);
            }
        }
        
        bool getValue(long index)
        {
            if (index < 8)
            {
                return getBit(bytes[0], (byte)index);
            }

            long byteIndex = (index & 7) == 0 ? ((index >> 3) - 1) : index >> 3;
            byte bitIndex = (byte)(index & 7);
            return getBit(bytes[byteIndex], bitIndex);
        }
        void setValue(long index, bool value)
        {
            if (index < 8)
            {
                bytes[0] = setBit(bytes[0], (byte)index, value);
                return;
            }

            long byteIndex = (index & 7) == 0 ? (index >> 3) - 1 : index >> 3;
            byte bitIndex = (byte)(index & 7);
            
            bytes[byteIndex] = setBit(bytes[byteIndex], bitIndex, value);
        }

        bool getBit(byte byt, byte index)
        {
            return ((byt & (1 << index)) >> index) == 1;
        }

        byte setBit(byte byt, byte index, bool value)
        {
            return (byte)((byt & ~(1 << index)) + (value ? 1 << index : 0));
        }

        public void Dispose()
        {
            GC.Collect(2, GCCollectionMode.Optimized);
        }

        private long limit;
        public long Limit { get { return limit; } }
        public byte[] Bytes { get { return this.bytes; } } 
    }

However, I would suggest you a much better method for prime number testing. For 64 bit numbers, no matter how large the number is, it gives the exact result in milliseconds.

public static bool IsPrime(ulong number)
{
    return number == 2 
        ? true 
        : (BigInterger.ModPow(2, number, number) == 2 
            ? (number & 1 != 0 && BinarySearchInA001567(number) == false) 
            : false)
}

public static bool BinarySearchInA001567(ulong number)
{
    // Is number in list?
    // todo: Binary Search in A001567 (https://oeis.org/A001567) below 2 ^ 64
    // Only 2.35 Gigabytes as a text file http://www.cecm.sfu.ca/Pseudoprimes/index-2-to-64.html
}
0

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.

2
  • ceiling will give you a false positive on a prime squared I think :)
    – leppie
    Mar 9, 2009 at 18:50
  • Not if you take up to the ceiling.
    – dirkgently
    Mar 9, 2009 at 18:59
0
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...

1
  • Nice try. Check my answer for an example on how to make it faster.
    – Hosam Aly
    Mar 10, 2009 at 10:52
0

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.

0

This works very fast for testing primes (vb.net)

Dim rnd As New Random()
Const one = 1UL

    Function IsPrime(ByVal n As ULong) As Boolean
        If n Mod 3 = 0 OrElse n Mod 5 = 0 OrElse n Mod 7 = 0 OrElse n Mod 11 = 0 OrElse n Mod 13 = 0 OrElse n Mod 17 = 0 OrElse n Mod 19 = 0 OrElse n Mod 23 = 0 Then
           return false
        End If

        Dim s = n - one

        While s Mod 2 = 0
            s >>= one
        End While

        For i = 0 To 10 - 1
            Dim a = CULng(rnd.NextDouble * n + 1)
            Dim temp = s
            Dim m = Numerics.BigInteger.ModPow(a, s, n)

            While temp <> n - one AndAlso m <> one AndAlso m <> n - one
                m = (m * m) Mod n
                temp = temp * 2UL
            End While

            If m <> n - one AndAlso temp Mod 2 = 0 Then
                Return False
            End If
        Next i

        Return True
    End Function
0

In case anyone else is interested, here's my conversion of Mohammad's procedure above to VBA. I added a check to exclude 1, 0, and negative numbers as they are all defined as not prime.

I have only tested this in Excel VBA:

Function IsPrime(input_num As Long) As Boolean
    Dim i As Long
    If input_num < 2 Then '1, 0, and negative numbers are all defined as not prime.
        IsPrime = False: Exit Function
    ElseIf input_num = 2 Then
        IsPrime = True: Exit Function '2 is a prime
    ElseIf input_num = 3 Then
        IsPrime = True: Exit Function '3 is a prime.
    ElseIf input_num Mod 2 = 0 Then
        IsPrime = False: Exit Function 'divisible by 2, so not a prime.
    ElseIf input_num Mod 3 = 0 Then
        IsPrime = False: Exit Function 'divisible by 3, so not a prime.
    Else
        'from here on, we only need to check for factors where
        '6k ± 1 = square root of input_num:
        i = 5
        Do While i * i <= input_num
            If input_num Mod i = 0 Then
                IsPrime = False: Exit Function
            ElseIf input_num Mod (i + 2) = 0 Then
                IsPrime = False: Exit Function
            End If
            i = i + 6
        Loop
        IsPrime = True
    End If
End Function
-3

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.

1
  • 1
    Marking this down because it is neither faster nor more accurate. The floor of the square root of an integer is exactly what he wants. IEEE floating-point arithmetic is required to produce the correct result on integers, i.e., sqrt(25) cannot be 4.999999. And it's slower because you've introduced a multiplication into the loop which wasn't there before. And finally, you've also introduced a bug, because idx * idx can overflow and produce negative values, causing an infinite loop. Consider test = 2147483647 and idx = 46340 and 46341. Jul 28, 2014 at 5:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.