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

FILE *in_file;
FILE *dat_file;
FILE *log_file;
int i;
unsigned long N;
mpz_t two;
mpz_t n;
mpz_t np1;
mpz_t nm1;
mpz_t nm1d2;
mpz_t y;
mpz_t s;
mpz_t t;
mpz_t temp;
mpz_t temp2;
mpz_t temp3;

int main ( int argc, char *argv[] ) {
  mpz_init_set_ui ( two, 2 );
  mpz_init ( n );
  mpz_init ( np1 );
  mpz_init ( nm1 );
  mpz_init ( nm1d2 );
  mpz_init ( y );
  mpz_init ( s );
  mpz_init ( t );
  mpz_init ( temp );
  mpz_init ( temp2 );
  mpz_init ( temp3 );
  // seek current line
  in_file = fopen ( "psp2_short", "r" );
  dat_file = fopen ( "x2myxp2.dat", "r" );
    fscanf ( dat_file, "%llu", &N );
    mpz_inp_str ( n, in_file, 10 );
    while( mpz_cmp_ui ( n, N ) != 0 ) mpz_inp_str ( n, in_file, 10 );
  fclose ( dat_file );
  // loop over file
  while ( 1 ) {
    mpz_inp_str ( n, in_file, 10 );
    // reject even numbers
    mpz_mod_ui ( temp, n, 2 );
    if ( mpz_cmp_ui ( temp, 0 ) == 0 ) continue;
    // reject square numbers
    if ( mpz_perfect_square_p ( n ) ) continue;
    // precomputes
      mpz_add_ui ( np1, n, 1 );
      mpz_sub_ui ( nm1, n, 1 );
      mpz_div_ui ( nm1d2, nm1, 2 );
    // reject Mod(2,n)^((n-1)/2)!=jacobi(2,n)
    mpz_powm ( temp, two, nm1d2, n );
    mpz_set_si ( temp2, mpz_jacobi ( two, n ) );
    mpz_mod ( temp2, temp2, n );
    if ( mpz_cmp ( temp, temp2 ) != 0 ) continue;
    // test y=2^r
    mpz_set_ui ( y, 0 );
    while ( 1 ) {
      mpz_add_ui ( y, y, 1 );
      mpz_mul_ui ( y, y, 2 );
      mpz_sub_ui ( y, y, 1 );
      mpz_mod ( y, y, n );
      // check for loop end
      if ( mpz_cmp_ui ( y, 0 ) == 0 ) break;
      // check jacobi(y,n)==-1
      if ( mpz_jacobi ( y, n ) != -1 ) continue;
      // check Mod(y,n)^((n-1)/2)==-1
      mpz_powm ( temp, y, nm1d2, n );
      mpz_add_ui ( temp, temp, 1 );
      mpz_mod ( temp, temp, n );
      if ( mpz_cmp_ui ( temp, 0 ) != 0 ) continue;
      // start Frobenius test (x+2)^(n+1)==-y+4
      mpz_set_ui ( s, 1 );
      mpz_set_ui ( t, 2 );
      for ( i = mpz_sizeinbase ( np1, 2 )-2; i >= 0; i-- ) {
        mpz_mul ( temp, s, s );
        mpz_mul ( temp2, t, t );
        mpz_mul ( temp3, s, t );
        mpz_mul ( temp, temp, y );
        mpz_add ( temp, temp, temp2 );
        mpz_add ( s, temp3, temp3 );
        mpz_mod ( s, s, n );
        mpz_mod ( t, temp, n );
        if ( mpz_tstbit ( np1, i ) ) {
          mpz_mul ( temp, s, y );
          mpz_add ( temp, temp, t );
          mpz_add ( temp, temp, t );
          mpz_add ( s, s, s );
          mpz_add ( s, s, t );
          mpz_mod ( s, s, n );
          mpz_mod ( t, temp, n );
        }
      }
      mpz_add ( t, t, y );
      mpz_sub_ui ( t, t, 4 );
      mpz_mod ( s, s, n );
      mpz_mod ( t, t, n );
    // output and log what is passed
      if ( ( mpz_cmp_ui ( s, 0 ) == 0 ) && ( mpz_cmp_ui ( t, 0 ) == 0 ) ) {
        printf ( "[" );
        mpz_out_str ( NULL, 10, n );
        printf ( ", " );
        mpz_out_str ( NULL, 10, y );
        printf ( "]\n" );
        log_file = fopen ( "x2myxp2.log", "a" );
          fprintf ( log_file, "[");
          mpz_out_str ( log_file, 10, n );
          fprintf ( log_file, ", " );
          mpz_out_str ( log_file, 10, y );
          fprintf ( log_file, "]\n" );
        fclose ( log_file );
      }
    }
    // log current line
    dat_file = fopen ( "x2myxp2.dat", "w" );
    mpz_out_str ( dat_file, 10, n );
    fclose ( dat_file );
  }
}
//Compile: gcc -O2 -o x2myxp2 x2myxp2.c -lgmp