/* http://www.mis.hiroshima-u.ac.jp/~hiroki/major/galois1-e.html */
/* by Hiroki Sumida-Takahashi  2015, November 2 */


/*   f0, 2*deg*gap  <  2^31   */

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

main()
 {
  int p,f,fn,pn1,pn2,k,i1,i2,i3,l,lz,lh,hnum,tmp;
  mpz_t pm,lm,alm,tmp1,tmp2,f0m; 
  FILE *file0,*file1,*file2,*file3,*file4,*file5;

  mpz_init(pm);
  mpz_init(lm);
  mpz_init(alm);
  mpz_init(tmp1);
  mpz_init(tmp2);
  mpz_init(f0m);

  pn1=115000;
  pn2=120000; 
  for(fn=0;fn<60;fn++)
    {	
    if (fn==0){f=5;}
    if (fn==1){f=8;}
    if (fn==2){f=12;}
    if (fn==3){f=13;}
    if (fn==4){f=17;}
    if (fn==5){f=21;}
    if (fn==6){f=24;}
    if (fn==7){f=28;}
    if (fn==8){f=29;}
    if (fn==9){f=33;}
    if (fn==10){f=37;}
    if (fn==11){f=40;}
    if (fn==12){f=41;}
    if (fn==13){f=44;}
    if (fn==14){f=53;}
    if (fn==15){f=56;}
    if (fn==16){f=57;}
    if (fn==17){f=60;}
    if (fn==18){f=61;}
    if (fn==19){f=65;}
    if (fn==20){f=69;}
    if (fn==21){f=73;}
    if (fn==22){f=76;}
    if (fn==23){f=77;}
    if (fn==24){f=85;}
    if (fn==25){f=88;}
    if (fn==26){f=89;}
    if (fn==27){f=92;}
    if (fn==28){f=93;}
    if (fn==29){f=97;}
    if (fn==30){f=101;}
    if (fn==31){f=104;}
    if (fn==32){f=105;}
    if (fn==33){f=109;}
    if (fn==34){f=113;}
    if (fn==35){f=120;}
    if (fn==36){f=124;}
    if (fn==37){f=129;}
    if (fn==38){f=133;}
    if (fn==39){f=136;}
    if (fn==40){f=137;}
    if (fn==41){f=140;}
    if (fn==42){f=141;}
    if (fn==43){f=145;}
    if (fn==44){f=149;}
    if (fn==45){f=152;}
    if (fn==46){f=156;}
    if (fn==47){f=157;}
    if (fn==48){f=161;}
    if (fn==49){f=165;}
    if (fn==50){f=168;}
    if (fn==51){f=172;}
    if (fn==52){f=173;}
    if (fn==53){f=177;}
    if (fn==54){f=181;}
    if (fn==55){f=184;}
    if (fn==56){f=185;}
    if (fn==57){f=188;}
    if (fn==58){f=193;}
    if (fn==59){f=197;}
    for(p=pn1;p<pn2;p++)
      {

      if ((fn==0)&&(p==pn1))
      {
        fn=44;
        f=149;
        p=117167;
      }
   
      mpz_set_ui(pm,p);
      if (mpz_probab_prime_p(pm,20)>0)
        {
        mp_bernc(p,f);
	file0=fopen("tmpbernc.txt","r");
	while(fscanf(file0,"%d\n %d\n %d\n",&p,&k,&f)!=EOF)
	  {
	    mp_padic(p,k,f,&alm,&i1);
            if (i1>1){
	      file1=fopen("bcn1.txt","a+");
	      fprintf(file1,"%d,%d,%d,%d\n",p,k,f,i1);
	      fclose(file1);
	    }
	    if (i1==0){
	      file2=fopen("bcn2.txt","a+");
	      fprintf(file2,"%d,%d,%d,%d\n",p,k,f,i1);
	      fclose(file2);
	    }
	    if (i1>0){
	      file3=fopen("bcn3.txt","a+");
              mpz_divexact_ui(tmp1,alm,p);
              tmp=mpz_get_ui(tmp1);
	      fprintf(file3,"%d,%d,%d,%d,%d\n",p,k,f,tmp,i1);
	      fclose(file3); 
	    }
	    if (i1>0)
	      {
		if (f%p==0){mpz_set_ui(f0m,f);}
                    else{mpz_set_ui(f0m,f*p);}
		lz=1; i2=-1; hnum=0;
		while((i2!=1)&&(hnum<5))
		  {
		    lh=0;
		    while(lh!=1)
		      {
                        mpz_set_ui(lm,1);
                        mpz_addmul_ui(lm,f0m,lz);
			if(mpz_probab_prime_p(lm,20)>0){lh=1;} 
			lz++;
		      }
    		    file5=fopen("prime.txt","a+");
                    mpz_out_str(file5,10,lm);
 		    fprintf(file5,"\n");
		    fclose(file5);
   	            i2=mp_hantei1(p,k,f,lm,alm,0);
                    hnum++;
		  }
  		if(i2!=1){
		  file4=fopen("bcn1.txt","a+");
		  fprintf(file4,"%d,%d,%d,%d\n",p,k,f,i2);
		  fclose(file4);
  /*    	  lz=1; i3=-1; hnum=0;
	  	  while((i3!=1)&&(hnum<5))
		    {
		      lh=0;
		      while(lh!=1)
		        {
                          mpz_set_ui(lm,1);
                          mpz_mul_ui(tmp1,f0m,p);
                          mpz_addmul_ui(lm,tmp1,lz);
			  if(mpz_probab_prime_p(lm,20)>0){lh=1;} 
			  lz++;
		        }
   	              i3=mp_hantei1(p,k,f,lm,alm,1);
                      hnum++;
		    }
	          if(i3!=1){
		    file4=fopen("bcn1k.txt","a+");
		    fprintf(file4,"%d,%d,%d,1,%d\n",p,k,f,i3);
		    fclose(file4);
	          }
  */
		} /* if i2!=1 */
	      } /* if i0>0 */
   	   } /* while */
         fclose(file0);
         } /* prime */
      } /* p */
   } /* fn */
   file5=fopen("end.txt","a+");
   fprintf(file5,"end \n");
   fclose(file5);
   mpz_clear(lm);
   mpz_clear(pm);
   mpz_clear(alm);
   mpz_clear(tmp1);
   mpz_clear(tmp2);
   mpz_clear(f0m);

 } /* main */



mp_bernc(int p, int f)
 {
  int fn1,i,gap,deg,a,kro,fn,ct1=400,ct2=20000,ct3=800000;
  unsigned long long int ii;
  mpz_t im,im1,im2,fm,pm,gcm,tmpm;

  mpz_init(im);
  mpz_init(im1);
  mpz_init(im2);
  mpz_init(pm);
  mpz_init(fm);
  mpz_init(gcm);  
  mpz_init(tmpm);

  mpz_set_ui(pm,p);

  mpz_mul(tmpm,pm,pm);
  mpz_mul(tmpm,tmpm,pm);
  mpz_mul_ui(tmpm,tmpm,2);  

  deg=p-3;
  gap=mpz_sizeinbase(tmpm,16)*4;
 
/**********  im=t/(e^{ft}-1)   **********/
  mpz_set_ui(im,0);
  mp_expmodp(f,p,deg+1,gap,ct1,ct2,ct3,&im);
  mpz_tdiv_q_2exp(im,im,gap);
  mp_polymodp(im,p,deg,gap,ct1,ct2,ct3,&im1);
  mpz_set(im,im1);
  mp_polyinvp(im,p,deg,gap,ct1,ct2,ct3,&im1);
  mpz_set(im,im1);
/********** im1=sum kro(f,a) e(ax) *********/
  mpz_set_ui(fm,f);
  mpz_set_ui(im1,0);
  if (f==1){mpz_set_ui(im1,1);}
  else
  {
   for (a=1;a<f;a++)
     {
     kro=mpz_kronecker_ui(fm,a);
     if (kro==-1){kro=p-1;} 
     if (kro!=0)
       {
       mp_expmodp(a,p,deg+1,gap,ct1,ct2,ct3,&im2);
       mpz_mul_ui(im2,im2,kro);
       mpz_add(im1,im1,im2);
       }
     if (f>2*p)
       {
       mp_polymodp(im1,p,deg,gap,ct1,ct2,ct3,&im2);
       mpz_set(im1,im2);
       }
     }
     mp_polymodp(im1,p,deg,gap,ct1,ct2,ct3,&im2);
     mpz_set(im1,im2);
  }

  mpz_mul(im,im,im1);
  mp_polymodp(im,p,deg,gap,ct1,ct2,ct3,&im1);
  mpz_set(im,im1);
  mp_polymodp0f(im,p,deg,gap,f);

  mpz_clear(im);
  mpz_clear(im1);
  mpz_clear(im2);
  mpz_clear(pm);
  mpz_clear(fm);
  mpz_clear(gcm);
  mpz_clear(tmpm);
 }
 
/*     *jm=exp(ax) mod p,  [deg,  gap]  */

mp_expmodp(int a, int p, int deg, int gap, int ct1, int ct2, int ct3, mpz_t *jm)
  { 
    int j,m,n,gap4=gap/4,ctnum,cttnum,ctttnum;
    mpz_t pm,am,rm,sm,tm;
    char temp1[(deg+1)*gap4+1]; 
    char tmp[100];
    char ttmp[ct1+gap4+1];
    char tttmp[ct2+ct1+1];
    char ttttmp[ct3+ct2+1]; 

    mpz_init(pm);
    mpz_init(am);
    mpz_init(rm);
    mpz_init(sm);
    mpz_init(tm);
    mpz_set_ui(pm,p);
    mpz_set_ui(am,a);
    mpz_invert(am,am,pm); 
    mpz_set_ui(rm,1);

    strcpy(temp1,"");
    strcpy(ttmp,"");
    strcpy(tttmp,"");
    strcpy(ttttmp,"");
    mpz_get_str(tmp,16,rm);
    strcat(ttmp,tmp);

    ctnum=0;
    cttnum=0;
    ctttnum=0;
    j=1;
    while(j<deg+1)
    {
      mpz_set_ui(sm,p-1-j);
      mpz_mul(rm,rm,sm);
      mpz_mul(rm,rm,am); 
      mpz_mod(rm,rm,pm);
      n=gap4-mpz_sizeinbase(rm,16);
      for (m=1;m<n+1;m++)
         {strcat(ttmp,"0");}
      mpz_get_str(tmp,16,rm);
      strcat(ttmp,tmp);
      j++;
      ctnum=ctnum+gap4;
      if (ctnum>ct1)
      {
        strcat(tttmp,ttmp);
	cttnum=cttnum+ctnum;
        strcpy(ttmp,"");
	ctnum=0;
        if (cttnum>ct2)
	{
          strcat(ttttmp,tttmp);
          ctttnum=ctttnum+cttnum;
	  strcpy(tttmp,"");
	  cttnum=0;
	  if (ctttnum>ct3)
          {
            strcat(temp1,ttttmp);
            strcpy(ttttmp,"");
            ctttnum=0;
          }
        }
      }

     }
     strcat(tttmp,ttmp);
     strcat(ttttmp,tttmp);
     strcat(temp1,ttttmp);
     mpz_clear(pm);
     mpz_clear(am);
     mpz_clear(rm);
     mpz_clear(sm);
     mpz_clear(tm);

     mpz_set_str(*jm,temp1,16);


 }


/* jm=f(x) -> *jm1=f^{-1}(x) mod p +O(x^{deg+1}) mod p,  [deg,  gap]  */

mp_polyinvp(mpz_t jm, int p, int deg, int gap, int ct1, int ct2, int ct3, mpz_t *jm1)
  {
    int i,j;
    mpz_t pm,rm,sm;

    mpz_init(pm);
    mpz_init(rm);
    mpz_init(sm);
    mpz_set_ui(pm,p);

    i=1;	    
    mpz_tdiv_r_2exp(*jm1,jm,gap);
    mpz_invert(*jm1,*jm1,pm);
    while(i<deg+2)
      {
         mpz_tdiv_r_2exp(sm,jm,2*i*gap);
         mpz_mul(rm,*jm1,sm);
         mpz_sub_ui(rm,rm,1);
	     mp_polymodp(rm,p,2*i-1,gap,ct1,ct2,ct3,&sm);
         mpz_mul(rm,*jm1,sm);
	     mp_polymodp(rm,p,2*i-1,gap,ct1,ct2,ct3,&sm);
         mpz_mul_ui(rm,sm,p-1);
	     mp_polymodp(rm,p,2*i-1,gap,ct1,ct2,ct3,&sm);
         mpz_add(*jm1,*jm1,sm);
	     mp_polymodp(*jm1,p,2*i-1,gap,ct1,ct2,ct3,&sm);
         mpz_set(*jm1,sm);
         i=i*2;
      }
    mpz_tdiv_r_2exp(*jm1,*jm1,(deg+1)*gap);
    
    mpz_clear(pm);
    mpz_clear(rm);
    mpz_clear(sm);

  }


/* *jm1= *im f(x)+O(x^deg+1) -> f(x) mod p,  [deg,  gap]  */

mp_polymodp(mpz_t am, int p, int deg, int gap, int ct1, int ct2, int ct3, mpz_t *am1)
{
  int  n,m,nn,gap4=gap/4,ctnum,cttnum,ctttnum;
  long long int x,y[100];
  char temp1[(deg+1)*gap4+1];
  char temp2[(deg+1)*gap4+1];
  char tmp[100];
  char ttmp[ct1+gap4+1];
  char tttmp[ct2+ct1+1];
  char ttttmp[ct3+ct2+1];
  mpz_t xm;

  strcpy(temp1,"");
  strcpy(temp2,"");

  mpz_init(xm);
  y[1]=1;
  for (m=2;m<gap4+1;m++)
    {y[m]=(y[m-1]*16)%p;}
 
  mpz_tdiv_r_2exp(am,am,(deg+1)*gap);
  mpz_get_str(temp2,16,am);
  nn=gap4*(deg+1)-mpz_sizeinbase(am,16);
  for (m=1;m<nn+1;m++)
    {strcat(temp1,"0");}
  strcat(temp1,temp2);

  strcpy(temp2,"");
  strcpy(ttmp,"");
  strcpy(tttmp,"");
  strcpy(ttttmp,"");

  ctnum=0;
  cttnum=0;
  ctttnum=0;
  for (n=0;n<deg+1;n++)
    {
        x=0;
	for (m=1;m<gap4+1;m++){
	  if (temp1[gap4*n+m-1]=='1'){x=x+y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='2'){x=x+2*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='3'){x=x+3*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='4'){x=x+4*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='5'){x=x+5*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='6'){x=x+6*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='7'){x=x+7*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='8'){x=x+8*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='9'){x=x+9*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='a'){x=x+10*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='b'){x=x+11*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='c'){x=x+12*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='d'){x=x+13*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='e'){x=x+14*y[gap4-m+1];}
	  if (temp1[gap4*n+m-1]=='f'){x=x+15*y[gap4-m+1];}
	x=x%p;
	}

        mpz_set_ui(xm,x);
	nn=gap4-mpz_sizeinbase(xm,16);
	for (m=1;m<nn+1;m++)
	  {strcat(ttmp,"0");}
	mpz_get_str(tmp,16,xm);
	strcat(ttmp,tmp);
        ctnum=ctnum+gap4;
	if (ctnum>ct1)
        {
          strcat(tttmp,ttmp);
	  cttnum=cttnum+ctnum;
	  strcpy(ttmp,"");
	  ctnum=0;
	  if (cttnum>ct2)
	  {
            strcat(ttttmp,tttmp);
	    ctttnum=ctttnum+cttnum;           
	    strcpy(tttmp,"");
	    cttnum=0;
	    if (ctttnum>ct3)
            {
              strcat(temp2,ttttmp);
              strcpy(ttttmp,"");
              ctttnum=0;
            }
          }
        }
    }

    strcat(tttmp,ttmp);
    strcat(ttttmp,tttmp);
    strcat(temp2,ttttmp);
    mpz_set_str(*am1,temp2,16);

    mpz_clear(xm);
}







/* am=\sum a_i x^i -> i (a_i mod p=0),  [deg, gap] */  





mp_polymodp0f(mpz_t am, int p, int deg, int gap, int f)
  { 

    int  n,m,l,nn,gap4=gap/4;
    long long int x,y[100];
    char temp1[(deg+1)*gap4+1];
    char temp2[(deg+1)*gap4+1];
    FILE *file;

    strcpy(temp1,"");
    strcpy(temp2,"");

    file=fopen("tmpbernc.txt","w");
    y[1]=1;
    for (m=2;m<gap4+1;m++)
      {y[m]=(y[m-1]*16)%p;}

    mpz_tdiv_r_2exp(am,am,(deg+1)*gap);
    mpz_get_str(temp2,16,am);
    nn=gap4*(deg+1)-mpz_sizeinbase(am,16);
    for (m=1;m<nn+1;m++)
      {strcat(temp1,"0");}
    strcat(temp1,temp2);

    for (n=0;n<deg+1;n++)
      {
        x=0;
        for (m=1;m<gap4+1;m++){
          if (temp1[gap4*n+m-1]=='1'){x=x+y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='2'){x=x+2*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='3'){x=x+3*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='4'){x=x+4*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='5'){x=x+5*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='6'){x=x+6*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='7'){x=x+7*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='8'){x=x+8*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='9'){x=x+9*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='a'){x=x+10*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='b'){x=x+11*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='c'){x=x+12*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='d'){x=x+13*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='e'){x=x+14*y[gap4-m+1];}
          if (temp1[gap4*n+m-1]=='f'){x=x+15*y[gap4-m+1];}
          x=x%p;
        }
    l=deg-n;

    if ((x%p==0)&&(l%2==0)&&(l!=0))
      {
        fprintf(file,"%d\n",p);
        fprintf(file,"%d\n",l);
        fprintf(file,"%d\n",f);
      } 
   }
   fclose(file);

}


mp_padic(int p, int k, int f, mpz_t *alm, int *i1)
{

  mpz_t pm, fm, am, f0m, p2m, p3m, fapm, c0m, c1m, ram, qam;
  mpz_t tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
  int kro;

  mpz_init(pm);
  mpz_init(fm);
  mpz_init(am);
  mpz_init(f0m);
  mpz_init(p2m);
  mpz_init(p3m);
  mpz_init(fapm);
  mpz_init(c0m);
  mpz_init(c1m);
  mpz_init(ram);
  mpz_init(qam);
  mpz_init(tmp1);
  mpz_init(tmp2);
  mpz_init(tmp3);
  mpz_init(tmp4);
  mpz_init(tmp5);
  mpz_init(tmp6);

  mpz_set_ui(pm,p);
  mpz_set_ui(fm,f);

  if (f%p==0){mpz_set_ui(f0m,f);} 
  else {mpz_set_ui(f0m,f*p);}
  mpz_mul(p2m,pm,pm);
  mpz_mul(p3m,p2m,pm);

  mpz_set_ui(c0m,0);
  mpz_set_ui(c1m,0);

  mpz_set_ui(am,1);
  while(mpz_cmp(am,f0m)<0)
    {
      mpz_gcd(tmp1,f0m,am);
      if (mpz_cmp_ui(tmp1,1)==0)
	{
	  mpz_invert(tmp1,am,pm);
	  mp_omega(tmp1,pm,3,&fapm);
	  mpz_mul(tmp1,am,fapm);
	  mpz_sub_ui(tmp1,tmp1,1);
	  mpz_divexact(ram,tmp1,pm);
	  mpz_invert(tmp1,am,p3m);
	  mpz_mul(tmp1,tmp1,f0m);
	  mpz_divexact_ui(qam,tmp1,p);
	  mpz_mod(ram,ram,p2m);
	  mpz_neg(qam,qam);
	  mpz_mod(qam,qam,p2m);
	  kro=mpz_kronecker(fm,am);
	  mpz_mul_ui(tmp1,p2m,p-1);
	  mpz_sub_ui(tmp1,tmp1,k);
	  mpz_powm(tmp1,fapm,tmp1,p3m);
	  mpz_mul_si(tmp1,tmp1,kro);

	  mpz_mul(tmp2,ram,qam);
	  mpz_mul(tmp3,ram,ram);
	  mpz_mul(tmp4,qam,qam);

	  mpz_mul_ui(tmp5,ram,2);
	  mpz_add(tmp5,tmp5,qam);
	  mpz_addmul_ui(tmp5,tmp2,p);

	  mpz_mul_ui(tmp6,tmp3,12);
	  mpz_add(tmp6,tmp6,tmp4);
	  mpz_submul_ui(tmp6,tmp2,12);
	  mpz_addmul(c0m,tmp1,tmp5);
	  mpz_mod(c0m,c0m,p2m);
	  mpz_addmul(c1m,tmp1,tmp6);
	  mpz_mod(c1m,c1m,p2m);
	}
      mpz_add_ui(am,am,1);
    }

  mpz_set_ui(tmp1,2);
  mpz_invert(tmp1,tmp1,p2m);
  mpz_mul(c0m,c0m,tmp1);
  mpz_mod(c0m,c0m,p2m);

  mpz_set_si(tmp1,-24);
  mpz_invert(tmp1,tmp1,p2m);
  mpz_mul(c1m,c1m,tmp1);
  mpz_mod(c1m,c1m,p2m);

  mpz_mul_ui(*alm,c1m,p);
  mpz_add(*alm,*alm,c0m);
  mpz_invert(tmp1,c1m,pm);
  mpz_mul(*alm,*alm,tmp1);
  mpz_divexact_ui(tmp1,f0m,p);
  mpz_mul(*alm,*alm,tmp1);
  mpz_mod(*alm,*alm,p2m);

  if (mpz_divisible_p(c0m,pm)==0)
    {*i1=0;}
  else
  {
    *i1=1;
    mpz_divexact(tmp1,f0m,pm);
    mpz_invert(tmp1,tmp1,pm);
    mpz_mul(tmp1,*alm,tmp1);
    mpz_mod(tmp1,tmp1,p2m);
    mpz_divexact(tmp1,tmp1,pm);
    if (mpz_cmp_ui(tmp1,0)==0)
      {*i1=*i1*2;}
    if (mpz_cmp_ui(tmp1,1)==0)
      {*i1=*i1*3;}
    if (mpz_cmp_ui(tmp1,k+1)==0)
      {*i1=*i1*5;}
    if (mpz_cmp_ui(tmp1,k)==0)
      {*i1=*i1*6;}      
    if(mpz_divisible_p(c1m,pm)!=0)
      {*i1=*i1*4;}
    if((f%p==0)&&((2*k)%(p-1)==0))
      {*i1=0;}
  }    

  mpz_clear(pm);
  mpz_clear(fm);
  mpz_clear(am);
  mpz_clear(f0m);
  mpz_clear(p2m);
  mpz_clear(p3m);
  mpz_clear(fapm);
  mpz_clear(c0m);
  mpz_clear(c1m);
  mpz_clear(ram);
  mpz_clear(qam);
  mpz_clear(tmp1);
  mpz_clear(tmp2);
  mpz_clear(tmp3);
  mpz_clear(tmp4);
  mpz_clear(tmp5);
  mpz_clear(tmp6);

}


int mp_hantei(int p, int k, int f, mpz_t lm, mpz_t alm, int n)
{
  mpz_t pm,fm,am,jm,gm,ym,yyym,lm1,lm2;
  mpz_t hm,hhm,cm,tmp0,tmp1,tmp2,tmp3,tmp4,tmp5;
  mpz_t pnm,pn1m,fpnm,kkm,aam,km,ccm,lzm,klastm;
  int i,g,w,f0,ff,rt,kro,kt;
  FILE *file; 

  mpz_init(pm);
  mpz_init(fm);
  mpz_init(gm);
  mpz_init(am);
  mpz_init(jm);
  mpz_init(ym);
  mpz_init(yyym);
  mpz_init(lm1);
  mpz_init(lm2);
  mpz_init(hm);
  mpz_init(hhm);
  mpz_init(cm);
  mpz_init(tmp0);
  mpz_init(tmp1);
  mpz_init(tmp2);
  mpz_init(tmp3);
  mpz_init(tmp4);
  mpz_init(tmp5);
  mpz_init(pnm);
  mpz_init(pn1m);
  mpz_init(fpnm);
  mpz_init(kkm);
  mpz_init(aam);
  mpz_init(km);
  mpz_init(ccm);
  mpz_init(lzm);
  mpz_init(klastm);

  mpz_set_ui(pm,p);
  mpz_set_ui(fm,f);

  if (f%p==0){f0=f;} 
  else {f0=f*p;}
  ff=f0/p;

  mpz_set_ui(pnm,1);
  for (i=1;i<n+1;i++)
    {mpz_mul(pnm,pnm,pm);}
  mpz_mul(pn1m,pnm,pm);
  mpz_mul_ui(fpnm,pnm,f0);
  mpz_sub_ui(tmp1,lm,1);
  mpz_divexact_ui(tmp1,tmp1,f0);
  mpz_divexact(lzm,tmp1,pnm);

  g=prmroot(lm,p);
  mpz_set_ui(gm,g);
  mpz_powm(ym,gm,lzm,lm);
  mpz_sub_ui(lm1,lm,1);
  mpz_sub_ui(lm2,lm,2);

  mpz_set_ui(cm,1);
 
  kt=0;
  mpz_set_ui(km,0);
  mpz_set(klastm,pnm);

  while(mpz_cmp(km,klastm)<0)
  {
    mpz_set_ui(ccm,1);
    mpz_set_ui(tmp1,1+f0);
    mpz_powm(kkm,tmp1,km,fpnm);
    mpz_add_ui(tmp1,alm,1);
    mpz_sub(tmp2,pnm,km);
    mpz_sub_ui(tmp2,tmp2,1);
    mpz_powm(aam,tmp1,tmp2,lm1);
    mpz_set_ui(am,1);
    while(mpz_cmp(am,pm)<0)
    {
      w=0;
      mpz_invert(tmp0,am,pm);
      mp_omega(tmp0,pm,n+3,&tmp1);
      mpz_powm_ui(tmp4,tmp1,p-1-k,pn1m);
      mpz_mul(tmp5,tmp4,lm2);
      mpz_mod(tmp5,tmp5,pn1m);
      while(w<ff)
      {
        mpz_set(hm,tmp1);
        mpz_addmul_ui(hm,pn1m,w);
        mpz_mul(hhm,hm,kkm);
        mpz_mod(hhm,hhm,fpnm);
        kro=mpz_kronecker(fm,hm);
        if (kro!=0)
        {
          mpz_powm(yyym,ym,hhm,lm);
          mpz_set_ui(tmp3,1);
          mpz_sub(yyym,tmp3,yyym);
          if (kro==1)
          {
            mpz_powm(yyym,yyym,tmp4,lm);
            mpz_mul(ccm,ccm,yyym);
            mpz_mod(ccm,ccm,lm);
          }
          else
          {
            mpz_powm(yyym,yyym,tmp5,lm);
            mpz_mul(ccm,ccm,yyym);
            mpz_mod(ccm,ccm,lm);
          }
        }
        w++;
      }
      mpz_add_ui(am,am,1);
    }
    mpz_powm(tmp1,ccm,aam,lm);
    mpz_mul(cm,cm,tmp1);
    mpz_mod(cm,cm,lm);
 /*
    file=fopen("c-result.txt","a+");
    mpz_out_str(file,10,lm);
    fprintf(file,",%d\n",kt);
    mpz_out_str(file,10,cm);
    fprintf(file,",%d\n",kt);
    fclose(file);
 */
    mpz_add_ui(km,km,1);
    kt++;

  }

  mpz_mul_ui(tmp1,lzm,f0);
  mpz_divexact(tmp1,tmp1,pm);
  mpz_powm(cm,cm,tmp1,lm);
  if (mpz_cmp_ui(cm,1)!=0)
    {rt=1;} 
    else 
    {rt=-1;}

  mpz_clear(fm);
  mpz_clear(gm);
  mpz_clear(am);
  mpz_clear(jm);
  mpz_clear(ym);
  mpz_clear(yyym);
  mpz_clear(lm1);
  mpz_clear(lm2);
  mpz_clear(hm);
  mpz_clear(hhm);
  mpz_clear(cm);
  mpz_clear(tmp0);
  mpz_clear(tmp1);
  mpz_clear(tmp2);
  mpz_clear(tmp3);
  mpz_clear(tmp4);
  mpz_clear(tmp5);
  mpz_clear(pnm);
  mpz_clear(pn1m);
  mpz_clear(fpnm);
  mpz_clear(kkm);
  mpz_clear(aam);
  mpz_clear(km);
  mpz_clear(ccm);
  mpz_clear(lzm);
  mpz_clear(klastm);

  return(rt);
  }



mp_omega(mpz_t am, mpz_t pm, int k, mpz_t *jm)
 {
  int i;
  mpz_t tmp,tmp1;

  mpz_init(tmp);
  mpz_init(tmp1);

  mpz_mod(*jm,am,pm);
  mpz_set_ui(tmp,1);
  for(i=1;i<k;i++)
    {
      mpz_mul(tmp,tmp,pm);
    }
  mpz_mul(tmp1,tmp,pm);
  mpz_powm(*jm,*jm,tmp,tmp1);

  mpz_clear(tmp);
  mpz_clear(tmp1);
 }



int prmroot(mpz_t lm, int p)   /*l-1 is divided by p*/
{
  mpz_t gm,hm,ltm,tmpm;
  int i,j,k,g,gg,p1,m1,hantei,pd[100];

  mpz_init(ltm);
  mpz_init(gm);
  mpz_init(hm);
  mpz_init(tmpm);
  
  pd[1]=p;
  mpz_sub_ui(ltm,lm,1);
  while(mpz_divisible_ui_p(ltm,p)!=0)
    {
      mpz_divexact_ui(ltm,ltm,p);
    }
  k=2;
  while(mpz_cmp_ui(ltm,1)>0)
    {
      prmdiv(ltm,&p1,&m1);
      pd[k]=p1;
      for(j=1;j<m1+1;j++)
	{mpz_divexact_ui(ltm,ltm,p1);}
      k++;
    }
  k=k-1; 

  mpz_set_ui(gm,2);
  gg=2;
  g=1;
  while(g==1)
    {
      hantei=0;
      mpz_sub_ui(ltm,lm,1);
      for(i=1;i<k+1;i++)
	{ 
          mpz_divexact_ui(tmpm,ltm,pd[i]);
	  mpz_powm(hm,gm,tmpm,lm);
	  if (mpz_cmp_ui(hm,1)==0){hantei=-1;}
	} 
      if (hantei==0){g=gg;}
      mpz_add_ui(gm,gm,1);
      gg++;
    }

  mpz_clear(ltm);
  mpz_clear(gm);
  mpz_clear(hm);
  mpz_clear(tmpm);
  return(g);
}



prmdiv(mpz_t lm, int *p, int *m)
  {
  int i;
  mpz_t nm;

  mpz_init(nm);
  mpz_set(nm,lm);
 
  if (mpz_probab_prime_p(nm,20)>0)
    {
      *p=mpz_get_ui(nm);
      *m=1;
    }
  else
    {
      i=2;
      *p=1;
      while(*p==1)
	{
	  if(mpz_divisible_ui_p(nm,i)!=0)
	    {*p=i;}
	  i++;
	}

      *m=1;
      mpz_divexact_ui(nm,nm,*p);
      while(mpz_divisible_ui_p(nm,*p)!=0)
	{
	 (*m)++;
	  mpz_divexact_ui(nm,nm,*p);
	}
    }
  mpz_clear(nm);

  } 


/*  p^(n+1) <tmp^3  */

int mp_hantei1(int p, int k, int f, mpz_t lm, mpz_t alm, int n)
{
  mpz_t pm,fm,am,jm,gm,ym,yyym,lm1,lm2;
  mpz_t hm,hhm,cm,tmp0,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
  mpz_t pnm,pn1m,fpnm,kkm,aam,km,ccm,lzm,klastm;
  unsigned long long int i,g,w,f0,ff,rt,kro,kt,tmp,i1,i2,i3,j1,j2,j3,l,tp;
  mpz_t cs1[5000],cs2[5000],cs3[5000];  
  FILE *file; 

  tmp=4000;

  l=mpz_get_ui(lm);

  mpz_init(pm);
  mpz_init(fm);
  mpz_init(gm);
  mpz_init(am);
  mpz_init(jm);
  mpz_init(ym);
  mpz_init(yyym);
  mpz_init(lm1);
  mpz_init(lm2);
  mpz_init(hm);
  mpz_init(hhm);
  mpz_init(cm);
  mpz_init(tmp0);
  mpz_init(tmp1);
  mpz_init(tmp2);
  mpz_init(tmp3);
  mpz_init(tmp4);
  mpz_init(tmp5);
  mpz_init(tmp6);
  mpz_init(pnm);
  mpz_init(pn1m);
  mpz_init(fpnm);
  mpz_init(kkm);
  mpz_init(aam);
  mpz_init(km);
  mpz_init(ccm);
  mpz_init(lzm);
  mpz_init(klastm);
  for(i=0;i<tmp;i++)
    {
      mpz_init(cs1[i]);
      mpz_init(cs2[i]);
      mpz_init(cs3[i]);
    } 


  mpz_set_ui(pm,p);
  mpz_set_ui(fm,f);

  if (f%p==0){f0=f;} 
  else {f0=f*p;}
  ff=f0/p;

  mpz_set_ui(pnm,1);
  for (i=1;i<n+1;i++)
    {mpz_mul(pnm,pnm,pm);}
  mpz_mul(pn1m,pnm,pm);
  mpz_mul_ui(fpnm,pnm,f0);
  mpz_sub_ui(tmp1,lm,1);
  mpz_divexact_ui(tmp1,tmp1,f0);
  mpz_divexact(lzm,tmp1,pnm);

  g=prmroot(lm,p);
  mpz_set_ui(gm,g);
  mpz_powm(ym,gm,lzm,lm);
  mpz_sub_ui(lm1,lm,1);
  mpz_sub_ui(lm2,lm,2);

  mpz_set_ui(cm,1);
 
  kt=0;
  mpz_set_ui(km,0);
  mpz_set(klastm,pnm);
  
/*
  if (n==1)
   {
    mpz_set_ui(cm,0);
    mpz_mul_ui(cm,cm,1000000);
    mpz_add_ui(cm,cm,1);
    kt=0;
    mpz_set_ui(km,kt);
   }       
*/

  while(mpz_cmp(km,klastm)<0)
  {
    mpz_set_ui(ccm,1);
    mpz_set_ui(tmp1,1+f0);
    mpz_powm(kkm,tmp1,km,fpnm);
    mpz_add_ui(tmp1,alm,1);
    mpz_sub(tmp2,pnm,km);
    mpz_sub_ui(tmp2,tmp2,1);
    mpz_powm(aam,tmp1,tmp2,lm1);
    mpz_set_ui(am,1);
    for(i=0;i<tmp;i++)
     {
      mpz_set_ui(cs1[i],1);
      mpz_set_ui(cs2[i],1);
      mpz_set_ui(cs3[i],1);
     }
    while(mpz_cmp(am,pm)<0)
    {
      w=0;
      mpz_invert(tmp0,am,pm);
      mp_omega(tmp0,pm,n+3,&tmp1);
      mpz_powm_ui(tmp4,tmp1,p-1-k,pn1m);
      
/**/
      mpz_set(tmp6,tmp4);
      mpz_mod_ui(tmp5,tmp6,tmp);
      i1=mpz_get_ui(tmp5);
      mpz_div_ui(tmp6,tmp6,tmp);
      mpz_mod_ui(tmp5,tmp6,tmp);
      i2=mpz_get_ui(tmp5);
      mpz_div_ui(tmp6,tmp6,tmp);
      mpz_mod_ui(tmp5,tmp6,tmp);
      i3=mpz_get_ui(tmp5);

      mpz_mul(tmp5,tmp4,lm2);
      mpz_mod(tmp5,tmp5,pn1m);

      mpz_set(tmp6,tmp5);
      mpz_mod_ui(tmp5,tmp6,tmp);
      j1=mpz_get_ui(tmp5);
      mpz_div_ui(tmp6,tmp6,tmp);
      mpz_mod_ui(tmp5,tmp6,tmp);
      j2=mpz_get_ui(tmp5);
      mpz_div_ui(tmp6,tmp6,tmp);
      mpz_mod_ui(tmp5,tmp6,tmp);
      j3=mpz_get_ui(tmp5);
      
/**/

      while(w<ff)
      {
        mpz_set(hm,tmp1);
        mpz_addmul_ui(hm,pn1m,w);
        mpz_mul(hhm,hm,kkm);
        mpz_mod(hhm,hhm,fpnm);
        kro=mpz_kronecker(fm,hm);
        if (kro!=0)
        {
          mpz_powm(yyym,ym,hhm,lm);
          mpz_set_ui(tmp3,1);
          mpz_sub(yyym,tmp3,yyym);
          if (kro==1)
          {
           mpz_mul(cs1[i1],cs1[i1],yyym);
           mpz_mod(cs1[i1],cs1[i1],lm);

           mpz_mul(cs2[i2],cs2[i2],yyym);
           mpz_mod(cs2[i2],cs2[i2],lm);

           mpz_mul(cs3[i3],cs3[i3],yyym);
           mpz_mod(cs3[i3],cs3[i3],lm);
         }
          else
          {
	    mpz_mul(cs1[j1],cs1[j1],yyym);
	    mpz_mod(cs1[j1],cs1[j1],lm);

	    mpz_mul(cs2[j2],cs2[j2],yyym);
	    mpz_mod(cs2[j2],cs2[j2],lm);

	    mpz_mul(cs3[j3],cs3[j3],yyym);
	    mpz_mod(cs3[j3],cs3[j3],lm);   
          }
       }
        w++;
      }
      mpz_add_ui(am,am,1);
    }
    mpz_set_ui(ccm,1);
    for(i=1;i<tmp;i++)
    {
/**/

      mpz_powm_ui(tmp6,cs1[i],i,lm);
      mpz_mul(ccm,ccm,tmp6);
      mpz_powm_ui(tmp6,cs2[i],i,lm);
      mpz_powm_ui(tmp6,tmp6,tmp,lm);
      mpz_mul(ccm,ccm,tmp6);
      mpz_powm_ui(tmp6,cs3[i],i,lm);
      mpz_powm_ui(tmp6,tmp6,tmp,lm);
      mpz_powm_ui(tmp6,tmp6,tmp,lm);
      mpz_mul(ccm,ccm,tmp6);

/**/

    }
    mpz_powm(tmp6,ccm,aam,lm);
    mpz_mul(cm,cm,tmp6);
    mpz_mod(cm,cm,lm);
/*
    file=fopen("c-result.txt","a+");
    mpz_out_str(file,10,lm);
    fprintf(file,",%d\n",kt);
    mpz_out_str(file,10,cm);
    fprintf(file,",%d\n",kt);
    fclose(file);
*/
    mpz_add_ui(km,km,1);
    kt++;

  }

  mpz_mul_ui(tmp1,lzm,f0);
  mpz_divexact(tmp1,tmp1,pm);
  mpz_powm(cm,cm,tmp1,lm);
  if (mpz_cmp_ui(cm,1)!=0)
    {rt=1;} 
    else 
    {rt=-1;}

  mpz_clear(fm);
  mpz_clear(gm);
  mpz_clear(am);
  mpz_clear(jm);
  mpz_clear(ym);
  mpz_clear(yyym);
  mpz_clear(lm1);
  mpz_clear(lm2);
  mpz_clear(hm);
  mpz_clear(hhm);
  mpz_clear(cm);
  mpz_clear(tmp0);
  mpz_clear(tmp1);
  mpz_clear(tmp2);
  mpz_clear(tmp3);
  mpz_clear(tmp4);
  mpz_clear(tmp5);
  mpz_clear(tmp6);
  mpz_clear(pnm);
  mpz_clear(pn1m);
  mpz_clear(fpnm);
  mpz_clear(kkm);
  mpz_clear(aam);
  mpz_clear(km);
  mpz_clear(ccm);
  mpz_clear(lzm);
  mpz_clear(klastm);
  for(i=0;i<tmp;i++)
   {
      mpz_clear(cs1[i]);
      mpz_clear(cs2[i]);
      mpz_clear(cs3[i]);
    }


  return(rt);
  }