/******************************************************************************/
/*									     */
/*  COMPUTE FOURTH POWERS OF GAUSSIAN SUMS				     */
/*  06/16/06 (dkc)							     */
/*									     */
/*  This C program computes fourth powers of generalized Gaussian sums.      */
/*									     */
/*****************************************************************************/
#include <stdio.h>
int main ()
{
/**************************************************************************/
/*  p is a prime, n is a divisor of p-1, and r is a primitive root of p   */
/*  the cyclotomic cosets are stored in c[n][(p-1)/n]			  */
/**************************************************************************/
//unsigned int n=4, p=5, r=2;
//unsigned int n=4, p=13, r=2;
//unsigned int n=4, p=17, r=3;
//unsigned int n=4, p=29, r=2;
//unsigned int n=4, p=37, r=2;
  unsigned int n=4, p=41, r=6;
//unsigned int n=4, p=53, r=2;
//unsigned int n=4, p=61, r=2;
//unsigned int n=4, p=73, r=5;
//unsigned int n=4, p=89, r=3;
//unsigned int n=4, p=97, r=5;
//unsigned int n=4, p=101, r=2;
//unsigned int n=4, p=109, r=6;
//unsigned int n=4, p=197, r=2;
//unsigned int n=4, p=677, r=2;
//unsigned int n=4, p=2917, r=5;
//unsigned int n=4, p=4357, r=2;
//unsigned int n=4, p=5477, r=2;
unsigned int c[4][1500];	      // dimensions are [n][(p-1)/n]
unsigned int sum[6000][2],temp[6000][2],save[6000][2];
unsigned int k[6000];
unsigned int h, i, j, l;
FILE *Outfp;
Outfp = fopen("gauss4.dat","w");
if ((p-1) != ((p-1)/n)*n) {
   fprintf(Outfp,"error:  n does not divide p-1 \n");
   goto bskip;
   }
/******************************************************************/
/*  generate permutation of 1,2,3...,(p-1) (by Fermat's theorem)  */
/******************************************************************/
k[0] = r;
for (i=1; i<p-1; i++) {
   k[i] = k[i-1]*r - ((k[i-1]*r)/p)*p;
   }
for (i=0; i<p-2; i++) {
   if (k[i] == 1) {
      fprintf(Outfp,"error:  r is not a primitive root of p \n");
      goto bskip;
      }
   }
/**********************************************/
/*  sort permutation into cyclotomic cosets   */
/**********************************************/
for (h=0; h<n; h++) {
   j=0;
   for (i=h; i<p-1; i+=n) {
      c[h][j] = k[i];
      j=j+1;
      }
   }
/*****************************/
/*  compute Gaussian sum   */
/*****************************/
for (i=0; i<(p-1)/n; i++) {
   temp[c[0][i]][0]=1;
   temp[c[0][i]][1]=0;
   printf(" %d",c[0][i]);
   }
printf("\n");
for (i=0; i<(p-1)/n; i++) {
   temp[c[1][i]][0]=0;
   temp[c[1][i]][1]=1;
   printf(" %d",c[1][i]);
   }
printf("\n");
for (i=0; i<(p-1)/n; i++) {
   temp[c[2][i]][0]=-1;
   temp[c[2][i]][1]=0;
   printf(" %d",c[2][i]);
   }
printf("\n");
for (i=0; i<(p-1)/n; i++) {
   temp[c[3][i]][0]=0;
   temp[c[3][i]][1]=-1;
   printf(" %d",c[3][i]);
   }
temp[0][0]=0;
temp[0][1]=0;
printf("\n");
printf("\n");
/**************************************/
/*  compute square of Gaussian sum  */
/**************************************/
/****************************/
/*  initialize sum squared  */
/****************************/
for (h=1; h<=p-1; h++) {
   sum[h+1][0]=temp[1][0]*temp[h][0]-temp[1][1]*temp[h][1];
   sum[h+1][1]=temp[1][0]*temp[h][1]+temp[1][1]*temp[h][0];
   save[h][0]=temp[h][0];
   save[h][1]=temp[h][1];
   }
sum[1][0]=0;
sum[1][1]=0;
/*******************/
/*  rotate array   */
/*******************/
for (h=p; h>1; h--) {
   temp[h][0]=temp[h-1][0];
   temp[h][1]=temp[h-1][1];
   }
temp[1][0]=0;
temp[1][1]=0;
/***************************/
/*  compute partial sums   */
/***************************/
for (h=1; h<=p-2; h++) {
   j=temp[p][0];
   l=temp[p][1];
   for (i=p; i>1; i--) {
      temp[i][0]=temp[i-1][0];
      temp[i][1]=temp[i-1][1];
      }
   temp[1][0]=j;
   temp[1][1]=l;
   for (i=1; i<=p; i++) {
      sum[i][0]=sum[i][0]+temp[i][0]*save[h+1][0]-temp[i][1]*save[h+1][1];
      sum[i][1]=sum[i][1]+temp[i][0]*save[h+1][1]+temp[i][1]*save[h+1][0];
      }
   }
/************************************/
/*  write square of Gaussian sum  */
/************************************/
for (i=1; i<=p; i++) {
   fprintf(Outfp," %d %d \n",sum[i][0],sum[i][1]);
   printf(" %d %d \n",sum[i][0],sum[i][1]);
   }
fprintf(Outfp,"\n");
printf("\n");
/****************************/
/*  initialize sum cubed    */
/****************************/
for (h=1; h<=p-1; h++) {
   temp[h+1][0]=0;
   temp[h+1][1]=0;
   }
temp[1][0]=0;
temp[1][1]=0;
/***************************/
/*  compute partial sums   */
/***************************/
for (h=1; h<=p-1; h++) {
   j=sum[p][0];
   l=sum[p][1];
   for (i=p; i>1; i--) {
      sum[i][0]=sum[i-1][0];
      sum[i][1]=sum[i-1][1];
      }
   sum[1][0]=j;
   sum[1][1]=l;
   for (i=1; i<=p; i++) {
      temp[i][0]+=sum[i][0]*save[h][0]-sum[i][1]*save[h][1];
      temp[i][1]+=sum[i][0]*save[h][1]+sum[i][1]*save[h][0];
      }
   }
/********************************/
/*  write cube of Gaussian sum	*/
/********************************/
for (i=1; i<=p; i++) {
   fprintf(Outfp," %d %d \n",temp[i][0],temp[i][1]);
   printf(" %d %d \n",temp[i][0],temp[i][1]);
   }
fprintf(Outfp,"\n");
printf("\n");

/************************************/
/*  initialize fourth power of sum  */
/************************************/
for (h=1; h<=p-1; h++) {
   sum[h+1][0]=0;
   sum[h+1][1]=0;
   }
sum[1][0]=0;
sum[1][1]=0;
/***************************/
/*  compute partial sums   */
/***************************/
for (h=1; h<=p-1; h++) {
   j=temp[p][0];
   l=temp[p][1];
   for (i=p; i>1; i--) {
      temp[i][0]=temp[i-1][0];
      temp[i][1]=temp[i-1][1];
      }
   temp[1][0]=j;
   temp[1][1]=l;
   for (i=1; i<=p; i++) {
      sum[i][0]+=temp[i][0]*save[h][0]-temp[i][1]*save[h][1];
      sum[i][1]+=temp[i][0]*save[h][1]+temp[i][1]*save[h][0];
      }
   }
/****************************************/
/*  write fourth power of Gaussian sum	*/
/****************************************/
for (i=1; i<=p; i++) {
   fprintf(Outfp," %d %d \n",sum[i][0],sum[i][1]);
   printf(" %d %d \n",sum[i][0],sum[i][1]);
   }
fprintf(Outfp,"\n");
printf("\n");

bskip:
fclose(Outfp);
return(0);
}