/*  gauss.c for f_n=f0 p */

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

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


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

   mpz_init(lm);
   mpz_init(ltm);
   mpz_init(lam);
   mpz_init(lbm);
 

   p=5;
   ds=58424;
   kw=1;
   lz0=1;
   lz1=1;
   l=16066601;
   ll=20448401; 

   mpz_set_ui(lam,l);
   mpz_set_ui(lbm,ll);
 
   if (ds%p!=0)
      {f0=p*p*ds;} 
     else 
      {f0=p*ds;}

   mpz_set(lm,lam);

   mpz_set_ui(ltm,4);
   mpz_mul_ui(ltm,ltm,2*f0);
   mpz_mul(ltm,ltm,lam);
   mpz_mul(ltm,ltm,lbm);
   mpz_add_ui(ltm,ltm,1);
  
     file1=fopen("g-num1.txt","w");
     file2=fopen("g-num2.txt","w");
     gauss1(file1,file2,p,ds,l,ll,ltm);
     fclose(file1);
     fclose(file2);

   printf("111\n");

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

   printf("222\n");


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

   printf("333\n");



     mpz_clear(lm);
     mpz_clear(ltm);
     mpz_clear(lam);
     mpz_clear(lbm);

  }


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

gauss1(FILE *file1, FILE *file2, int p, int ds, int l, int ll, mpz_t ltm)
 {
   int f0, f1, lh, g, gt, lz;
   int n, i,k;
   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); 
 
   if (ds%p!=0)
      {f0=p*p*ds;} 
     else 
      {f0=p*ds;}

   mpz_set_ui(lm,l);

   g=prmroot(lm,p);
 printf("%d\n",g);  
   mpz_set_ui(gm,g);

   gt=prmroot1(ltm,l,ll);
 printf("%d\n",gt);  
   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 */ 


gauss2(FILE *file1, FILE *file2, FILE *file3, int deg, mpz_t ltm)
{
   FILE *file1t, *file2t, *file3t, *file4t;
   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    */


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

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


   if (ds%p!=0)
      {f0=p*p*ds;} 
     else 
      {f0=p*ds;}

   mpz_set_ui(f0m,f0);
   mpz_set_ui(fm,-ds);
   mpz_set_ui(lm,l);
   mpz_set_ui(p2m,p*p);

   gt=prmroot1(ltm,l,ll);
   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,10);
     }

   for(n=f0-1;n>-1;n--)
     {
      mpz_inp_str(gsum,file1,10);
      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);
         mpz_powm_ui(tmpm,tmpm,p-1,p2m);
         mpz_sub_ui(tmpm,tmpm,1);
         mpz_divexact_ui(tmpm,tmpm,p);
         mpz_mul_ui(tmpm,tmpm,n);
         mpz_mul_ui(tmpm,tmpm,p);
         mpz_add_ui(tmpm,tmpm,n);
         mpz_mod(tmpm,tmpm,p2m);

         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,"%d,%d,%d \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(gtm);
   mpz_clear(tmpm);
   mpz_clear(sm);
   mpz_clear(gsum);
   mpz_clear(ggsum);
   mpz_clear(fm);
   mpz_clear(f0m);
   mpz_clear(omegam);
   mpz_clear(p2m);   

}



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



int prmroot1(mpz_t lm, int p, int pp)   /*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);
    }
  pd[2]=pp;
  while(mpz_divisible_ui_p(ltm,pp)!=0)
    {
      mpz_divexact_ui(ltm,ltm,pp);
    }
  k=3;
  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);

  } 


henkan1(FILE *file1, FILE *file1t, int deg, int gap)
  {
   mpz_t im1;
   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);
 }

henkan2(FILE *file1, FILE *file1t, int deg, int gap)
  {
   FILE *ftmp;
   mpz_t im1;
   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);
}


henkan3(FILE *file1, FILE *file1t, int deg, mpz_t lm)
  {
   mpz_t im1;
   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,10,im1);
      fprintf(file1t,"\n");
    }
    mpz_clear(im1);
 }




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

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

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

printf("4444\n");  
   mpz_mul(im3,im1,im2);
printf("5555\n");
  
   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);


   }


/*
  FILE *file2;
 
  file2=fopen("g-sum1.txt","a+");
  mpz_out_str(file2,10,lm);
  fprintf(file2,"\n");
  fclose(file2);
*/