﻿ find consecutive roots
```/*****************************************************************************/
/*									     */
/*  COUNT PAIRS OF CONSECUTIVE ROOTS OF X**(P-1)=1(MOD P**2)		     */
/*  02/13/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,S[2],T[2];
unsigned int js,jsh,jj,m[50000],p[50000];
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)
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;
m[0]=jj;
/*****************/
/* compute roots */
/*****************/
for (i=1; i<jh; i++) {
mul32_32(m[i-1], jj, S);
div64_32(S, T, js);
mul32_32(T[1], js, T);
sub64(S, T);
m[i]=T[1];
if (m[i]>jsh)
m[i]=js-m[i];
}
if (m[jh-1]!=1) {
printf("error: p=%d r=%d \n",j,k);
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+2;
}
s[nq/2]=s[nq/2]+1;
printf(" p=%d, %d \n",j,nq);
}
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);
}
```