How does this primality test make the code 5000x faster?
Warning: This code is a solution for Project Euler problem #50. If you
don't want it spoiled, don't look here.
Here I have code that searches for a long sequence of consecutive prime
numbers, which summed together are also a prime. At one point, I need to
test whether a sum is prime.
I have two tests, which are ifdef'd in the function computeMaxPrime. The
first checks the sum against a std::set of prime numbers. The second uses
a Miller-Rabin test implemented by GMP. The function only gets called 6
times. When I use the first test, the function computeMaxPrime takes .12
seconds. When I use the second test, it only takes ~.00002 seconds. Can
someone explain how that's possible? I wouldn't think 6 calls to check
whether a number is in a set would take 100 ms. I also tried using an
unordered_set, and it performs the same.
I thought that maybe it was a timing issue, but I've verified it via
timing the whole program execution from Terminal (on OSX). I've also
verified that if I change the test to use the Miller-Rabin test first, and
then confirm using the set, it makes a single call to the set and the
clock reports .02 seconds, exactly what I would expect (1/6th the total
time of only using the set test).
#include "PrimeGenerator2.h"
#include <set>
#include <stdio.h>
#include <time.h>
#include <gmp.h>
typedef std::set<u_int64t> intSet;
bool
isInIntSet (
intSet set,
u_int64t key)
{
return (set.count(key) > 0);
}
bool
isPrime (
u_int64t key)
{
mpz_t integ;
mpz_init (integ);
mpz_set_ui (integ, key);
return (mpz_probab_prime_p (integ, 25) > 0);
}
void
computeInitialData (
const u_int64t limit,
intSet *primeSet,
intList *sumList,
u_int64t *maxCountUpperBound)
{
PrimeSieve sieve;
u_int64t cumSum = 0;
u_int64t pastUpperBound = 0;
*maxCountUpperBound = 0;
for (u_int64t prime = sieve.NextPrime(); prime < limit; prime =
sieve.NextPrime()) {
primeSet->insert(prime);
cumSum += prime;
sumList->push_back(cumSum);
if (cumSum < limit)
(*maxCountUpperBound)++;
else
pastUpperBound++;
}
}
u_int64t
computeMaxPrime (
const u_int64t limit,
const intSet &primeSet,
const intList &sumList,
const u_int64t maxCountUpperBound)
{
for (int maxCount = maxCountUpperBound; ; maxCount--) {
for (int i = 0; i + maxCount < sumList.size(); i++) {
u_int64t sum;
sum = sumList[maxCount + i] - sumList[i];
if (sum > limit)
break;
#if 0
if (isInIntSet (primeSet, sum))
return sum;
#else
if (isPrime (sum))
return sum;
No comments:
Post a Comment