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 <gmpxx.h>
#include <vector>
#include <cmath>
#include <iostream>

// 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> int_vector;
typedef std::vector<int_vector> matrix;
typedef std::vector<mpz_class> mpz_vector;

template <typename T> // Takes int_vector or mpz_vector
void print_vector(const T &x)
{
for(size_t i=0; i<x.size(); i++)
std::cout << x[i] << ", ";

std::cout << '\n';
}

// Sloppy coding
// Return a list of primes
int_vector eratosthenes(int bound)
{
int_vector primes;

std::vector<bool> A (bound, 1);
A[0] = 0; A[1] = 0; // 0 and 1 aren't prime

for(int i=2; i<sqrt(bound); i++)
{
if (A[i])
{
for(int j = i*i; j<=bound; j+=i)
A[j] = 0;
}
}

for(int i=0; i<bound; i++)
{
if (A[i])
primes.push_back(i);
}
return primes;
}

// Return a vector of a number's factors (ex. [0, 1, 2, 0]) and a boolean of
// whether it's smooth or not
typedef std::pair<int_vector, bool> 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<factor_base.size(); i++)
{
mpz_class factor = factor_base[i];
while (n % factor == 0)
{
n /= factor;
factors[i] ^= 1; // + 1 (mod 2) matrices
}
}
bool is_smooth = (n==1);
vb_pair return_pair(factors, is_smooth);
return return_pair;
}

// ------------------------------------------------

int main()
{
// Test numbers: 502560280658509, 90283
const mpz_class n("502560280658509");

int_vector primes = eratosthenes(TRIAL_BOUND);
mpz_vector factor_base;

// Create factor base
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);
}
}

if (DEBUG)
{
std::cout << "Factor base: ";
print_vector(factor_base);
}

// Find smooth numbers with x = sqrt(n) + j
mpz_class j = 1;
mpz_class sqrt_n = sqrt(n);

// Allocate factor_base+1 size
mpz_vector smooth_numbers(factor_base.size()+1, 0);
mpz_vector smooth_x(factor_base.size()+1, 0);

// Corresponds to smooth numbers
matrix smooth_factors(factor_base.size()+1, {0});
int smooth_count = 0;
bool not_done = true;

while (not_done) // pi(B) + 1
{
mpz_vector current_chunk(SIEVE_CHUNK);
mpz_vector current_x(SIEVE_CHUNK);
for(int i=0; i<SIEVE_CHUNK; i++)
{
mpz_class current;
mpz_class x = sqrt_n + j + i; // Current addition to x
// current = (j+i)^2 mod n
mpz_powm_ui(current.get_mpz_t(), x.get_mpz_t(), 2, n.get_mpz_t());

current_chunk[i] = current;
current_x[i] = x;
}
j += SIEVE_CHUNK;

// Actual factoring
for(size_t i=0; i<current_chunk.size(); i++)
{
vb_pair factored = factor_smooth(current_chunk[i], factor_base);
if (factored.second) // Is smooth
{
if (smooth_count > 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 << "Smooth x: ";
print_vector(smooth_x);
std::cout << "Smooth numbers: ";
print_vector(smooth_numbers);

std::cout << "Smooth factors:\n";
for(size_t i=0; i<smooth_factors.size(); i++)
print_vector(smooth_factors[i]);

std::cout << '\n';
}

// Gaussian Elimination -----------------------------------
// Transpose the matrix
int Ai = smooth_factors[0].size(); // row
int Aj = smooth_factors.size(); // column
matrix A(Ai, int_vector(Aj, 0));

for(int i=0; i<Ai; i++)
{
for(int j=0; j<Aj; j++)
{
A[i][j] = smooth_factors[j][i];
}
}

if (DEBUG)
{
std::cout << "Transposed matrix A:\n";
for(size_t i=0; i<A.size(); i++)
print_vector(A[i]);

std::cout << '\n';
}

for(int k=0; k<Ai; k++)
{
// Swap with pivot if current diagonal is 0
if (A[k][k] == 0)
{
for(int l=k; l<Ai; l++)
{
if (A[l][k]==1)
{
A[l].swap(A[k]);
break;
}
}
}
// For rows below pivot
for(int i=k+1; i<Ai; i++)
{
// If row can be subtracted, subtract every element (using xor)
if (A[i][k])
{
for(int j=0; j<Aj; j++)
A[i][j] ^= A[k][j];
//for(size_t i=0; i<A.size(); i++)
//   print_vector(A[i]);
//std::cout << '\n';
}
}
}

// Find line between free and pivot variables
int f;
for(f=0; f<Aj; f++)
{
if (A[f][f] != 1)
break;
}

// Back substitution on upper triangular matrix
for(int k=f-1; k>=0; k--)
{
for(int i=k-1; i>=0; i--)
{
if (A[i][k])
{
for(int j=0; j<Aj; j++)
A[i][j] ^= A[k][j];
}
}
}

if (DEBUG)
{
std::cout << "Fully reduced matrix:\n";
for(size_t i=0; i<A.size(); i++)
print_vector(A[i]);

std::cout << '\n';
}

int_vector null_space(Aj, 0);
// Subject to change
// First free variable is 1, rest are 0
null_space[f] = 1;

for(int i=0; i<f; i++)
null_space[i] = A[i][f];

mpz_class x_square = 1;
mpz_class y = 1;
for(size_t i=0; i<null_space.size(); i++)
if (null_space[i])
{
x_square *= smooth_numbers[i];
y *= smooth_x[i];
}

if (DEBUG)
{
std::cout << "Null space: ";
print_vector(null_space);
std::cout << "Square: " << x_square << std::endl;

}
mpz_class x, rem;
mpz_sqrtrem(x.get_mpz_t(), rem.get_mpz_t(), x_square.get_mpz_t());

if (DEBUG)
{
if (rem==0)
std::cout << "Remainder 0\n";

std::cout << "x: " << x << '\n' << "y: " << y << "\n\n";
}

mpz_class factor_1;
mpz_class dif = y - x;
mpz_gcd(factor_1.get_mpz_t(), n.get_mpz_t(), dif.get_mpz_t());
if (factor_1 == 1 || factor_1 == n)
std::cout << "Factoring failure: try again with different parameters\n";

mpz_class factor_2 = n / factor_1;

std::cout << "Factor 1: " << factor_1 << '\n';
std::cout << "Factor 2: " << factor_2 << '\n';

return 0;
}


1. Consider allowing your build-system to override the interesting parameters for automating testing and being able to deploy separate release- and debug-versions without manual fiddling with the sources.

#if !SMOOTH_BOUND
#define SMOOTH_BOUND 500
#endif
#if !TRIAL_BOUND
#define TRIAL_BOUND 400
#endif
#if !SIEVE_CHUNK
#define SIEVE_CHUNK 60
#endif


For DEBUG, you might also want to lean on NDEBUG like assert from the standard library.

#ifndef DEBUG
#ifdef NDEBUG
#define DEBUG false
#else
#define DEBUG true
#endif
#endif

2. This code made me wonder:

void *_Unwind_Resume;
void *__gxx_personality_v0;


Why doesn't your standard library provide them, if needed? Did you mess up linking, perhaps by using gcc instead of g++ for that?
See What is __gxx_personality_v0 for?

3. I'm not a fan of hiding perfectly good types behind a typedef. As long as they aren't too long and unwieldy, or the typedef adds semantic information.
Neither seems to be the case here...

4. Do you know for-range-loops? No need to fiddle with indices, or the often superior iterators:

for(const auto& v : x)
std::cout << v << ", ";

5. There's no reason not to generalize print_vector further:

Let it accept any range of printable values, and allow printing to a different output-stream.

6. You might want to test whether using a std::vector<char> and thus avoiding the specialization for std::vector<bool> gives you a performance boost.
And you probably want to cache the square-root, instead of relying on the compiler to do it for you.

Also, you should be more consistent in whether bound is the last or just after the last element you want to consider.
Also, there's no reason you rely on sqrt being completely accurate, just check by squaring the successor. See C++ sqrt function precision for full squares

7. In factor_smooth you are making a useless copy of the vector you return, not even a gratuitious move. Don't do that.

You are also making useless copies of the elements of factor_base.
Also, use in-place construction:

std::pair<std::vector<int>, bool>
factor_smooth(mpz_class n, const mpz_vector &factor_base) {
std::pair<std::vector<int>, bool> r(
std::piecewise_construct,
std::make_tuple(factor_base.size()), // Each item corresponds to number in factor base
std::make_tuple());
auto&& factors = r.first;

for(size_t i=0; i<factor_base.size(); i++) {
auto&& factor = factor_base[i];
while (n % factor == 0) {
n /= factor;
factors[i] ^= 1; // + 1 (mod 2) matrices
}
}
r.second = (n==1);
return r;
}

8. <cmath> is not guaranteed to introduce its names into the global namespace, though it may. Either use <math.h> or preferably qualify with the namespace.

Style:

1. Try being a bit more consistent in your indentation. For whatever reason, some blocks are as a whole indented some amount extra. Perhaps you removed some nesting but didn't correct the indentation?

2. Consider consistently using a single space on each side of (most) binary operators.
Also, decide: Either there's a space between if/while/... and their parentheses, or there isn't, but be consistent.

3. main is the proverbial megamoth: It should be broken down into lots of smaller functions.

Regarding avoidance of returning big types: As long as return-value-optimization applies, that's efficient, don't worry about it.
Also, moving types, as long as most of the bulk isn't in the type itself is efficient.

On coliru, incorporating most of the advice: http://coliru.stacked-crooked.com/a/6a3bb97b48e67c2e

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

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


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

• I've done some research, and push_back is not too bad for performance, since the vector size is automatically doubled. I do think reserve should be used. What do you think?
– qwr
Nov 27, 2015 at 23:05
• Of course, you can use reserve. It depends mostly on how OCD you want to get on performance issues :). Nov 28, 2015 at 4:51
• Are you sure gcc won't optimize the first suggestion anyway?
– qwr
Nov 29, 2015 at 3:56
• My programs have to run 3d math at 30fps at least, and most of the time I like to run them in debug build with the debugger attached, where a lot of optimizations are disabled and another amount of additional runtime checks are performed. So yeah, gcc may optimize a lot in many cases and it depends on the performance requirements of your program. Dec 11, 2015 at 3:44