HiveBrain v1.2.0
Get Started
← Back to all entries
patterncppMinor

Quadratic sieve

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
quadraticsievestackoverflow

Problem

Here is an implementation I pieced together from various papers on the subject. It took a lot of fiddling to get it to run at all.

My ideas for better performance and readability, based on looking at other GMP examples:

  • Always pass by reference and never return any "large type"



  • Separate matrix Gaussian elimination from main



  • Avoid pushback() and all dynamic allocation. Avoid copying whenever possible.



Please tell me if my ideas are in the right direction.

```
#include
#include
#include
#include

// Constants
// The optimal smoothness bound is exp((0.5 + o(1)) sqrt(log(n)log(log(n)))).
const int SMOOTH_BOUND = 500;
const int TRIAL_BOUND = 400;
const int SIEVE_CHUNK = 60;

const bool DEBUG = true;

void *_Unwind_Resume;
void *__gxx_personality_v0;

typedef std::vector int_vector;
typedef std::vector matrix;
typedef std::vector mpz_vector;

template // Takes int_vector or mpz_vector
void print_vector(const T &x)
{
for(size_t i=0; i A (bound, 1);
A[0] = 0; A[1] = 0; // 0 and 1 aren't prime

for(int i=2; i vb_pair;

vb_pair factor_smooth(mpz_class n, const mpz_vector &factor_base)
{
// Each item in factors corresponds to number in factor base
int_vector factors(factor_base.size(), 0);

for(size_t i=0; i SMOOTH_BOUND) // Up to smooth limit
break;
mpz_class p_mpz = p;
// Use only primes that match (n|p) = 1
if (mpz_legendre(n.get_mpz_t(), p_mpz.get_mpz_t()) == 1)
{
factor_base.push_back(p);
}
}

if (DEBUG)
{
std::cout factor_base.size())
{
not_done = false;
break;
}
smooth_x[smooth_count] = current_x[i];
smooth_numbers[smooth_count] = current_chunk[i];
smooth_factors[smooth_count] = factored.first;
smooth_count++;
}
}
}

if (DEBUG)
{
std::cout =0; k--)
{
for

Solution

-
One thing you can do is to remove the .size() and sqrt() calls in loops. You should avoid calculating redundant square roots because they're too heavy.

Instead of:

for(size_t i=0; i<null_space.size(); i++)

for(int i=2; i<sqrt(bound); i++)


Do:

for(size_t i=0, nSize = null_space.size(); i<nSize; i++)

for(int i=2, boundsqrt=sqrt(bound); i<boundsqrt; i++)


-
It's often better if you don't use push_back(), even if you have to loop twice. Count the items, call resize(), then loop again for storing the data. If you care about performance but you still want to continue using std::vector, you should be careful on not resizing the vector all the time inside a loop. Looping twice can be faster than the way you're currently doing because you avoid reallocating the vector.

-
Smaller functions improve both speed and readability. You can move more code into functions, like:

void createFactorBase(mpz_vector& factor_base, mpz_vector& primes, const mpz_class& n)
{
    mpz_class two = 2;
    factor_base.push_back(two);
    for(size_t i=0; i SMOOTH_BOUND) // Up to smooth limit
            break;
        mpz_class p_mpz = p;
        // Use only primes that match (n|p) = 1
        if (mpz_legendre(n.get_mpz_t(), p_mpz.get_mpz_t()) == 1)
        {
            factor_base.push_back(p); // without push_back() tho
        }
    }
}


[Update] You might be interested in watching some stuff like this: CppCon 2014: Mike Acton "Data-Oriented Design and C++"

Code Snippets

for(size_t i=0; i<null_space.size(); i++)

for(int i=2; i<sqrt(bound); i++)
for(size_t i=0, nSize = null_space.size(); i<nSize; i++)

for(int i=2, boundsqrt=sqrt(bound); i<boundsqrt; i++)
void createFactorBase(mpz_vector& factor_base, mpz_vector& primes, const mpz_class& n)
{
    mpz_class two = 2;
    factor_base.push_back(two);
    for(size_t i=0; i<primes.size(); i++)
    {
        int p = primes[i];
        if (p > SMOOTH_BOUND) // Up to smooth limit
            break;
        mpz_class p_mpz = p;
        // Use only primes that match (n|p) = 1
        if (mpz_legendre(n.get_mpz_t(), p_mpz.get_mpz_t()) == 1)
        {
            factor_base.push_back(p); // without push_back() tho
        }
    }
}

Context

StackExchange Code Review Q#111984, answer score: 4

Revisions (0)

No revisions yet.