/*****************************************************************************/
/*									     */
/*  COUNT PAIRS OF CONSECUTIVE ROOTS OF X**(P-1)=1(MOD P**2)		     */
/*  02/13/07 (dkc)							     */
/*									     */
/*  Note:  This program is the same as "test1s.c" except that larger p values*/
/*  can be processed.							     */
/*									     */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
#include "input.h"
void mul64_64(unsigned int a0, unsigned int a2, unsigned int m0,
	      unsigned int m2, unsigned int *product);
void div128_64(unsigned int a0, unsigned int a1, unsigned int a2,
	       unsigned int a3, unsigned int *quotient, unsigned int d2,
	       unsigned int d3);
void sub128(unsigned int *a, unsigned int *b);
void mul32_32(unsigned int a0, unsigned int m0, unsigned int *product);
void div64_32(unsigned int *dividend, unsigned int *quotient,
	      unsigned int divisor);
void sub64(unsigned int *a, unsigned int *b);
int main () {
unsigned int g,j,k,j1,h,jh,i,nq,S[4],T[4],U[2],V[2],W[2],X[2];
unsigned int m[50000*2];
int q,f,ilb[50001],irb[50001];
unsigned int s[101];
FILE *Outfp;
Outfp = fopen("output.dat","w");
/**********************/
/*  clear histogram   */
/**********************/
for (i=0; i<101; i++)
   s[i]=0;
/*************************************/
/*  load primes and primitive roots  */
/*************************************/
for (g=0; g<insize; g++) {
   j=input[2*g];		    // p (prime)
   if (j<60000)
      continue;
   k=input[2*g+1];		    // primitive root
   j1=j-1;			    // p-1
   jh=j1/2;			    // (p-1)/2
   mul32_32(j, j, U);		    // p*p
   V[0]=0;
   V[1]=1;
   sub64(U, V); 		    // p*p-1
   div64_32(V, V, 2);		    // (p*p-1)/2
/**********************************/
/* compute least residue of k**p  */
/**********************************/
   T[0]=0;
   T[1]=k;
   for (i=1; i<j; i++) {
      mul64_64(T[0], T[1], 0, k, S);
      div128_64(S[0], S[1], S[2], S[3], T, U[0], U[1]);
      mul64_64(T[2], T[3], U[0], U[1], T);
      sub128(S, T);
      T[0]=T[2];
      T[1]=T[3];
      }
   m[0]=T[0];
   m[1]=T[1];
   W[0]=V[0];
   W[1]=V[1];
   sub64(m, W);
   if ((W[0]&0x80000000)==0)
      sub64(U, m);
/*****************/
/* compute roots */
/*****************/
   for (i=1; i<jh; i++) {
      mul64_64(m[2*i-2], m[2*i-1], m[0], m[1], S);
      div128_64(S[0], S[1], S[2], S[3], T, U[0], U[1]);
      mul64_64(T[2], T[3], U[0], U[1], T);
      sub128(S, T);
      m[2*i]=T[2];
      m[2*i+1]=T[3];
      W[0]=V[0];
      W[1]=V[1];
      sub64(m+2*i, W);
      if ((W[0]&0x80000000)==0)
	 sub64(U, m+2*i);
      }
   if ((m[2*jh-2]!=0)||(m[2*jh-1]!=1)) {
      printf("error: p=%d r=%d %d %d \n",j,k,m[2*jh-2],m[2*jh-1]);
      goto zskip;
      }
/*****************************************/
/* sort roots using monkey-puzzle sort	 */
/*****************************************/
   nq=0;
   ilb[1]=0;
   irb[1]=0;
   for (h=2; h<=jh; h++) {
      ilb[h]=0;
      irb[h]=0;
      q=1;
lb:   W[0]=m[2*q-2];
      W[1]=m[2*q-1];
      sub64(m+2*(h-1), W);
      if ((W[0]&0x80000000)==0)
	 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;
   if (f!=0) {
      if (((m[2*q-2]-X[0])==0)&&((m[2*q-1]-X[1])==1))
	 nq=nq+1;
      }
   X[0]=m[2*q-2];
   X[1]=m[2*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:s[nq]=s[nq]+1;
   printf(" p=%d, %d \n",j,nq*2);
   }
fprintf(Outfp,"\n");
printf("\n");
for (i=0; i<101; i++) {
   fprintf(Outfp," i=%d, %d \n",i,s[i]);
   printf(" i=%d, %d \n",i,s[i]);
   }
zskip:
fclose(Outfp);
return(0);
}