#include #include #include #include #include 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_vector; typedef std::vector matrix; typedef std::vector mpz_vector; template // 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 A(bound, 1); A[0] = 0; A[1] = 0; // 0 and 1 aren't prime for (int i = 2; i 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 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= 0; k--) { for (int i = k - 1; i >= 0; i--) { if (A[i][k]) { for (int j = 0; j= 0) // exit gcd with value m return m; else return -m; }