/*****************************************************************************/
/*									     */
/*  COUNT PAIRS OF CONSECUTIVE ROOTS OF THE CONGRUENCE X**((P-1)/N)=Y(MOD P) */
/*  06/01/07 (dkc)							     */
/*									     */
/*  This C program tests propositions for n=3.				     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "input.h"
unsigned int croots(unsigned int *input, unsigned int index, unsigned int n,
		    unsigned int *output);
void res64_32(unsigned int exp0, unsigned int exp1, unsigned int q0,
	      unsigned int q1, unsigned int *output, unsigned int base);
unsigned int r[100000];
int main () {
unsigned int n=3;		 // n value
unsigned int h,i,sum,p,p1,pn,temp,flag0,flag1,flag2;
unsigned int s[14];
unsigned int U[2];
FILE *Outfp;
Outfp = fopen("out3.dat","w");
for (h=0; h<insize; h++) {
   p=input[2*h];		     // load p
   if (p>40000) 		     // reduce execution time
      break;
   i=croots(input, h, n, s);	     // count consecutive roots of congruences
   if (i==0)			     // continue if n does not divide p-1
      continue;
   if (i==2) {
      printf(" error: bad primitive root \n");
      break;
      }
   printf(" p=%d ",input[2*h]);
/***************************************************/
/* check if the sum of S[i] is equal to (p-1)/n-1  */
/***************************************************/
   sum=0;
   for (i=0; i<n; i++) {
      sum=sum+s[i];
      printf(" %d ",s[i]);
      }
   printf("\n");
   p1=p-1;
   pn=p1/n;
   if (sum!=pn-1) {
      printf(" error: incorrect sum \n");
      break;
      }
/***********************************/
/*  check if s[1]-s[2]=s[2]-s[3]   */
/***********************************/
   flag0=0;
   if ((s[0]-s[1])==(s[1]-s[2]))
      flag0=1;
   if ((s[1]-s[0])==(s[0]-s[2]))
      flag0=1;
/****************************/
/*  check if r is a square  */
/****************************/
   flag1=0;
   temp=(unsigned int)(sqrt((double)pn)+0.01);
   if (temp*temp==pn)
      flag1=1;
/************************************/
/*  check if s[3]-s[2]+2=s[1]-s[3]  */
/************************************/
   flag2=0;
   if ((2*s[2]+2)==(s[0]+s[1]))
      flag2=1;
/*****************************************************/
/* check if r is a square when s[1]-s[2]=s[2]-s[3]   */
/*****************************************************/
   if (flag0==1) {
      fprintf(Outfp," p=%d ",p);
      for (i=0; i<n; i++)
	 fprintf(Outfp," %d ",s[i]);
      fprintf(Outfp,"\n");
      if (flag1==0) {
	 fprintf(Outfp,"p=%d, error: (p-1)/n is not a square \n",p);
	 printf(" error: (p-1)/n is not a square \n");
	 }
      }
/*******************************************************/
/* check if r is a square when s[3]-s[2]+2=s[1]-s[3]   */
/*******************************************************/
   if (flag2==1) {
      if (flag1==0) {
	 fprintf(Outfp,"p=%d, error: (p-1)/n is not a square \n",p);
	 printf(" error: (p-1)/n is not a square \n");
	 }
      }
/***********************************************************************/
/* check if s[1]-s[2]=s[2]-s[3] when r is a square and 2**r<>1(mod p)  */
/***********************************************************************/
   if (flag1==1) {
      res64_32(0, pn, 0, p, U, 2);
      if ((U[0]!=0)||(U[1]!=1)) {
	 if (flag0==0) {
	    fprintf(Outfp,"p=%d, error: differences not equal \n",p);
	    printf(" error: differences not equal \n");
	    }
	 }
/*************************************************************************/
/* check if s[3]-s[2]+2=s[1]-s[3] when r is a square and 2**r==1(mod p)  */
/*************************************************************************/
      else {
	 if (flag2==0) {
	    fprintf(Outfp,"p=%d, error: differences not equal \n",p);
	    printf(" error: differences not equal \n");
	    }
	 }
      }
/****************************************/
/*  check if r/2 is odd when s[1]=s[3]	*/
/****************************************/
      if ((s[0]==s[2])||(s[1]==s[2])) {
	 if ((pn/4)*4==pn) {
	    fprintf(Outfp,"p=%d, error: r/2 is even \n",p);
	    printf(" error: r/2 is even \n");
	    }
	 }
/*************************/
/*  check if s[1]=s[2]	 */
/*************************/
      if (s[0]==s[1]) {
	 fprintf(Outfp,"p=%d, error: s[1] equals s[2] \n",p);
	 printf(" error: s[1] equals s[2] \n");
	 }
   }
fclose(Outfp);
return(0);
}