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