// exit(0) iff is 2-Euler PRP and Lucas PRP
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>

int j, j2=0, i;
mpz_t n, np1, nm1, nm1d2, t, r, d, g, a, u, v, ures, vres;

int main(int argc, char *argv[]) {
    if (argc==2) mpz_init_set_str(n, argv[1], 10);
    else mpz_inp_str(n, stdin, 10);
    // trivialities
    if (mpz_cmp_ui(n, 0)==-1) exit(1);
    if (mpz_perfect_square_p(n)) exit(1);
    if (mpz_cmp_ui(n, 2)==0 || mpz_cmp_ui(n, 5)==0) exit(0);
    // Trial division
    mpz_init(r);
    mpz_init(d);
    mpz_sqrt(r, n);
    if (mpz_cmp_ui(r,100000)==1) j=100000; else j=mpz_get_ui(r);
    for(i=2; i<=j; i++) {
        mpz_mod_ui(d, n, i);
        if (mpz_cmp_ui(d, 0)==0) exit(1);
    }
    if (mpz_cmp_ui(r, --i)==0) exit(0);
    // Euler PRP test
    mpz_init(nm1);
    mpz_init(nm1d2);
    mpz_sub_ui(nm1, n, 1);
    mpz_div_ui(nm1d2, nm1, 2);
    mpz_init(t);
    mpz_set_ui(t, 2);
    j = mpz_jacobi(t, n);
    mpz_powm(r, t, nm1d2, n);
    if (mpz_cmp_ui(r, 1)==0) j2=1;
    else if (mpz_cmp(r, nm1)==0) j2=-1;
    if (j!=j2) exit(1);
    // find strong discriminant and appropriate gcd
    mpz_init(g);
    i=2;
    mpz_set_ui(d, 8);
    mpz_set_ui(t, 4);
    mpz_gcd_ui(g, nm1, 3);
    while (mpz_jacobi(d, n)!=-1 || mpz_cmp_ui(g, 2)==-1) {
        i++;
        mpz_pow_ui(r, t, i);
        mpz_sub_ui(d, r, 8);
        mpz_gcd_ui(g, nm1, (i-1)*(2*i-1));

    }
    // Lucas chain for (n+1)/2
    mpz_init(a);
    mpz_init(u);
    mpz_init(v);
    mpz_init(ures);
    mpz_init(vres);
    mpz_init(np1);
    mpz_add_ui(np1, n, 1);
    mpz_set_ui(t, 2);
    mpz_pow_ui(r, t, 2*i-1);
    mpz_sub_ui(a, r, 2);
    mpz_set_ui(u, 2);
    mpz_set(v, a);
    for (i = mpz_sizeinbase(np1, 2)-1; i > 0; i--) {
        if (mpz_tstbit( np1, i )) { 
            mpz_mul( u, u, v);
            mpz_sub(u, u, a);
            mpz_mul(v, v, v);
            mpz_sub_ui(v, v, 2);
        } else {
            mpz_mul(v, u, v);
            mpz_sub(v, v, a);
            mpz_mul(u, u, u);
            mpz_sub_ui(u, u, 2);
        }
        mpz_mod(u, u, n);
        mpz_mod(v, v, n);
    }
    mpz_set_ui(t, 2);
    mpz_mul_si(ures, t, j);
    mpz_mul_si(vres, a, j);
    mpz_mod(ures, ures, n); 
    mpz_mod(vres, vres, n);
    if (mpz_cmp( u, ures )!=0 || mpz_cmp(v, vres)!=0) exit(1);  
    exit(0);
}

// gcc -o prp_v2 prp_v2.c -lgmp
// bash use case:-
/*
for i in {1..2281}; do 
    if ./prp_v2 $i; then 
        if echo "print(2^$i-1)" | gp -q | ./prp_v2; then
            echo "M$i Is 2-Euler and Lucas PRP";
        fi
    fi
done
*/