Efficient prime factorization for large numbers

I've been working on a little problem where I need to compute 18-digit numbers into their respective prime factorization. Everything compiles and it runs just fine, considering that it actually works, but I am looking to reduce the run time of the prime factorization. I have implemented recursion and threading but I think I might need some help in understanding possible algorithms for large number computation.

Every time I run this on the 4 numbers I have pre-made, it takes about 10 seconds. I would like to reduce this to possibly 0.06 seconds if there are any ideas out there.

I noticed a few algorithms like Sieve of Eratosthenes and producing a list of all the prime numbers prior to computing. I'm just wondering if someone could elaborate on it. For instance, I'm having issues understanding how to implement Sieve of Eratosthenes into my program or if it would even be a good idea. Any and all pointers on how to approach this better would be really helpful!

#include <iostream>
#include <vector>
#include <chrono>

using namespace std;
using namespace std::chrono;

vector<long long> inputVector;
bool developer = false;
vector<long long> primeVector;

{
long long initValue;        // the number being prime factored
vector<long long> factors;  // all of the factor values
public:
void setInitValue(long long n)
{
initValue = n;
}
{
factors.push_back(m);
}
void setVector(vector<long long> m)
{
factors = m;
}
long long getInitValue()
{
return initValue;
}
vector<long long> getVector()
{
return factors;
}
};

// find primes recursively and have them returned in vectors
vector<long long> getPrimes(long long n, vector<long long> vec)
{
double sqrt_of_n = sqrt(n);

for (int i = 2; i <= sqrt_of_n; i++)
{
if (n % i == 0)
{
return vec.push_back(i), getPrimes(n / i, vec); //cause recursion
}
}

// pick up the last prime factorization number
vec.push_back(n);

//return the finished vector
return vec;
}

void getUserInput()
{
long long input = -1;
cout << "Enter all of the numbers to find their prime factors. Enter 0 to compute" << endl;
do
{
cin >> input;
if (input == 0)
{
break;
}
inputVector.push_back(input);
} while (input != 0);
}

int main()
{
vector<long long> temp1;   // empty vector
vector<long long> result1; // temp vector

if (developer == false)
{
getUserInput();
}
else
{
cout << "developer mode active" << endl;
long long a1 = 771895004973090566;
long long b1 = 788380500764597944;
long long a2 = 100020000004324000;
long long b2 = 200023423420000000;
inputVector.push_back(a1);
inputVector.push_back(b2);
inputVector.push_back(b1);
inputVector.push_back(a2);
}

high_resolution_clock::time_point time1 = high_resolution_clock::now();

// give each thread a number to comput within the recursive function
for (int i = 0; i < inputVector.size(); i++)
{
prime.setInitValue(inputVector.at(i));
prime.setVector(result1 = getPrimes(inputVector.at(i), temp1));
primes.push_back(prime);
}));
}

// allow all of the threads to join back together.
{
th.join();
}

high_resolution_clock::time_point time2 = high_resolution_clock::now();

// print all of the information
for (int i = 0; i < primes.size(); i++)
{
vector<long long> temp = primes.at(i).getVector();

for (int m = 0; m < temp.size(); m++)
{
cout << temp.at(m) << " ";
}
cout << endl;
}

cout << endl;

// so the running time
auto duration = duration_cast<microseconds>(time2 - time1).count();

cout << "Duration: " << (duration / 1000000.0) << endl;

return 0;
}

• A simple improvement would be to change your for (int i = 2; i <= sqrt_of_n; i++) loop to for (int i = 3; i <= sqrt_of_n; i = i +2) Oct 13 '14 at 16:07
• wow yeah haha that shaved off about 3 seconds off of the program. Thanks! Oct 13 '14 at 16:08
• It looks like you are doing trial division. Have you looked at more efficient algorithms? Oct 14 '14 at 16:33

Obviously, the best way to speed up your code woud be use another algorithm. However, the other posts already do a great job covering that part. Therefore, I will review the other things in your code that could have been done better:

• It seems that you forgot to include the header <cmath> for std::sqrt. Your implementation may include it from another header, but that's not guaranteed by the standard. My implementation produced a compile-time error because of the lack of <cmath>.

• threads.push_back(thread([&]{ /* ... */ })); is redundant: we already know that threads is a std::vector<std::thread>, so let's write threads.emplace_back([&]{ /* ... */ }); instead.

• You have problems due to paralellism. When I try to run your program, I get an std::out_of_range exception thrown which says that I tried to access the nth element of a std::vector of size n. I did not understand all the details, but it comes from this loop:

for (int i = 0; i < inputVector.size(); i++)
{
prime.setInitValue(inputVector.at(i));
prime.setVector(result1 = getPrimes(inputVector.at(i), temp1));
primes.push_back(prime);
});
}


In this loop, the lambda captures all the variables by reference and that's a bad idea: once the thread is spawn, i might be incremented before inputVector.at(i) is called. To avoid this problem, you should take i by value in the capture instead of taking it by reference:

threads.emplace_back([&, i]{
/* ... */
});


Capturing i by value ensures that it won't be modified by another thread meanwhile. In multithreaded code, it is really important to think about what should be taken by value and what should be taken by reference.

• Moreover, on the same piece of code, you have concurrent stores on primes through the method push_back. While I didn't get any error - it seems that you didn't either -, this may be a data race. Therefore, you should use an std::mutex to lock primes before pushing values to it (that comment is also valid for result1):

std::mutex primes_mutex;
for (int i = 0; i < inputVector.size(); i++)
{
prime.setInitValue(inputVector.at(i));
std::lock_guard<std::mutex> lock(primes_mutex);
prime.setVector(result1 = getPrimes(inputVector.at(i), temp1));
primes.push_back(prime);
});
}


Note that I used a, std::lock_guard to guard the std::mutex. It avoids to forget a potential call to unlock since it invokes the unlocking function when its destructor is called.

• I love to focus on small pieces of code, so there we go: the loop above can still be improved thanks to the C++11 range-based for loop:

std::mutex primes_mutex;
for (long long val: inputVector)
{
prime.setInitValue(val);
std::lock_guard<std::mutex> lock(primes_mutex);
prime.setVector(result1 = getPrimes(val, temp1));
primes.push_back(prime);
});
}


This iteration style allows to avoid the double concurrent lookup of inputVector at the same position. While it was already data-race free (the standard guarantees that concurrent loads do not induce data races), the less variables to read and write, the less there is to think about.

• Also (yeah, still the same loop), I don't understand why you left the initialization of prime out of the thread: there is no point in leaving it out, and it forces the lambda to capture yet another parameter by reference, which would have been avoided had you declared prime directly into the new thread.

using namespace std

Naming

• PrimeNumber

This class doesn't deserve this name, because it doesn't represent a prime number. It represents a big big number of which you want to retrieve the prime factors. This leads to the getVector() method which should be better named getPrimeFactors(). In addition to this, vector<PrimeNumber> primes; should be renamed also.

• getPrimes()

This method does not get primes, but prime factors, so a rename should be done to e.g getPrimeFactors(). As you see, we would now have two times the same method name, so you should consider to do the calculation of the prime factors inside the to be renamed PrimeNumber class.

Comments should be comment on why is something done instead of what is done. See

//return the finished vector
return vec;


Because comments can become misleading fast, use them with care and update them also if your code changes. See

vector<long long> temp1;   // empty vector
vector<long long> result1; // temp vector


In addition, if you need a comment to describe what a variable is for, or what it should represent, then you didn't name your variable good.

Single responsibility principle

Your main() method does a way to much.

• processing input from user
• calling the calculation
• printing the results
• measuring the time

You should extract some of these to separate methods.

What I really like: You are using braces {} for single if statements also.

Algorithm

To speed up your algorithm you should

1. precalculate an amount of prime numbers up to a defined border, by using e.g Sieve of Eratosthenes.
2. Check if these primen umbers are sufficient to get all prime factors.
3. If they are not sufficient, get the next prime number, store it and go back to 2. to check, if the new prime number is a prime factor of your number.

Some ideas:

• use a mod-6 wheel (check div by 2,3; set incr=2, limit=sqrt(n); then loop from 5 to limit with i += incr, incr ^= 6). That skips multiples of 2 and 3 by alternating between incrementing by 2 and by 4.

• The recursion seems awkward. We can just recalculate the limit as we find a factor, and push the result.

E.g. something like (not tested):

vector<unsigned long long> getFactors(unsigned long long n) {
unsigned long long i, inc, limit;
vector<unsigned long long> vec;

if (n <= 1) {
if (n == 1)  vec.push_back(1);
return vec;
}
while (n % 2 == 0) {
n /= 2;
vec.push_back(2);
}
while (n % 3 == 0)
n /= 3;
vec.push_back(3);
}
limit = (unsigned long long) sqrt(n);
for (i = 5, inc = 2;  i <= limit;  i += inc, inc ^= 6) {
if (n % i == 0) {
do {
n /= i;
vec.push_back(i);
} while (n % i == 0);
limit = (unsigned long long) sqrt(n);
}
}
if (n > 1)
vec.push_back(n);
return vec;
}


The problem is that for arbitrary 18 digit numbers this is still going to be slow and may or may not hit your target time (it's about 0.11s each for 1000 numbers on my computer with 10 other computations running). I can speed that up about 3x using primes to 2000 then a mod-30 wheel. Going to Pollard-Rho makes it run thousands of times faster.

If you really need it fast, you're going to want something like Pollard's Rho, Brent's improvements to Rho, P-1, or SQUFOF. Trial division (even efficiently done) is just too slow for general numbers this size. The basic Pollard's Rho is pretty simple -- the only real trick is doing an efficient mulmod. Brent isn't too hard -- more or less just running the gcds in batches with backtracking when needed. For this size, P-1 is probably going to be slower than Rho/Brent. SQUFOF is more complicated to implement. An additional complication is that now you need a primality test to know when you're done (deterministic M-R with at most 7 selected bases, or BPSW). It's possible to skip that by just using Rho etc. opportunistically, but that doesn't fully utilize the idea.

For the point of view of an algorithm, you can generate all primes up to a number N as follows:

• For the next prime 'p', mark powers of p and p multiplied by any subsequent "primes" as non-prime. I put "primes" here in quotes because they may well not be prime.

Thus for 2, you mark 4 and up.

For 3, start with 9 and all the powers and also 3 multiplied by all odd numbers. You "mark" a number with the prime currently being used. So if a number if marked with the current one you use it anyway. For example, you have marked 15 with a "3" but you still use it for now 3*15 = 45.

The numbers you have "marked" then get removed from a "list" so that the next iteration only runs through unmarked numbers.

Once you have done this you have all primes from 2 to N.

(Note: this is Eratosthenes's algorithm, that you suggested you might want to use).

For your huge primes, you have a smaller list of numbers to test against.

For the actual code itself, I do not like the class being called "PrimeNumber" and then having a vector of its "factors". Because a prime number doesn't have factors other than 1 and itself. So it isn't actually a prime number and it's a misnomer.

I would typedef the long long and not put using namespace std in the header file.

Also make the code const-correct.