#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;
}