Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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");
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement