#include <iostream>
#include <mpir.h>
#include <vector>
#include <stdint.h>
using namespace std;
vector<unsigned long long> factor_base;
vector<unsigned long long> start;
vector<unsigned long long> factor_base_index;
//ms == 'matrix size'
const int ms = 21000;
uint8_t vec[ms][ms];
uint8_t ident[ms][ms];
//dcs == 'data collection size'
unsigned long long dcs;
unsigned long long pointer = 0;
unsigned long long closenuf;
//solution found boolean
bool sol_found = false;
//T is used to calculate the threshold when sieving
double T;
void sieve(unsigned long long b1);
void refine(const mpz_t n);
void data_collection(const mpz_t n, mpz_t *tX, mpz_t *tY, mpz_t *zero, mpz_t xprog);
void solve_congruences(const mpz_t n);
void quadratic_sieve(mpz_t *zero, mpz_t *X, mpz_t *tX, mpz_t *Y, mpz_t *tY);
void build_matrix(mpz_t *B);
void forward_elimination();
void backwards_elimination();
void find_solution(const mpz_t n, mpz_t *Y, mpz_t *X);
void update_congruences();
void calculate_threshold(const mpz_t n);
int main(){
//get n, the number to be factored
mpz_t n;
mpz_init(n);
char input[1024];
cout << "n = ";
cin >> input;
mpz_set_str(n, input, 10);
//set xprog, used to save the progress of X, used for data collection
mpz_t xprog;
mpz_init(xprog);
mpz_sqrt(xprog, n);
mpz_add_ui(xprog, xprog, 1);
//get b1, the upper boundary for the factor base
unsigned long long b1;
cout << "b1 = ";
cin >> b1;
//get T, used for calculating the threshold for sieving
cout << "T = ";
cin >> T;
//sieve out all the composite numbers using the sieve of Eratosthenes
cout << "Sieving for prime factors...\n";
sieve(b1);
cout << "Complete!\n";
//refine the factor base, leaving only quadratic residues
cout << "Searching for quadratic residues...\n";
refine(n);
cout << "Complete!\n";
//set the data collection size to 1 greater than the size of the factor base
dcs = factor_base.size() + 1;
//use the values of T and n to calculate the sieving threshold
calculate_threshold(n);
//initialise arrays to be used for data collection
unsigned long long i;
mpz_t *X = new mpz_t[250000];
mpz_t *Y = new mpz_t[250000];
mpz_t *tX = new mpz_t[250000];
mpz_t *tY = new mpz_t[250000];
mpz_t *zero = new mpz_t[250000];
mpz_t *B = new mpz_t[250000];
for(i = 0; i < 250000; i++){
mpz_init_set_str(X[i], "0", 10);
mpz_init_set_str(Y[i], "0", 10);
mpz_init_set_str(B[i], "0", 10);
mpz_init_set_str(tX[i], "0", 10);
mpz_init_set_str(tY[i], "0", 10);
mpz_init_set_str(zero[i], "0", 10);
}
//the congruences are stored in the arrays 'start' and 'factor_base_index'
//'start' contains the index value in the tY array to start the loop on
//'factor_base_index' contains the value of p, the prime which the indexed value in tY is divisible by
solve_congruences(n);
cout << "Collecting data...\n";
while(pointer <= dcs){
data_collection(n, tX, tY, zero, xprog);
quadratic_sieve(zero, X, tX, Y, tY);
update_congruences();
}
for(i = 0; i <= dcs; i++){
mpz_set(B[i], Y[i]);
}
//use the values from the B array for trial division to build the matrix
cout << "Building matrix...\n";
build_matrix(B);
cout << "Complete!\n";
mpz_clear(*tX);
mpz_clear(*tY);
mpz_clear(*B);
while(sol_found == false){
//call the gaussian elimination functions
cout << "Reducing matrix...\n";
forward_elimination();
backwards_elimination();
cout << "Looking for a solution...\n";
find_solution(n, Y, X);
cout << "fb: " << factor_base.size() << "\ndcs:" << dcs << "\n";
}
system("pause");
return 0;
}
void sieve(unsigned long long b1){
//sieve of Eratosthenes
vector<int> consecutive;
unsigned long long p = 0;
unsigned long long i;
cout << "Building Prime Base...\n";
consecutive.push_back(0);
consecutive.push_back(0);
consecutive.push_back(2);
for(i = 3; i <= b1; i+=2) {
consecutive.push_back(i);
consecutive.push_back(0);
}
while(p <= b1) {
if(consecutive[p] != 0) {
p = consecutive[p];
for(i = p; i <= b1; i+=p) {
if(i > p && consecutive[i] != 0) {
consecutive[i] = 0;
}
}
}
p++;
}
for(i = 0; i <= b1; i++) {
if(consecutive[i] != 0) {
factor_base.push_back(consecutive[i]);
}
}
consecutive.erase(consecutive.begin(), consecutive.end());
}
void refine(const mpz_t n){
mpz_t temp;
mpz_init(temp);
unsigned long long i;
//loop through the factor base, if the legendre symbol isn't +1, discard that factor
for(i = 1; i < factor_base.size(); i++){
mpz_set_ui(temp, factor_base[i]);
if(mpz_legendre(n, temp) != 1){
factor_base.erase(factor_base.begin()+i);
i--;
}
}
mpz_clear(temp);
}
void data_collection(const mpz_t n, mpz_t *tX, mpz_t *tY, mpz_t *zero, mpz_t xprog){
mpz_t temp;
mpz_init(temp);
unsigned long long cntr = 0;
//xprog saves the last value that X was set to, so we can start from where we left off
//calculate x^2 - n and save value in the tX array if it is greater than 1
//reset the zero array to... zero
while(cntr < 250000){
mpz_set(tX[cntr], xprog);
mpz_mul(temp, xprog, xprog);
mpz_sub(temp, temp, n);
mpz_set_ui(zero[cntr], 0);
if(mpz_cmp_ui(temp, 1) > 0){
mpz_set(tY[cntr], temp);
cntr++;
}
mpz_add_ui(xprog, xprog, 1);
}
mpz_clear(temp);
}
void solve_congruences(const mpz_t n){
//calculate ceiling of square root of n
mpz_t sqrt;
mpz_init(sqrt);
mpz_sqrt(sqrt, n);
mpz_add_ui(sqrt, sqrt, 1);
//solve the case of (X+sqrt)^2 - n = 0 (mod 2)
unsigned long long X = 0;
mpz_t temp;
mpz_init(temp);
//increment X until a solution is found
while(1){
mpz_add_ui(temp, sqrt, X);
mpz_mul(temp, temp, temp);
mpz_sub(temp, temp, n);
mpz_mod_ui(temp, temp, 2);
if(mpz_cmp_ui(temp, 0) == 0){
start.push_back(X);
factor_base_index.push_back(2);
break;
}
X++;
}
//now iterate through the factor base, and solve (X+sqrt)^2 - n = 0 (mod p)
//once solved, push X onto start array, and p onto factor_base_index. Repeat twice for both solutions.
unsigned long long i;
unsigned long long solved;
for(i = 1; i < factor_base.size(); i++){
X = 0;
solved = 0;
while(solved < 2){
mpz_add_ui(temp, sqrt, X);
mpz_mul(temp, temp, temp);
mpz_sub(temp, temp, n);
mpz_mod_ui(temp, temp, factor_base[i]);
if(mpz_cmp_ui(temp, 0) == 0){
start.push_back(X);
factor_base_index.push_back(factor_base[i]);
solved++;
}
X++;
}
cout << "Solved " << i << " of " << factor_base.size() << " congruences\r";
}
mpz_clear(sqrt);
mpz_clear(temp);
}
void quadratic_sieve(mpz_t *zero, mpz_t *X, mpz_t *tX, mpz_t *Y, mpz_t *tY){
mpz_t temp;
mpz_init(temp);
unsigned long long s;
unsigned long long st;
unsigned long long inc;
unsigned long long cntr;
mpz_t base2;
mpz_init(base2);
size_t sz;
for(s = 0; s < start.size(); s++){
//starting at st, and incrementing by inc, divide each value by inc
inc = factor_base_index[s];
mpz_set_ui(base2, inc);
sz = mpz_sizeinbase(base2, 2);
for(st = start[s]; st < 250000; st += inc){
//increment the zero array by ~log2(inc)
mpz_add_ui(zero[st], zero[st], sz);
}
}
for(s = 0; s < 250000; s++){
//values which are at or above the threshold (closenuf) enter the trial division stage
if(mpz_cmp_ui(zero[s], closenuf) >= 0){
//the potential smooth number is backed up in the X and Y arrays
mpz_set(X[pointer], tX[s]);
mpz_set(Y[pointer], tY[s]);
//the potential smooth number is divided by each factor in the factor base
for(cntr = 0; cntr < factor_base.size(); cntr++){
mpz_mod_ui(temp, tY[s], factor_base[cntr]);
while(mpz_cmp_ui(temp, 0) == 0){
mpz_div_ui(tY[s], tY[s], factor_base[cntr]);
mpz_mod_ui(temp, tY[s], factor_base[cntr]);
}
//if the number is smooth increment the pointer, otherwise the value in the X/Y arrays will be
//overwritten on the next iteration of the loop
if(mpz_cmp_ui(tY[s], 1) == 0){
pointer++;
cout << pointer << "/" << dcs << " \r";
break;
}
}
}
}
mpz_clear(temp);
mpz_clear(base2);
}
void update_congruences(){
unsigned long long difference;
unsigned long long overlap;
unsigned long long offset;
unsigned long long i;
for(i = 0; i < start.size(); i++){
difference = 250000 - start[i];
overlap = difference % factor_base_index[i];
offset = factor_base_index[i] - overlap;
start[i] = offset;
}
}
void build_matrix(mpz_t *B){
//initialise arrays
unsigned long long outer;
unsigned long long inner;
for(outer = 0; outer < ms; outer++){
for(inner = 0; inner < ms; inner++){
vec[outer][inner] = 0;
ident[outer][inner] = 0;
}
}
for(outer = 0; outer < ms; outer++){
ident[outer][outer] = 1;
}
//build matrix
mpz_t temp;
mpz_init(temp);
uint8_t flipper = 0;
unsigned long long row = 0;
unsigned long long col = 0;
for(outer = 0; outer < dcs; outer++){
while(col < factor_base.size()){
mpz_mod_ui(temp, B[outer], factor_base[col]);
if(mpz_cmp_ui(temp, 0) == 0){
flipper++;
flipper = flipper % 2;
vec[row][col] = flipper;
mpz_div_ui(B[outer], B[outer], factor_base[col]);
}
else
{
col++;
flipper = 0;
}
}
if(mpz_cmp_ui(B[outer], 1) == 0){
row++;
}
col = 0;
flipper = 0;
if(row > dcs){
break;
}
}
dcs = row;
mpz_clear(temp);
}
void forward_elimination(){
long int pivot;
long int row;
long int i;
long int row_iterator;
long int size;
size = factor_base.size();
for(pivot = 0; pivot < size; pivot++){
cout << pivot << "/" << size << "\r";
for(row = pivot; row < dcs-1; row++){
if(vec[row][pivot] == 1){
for(row_iterator = row+1; row_iterator < dcs; row_iterator++){
if(vec[row_iterator][pivot] == 1){
for(i = 0; i < size; i++){
if(vec[row][i] == 1){
vec[row_iterator][i] += vec[row][i];
vec[row_iterator][i] = vec[row_iterator][i] % 2;
}
}
for(i = 0; i < dcs; i++){
if(ident[row][i] == 1){
ident[row_iterator][i] += ident[row][i];
ident[row_iterator][i] = ident[row_iterator][i] % 2;
}
}
}
}
}
}
}
}
void backwards_elimination(){
long int pivot;
long int row;
long int i;
long int row_iterator;
long int size;
size = factor_base.size();
for(pivot = size - 1; pivot >= 0; pivot--){
cout << pivot << "/" << size << "\r";
for(row = pivot; row > 0; row--){
if(vec[row][pivot] == 1){
for(row_iterator = row-1; row_iterator >= 0; row_iterator--){
if(vec[row_iterator][pivot] == 1){
for(i = 0; i < size; i++){
if(vec[row][i] == 1){
vec[row_iterator][i] += vec[row][i];
vec[row_iterator][i] = vec[row_iterator][i] % 2;
}
}
for(i = 0; i < dcs; i++){
if(ident[row][i] == 1){
ident[row_iterator][i] += ident[row][i];
ident[row_iterator][i] = ident[row_iterator][i] % 2;
}
}
}
}
}
}
}
}
void find_solution(const mpz_t n, mpz_t *Y, mpz_t *X){
unsigned long long i;
unsigned long long c;
unsigned long long indexer;
bool flag;
vector<unsigned long long> luckydip;
mpz_t Xt;
mpz_t Yt;
mpz_t temp;
mpz_init(Xt);
mpz_init(Yt);
mpz_init(temp);
for(i = 0; i < dcs; i++){
mpz_set_ui(Xt, 1);
mpz_set_ui(Yt, 1);
flag = true;
for(c = 0; c < factor_base.size(); c++){
if(vec[i][c] != 0){
flag = false;
break;
}
}
if(flag == true){
for(c = 0; c < dcs; c++){
if(ident[i][c] == 1){
luckydip.push_back(c);
}
}
for(indexer = 0; indexer < luckydip.size(); indexer++){
mpz_mul(Xt, Xt, X[luckydip[indexer]]);
mpz_mul(Yt, Yt, Y[luckydip[indexer]]);
}
mpz_sqrt(Yt, Yt);
mpz_sub(temp, Xt, Yt);
mpz_gcd(temp, temp, n);
if(mpz_cmp_ui(temp, 1) != 0 && mpz_cmp(temp, n) != 0){
mpz_out_str(stdout, 10, temp);
cout << "\n\a";
mpz_div(temp, n, temp);
mpz_out_str(stdout, 10, temp);
cout << "\n";
sol_found = true;
return;
}
}
luckydip.erase(luckydip.begin(), luckydip.end());
}
}
void calculate_threshold(const mpz_t n){
unsigned long long target;
size_t sz;
mpz_t base2;
mpz_init(base2);
sz = mpz_sizeinbase(n, 2);
target = sz / 2;
mpz_set_ui(base2, 125000);
sz = mpz_sizeinbase(base2, 2);
target += sz;
mpz_set_ui(base2, factor_base[factor_base.size() - 1]);
sz = mpz_sizeinbase(base2, 2);
closenuf = T * sz;
closenuf = target - closenuf;
cout << "TARGET = " << target << " bits\n";
cout << "THRESHOLD = " << closenuf << " bits\n";
cout << "T = " << T << endl;
system("pause");
}
