﻿ compute square of Gaussian sum
```/*****************************************************************************/
/*									     */
/*  COMPUTE SQUARES OF GAUSSIAN SUMS					     */
/*  11/14/97 (dkc)							     */
/*									     */
/*  This C program computes squares of 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[2][500]; 		     //  c[n][(p-1)/n]
unsigned int sum[1002],temp[1002],save[1002];
unsigned int k[1002];
unsigned int g, h, i, j;
FILE *Outfp;
Outfp = fopen("gauss2b.dat","w");
n=2;
for (g=0; g<100; g++) {
p=input[2*g];
r=input[2*g+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   */
/***************************/
j=1;
for (h=0; h<n; h++) {
j=j*-1;
for (i=0; i<(p-1)/n; i++) {
temp[c[h][i]]=j;
}
}
temp[0]=0;
/************************************/
/*  compute square of Gaussian sum  */
/************************************/
/****************************/
/*  initialize sum squared  */
/****************************/
for (h=1; h<=p-1; h++) {
sum[h+1]=temp[h]*temp[1];
save[h]=temp[h];
}
sum[1]=0;
/*******************/
/*  rotate array   */
/*******************/
for (h=p; h>1; h--) temp[h]=temp[h-1];
temp[1]=0;
/***************************/
/*  compute partial sums   */
/***************************/
for (h=1; h<=p-2; h++) {
j=temp[p];
for (i=p; i>1; i--) {
temp[i]=temp[i-1];
}
temp[1]=j;
for (i=1; i<=p; i++) {
sum[i]=sum[i]+temp[i]*save[h+1];
}
}
/**********************************/
/*  write square of Gaussian sum  */
/**********************************/
fprintf(Outfp," p=%d r=%d sum=%d \n",p,r,sum[1]);
printf(" p=%d r=%d sum=%d \n",p,r,sum[1]);
for (i=2; i<p; i++) {
if (sum[i]!=sum[1]) {
fprintf(Outfp,"error \n");
printf("error \n");
goto bskip;
}
}
if (((sum[1]-sum[p])!=p)&&((sum[p]-sum[1])!=p)) {
fprintf(Outfp,"error \n");
printf("error \n");
goto bskip;
}
if ((p-1)==((p-1)/4)*4) {
if (sum[1]!=-1) {
fprintf(Outfp,"error \n");
printf("error \n");
goto bskip;
}
}
else {
if (sum[1]!=1) {
fprintf(Outfp,"error \n");
printf("error \n");
goto bskip;
}
}
}
bskip:
fclose(Outfp);
return(0);
}
```