How to speed up this primality test

Question

I would like to find the largest prime factor of a given number. After several attempts, I've enhanced the test to cope with rather big numbers (i.e. up to one billion in milliseconds). The problem is now if go beyond one billion, the execution time goes forever, so to speak. I wonder if I can do more improvements and reduce the execution time. I'm hoping for better execution time because in this link Prime Factors Calculator, the execution time is incredibly fast. My target number at this moment is 600851475143. The code is rather self-explanatory. Note: I've considered Sieve of Eratosthenes algorithm with no luck regarding the execution time.

#include <iostream>
#include <cmath>

bool isPrime(int n)
{
    if (n==2)
        return true;

    if (n%2==0)
        return false;

    for (int i(3);i<=sqrt(n);i+=2) // ignore even numbers and go up to sqrt(n)
        if (n%i==0)
            return false;

    return true;
}

int main()
{
    int max(0);
    long long target(600851475143);

    if( target%2 == 0 )
        max = 2;

    for ( int i(3); i<target; i+=2 ){ // loop through odd numbers. 
        if( target%i == 0 )  // check for common factor
            if( isPrime(i) ) // check for prime common factor
                max = i;
    }

    std::cout << "The greatest prime common factor is " << max << "\n";


    return 0;
}

Answer 1

One obvious optimization that I can see is:

for (int i(3);i<=sqrt(n);i+=2) // ignore even numbers and go up to sqrt(n)

instead of calculating sqrt everytime you can cache the result in a variable.

auto maxFactor = static_cast<int>sqrt(n);
for (int i(3); i <= maxFactor; i+=2);

The reason I believe this could lead to speed up is sqrt deals with floating point arithematic and compilers usually aren't generous in optimizing floating point arithematic. gcc has a special flag ffast-math to enable floating point optimizations explicitely.

For numbers upto the target range that you mentioned, you will need better algorithms. repeated divisioning should suffice.

Here is the code (http://ideone.com/RoAmHd) which hardly takes any time to finish:

int main() {
    long long input = 600851475143;
    long long mx = 0;
    for (int x = 2; x <= input/x; ++x){
        while(input%x==0) {input/=x; mx = x; }

    }
    if (input > 1){
        mx = input;
    }
    cout << mx << endl;
    return 0;
}

The idea behind repeated division is if a number is already a factor of p, it is also a factor of p^2, p^3, p^4..... So we keep eliminating factors so only prime factors remain that eventually get to divide the number.

Answer 2

You don't need a primality test. Try this algorithm:

function factors(n)
    f := 2
    while f * f <= n
        if n % f == 0
            output f
            n := n / f
        else
            f := f + 1
    output n

You don't need a primality test because the trial factors increase by 1 at each step, so any composite trial factors will have already been handled by their smaller constituent primes.

I'll leave it to you to implement in C++ with appropriate data types. This isn't the fastest way to factor integers, but it is sufficient for Project Euler 3.

Answer 3

Note that except for 2 and 3, all prime numbers are adjacent to multiples of 6.

The following code reduces the total number of iterations by:

• Leveraging the fact mentioned above
• Decreasing target every time a new prime factor is found

---

#include <iostream>


bool CheckFactor(long long& target,long long factor)
{
    if (target%factor == 0)
    {
        do target /= factor;
        while (target%factor == 0);
        return true;
    }
    return false;
}


long long GetMaxFactor(long long target)
{
    long long maxFactor = 1;

    if (CheckFactor(target,2))
        maxFactor = 2;

    if (CheckFactor(target,3))
        maxFactor = 3;

    // Check only factors that are adjacent to multiples of 6
    for (long long factor = 5, add = 2; factor*factor <= target; factor += add, add = 6-add)
    {
        if (CheckFactor(target,factor))
            maxFactor = factor;
    }

    if (target > 1)
        return target;
    return maxFactor;
}


int main()
{
    long long target = 600851475143;
    std::cout << "The greatest prime factor of " << target << " is " << GetMaxFactor(target) << std::endl;
    return 0;
}

Answer 4

for ( int i(3); i<target; i+=2 ){ // loop through odd numbers. 
    if( target%i == 0 )  // check for common factor
        if( isPrime(i) ) // check for prime common factor
            max = i;

It is the first two lines of this code, not primality checks, which take almost all time. You divide target to all numbers from 3 to target-1. This takes about target/2 divisions.

Besides, target is long long, while i is only int. It is possible that the size is too small, and you get an infinite loop.

Finally, this code does not calculate the greatest prime common factor. It calculate the greatest prime divisor of target, and does it very inefficiently. So what do you really need?

And it is a bad idea to call anything "max" in c++, because max is a standard function.

Answer 5

Here is my basic version:

int main() {
    long long input = 600851475143L;

    long long pMax = 0;

    // Deal with prime 2.
    while (input % 2 == 0) {
        input /= 2;
        pMax = 2;
    }

    // Deal with odd primes.
    for (long long x = 3; x * x <= input; x += 2) {
        while (input % x == 0) { 
            input /= x;
            pMax = x;
        }
    }

    // Check for unfactorised input - must be prime.
    if (input > 1) {
        pMax = input;
    }

    std::cout << "The greatest prime common factor is " << pMax << "\n";

    return 0;
}

It might be possible to speed things up further by using a Newton-Raphson integer square root method to set up a (mostly) fixed limit for the loop. If available that would need a rewrite of the main loop.

    long long limit = iSqrt(input)
    for (long long x = 3; x <= limit; x += 2) {
        if (input % x == 0) {
            pMax = x;
            do {
                input /= x;
            } while (input % x == 0);
            limit = iSqrt(input); // Value of input changed so reset limit.
        }
    }

The square root is only calculated when a new factor is found and the value of input has changed.