﻿ compute cube of Gaussian sum
```/*****************************************************************************/
/*									     */
/*  COMPUTE CUBES OF GAUSSIAN SUMS					     */
/*  11/14/97 (dkc)							     */
/*									     */
/*  This C program computes cubes of generalized Gaussian sums. 	     */
/*									     */
/*****************************************************************************/
#include "input.h"
#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, p, r;
unsigned int c[3][1000];	      // dimensions are [n][(p-1)/n]
unsigned int sum[3000][2],temp[3000][2],save[3000][2];
unsigned int k[3000];
unsigned int g, h, i, j, l, t, count;
int d, e;
FILE *Outfp;
Outfp = fopen("gauss3b.dat","w");
count=0;
n=3;
for (g=0; g<200; g++) {
p=input[2*g];
r=input[2*g+1];
if ((p-1)!=((p-1)/n)*n)
continue;
count+=1;
/******************************************************************/
/*  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;	     // set to 1
temp[c[0][i]][1]=0;
}
for (i=0; i<(p-1)/n; i++) {	     // set to rho
temp[c[1][i]][0]=0;
temp[c[1][i]][1]=1;
}
for (i=0; i<(p-1)/n; i++) {
temp[c[2][i]][0]=-1;	     // set to rho**2
temp[c[2][i]][1]=-1;
}
temp[0][0]=0;
temp[0][1]=0;
/**************************************/
/*  compute square of Gaussian sum  */
/**************************************/
/****************************/
/*  initialize sum squared  */
/****************************/
for (h=1; h<=p-1; h++) {
t=temp[1][1]*temp[h][1];
sum[h+1][0]=temp[1][0]*temp[h][0]-t;
sum[h+1][1]=temp[1][0]*temp[h][1]+temp[1][1]*temp[h][0]-t;
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++) {
t=temp[i][1]*save[h+1][1];
sum[i][0]+=temp[i][0]*save[h+1][0]-t;
sum[i][1]+=temp[i][0]*save[h+1][1]+temp[i][1]*save[h+1][0]-t;
}
}
/****************************/
/*  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++) {
t=sum[i][1]*save[h][1];
temp[i][0]+=sum[i][0]*save[h][0]-t;
temp[i][1]+=sum[i][0]*save[h][1]+sum[i][1]*save[h][0]-t;
}
}
/********************************/
/*  write cube of Gaussian sum	*/
/********************************/
fprintf(Outfp," %d, %d, %d, \n",p,temp[1][0],temp[1][1]);
printf(" %d, %d, %d, \n",p,temp[1][0],temp[1][1]);
d=(int)temp[1][0]-(int)temp[p][0];
e=(int)p;
if (d!=(d/e)*e) {
fprintf(Outfp," error \n");
printf(" error \n");
goto bskip;
}
}
bskip:
fprintf(Outfp," count=%d \n",count);
printf(" count=%d \n",count);
fclose(Outfp);
return(0);
}
```