/* gauss64.c   by Hiroki Sumida-Takahashi 2021, March 18*/
/*   2l=2*(f0*lz+1), 2*f0  4*deg *gap   < 2^63   */
/* if imaginary, change ds to -ds in the def of fm in gauss3 */


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


void gauss1(FILE *, FILE *, long long int, long long int, long long int, mpz_t);
void gauss2(FILE *, FILE *, FILE *, long long int, mpz_t);
void gauss3(FILE *, FILE *, long long int, long long int, long long int, long long int, mpz_t);
long long int prmroot(mpz_t, long long int);
void prmdiv(mpz_t, long long int *, long long int *);
void henkan1(FILE *, FILE *, long long int, long long int);
void henkan2(FILE *, FILE *, long long int, long long int);
void henkan3(FILE *, FILE *, long long int, mpz_t);
void multp(FILE *, FILE *, FILE *, long long int, long long int);


/* gauss64.c  [compute Gauss sum g(l)^{e_\chi_d omega^k} mod lt] ->gsum.txt */ 


main()
  {
    FILE *file1, *file2, *file3, *file4;
   long long int p,ds,kw,lh,lz,lz0,lz1,f0,l,deg;
   mpz_t lm,ltm;

   char filename4[]="excep2.txt";

   file4=fopen(filename4,"r");

 while(fscanf(file4,"%lld %lld %lld \n",&p,&ds,&kw)!=EOF)
  {

   printf("%lld %lld %lld\n",p,ds,kw);

   mpz_init(lm);
   mpz_init(ltm);

   lz0=1;
   lz1=1;

   f0=p*ds;

   lh=0;
   lz=lz0;
   while(lh!=1)
     {
     l=1+lz*f0;
     mpz_set_ui(lm,l);
     if(mpz_probab_prime_p(lm,20)>0){lh=1;} 
     lz++;
     }


   lh=0;
   lz=lz1;
   while(lh!=1)
     {
     mpz_set_ui(ltm,2);
     mpz_mul_ui(ltm,ltm,f0);
     mpz_mul_ui(ltm,ltm,lz);
     mpz_mul_ui(ltm,ltm,l);
     mpz_add_ui(ltm,ltm,1);
     if(mpz_probab_prime_p(ltm,20)>0){lh=1;} 
     lz++;
     }
   
     file1=fopen("g-num1.txt","w");
     file2=fopen("g-num2.txt","w");
     gauss1(file1,file2,p,ds,l,ltm);
     fclose(file1);
     fclose(file2);

     deg=2*f0-1; 
     file1=fopen("g-num1.txt","r");
     file2=fopen("g-num2.txt","r");
     file3=fopen("g-result.txt","w");
     gauss2(file1,file2,file3,deg,ltm);
     fclose(file1);
     fclose(file2);
     fclose(file3);

     file1=fopen("g-result.txt","r");
     file2=fopen("g-sum.txt","a+");
     gauss3(file1,file2,p,ds,kw,l,ltm);
     fclose(file1);
     fclose(file2);

     mpz_clear(lm);
     mpz_clear(ltm);

  }/* while */
 
 } /* main */


/*  gauss1 [set coeffcients of g(T) for f(T)*g(T )] -> file1, file2  */

void gauss1(FILE *file1, FILE *file2, long long int p, long long int ds, long long int l, mpz_t ltm)
 {
   long long int f0, f1, lh, g, gt, lz;
   long long int n, i;
   mpz_t lm,gm,gtm,tmpm,sm,sim,tm,sum;


   mpz_init(lm);
   mpz_init(gm); 
   mpz_init(gtm); 
   mpz_init(tmpm); 
   mpz_init(sm); 
   mpz_init(sim); 
   mpz_init(tm); 
   mpz_init(sum); 

   f0=p*ds;

   mpz_set_ui(lm,l);
   g=prmroot(lm,p);
   mpz_set_ui(gm,g);
   gt=prmroot(ltm,l);
   mpz_set_ui(gtm,gt);

   mpz_sub_ui(tmpm,ltm,1);
   mpz_divexact_ui(tmpm,tmpm,2*f0);
   mpz_powm(sm,gtm,tmpm,ltm);
   mpz_invert(sim,sm,ltm);
   
   mpz_sub_ui(tmpm,ltm,1);
   mpz_divexact(tmpm,tmpm,lm);
   mpz_powm(tm,gtm,tmpm,ltm);

   for(n=2*f0-1;n>-1;n--)
     {
     mpz_powm_ui(tmpm,sim,n,ltm);
     mpz_powm_ui(tmpm,tmpm,n,ltm);
     mpz_out_str(file1,16,tmpm);
     fprintf(file1,"\n");

     if (n<f0)
       {
       mpz_set_ui(sum,0);
       for(i=0;i<(l-1)/f0;i++)
         {
	  mpz_set_ui(tmpm,f0);
	  mpz_mul_ui(tmpm,tmpm,i);
          mpz_add_ui(tmpm,tmpm,n);
          mpz_powm(tmpm,gm,tmpm,lm);
          mpz_powm(tmpm,tm,tmpm,ltm);
          mpz_add(sum,sum,tmpm);
          mpz_mod(sum,sum,ltm);
         }
       mpz_powm_ui(tmpm,sm,n,ltm);
       mpz_powm_ui(tmpm,tmpm,n,ltm);
       mpz_mul(sum,sum,tmpm);
       mpz_mod(sum,sum,ltm);
       mpz_out_str(file2,16,sum);
       fprintf(file2,"\n");
      }
      else fprintf(file2,"0\n");
    }

   mpz_clear(lm);
   mpz_clear(gm);
   mpz_clear(gtm);
   mpz_clear(tmpm);
   mpz_clear(sm);
   mpz_clear(sim);
   mpz_clear(tm);
   mpz_clear(sum);

}

/* gauss2   [convolution f(T)*g(T)]   file1, file2 -> file3 */ 


void gauss2(FILE *file1, FILE *file2, FILE *file3, long long int deg, mpz_t ltm)
{
   FILE *file1t, *file2t, *file3t, *file4t;
   long long int gap;
   mpz_t im1; 

   mpz_init(im1);

   mpz_mul(im1,ltm,ltm);
   mpz_mul_ui(im1,im1,deg+1);
   gap=mpz_sizeinbase(im1,16)+1;

   file1t=fopen("g-num1t.txt","w");  
   henkan1(file1,file1t,deg,gap);
   fclose(file1t);

   file2t=fopen("g-num2t.txt","w");  
   henkan1(file2,file2t,deg,gap);
   fclose(file2t);

   file1t=fopen("g-num1t.txt","r");
   file2t=fopen("g-num2t.txt","r");  
   file3t=fopen("g-tmp1.txt","w");
   multp(file1t,file2t,file3t,deg,gap);
   fclose(file1t);
   fclose(file2t);
   fclose(file3t);

   file3t=fopen("g-tmp1.txt","r");
   file4t=fopen("g-tmp2.txt","w");  
   henkan2(file3t,file4t,deg,gap);
   fclose(file3t);
   fclose(file4t);

   file4t=fopen("g-tmp2.txt","r");
   henkan3(file4t,file3,deg,ltm);
   fclose(file4t);
   

   mpz_clear(im1);

}

/*  gauss3  [multiply gauss sums]  file1  ->  file2    */


void gauss3(FILE *file1, FILE *file2, long long int p, long long int ds, long long int kw, long long int l, mpz_t ltm)
 {
   long long int f0, f1, lh, g, gt, lz, kro;
   long long int n, i;
   mpz_t lm,lt1m,gm,gtm,tmpm,sm,gsum,ggsum,fm,f0m,omegam;

   mpz_init(lm); 
   mpz_init(lt1m); 
   mpz_init(gm); 
   mpz_init(gtm); 
   mpz_init(tmpm); 
   mpz_init(sm); 
   mpz_init(gsum); 
   mpz_init(ggsum); 
   mpz_init(fm);
   mpz_init(f0m);
   mpz_init(omegam);

   f0=p*ds;
   mpz_set_ui(f0m,f0);
   mpz_set_ui(fm,-ds);   /* real: ds, img: -ds */
   mpz_set_ui(lm,l);

   g=prmroot(lm,p);
   mpz_set_ui(gm,g);
   gt=prmroot(ltm,l);
   mpz_set_ui(gtm,gt);
   mpz_set(lt1m,ltm);
   mpz_sub_ui(lt1m,lt1m,1);

   mpz_sub_ui(tmpm,ltm,1);
   mpz_divexact_ui(tmpm,tmpm,2*f0);
   mpz_powm(sm,gtm,tmpm,ltm);
   
   mpz_set_ui(ggsum,1); 
   for(n=1;n<f0+1;n++)
     {
      mpz_inp_str(gsum,file1,16);
     }

   for(n=f0-1;n>-1;n--)
     {
      mpz_inp_str(gsum,file1,16);
      mpz_gcd_ui(tmpm,f0m,n);
      if (mpz_cmp_ui(tmpm,1)==0)
        {
         mpz_powm_ui(tmpm,sm,n,ltm);
         mpz_powm_ui(tmpm,tmpm,n,ltm);
         mpz_mul(gsum,tmpm,gsum);
         mpz_mod(gsum,gsum,ltm);
         mpz_set_ui(tmpm,n%p);
         mpz_powm_ui(omegam,tmpm,p-2+kw,lt1m);
         mpz_powm(tmpm,gsum,omegam,ltm);
         kro=mpz_kronecker_ui(fm,n);
         if (kro==1)
           {
            mpz_mul(ggsum,ggsum,tmpm);
            mpz_mod(ggsum,ggsum,ltm);
           }
         else
           {
            mpz_invert(tmpm,tmpm,ltm);
            mpz_mul(ggsum,ggsum,tmpm);
            mpz_mod(ggsum,ggsum,ltm);
           }
        }
     }

   fprintf(file2,"%lld,%lld,%lld \n",p,ds,kw);
   mpz_out_str(file2,10,lm);
   fprintf(file2,"\n");
   mpz_out_str(file2,10,ltm);
   fprintf(file2,"\n");
   mpz_out_str(file2,10,ggsum);
   fprintf(file2,"\n");
   mpz_divexact_ui(lt1m,lt1m,p);
   mpz_powm(ggsum,ggsum,lt1m,ltm);
   mpz_out_str(file2,10,ggsum);
   fprintf(file2,"\n");

   mpz_clear(lm);
   mpz_clear(lt1m);
   mpz_clear(gm);
   mpz_clear(gtm);
   mpz_clear(tmpm);
   mpz_clear(sm);
   mpz_clear(gsum);
   mpz_clear(ggsum);
   mpz_clear(fm);
   mpz_clear(f0m);
   mpz_clear(omegam);


}



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);

  } 


void henkan1(FILE *file1, FILE *file1t, long long int deg, long long int gap)
  {
   mpz_t im1;
   long long int i, n, x;

   mpz_init(im1);
   for(i=0;i<deg+1;i++)
    { 
      mpz_inp_str(im1,file1,16);
      x=gap-mpz_sizeinbase(im1,16); 
      for (n=1;n<x+1;n++)
  	{fprintf(file1t,"0");}
      mpz_out_str(file1t,16,im1);
    }
    mpz_clear(im1);
 }

void henkan2(FILE *file1, FILE *file1t, long long int deg, long long int gap)
  {
   FILE *ftmp;
   mpz_t im1;
   long long int n;
   char s[100];

  mpz_init(im1);
  ftmp=fopen("g-ftmp.txt","w");
  for (n=0;n<deg+1;n++){
    fgets(s,gap+1,file1);
    fprintf(ftmp,"%s\n",s);
   }
  fclose(ftmp);
  ftmp=fopen("g-ftmp.txt","r");
  for (n=0;n<deg+1;n++){
    mpz_inp_str(im1,ftmp,16);
    mpz_out_str(file1t,16,im1);
    fprintf(file1t,"\n");
   }
  fclose(ftmp);  
  mpz_clear(im1);
}


void henkan3(FILE *file1, FILE *file1t, long long int deg, mpz_t lm)
  {
   mpz_t im1;
   long long int i;

   mpz_init(im1);
   for(i=0;i<deg+1;i++)
    { 
      mpz_inp_str(im1,file1,16);
      mpz_mod(im1,im1,lm);
      mpz_out_str(file1t,16,im1);
      fprintf(file1t,"\n");
    }
    mpz_clear(im1);
 }




void multp(FILE *file1, FILE *file2, FILE *file3, long long int deg, long long int gap)
   {
   mpz_t im1,im2,im3;
   long long int n,x;

   mpz_init(im1);
   mpz_init(im2);
   mpz_init(im3);

   mpz_inp_str(im1,file1,16);
   mpz_inp_str(im2,file2,16);

   mpz_mul(im3,im1,im2);

   mpz_tdiv_q_2exp(im1,im3,4*(deg+1)*gap);
   mpz_tdiv_r_2exp(im2,im3,4*(deg+1)*gap);
   mpz_add(im3,im1,im2); 
  
   x=(deg+1)*gap-mpz_sizeinbase(im3,16);
   for (n=1;n<x+1;n++)
	{fprintf(file3,"0");}
	mpz_out_str(file3,16,im3);

   mpz_clear(im1);
   mpz_clear(im2);
   mpz_clear(im3);


   }


/*
   p=19;
   ds=5;
   kw=8;
   lz0=398;
   lz1=1;
*/

/* 
   p=34301;
   ds=8;
   kw=114;
   lz0=1;
   lz1=1;
*/

/*
   p=55621;
   ds=56;
   kw=9294;
   lz0=1;
   lz1=1;
*/

/* 
   p=62791;
   ds=193;
   kw=57100;
   lz0=1;
   lz1=1;
*/

/* 
   p=157229;
   ds=8;
   kw=140434;
   lz0=1;
   lz1=1;
*/

/* 
   p=384487;
   ds=161;
   kw=13724;
   lz0=1;
   lz1=1;
*/

/*
   p=937943;
   ds=167;
   kw=11057;
   lz0=1;
   lz1=1;
*/