﻿ compute fourth power of Gaussian sum
```/******************************************************************************/
/*									     */
/*  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);
}
```