Description

#include <vector>
#include <iostream>
#include <cmath>
#include <iostream>
#include <boost/multiprecision/cpp_int.hpp>
using namespace boost::multiprecision;
int1024_t legendre(int1024_t, int1024_t);
int1024_t gcd(int1024_t, int1024_t);
// 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<int1024_t> mpz_vector;
template <typename T> // Takes int_vector or mpz_vector
inline void print_vector(const T &x)
{
for (auto y : x)
std::cout << y << ", ";
std::cout << '\n';
}
inline void print_matrix(const matrix &m)
{
for (auto x : m)
print_vector(x);
}
// 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
void factor_smooth(int_vector &factors, bool &is_smooth, int1024_t n, const mpz_vector &factor_base) // n copy
{
for (size_t i = 0; i<factor_base.size(); i++)
{
int1024_t factor = factor_base[i];
while (n % factor == 0)
{
n /= factor;
factors[i] ^= 1; // + 1 (mod 2) matrices
}
}
is_smooth = (n == 1);
}
void create_factor_base(mpz_vector &factor_base, const int_vector &primes, const int1024_t &n)
{
int1024_t 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;
int1024_t p_mpz = p;
// Use only primes that match (n|p) = 1
if (legendre(p, n) == 1)
factor_base.push_back(p);
}
if (DEBUG)
{
std::cout << "Factor base: ";
print_vector(factor_base);
}
}
int1024_t legendre(int1024_t a, int1024_t m)
{
int1024_t prime = m;
int t = 1;
int c = 0;
int1024_t temp;
a = a%m;
while (a != 0)
{
c = 0;
while (a % 2 == 0)
{
a = a / 2;
c = 1 - c;
}
if (c == 1)
{
if (m % 8 == 3 || m % 8 == 5)
{
t = -t;
}
}
if (a % 4 == 3 & m % 4 == 3)
{
t = -t;
}
temp = m;
m = a;
a = temp;
a = a%m;
}
return m;
}
void sieve(mpz_vector &smooth_numbers, mpz_vector &smooth_x, matrix &smooth_factors,
const int1024_t &n, const mpz_vector &factor_base)
{
int smooth_count = 0;
int1024_t sqrt_n = sqrt(n);
bool not_done = true;
int1024_t j = 1;
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++)
{
int1024_t current;
int1024_t 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 = pow(x, 2) % n;
current_chunk[i] = current;
current_x[i] = x;
}
// To do: add Shanks-Tonelli
j += SIEVE_CHUNK;
// Actual factoring
for (size_t i = 0; i<current_chunk.size(); i++)
{
// Each item in factors corresponds to number in factor base
int_vector factors(factor_base.size(), 0);
bool is_smooth;
factor_smooth(factors, is_smooth, current_chunk[i], factor_base);
if (is_smooth) // Is smooth
{
if (smooth_count > int(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] = factors;
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";
print_matrix(smooth_factors);
std::cout << '\n';
}
}
void gaussian_elimination(int1024_t &y, const matrix &smooth_factors,
const mpz_vector &smooth_numbers, const mpz_vector &smooth_x, int1024_t &x)
{
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";
print_matrix(A);
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";
print_matrix(A);
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];
int1024_t x_square = 1;
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;
}
int1024_t rem;
//mpz_sqrtrem(x.get_mpz_t(), rem.get_mpz_t(), x_square.get_mpz_t());
x = (int1024_t)floor((int)sqrt(x_square));
rem = x_square - pow(x, 2);
if (DEBUG)
{
if (rem == 0)
std::cout << "Remainder 0\n";
std::cout << "x: " << x << '\n' << "y: " << y << "\n\n";
}
}
int main()
{
// Test numbers: 502560280658509, 90283
const int1024_t n("502560280658509");
int_vector primes = eratosthenes(TRIAL_BOUND);
mpz_vector factor_base;
create_factor_base(factor_base, primes, n);
// Find smooth numbers with x = sqrt(n) + j
// 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 });
sieve(smooth_numbers, smooth_x, smooth_factors, n, factor_base);
// Gaussian Elimination -----------------------------------
// Transpose the matrix
int1024_t y, x;
gaussian_elimination(y, smooth_factors, smooth_numbers, smooth_x, x);
// Final division
int1024_t factor_1;
int1024_t dif = y - x;
factor_1 = gcd(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";
int1024_t factor_2 = n / factor_1;
std::cout << "Factor 1: " << factor_1 << '\n';
std::cout << "Factor 2: " << factor_2 << '\n';
return 0;
}
int1024_t gcd(int1024_t m, int1024_t n) // function definition
{ // block begin
int1024_t r; // declaration of remainder
while (n != 0) { // not equal
r = m % n; // modulus operator
m = n; // assignment
n = r;
} // end while loop
if (m >= 0) // exit gcd with value m
return m;
else
return -m;
}