/*****************************************************************************/
/*									     */
/*  COUNT PAIRS OF CONSECUTIVE ROOTS OF X**(P-1)=Y(MOD P*P), Y**P=1(MOD P*P) */
/*  02/14/07 (dkc)							     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "input.h"
void div64_32(unsigned int *dividend, unsigned int *quotient,
	      unsigned int divisor);
void mul32_32(unsigned int a0, unsigned int m0, unsigned int *product);
void sub64(unsigned int *a, unsigned int *b);
int main () {
unsigned int g,j,k,j1,h,jh,i,nq,kk,sum,S[2],T[2];
unsigned int js,jsh,jj,m[5000],p[5000];
int q,f,ilb[5001],irb[5001];
unsigned int s[101],t[101];
unsigned int l[5000];
FILE *Outfp;
Outfp = fopen("out1.dat","w");
/***********************/
/*  clear histograms   */
/***********************/
for (i=0; i<101; i++) {
   s[i]=0;
   t[i]=0;
   }
/*************************************/
/*  load primes and primitive roots  */
/*************************************/
for (g=0; g<insize; g++) {
   j=input[2*g];		    // p (prime)
   if (j>10000)
      break;
   k=input[2*g+1];		    // primitive root
   j1=j-1;			    // p-1
   js=j*j;			    // p*p
   jsh=(js-1)/2;		    // (p*p-1)/2
   jh=j1/2;			    // (p-1)/2
/**********************************/
/* compute least residue of k**p  */
/**********************************/
   jj=k;
   for (i=1; i<j; i++) {
      mul32_32(jj, k, S);
      div64_32(S, T, js);
      mul32_32(T[1], js, T);
      sub64(S, T);
      jj=T[1];
      }
   if (jj>jsh)
      jj=js-jj;
   l[0]=jj;
/*****************************************/
/* compute roots of x**(p-1)=1(mod p*p)  */
/*****************************************/
   for (i=1; i<jh; i++) {
      mul32_32(l[i-1], jj, S);
      div64_32(S, T, js);
      mul32_32(T[1], js, T);
      sub64(S, T);
      l[i]=T[1];
      if (l[i]>jsh)
	 l[i]=js-l[i];
      }
   sum=0;
   for (i=0; i<j; i++) {
/****************************************/
/* compute roots of x**(p-1)=y(mod p*p) */
/****************************************/
      kk=i*j+1; 		    // y
      if (kk>jsh)
	 kk=js-kk;
      for (h=0; h<jh; h++) {
	 mul32_32(l[h], kk, S);
	 div64_32(S, T, js);
	 mul32_32(T[1], js, T);
	 sub64(S, T);
	 m[h]=T[1];
	 if (m[h]>jsh)
	    m[h]=js-m[h];
	 }
      if ((i==0)&&(m[jh-1]!=1)) {
	 printf("error p=%d r=%d %d \n",j,k,m[jh-1]);
	 goto zskip;
	 }
/*****************************************/
/* sort roots using monkey-puzzle sort	 */
/*****************************************/
      ilb[1]=0;
      irb[1]=0;
      for (h=2; h<=jh; h++) {
	 ilb[h]=0;
	 irb[h]=0;
	 q=1;
lb:	 if (m[h-1]>m[q-1])
	    goto le;
	 if(ilb[q]==0)
	    goto ld;
	 q=ilb[q];
	 goto lb;
ld:	 irb[h]=-q;
	 ilb[q]=(int)h;
	 continue;
le:	 if (irb[q]<0)
	    goto lf;
	 if (irb[q]==0)
	    goto lf;
	 q=irb[q];
	 goto lb;
lf:	 irb[h]=irb[q];
	 irb[q]=(int)h;
	 }
      f=-1;
      q=1;
      goto li;
lh:   q=ilb[q];
li:   if (ilb[q]>0)
	 goto lh;
lj:   f=f+1;
      p[f]=m[q-1];
      if (irb[q]<0)
	 goto lk;
      if (irb[q]==0)
	 goto ll;
      q=irb[q];
      goto li;
lk:   q=-irb[q];
      goto lj;
/*************************************/
/* count pairs of consecutive roots  */
/*************************************/
ll:   nq=0;
      for (h=1; h<jh; h++) {
	 if ((p[h]-p[h-1])==1)
	    nq=nq+1;
	 }
      if (i!=0)
	 s[nq]=s[nq]+1;
      else
	 t[nq]=t[nq]+1;
      sum=sum+nq;
      }
   sum=sum*2+1;
   if (sum!=j-2) {
      printf("error: sum=%d p-2=%d \n",sum,j-2);
      goto zskip;
      }
   printf("p=%d \n",j);
   }
for (i=0; i<101; i++) {
   fprintf(Outfp," i=%d, %d %d \n",i,s[i],t[i]);
   printf(" i=%d, %d %d \n",i,s[i],t[i]);
   }
zskip:
fclose(Outfp);
return(0);
}