/****** Written by: Allan Menezes July 2018 ****/
#include "stdio.h"
#include "memory.h"
#include "stdlib.h"
#include "string.h"
#include "math.h"
#include <gmp.h>
#include "time.h"

#define S1 2		
#define S2 2
#define N1 100000000
//define N2 100000000000000000
#define MAXPROCS 50
#define CHUNK 1
#define WORKTAG 1
#define DIETAG 2/* Local functions */
int
main(int argc, char **argv)
{
  
  unsigned long int i,mm=0,n,MM,NN,pexit=0,pexit1=0,jj=0,a,km,phi=0,mk=0,bb=0,pq[20],kk=0,insertp=0,pow=0,once,bmax=0;
 // char ch[100],chh,aa[100];
  FILE * fpw, * fpo, *fpi;
  time_t start,finish;
  time_t    tt;
struct tm *tm;
 char   bufi[255],buff[255];
   double e,t1,t2;
 unsigned long int j,k,s,pm,mn,b,b1,b2,ja;
 
mpz_t p,q,r;
mpz_init(p);
mpz_init(q);
mpz_init(r);
  char ch,c[100000],chm[10000];

//N2=strtol(argv[2],NULL,10);
fpw=fopen(argv[1],"w");
//fpi=fopen(argv[3],"r");
fpo=fopen(argv[2],"r");
s=0;
b=2;
while(fgets(c,100,fpo)!=NULL)
{

jj++;
 if(jj%1000==0)
  fprintf(stderr,"\n %ld ",jj);
a=strtol(c,NULL,10);
n=(a-1)/2;
mpz_set_ui(r,a);
mpz_set_ui(p,b);
ja=mpz_jacobi(p,r);
b1=2;
j=ja;
mpz_set_ui(q,n);
mpz_powm(q,p,q,r);
mpz_set_si(p,ja);
mpz_sub(p,q,p);
mpz_mod(q,p,r);

 if (mpz_cmp_ui(q,0)==0)
  {
b=3;
mpz_set_ui(p,b);
while(ja==j|ja==0)
{
b++;
mpz_set_ui(p,b);
ja=mpz_jacobi(p,r);
}

b2=b;
 if (b>bmax)
     bmax=b;

mpz_set_ui(q,n);
mpz_powm(q,p,q,r);
mpz_set_si(p,ja);
mpz_sub(p,q,p);
mpz_mod(q,p,r);
 if (mpz_cmp_ui(q,0)==0)
  {
  s++;
  fprintf(fpw,"%ld. %ld %ld %ld\n",s,b1,b2,a);
  fprintf(stderr,"%ld. %ld %ld %ld\n",s,b1,b2,a);
  }
}
}
fprintf(fpw,"\n The total number of spsps base 2 and bmax =%ld of <= 64 bits is %ld",bmax,s);

fprintf(stderr,"\n The total number of spsps base 2 and bmax =%ld  of <= 64 bits is %ld",bmax,s);
//fprintf(fpw,"\n The total number s.t. for least p, p-1|n-1 is %ld and exceptions are %ld",i-1,kk);
fclose(fpw);
fclose(fpo);
}