ECM

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