#include <iostream>
#include <mpir.h>
#include <vector>

using namespace std;

void NextValues(mpz_t X, mpz_t Z, mpz_t k, mpz_t n, mpz_t a, mpz_t b);
void x_sub_2i(mpz_t ret, mpz_t r, mpz_t s, mpz_t a, mpz_t b, mpz_t n);
void z_sub_2i(mpz_t ret, mpz_t r, mpz_t s, mpz_t a, mpz_t b, mpz_t n);
void x_sub_2i_plus_1(mpz_t ret, mpz_t r, mpz_t s, mpz_t u, mpz_t v, mpz_t a, mpz_t b, mpz_t n, mpz_t Z);
void z_sub_2i_plus_1(mpz_t ret, mpz_t r, mpz_t s, mpz_t u, mpz_t v, mpz_t n, mpz_t X);
void ECM(mpz_t n, mpz_t X, mpz_t Y, mpz_t a, mpz_t MAX);

int main(){
	mpz_t n;
	mpz_t X;
	mpz_t Y;
	mpz_t a;
	mpz_t MAX;
	mpz_t rn;
	mpz_init(n);
	mpz_init(X);
	mpz_init(Y);
	mpz_init(a);
	mpz_init(MAX);
	gmp_randstate_t r_state;
	
	char input[512];
    unsigned long int seed;
    
	cout << "ECM\n";
	cout << "n = ";
	cin >> input;
	mpz_set_str(n, input, 10);
	cout << "Random Seed: ";
	cin >> seed;
	cout << "max = ";
	cin >> input;
	mpz_set_str(MAX, input, 10);
	
	gmp_randinit_default (r_state);
    gmp_randseed_ui(r_state, seed);
	mpz_init(rn);

	while(1){
		mpz_urandomm(rn,r_state,n);
		mpz_set(X, rn);
		mpz_urandomm(rn,r_state,n);
		mpz_set(Y, rn);
		mpz_urandomm(rn,r_state,n);
		mpz_set(a, rn);
		ECM(n, X, Y, a, MAX);
	}
	
	return 0;
}

void ECM(mpz_t n, mpz_t X, mpz_t Y, mpz_t a, mpz_t MAX){
	mpz_t b;
	mpz_t Y2;
	mpz_t X3;
	mpz_t aX;
	mpz_t a3;
	mpz_t b2;
	mpz_t g;
	mpz_t Z;
	mpz_t k;
	mpz_t temp1;
	mpz_t temp2;
	mpz_init(b);
	mpz_init(Y2);
	mpz_init(X3);
	mpz_init(aX);
	mpz_init(a3);
	mpz_init(b2);
	mpz_init(g);
	mpz_init(Z);
	mpz_init(k);
	mpz_init(temp1);
	mpz_init(temp2);
	unsigned long long i;

	//b = Y^2 - X^3 - aX mod n
	mpz_mul(Y2, Y, Y);
	mpz_pow_ui(X3, X, 3);
	mpz_mul(aX, a, X);
	mpz_sub(temp1, Y2, X3);
	mpz_sub(temp1, temp1, aX);
	mpz_mod(b, temp1, n);

	//g = gcd(4a^3 + 27b^2, n)
	mpz_pow_ui(a3, a, 3);
	mpz_mul(b2, b, b);
	mpz_mul_si(temp1, a3, 4);
	mpz_mul_si(temp2, b2, 27);
	mpz_add(temp1, temp1, temp2);
	mpz_gcd(g, temp1, n);

	if(mpz_cmp_ui(g, 1) != 0){
		mpz_out_str(stdout, 10, g);
		cout << "\n\a";
		system("pause");
		return;
	}

	mpz_set_si(Z, 1);
	mpz_set_si(k, 2);

	while(mpz_cmp(k, MAX) <= 0){
		for(i = 1; i <= 10; i++){
			NextValues(X, Z, k, n, a, b);
			mpz_add_ui(k, k, 1);
		}

		mpz_gcd(g, Z, n);
		if(mpz_cmp_ui(g, 1) != 0){
			mpz_out_str(stdout, 10, g);
			cout << "\n\a";
			system("pause");
			return;
		}
	}

	mpz_clear(b);
	mpz_clear(Y2);
	mpz_clear(X3);
	mpz_clear(aX);
	mpz_clear(a3);
	mpz_clear(b2);
	mpz_clear(g);
	mpz_clear(Z);
	mpz_clear(k);
	mpz_clear(temp1);
	mpz_clear(temp2);
}

void NextValues(mpz_t X, mpz_t Z, mpz_t k, mpz_t n, mpz_t a, mpz_t b){
	unsigned int i = 0;
	unsigned int length;
	vector<int> C;
	mpz_t ik;
	mpz_init(ik);
	mpz_set(ik, k); //ik is used to preserve k
	mpz_t temp1;
	mpz_init(temp1);
	mpz_t X1;
	mpz_t Z1;
	mpz_t X2;
	mpz_t Z2;
	mpz_init(X1);
	mpz_init(Z1);
	mpz_init(X2);
	mpz_init(Z2);
	mpz_t U1;
	mpz_t U2;
	mpz_init(U1);
	mpz_init(U2);
	mpz_t temp;
	mpz_init(temp);

	//binary expansion of ik
	while(mpz_cmp_ui(ik, 0) > 0){
		mpz_mod_ui(temp1, ik, 2);
		if(mpz_cmp_ui(temp1, 0) == 0){
			C.push_back(0);
			mpz_div_ui(ik, ik, 2);
		}
		else
		{
			C.push_back(1);
			mpz_sub_ui(ik, ik, 1);
			mpz_div_ui(ik, ik, 2);
		}
	}

	length = C.size();

	mpz_set(X1, X);
	mpz_set(Z1, Z);
	
	x_sub_2i(X2, X, Z, a, b, n);
	z_sub_2i(Z2, X, Z, a, b, n);

	for(i = length - 1; i >= 1; i--){
		x_sub_2i_plus_1(U1, X1, Z1, X2, Z2, a, b, n, Z);
		z_sub_2i_plus_1(U2, X1, Z1, X2, Z2, n, X);

		if(C[i] == 0){
			x_sub_2i(temp, X1, Z1, a, b, n);
			z_sub_2i(Z1, X1, Z1, a, b, n);
			mpz_set(X1, temp);
			mpz_set(X2, U1);
			mpz_set(Z2, U2);
		}
		else
		{
			x_sub_2i(temp, X2, Z2, a, b, n);
			z_sub_2i(Z2, X2, Z2, a, b, n);
			mpz_set(X2, temp);
			mpz_set(X1, U1);
			mpz_set(Z1, U2);
		}
	}
	
	C.erase(C.begin(), C.end());
	mpz_set(X, X1);
	mpz_set(Z, Z1);

	mpz_clear(ik);
	mpz_clear(temp1);
	mpz_clear(X1);
	mpz_clear(Z1);
	mpz_clear(X2);
	mpz_clear(Z2);
	mpz_clear(U1);
	mpz_clear(U2);
	mpz_clear(temp);

	return;
}

void x_sub_2i(mpz_t ret, mpz_t r, mpz_t s, mpz_t a, mpz_t b, mpz_t n){
	mpz_t term;
	mpz_t r2;
	mpz_t s2;
	mpz_t value;
	mpz_t term2;
	mpz_t s3;
	mpz_t temp;
	mpz_init(term);
	mpz_init(r2);
	mpz_init(s2);
	mpz_init(value);
	mpz_init(term2);
	mpz_init(s3);
	mpz_init(temp);

	//term = r^2 - a * s^2 MOD n
	mpz_mul(r2, r, r);
	mpz_mul(s2, s, s);
	mpz_mul(temp, a, s2);
	mpz_sub(temp, r2, temp);
	mpz_mod(term, temp, n);

	//value = term^2 - 8 * b * r * s^3 MOD n
	mpz_mul(term2, term, term);
	mpz_pow_ui(s3, s, 3);
	mpz_mul_si(temp, b, 8);
	mpz_mul(temp, temp, r);
	mpz_mul(temp, temp, s3);
	mpz_sub(temp, term2, temp);
	mpz_mod(value, temp, n);

	mpz_set(ret, value);

	mpz_clear(term);
	mpz_clear(r2);
	mpz_clear(s2);
	mpz_clear(value);
	mpz_clear(term2);
	mpz_clear(s3);
	mpz_clear(temp);
	return;
}

void z_sub_2i(mpz_t ret, mpz_t r, mpz_t s, mpz_t a, mpz_t b, mpz_t n){
	mpz_t term;
	mpz_t r3;
	mpz_t s2;
	mpz_t s3;
	mpz_t value;
	mpz_t temp;
	mpz_init(term);
	mpz_init(r3);
	mpz_init(s2);
	mpz_init(s3);
	mpz_init(value);
	mpz_init(temp);
	mpz_t temp2;
	mpz_init(temp2);

	//term = r^3 + a * r * s^2 + b * s^3 MOD n
	mpz_pow_ui(r3, r, 3);
	mpz_pow_ui(s2, s, 2);
	mpz_pow_ui(s3, s, 3);
	mpz_mul(temp, a, r);
	mpz_mul(temp, temp, s2);
	mpz_mul(temp2, b, s3);
	mpz_add(temp, r3, temp);
	mpz_add(temp, temp, temp2);
	mpz_mod(term, temp, n);

	//value = 4 * s * term MOD n
	mpz_mul_ui(temp, s, 4);
	mpz_mul(temp, temp, term);
	mpz_mod(value, temp, n);
	mpz_set(ret, value);

	mpz_clear(term);
	mpz_clear(r3);
	mpz_clear(s2);
	mpz_clear(s3);
	mpz_clear(value);
	mpz_clear(temp);
	mpz_clear(temp2);
	return;
}

void x_sub_2i_plus_1(mpz_t ret, mpz_t r, mpz_t s, mpz_t u, mpz_t v, mpz_t a, mpz_t b, mpz_t n, mpz_t Z){
	mpz_t term1;
	mpz_t term2;
	mpz_t term12;
	mpz_t value;
	mpz_t ru;
	mpz_t asv;
	mpz_t bsv;
	mpz_t rv;
	mpz_t su;
	mpz_t temp1;
	mpz_t temp2;
	mpz_t temp3;
	mpz_t temp4;
	mpz_init(term1);
	mpz_init(term2);
	mpz_init(term12);
	mpz_init(value);
	mpz_init(ru);
	mpz_init(asv);
	mpz_init(bsv);
	mpz_init(rv);
	mpz_init(su);
	mpz_init(temp1);
	mpz_init(temp2);
	mpz_init(temp3);
	mpz_init(temp4);

	//term1 = r * u - a * s * v MOD n
	mpz_mul(ru, r, u);
	mpz_mul(asv, a, s);
	mpz_mul(asv, asv, v);
	mpz_sub(temp1, ru, asv);
	mpz_mod(term1, temp1, n);

	//term2 = b * s * v * (r * v + s * u) MOD n
	mpz_mul(bsv, b, s);
	mpz_mul(bsv, bsv, v);
	mpz_mul(rv, r, v);
	mpz_mul(su, s, u);
	mpz_add(temp1, rv, su);
	mpz_mul(temp1, bsv, temp1);
	mpz_mod(term2, temp1, n);

	//value = Z * ( term1^2 - 4 * term2 ) MOD n
	mpz_pow_ui(term12, term1, 2);
	mpz_mul_ui(temp1, term2, 4);
	mpz_sub(temp1, term12, temp1);
	mpz_mul(temp1, Z, temp1);
	mpz_mod(value, temp1, n);
	mpz_set(ret, value);

	mpz_clear(term1);
	mpz_clear(term2);
	mpz_clear(term12);
	mpz_clear(value);
	mpz_clear(ru);
	mpz_clear(asv);
	mpz_clear(bsv);
	mpz_clear(rv);
	mpz_clear(su);
	mpz_clear(temp1);
	mpz_clear(temp2);
	mpz_clear(temp3);
	mpz_clear(temp4);
	return;
}

void z_sub_2i_plus_1(mpz_t ret, mpz_t r, mpz_t s, mpz_t u, mpz_t v, mpz_t n, mpz_t X){
	mpz_t term;
	mpz_t us;
	mpz_t rv;
	mpz_t term2;
	mpz_t value;
	mpz_t temp;
	mpz_init(term);
	mpz_init(us);
	mpz_init(rv);
	mpz_init(term2);
	mpz_init(value);
	mpz_init(temp);
	
	//term = u * s - r * v MOD n
	mpz_mul(us, u, s);
	mpz_mul(rv, r, v);
	mpz_sub(temp, us, rv);
	mpz_mod(term, temp, n);

	//value = X * term^2 MOD n
	mpz_pow_ui(term2, term, 2);
	mpz_mul(temp, term2, X);
	mpz_mod(value, temp, n);
	mpz_set(ret, value);

	mpz_clear(term);
	mpz_clear(us);
	mpz_clear(rv);
	mpz_clear(term2);
	mpz_clear(value);
	mpz_clear(temp);
	return;
}