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