// Compile with: gcc -o gmp_optimised_lucas gmp_optimised_lucas.c -lm -lgmp

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <gmp.h>

int main(int argc, char* argv[]) {
    char str[1000000], *string=str;
    int i, r, r_times_2_minus_2;
    mpz_t a;
    mpz_t n;
    mpz_t align;
    mpz_t gcd;
    mpz_t discriminant;
    mpz_t s;
    mpz_t t;
    mpz_t temp1;
    mpz_t temp2;
    mpz_t temp3;
    mpz_t n_minus_1;
    mpz_t n_plus_1;

    if ( argc != 1 ) {
        printf( "Usage:                      echo n | ./gmp_optimised_lucas\n" );
        printf( "Or like:   echo \"print(n)\" | gp -q | ./gmp_optimised_lucas\n" );
        printf( "Or like: ./pfgw64 -k -f0 -od -q\"n\" | ./gmp_optimised_lucas\n" );
        exit ( -1 );
    }

    gets ( str );
    if ( ( string = strchr( str, ':' ) ) && string[1] == ' ' ) string += 2;
    else string = str;

    mpz_init_set_str ( n, str, 10 );

    if ( mpz_cmp_ui ( n, 2 )==0 || mpz_cmp_ui ( n, 3 )==0 || mpz_cmp_ui ( n, 7 )==0 ) {
        printf ( "Is prime.\n" );
        exit ( 1 );
    }

    mpz_init ( gcd );
    mpz_gcd_ui ( gcd, n, 42 );
    if ( mpz_cmp_ui( gcd, 1 ) != 0 ) {
         printf ( "Is composite.\n" );
         exit ( 1 );
    }

    if ( mpz_perfect_square_p ( n ) ) {
        printf ( "Is a square number.\n" );
        exit ( 1 );
    }

    mpz_init ( n_minus_1 );
    mpz_sub_ui ( n_minus_1, n, 1 );
    mpz_init_set_ui ( discriminant, 8 );
    r = 2; 
    mpz_gcd_ui ( gcd, n_minus_1, 3 );
    while ( mpz_jacobi ( discriminant, n ) != -1 || mpz_cmp_ui ( gcd, 1 ) != 0 ) {
        mpz_add_ui ( discriminant, discriminant, 8 ); 
        mpz_mul_2exp ( discriminant, discriminant, 2 );
        mpz_sub_ui ( discriminant, discriminant, 8 );
        r++;
        mpz_gcd_ui ( gcd, n_minus_1, ( r - 1 ) * ( ( r << 1 ) - 1 ) );
    }
    r_times_2_minus_2 = ( r << 1 ) - 2;
    mpz_clear ( n_minus_1 );
    mpz_clear ( discriminant );
    mpz_clear ( gcd );

    mpz_init(align);
    mpz_mul_2exp ( align, n, 32-mpz_mod_ui ( align, n, 32 ) );

    mpz_init_set_ui ( s, 2 );
    mpz_init_set_ui ( t, 0 );
    mpz_init ( temp1 );
    mpz_init ( temp2 );
    mpz_init ( temp3 );
    mpz_init ( n_plus_1 );
    mpz_add_ui ( n_plus_1, n, 1 );
    for ( i = mpz_sizeinbase ( n_plus_1, 2 ) - 2; i>0; i-- ) {
        mpz_sub ( temp2, t, s );
        mpz_mul_2exp ( temp1, s, r_times_2_minus_2 );
        mpz_add ( temp1, temp1, temp2 );
        mpz_mul ( temp1, temp1, s );
        mpz_mul_2exp ( temp1, temp1, 1 );
        mpz_add ( temp3, t, s );
        mpz_mul ( t, temp2, temp3 );
        mpz_mod ( t, t, align );
        mpz_mod ( s, temp1, align );
        if ( mpz_tstbit ( n_plus_1, i ) ) {
            mpz_mul_2exp ( temp1, s, r_times_2_minus_2 );
            mpz_sub ( temp1, temp1, s );
            mpz_mul_2exp ( temp1, temp1, 1 );
            mpz_add ( temp1, temp1, t );
            mpz_mul_2exp ( temp1, temp1, 1 );
            mpz_neg ( t, s );
            mpz_mul_2exp ( t, t, 1 );
            mpz_mod ( s, temp1, align );
        }
    }
    mpz_clear ( temp1 );
    mpz_clear ( temp2 );
    mpz_clear ( temp3 );
    mpz_clear ( n_plus_1 );
    mpz_clear ( align );
    mpz_mod( s, s, n );
    mpz_mod( t, t, n );
    if ( mpz_cmp_ui ( s, 0 ) == 0 && mpz_cmp_ui ( t, 2 ) == 0 ) {
        mpz_clear ( s );
        mpz_clear ( t );
        printf ( "Is PRP!\n" );
        return ( 0 );
    } else {
        mpz_clear ( s );
        mpz_clear ( t );
        printf ( "Is composite\n" );
        return ( 1 );
    }
 }