/* http://www.mis.hiroshima-u.ac.jp/~hiroki/major/galois1-e.html */
/* by Hiroki Sumida-Takahashi 2020, March 10 */


/*   p<160,000,000  */

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

void mp_bernc(long long int, long long int, long long int);
void mp_expmodp(long long int, long long int, long long int, long long int, long long int, long long int, long long int, mpz_t *, long long int);
void mp_polyinvp(mpz_t, long long int, long long int, long long int, long long int, long long int, long long int, mpz_t *, long long int);
void mp_mul(mpz_t *, mpz_t, mpz_t);
void mp_polymodp(mpz_t, long long int, long long int, long long int, long long int, long long int, long long int, mpz_t *, long long int);
void mp_polymodp0f(mpz_t, long long int, long long int, long long int, long long int,long long int);
void mp_padic(long long int, long long int, long long int, mpz_t *, long long int *);
long long int mp_hantei(long long int, long long int, long long int, mpz_t, mpz_t, long long int);
void mp_omega(mpz_t, mpz_t, long long int, mpz_t *);
long long int prmroot(mpz_t, long long int);
void prmdiv(mpz_t, long long int *, long long int *);
long long int mp_hantei1(long long int, long long int, long long int, mpz_t, mpz_t, long long int);

void main()
 {
  long long int p,f,fn,pn1,pn2,k,i1,i2,i3,l,lz,lh,hnum,tmp,base=32;
  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=710000;
  pn2=715000;
  for(fn=0;fn<60;fn++)
    {
    if (fn==-1){f=1;}
    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=711463;
      }
   
      mpz_set_ui(pm,p);
      if (mpz_probab_prime_p(pm,20)>0)
        {
	mp_bernc(p,f,base);
	file0=fopen("tmpbernc.txt","r");
	while(fscanf(file0,"%lld\n %lld\n %lld\n",&p,&k,&f)!=EOF)
	  {
	    mp_padic(p,k,f,&alm,&i1);
            if (i1>1){
	      file1=fopen("bcn1.txt","a+");
	      fprintf(file1,"%lld,%lld,%lld,%lld\n",p,k,f,i1);
	      fclose(file1);
	    }
	    if (i1==0){
	      file2=fopen("bcn2.txt","a+");
	      fprintf(file2,"%lld,%lld,%lld,%lld\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,"%lld,%lld,%lld,%lld,%lld\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);
   	          /*  i2=mp_hantei(p,k,f,lm,alm,0);*/
                    hnum++;
		  }
  		if(i2!=1){
		  file4=fopen("bcn1.txt","a+");
		  fprintf(file4,"%lld,%lld,%lld,%lld\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,"%lld,%lld,%lld,1,%lld\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 */



void mp_bernc(long long int p, long long int f, long long int base)
 {
   long long int fn1,i,gap,deg,a,kro,fn,ct1=400,ct2=20000,ct3=800000,x;
  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,base)*5;
 
/**********  im=t/(e^{ft}-1)   **********/

  mpz_set_ui(im,0);
  mp_expmodp(f,p,deg+1,gap,ct1,ct2,ct3,&im,base);
  mpz_tdiv_q_2exp(im,im,gap);
  mp_polymodp(im,p,deg,gap,ct1,ct2,ct3,&im1,base);
  mpz_set(im,im1);
  mp_polyinvp(im,p,deg,gap,ct1,ct2,ct3,&im1,base);
  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,base);
       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,base);
       mpz_set(im1,im2);
       }
     }
   mp_polymodp(im1,p,deg,gap,ct1,ct2,ct3,&im2,base);
     mpz_set(im1,im2);
  }
  mp_mul(&im,im,im1);
  mp_polymodp(im,p,deg,gap,ct1,ct2,ct3,&im1,base);
  mpz_set(im,im1);
  mp_polymodp0f(im,p,deg,gap,f,base);

  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]  */

void mp_expmodp(long long int a, long long int p, long long int deg, long long int gap, long long int ct1, long long int ct2, long long int ct3, mpz_t *jm, long long int base)
  { 
    long long int j,m,n,gap5=gap/5,ctnum,cttnum,ctttnum,bb,cc,dgg=5000000;
    mpz_t pm,am,rm,sm,tm;
    char temp1[dgg*gap5+1]; 
    char tmp[100];
    char ttmp[ct1+gap5+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); 
 
    ctnum=0;
    cttnum=0;
    ctttnum=0;
    
    strcpy(temp1,"");
    strcpy(tmp,"");
    strcpy(ttmp,"");
    strcpy(tttmp,"");
    strcpy(ttttmp,"");

    mpz_set_ui(*jm,0);
    mpz_set_ui(rm,1);
    mpz_get_str(tmp,base,rm);
    strcat(ttmp,tmp);

    j=1;
    while(deg>0)
      {  
      if (deg>=dgg)
          {bb=dgg;}   
      if (deg<dgg)
          {bb=deg;}
      deg=deg-bb;
      cc=0;
      while(cc<bb)
      {
        mpz_set_ui(sm,p-1-j);
        mpz_mul(rm,rm,sm);
        mpz_mul(rm,rm,am); 
        mpz_mod(rm,rm,pm);
        n=gap5-mpz_sizeinbase(rm,base);
        for (m=1;m<n+1;m++)
           {strcat(ttmp,"0");}
        mpz_get_str(tmp,base,rm);
        strcat(ttmp,tmp);
        j++;
        cc++;
        ctnum=ctnum+gap5;
        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_set_str(sm,temp1,base);
      mpz_mul_2exp(*jm,*jm,bb*gap);
      mpz_add(*jm,*jm,sm);

      strcpy(temp1,"");

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

      ctnum=0;
      cttnum=0;
      ctttnum=0;

   }

   mpz_clear(pm);
   mpz_clear(am);
   mpz_clear(rm);
   mpz_clear(sm);
   mpz_clear(tm);

 }


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

void mp_polyinvp(mpz_t jm, long long int p, long long int deg, long long int gap, long long int ct1, long long int ct2, long long int ct3, mpz_t *jm1,long long int base)
  {
    long long int i,j,k,ii,gap5=gap/5,x;
    mpz_t pm,ppm,rm,sm;

    mpz_init(pm);
    mpz_init(ppm);
    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_set(rm,jm);
         if (2*i<deg){ii=2*i;}
         if (2*i>=deg){ii=deg;}
         mpz_set_ui(sm,0);
         mpz_tdiv_r_2exp(sm,rm,ii*gap);
         mp_mul(&rm,*jm1,sm);
         mpz_sub_ui(rm,rm,1);
         mp_polymodp(rm,p,ii,gap,ct1,ct2,ct3,&sm,base);
         mp_mul(&rm,*jm1,sm);
         mp_polymodp(rm,p,ii,gap,ct1,ct2,ct3,&sm,base);
         mpz_mul_ui(rm,sm,p-1);
         mp_polymodp(rm,p,ii,gap,ct1,ct2,ct3,&sm,base);
         mpz_add(*jm1,*jm1,sm);
         mp_polymodp(*jm1,p,ii,gap,ct1,ct2,ct3,&sm,base);
         mpz_set(*jm1,sm);
         i=i*2;
        }
    mpz_set(sm,*jm1);
    mpz_set_ui(*jm1,0);
    for(i=0;i<gap5;i++) 
      {
    mpz_tdiv_r_2exp(pm,sm,(deg+1)*5);
    mpz_tdiv_q_2exp(sm,sm,(deg+1)*5);
       for(j=0;j<i;j++)
         {mpz_mul_2exp(pm,pm,(deg+1)*5);}
          mpz_add(*jm1,*jm1,pm);
        }
    mpz_clear(pm);
    mpz_clear(ppm);
    mpz_clear(rm);
    mpz_clear(sm);

  }



void mp_mul(mpz_t *qm, mpz_t xm, mpz_t ym)
{
  long long int x=7000000000;
  mpz_t am,bm,cm,dm;

  mpz_init(am);
  mpz_init(bm);
  mpz_init(cm);
  mpz_init(dm);  

  mpz_set(am,xm);
  mpz_set(bm,ym);

  mpz_tdiv_r_2exp(cm,am,x);
  mpz_tdiv_q_2exp(am,am,x);
    
  mpz_tdiv_r_2exp(dm,bm,x);
  mpz_tdiv_q_2exp(bm,bm,x);

  mpz_mul(*qm,cm,dm);
  mpz_mul(dm,am,dm);
  mpz_mul_2exp(dm,dm,x);
  mpz_add(*qm,*qm,dm);
  mpz_mul(dm,bm,cm);
  mpz_mul_2exp(dm,dm,x);
  mpz_add(*qm,*qm,dm);
  mpz_mul(dm,am,bm);
  mpz_mul_2exp(dm,dm,2*x);
  mpz_add(*qm,*qm,dm);


  mpz_clear(am);
  mpz_clear(bm);
  mpz_clear(cm);
  mpz_clear(dm);

}



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

void mp_polymodp(mpz_t am, long long int p, long long int deg, long long int gap, long long int ct1, long long int ct2, long long int ct3, mpz_t *am1, long long int base)
{
  long long int  n,m,nn,gap5=gap/5,ctnum,cttnum,ctttnum,degap,dgg=5000000,
    bb=dgg*gap,cc,hg=100000;
  long long int x,y[100];
  char temp1[dgg*gap5+1];
  char temp2[dgg*gap5+1];
  char tmp[100];
  char ttmp[ct1+gap5+1];
  char tttmp[ct2+ct1+1];
  char ttttmp[ct3+ct2+1];
  mpz_t xm;


  mpz_init(xm);
  mpz_set_ui(*am1,0);
  y[1]=1;
  for (m=2;m<gap5+1;m++){y[m]=(y[m-1]*base)%p;}
 
  degap=(deg+1)*gap;
  cc=0;
  while(degap>0)
  {

    if (degap<bb)
    {
       bb=degap;
       dgg=bb/gap;
    }
    mpz_tdiv_r_2exp(xm,am,bb);
    mpz_tdiv_q_2exp(am,am,bb);
    nn=bb/5-mpz_sizeinbase(xm,base);
    strcpy(temp1,"");
    strcpy(temp2,"");
    for (m=0;m<hg/100;m++)
      {strcat(temp2,"0");}
    for (m=0;m<100;m++)
      {strcat(temp1,temp2);}
    strcpy(temp2,"");   
    n=nn/hg; 
    for (m=0;m<n;m++)
      {strcat(temp2,temp1);}
    strcpy(temp1,"");
    n=nn%hg;
    for (m=0;m<n;m++)
      {strcat(temp1,"0");}
    strcat(temp1,temp2);

    mpz_get_str(temp2,base,xm);
    strcat(temp1,temp2);

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

    ctnum=0;
    cttnum=0;
    ctttnum=0;

    for (n=0;n<dgg;n++)
      {
	x=0;
	for (m=1;m<gap5+1;m++){
	  if (temp1[gap5*n+m-1]=='1'){x=x+y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='2'){x=x+2*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='3'){x=x+3*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='4'){x=x+4*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='5'){x=x+5*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='6'){x=x+6*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='7'){x=x+7*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='8'){x=x+8*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='9'){x=x+9*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='a'){x=x+10*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='b'){x=x+11*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='c'){x=x+12*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='d'){x=x+13*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='e'){x=x+14*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='f'){x=x+15*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='g'){x=x+16*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='h'){x=x+17*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='i'){x=x+18*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='j'){x=x+19*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='k'){x=x+20*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='l'){x=x+21*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='m'){x=x+22*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='n'){x=x+23*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='o'){x=x+24*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='p'){x=x+25*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='q'){x=x+26*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='r'){x=x+27*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='s'){x=x+28*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='t'){x=x+29*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='u'){x=x+30*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='v'){x=x+31*y[gap5-m+1];}
	  x=x%p;
	}
        mpz_set_ui(xm,x);
	nn=gap5-mpz_sizeinbase(xm,base);
        for (m=0;m<nn;m++)
	  {strcat(ttmp,"0");}
	mpz_get_str(tmp,base,xm);
	strcat(ttmp,tmp);
        ctnum=ctnum+gap5;
	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(xm,temp2,base);
    mpz_mul_2exp(xm,xm,cc);
    mpz_add(*am1,*am1,xm);
    cc=cc+bb; 
    degap=degap-bb;
  } 
  mpz_clear(xm);
}





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


void mp_polymodp0f(mpz_t am, long long int p, long long int deg, long long int gap, long long int f,long long int base)
{
  long long int  n,m,nn,gap5=gap/5,degap,dgg=5000000,
    bb=dgg*gap,dggg,l,hg=100000;
  long long int x,y[100];
  char temp1[dgg*gap5+1];
  char temp2[dgg*gap5+1];
  mpz_t xm;
  FILE *file;

 
  file=fopen("tmpbernc.txt","w");
  mpz_init(xm);

  y[1]=1;
  for (m=2;m<gap5+1;m++){y[m]=(y[m-1]*base)%p;}
 
  degap=(deg+1)*gap;
  dggg=0;
  while(degap>0)
  {
    if (degap<bb)
    {
       bb=degap;
       dgg=bb/gap;
    }
    dggg=dggg+dgg;
    mpz_tdiv_r_2exp(xm,am,bb);
    mpz_tdiv_q_2exp(am,am,bb);
    nn=bb/5-mpz_sizeinbase(xm,base);

    strcpy(temp1,"");
    strcpy(temp2,"");
    for (m=0;m<hg/100;m++)
      {strcat(temp2,"0");}
    for (m=0;m<100;m++)
      {strcat(temp1,temp2);}
    strcpy(temp2,"");   
    n=nn/hg; 
    for (m=0;m<n;m++)
      {strcat(temp2,temp1);}
    strcpy(temp1,"");
    n=nn%hg;
    for (m=0;m<n;m++)
      {strcat(temp1,"0");}
    strcat(temp1,temp2);

    mpz_get_str(temp2,base,xm);
    strcat(temp1,temp2);
    for (n=0;n<dgg;n++)
      {
	x=0;
	for (m=1;m<gap5+1;m++){
	  if (temp1[gap5*n+m-1]=='1'){x=x+y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='2'){x=x+2*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='3'){x=x+3*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='4'){x=x+4*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='5'){x=x+5*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='6'){x=x+6*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='7'){x=x+7*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='8'){x=x+8*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='9'){x=x+9*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='a'){x=x+10*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='b'){x=x+11*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='c'){x=x+12*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='d'){x=x+13*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='e'){x=x+14*y[gap5-m+1];}
	  if (temp1[gap5*n+m-1]=='f'){x=x+15*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='g'){x=x+16*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='h'){x=x+17*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='i'){x=x+18*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='j'){x=x+19*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='k'){x=x+20*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='l'){x=x+21*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='m'){x=x+22*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='n'){x=x+23*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='o'){x=x+24*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='p'){x=x+25*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='q'){x=x+26*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='r'){x=x+27*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='s'){x=x+28*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='t'){x=x+29*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='u'){x=x+30*y[gap5-m+1];}
          if (temp1[gap5*n+m-1]=='v'){x=x+31*y[gap5-m+1];}
	  x=x%p;
	}
       l=dggg-n-1;
 
      if ((x%p==0)&&(l%2==0)&&(l!=0))
        {
          fprintf(file,"%lld\n",p);
          fprintf(file,"%lld\n",l);
          fprintf(file,"%lld\n",f);
         } 
     }
    degap=degap-bb;
  } 
  mpz_clear(xm);
  fclose(file);

}


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

  mpz_t pm, fm, am, f0m, p2m, p3m, fapm, c0m, c1m, ram, qam;
  mpz_t tmp1,tmp2,tmp3,tmp4,tmp5,tmp6;
  long long 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);

}


long long int mp_hantei(long long int p, long long int k, long long int f, mpz_t lm, mpz_t alm, long long 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;
  long long 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,",%lld\n",kt);
    mpz_out_str(file,10,cm);
    fprintf(file,",%lld\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);
  }



void mp_omega(mpz_t am, mpz_t pm, long long int k, mpz_t *jm)
 {
  long long 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);
 }



long long int prmroot(mpz_t lm, long long int p)   /*l-1 is divided by p*/
{
  mpz_t gm,hm,ltm,tmpm;
  long long 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);
}



void prmdiv(mpz_t lm, long long int *p, long long int *m)
  {
  long long 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  */

long long int mp_hantei1(long long int p, long long int k, long long int f, mpz_t lm, mpz_t alm, long long 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,",%lld\n",kt);
    mpz_out_str(file,10,cm);
    fprintf(file,",%lld\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);
  }