/******************************************************************************/
/*									      */
/*  FIND ASSOCIATED CYCLES (where (L,K) values are in arithmetic progression) */
/*  06/09/13 (dkc)							      */
/*									      */
/*  This C program finds associated 3n+c cycles for c<200000.  Also, whether  */
/*  c divides 2^(K+L)-3^K is determined.				      */
/*									      */
/******************************************************************************/
#include <math.h>
#include <stdio.h>
#include "newlk.h"
#include "newlk1.h"
#include "newlk2.h"
#include "newlk3.h"

unsigned int divn_32(unsigned int *a0, unsigned int *quotient, unsigned int dn,
	     unsigned int n);
void lshiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void addn(unsigned int *a, unsigned int *b, unsigned int n);
void subn(unsigned int *c, unsigned int *d, unsigned int n);
void copyn(unsigned int *a, unsigned int *b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);
unsigned int orn(unsigned int *a, unsigned int n);

unsigned int test(unsigned int c, unsigned int L, unsigned int K, unsigned int m,
		  unsigned int flag) {
unsigned int A[2048],B[2048],C[2048],Q[2048],Z[2048],o,h;
setn(A, 1, m);
for (h=0; h<K; h++) {	  // 3**K
   copyn(A, B, m);
   addn(A, A, m);
   addn(B, A, m);
   }
setn(B, 1, m);
lshiftn(B, C, K+L, m);	  // 2**(K+L)
subn(C, A, m);		  // 2**(K+L)-3**K
if ((A[0]&0x80000000)!=0) {  // |2**(K+L)-3**K|
   setn(Z, 0, m);
   subn(Z, A, m);
   }
o=divn_32(A, Q, c, m);	  // |2**(K+L)-3**K|/c
if (o!=1)
   return(2);
if (flag!=0) {
   setn(B, 0, m);
   for (h=0; h<c; h++)	     // (|2**(K+L)-3**K|/c)*c
      addn(Q, B, m);
   subn(A, B, m);	     // compare to |2**(K+L)-3**K|
   }
else {
   copyn(Q, B, m);
   if (B[m-1]==1)
      B[m-1]=0;
   else
      B[m-1]=0xffffffff;
   }
o=orn(B, m);
if (o==0)
   return(1);
else
   return(0);
}

int main () {
unsigned int tflag=0;  // set to 1 to not count multiples of (K,L) values
int limit=1;  // maximum multiple of L1 (that divides L3-L2), set to 1
unsigned int flag4=1;  // output flag (normally set to 1)
unsigned int flag5=1;  // L-K flag (normally set to 1)
unsigned int flag6=0;  // divide by 6 flag (normally set to 0)
unsigned int flag7=0;  // rounding flag (normally set to 0)
unsigned int flag8=1;  // divide test (normally set to 0)
unsigned int flag9=0;  // non-associated cycle flag (normally set to 0)
unsigned int flag10=0; // minimum |L-K| flag (normally set to 0)
unsigned int flag11=1; // compare first sorted (L,K) value (normally set to 1)
unsigned int swap=0;   // (L,K) swap index (normally set to 0), used with "flag10"
int bound=3;  // |L-K| bound
unsigned int select=4;	// cycle count for L-K histograms
unsigned int climit=200000; // maximum c value
unsigned int owrite=0;	// output flag (normally set to 1)
unsigned int klcount,total,total1,ncyc[20],ncyco[20],index,flag[100],offset4;
unsigned int c,flag1,flag2,flag3,first,kount,cnt,iters,mismatch,total2,total3;
int L1,K1,L2,K2,L3,K3,L4,K4,out[100*4],min,mink,dL,dK,mindif,delta,mindif1,offset;
unsigned int i,k,l,smallcnt,parity1,parity2,parity3,parity4,start,histo17[200];
unsigned int histo1[900],histo2[900],histo3[900],histo4[900],histo5[900],histo6[900];
unsigned int histo7[900],histo8[900],histo9[900],histo10[200],histo11[200];
unsigned int histo12[200],histo13[900],swit,histo14[200],histo15[100],histo16[300];
int sum,sum2,j,offset1,L5,K5,delta1,delta2,L6,K6,offset2,offset3;
double lambda,ratio,mean,var;
FILE *Outfp;
Outfp = fopen("test23r.dat","w");
if ((flag9!=0)&&(flag5==0)) {
   printf("error: flag5 must be set when flag9 is set \n");
   goto zskip;
   }
for (i=0; i<900; i++) {
   histo1[i]=0; 		      // clear histograms of L-K values
   histo2[i]=0; 		      // (for different numbers of cycles)
   histo3[i]=0;
   histo4[i]=0;
   histo5[i]=0;
   histo6[i]=0;
   histo7[i]=0;
   histo8[i]=0;
   histo9[i]=0;
   }
for (i=0; i<200; i++) {
   histo10[i]=0;		      // smallest L
   histo11[i]=0;		      // next smallest L
   histo12[i]=0;		      // largest L
   histo14[i]=0;		      // smallest |L-K| in associated cycles
   }
for (i=0; i<900; i++)		      // clear histogram of L-K values
   histo13[i]=0;
for (i=0; i<100; i++)
   histo15[i]=0;
for (i=0; i<200; i++)
   histo17[i]=0;
for (i=0; i<300; i++)
   histo16[i]=0;
offset=100;
offset1=275;
offset2=50;
offset3=75;
offset4=300;
swit=0; 			      // clear switched (L1,K1) (L2,K2) count
mismatch=0;			      // clear L1>K1, L2>K2 or L1<K1,L2<K2 count
iters=0;			      // clear associated cycle count
for (i=0; i<20; i++) {		      // clear histogram of cycle counts
   ncyc[i]=0;
   ncyco[i]=0;
   }
total=0;
total1=0;
total2=0;
total3=0;
smallcnt=0;			      // clear small |L-K| count
parity1=0;			      // clear rounding parity flags
parity2=0;
parity3=0;
parity4=0;
c=1;
for (l=0; l<(33333+6667+6667+6666+6667+6667); l++) {
   if (c>climit)
      break;
   flag2=0;			     // associated cycle flag
   if (c<70000)
      klcount=count[l]; 	     // cycle count
   else {
      if (c<100000)
	 klcount=count1[l-23333];	// cycle count
      else {
	 if (c<150000)
	    klcount=count2[l-33333];   // cycle count
	 else
	    klcount=count3[l-50000];   // cycle count
	 }
      }
   for (k=0; k<klcount; k++) {	     // copy (L,K) values
      if (c<70000) {
	 L1=(int)kl[2*(k+total)];
	 K1=(int)kl[2*(k+total)+1];
	 }
      else {
	 if (c<100000) {
	    L1=(int)kl1[2*(k+total1)];
	    K1=(int)kl1[2*(k+total1)+1];
	    }
	 else {
	    if (c<150000) {
	       L1=(int)kl2[2*(k+total2)];
	       K1=(int)kl2[2*(k+total2)+1];
	       }
	    else {
	       L1=(int)kl3[2*(k+total3)];
	       K1=(int)kl3[2*(k+total3)+1];
	       }
	    }
	 }
      out[2*k]=L1;
      out[2*k+1]=K1;
      }
//
// cycles without attachment points and no corresponding interrelated cycles
// with attachment points
//
   if (c==1) {
      out[2*klcount]=0;
      out[2*klcount+1]=1;
      klcount=klcount+1;
      out[2*klcount]=1;
      out[2*klcount+1]=1;
      klcount=klcount+1;
      out[2*klcount]=1;
      out[2*klcount+1]=2;
      klcount=klcount+1;
      }
   if (c==11) {
      out[2*klcount]=1;
      out[2*klcount+1]=3;
      klcount=klcount+1;
      }
   if (c==49) {
      out[2*klcount]=1;
      out[2*klcount+1]=4;
      klcount=klcount+1;
      }
   if (c==179) {
      out[2*klcount]=1;
      out[2*klcount+1]=5;
      klcount=klcount+1;
      }
   if (c==601) {
      out[2*klcount]=1;
      out[2*klcount+1]=6;
      klcount=klcount+1;
      }
   if (c==1931) {
      out[2*klcount]=1;
      out[2*klcount+1]=7;
      klcount=klcount+1;
      }
   if (c==6049) {
      out[2*klcount]=1;
      out[2*klcount+1]=8;
      klcount=klcount+1;
      }
   if (c==18659) {
      out[2*klcount]=1;
      out[2*klcount+1]=9;
      klcount=klcount+1;
      }
   if (c==57001) {
      out[2*klcount]=1;
      out[2*klcount+1]=10;
      klcount=klcount+1;
      }
   if (c==173051) {
      out[2*klcount]=1;
      out[2*klcount+1]=11;
      klcount=klcount+1;
      }
//
   if (c==791) {
      out[2*klcount]=2;
      out[2*klcount+1]=8;
      klcount=klcount+1;
      }
//
   if (c==85) {
      out[2*klcount]=4;
      out[2*klcount+1]=8;
      klcount=klcount+1;
      }
   if (c==145) {
      out[2*klcount]=4;
      out[2*klcount+1]=8;
      klcount=klcount+1;
      }
   if (c==2167) {
      out[2*klcount]=4;
      out[2*klcount+1]=12;
      klcount=klcount+1;
      }
   if (c==8497) {
      out[2*klcount]=4;
      out[2*klcount+1]=15;
      klcount=klcount+1;
      }
   if (c==16133) {
      out[2*klcount]=12;
      out[2*klcount+1]=24;
      klcount=klcount+1;
      }
   if (c==29267) {
      out[2*klcount]=4;
      out[2*klcount+1]=16;
      klcount=klcount+1;
      }
   if (c==53095) {
      out[2*klcount]=4;
      out[2*klcount+1]=16;
      klcount=klcount+1;
      }
   if (c==66469) {
      out[2*klcount]=3;
      out[2*klcount+1]=13;
      klcount=klcount+1;
      }
   if (c==78313) {
      out[2*klcount]=5;
      out[2*klcount+1]=19;
      klcount=klcount+1;
      }
   if (c==117599) {
      out[2*klcount]=3;
      out[2*klcount+1]=13;
      klcount=klcount+1;
      }
   if (c==175559) {
      out[2*klcount]=5;
      out[2*klcount+1]=18;
      klcount=klcount+1;
      }
//
   ncyco[klcount]=ncyco[klcount]+1;
   first=0;
   for (k=0; k<klcount; k++) {	     // histogram L-K values
      L1=out[2*k];
      K1=out[2*k+1];
//
// check for L1=1 or K1=1
//
      if ((L1==1)||(K1==1)) {
	 i=test(c, (unsigned int)L1, (unsigned int)K1, 1024, 0);
	 if (i==2) {
	    printf("division error \n");
	    goto zskip;
	    }
	 if ((i==0)&&(flag10==0)) {
	    printf("error: c not equal to 2^(K+L)-3^K, c=%d, L=%d, K=%d \n",c,L1,K1);
	    goto zskip;
	    }
	 }
      if (flag8!=0) {		     // division test flag
	 if (c<200000)		     // temporary!
	    goto xjump;
	 if ((c==111611)||(c==176123)||(c==180833))
	    i=test(c, (unsigned int)L1, (unsigned int)K1, 2048, 1);
	 else {
	    if ((c==99391)||(c==110273)||(c==121189)||(c==122323)||(c==176023)||(c==177019)||(c==177131)||(c==185357)||(c==186689)||(c==194839)||(c==196547))
	       i=test(c, (unsigned int)L1, (unsigned int)K1, 1024, 1);
	    else
	       i=test(c, (unsigned int)L1, (unsigned int)K1, 512, 1);
	    }
	 if (i==0) {
	    printf("error: c does not divide 2^(K+L)-3^K: c=%d L=%d K=%d \n",c,L1,K1);
	    fprintf(Outfp,"error: c does not divide 2^(K+L)-3^K: c=%d L=%d K=%d \n",c,L1,K1);
	    goto zskip;
	    }
	 }
xjump:
      delta=L1-K1;
      if (flag6!=0) {		     // divide by 6 flag
	 delta=delta/6;
	 if (flag7!=0) {	     // rounding flag
	    j=(L1-K1)-delta*6;
	    if (j<0)
	       j=-j;
	    if (delta<0) {
	       if (j>3)
		  delta=delta-1;
	       if (j==3) {
		  if (parity4==1)
		     delta=delta-1;
		  parity4=parity4^1;
		  }
	       }
	    if (delta>0) {
	       if (j>3)
		  delta=delta+1;
	       if (j==3) {
		  if (parity4==1)
		     delta=delta+1;
		  parity4=parity4^1;
		  }
	       }
	    }
	 }
      if (((offset4+delta)<0)||((offset4+delta)>899)) {
	 printf("array not big enough (13), delta+offset4=%d \n",delta+offset4);
	 goto zskip;
	 }
      else
	 histo13[delta+offset4]=histo13[delta+offset4]+1;
//
      delta=L1-K1;		     // count small |L-K| values
      if (delta<0)
	 delta=-delta;
      if ((first==0)&&(delta<=bound)) {
	 smallcnt=smallcnt+1;
	 first=1;
	 }
      }
//
   for (k=0; k<klcount; k++) {	    // sort (L,K) values
      L1=out[2*k];
      K1=out[2*k+1];
      index=k;
      min=L1;
      mink=K1;
      for (i=k+1; i<klcount; i++) {
	 L2=out[2*i];
	 K2=out[2*i+1];
	 if (L2==min) {
	    if (K2<mink) {
	       index=i;
	       mink=K2;
	       }
	    }
	 if (L2<min) {
	    index=i;
	    min=L2;
	    mink=K2;
	    }
	 }
      if (index!=k) {
	 out[2*k]=out[2*index];
	 out[2*k+1]=out[2*index+1];
	 out[2*index]=L1;
	 out[2*index+1]=K1;
	 }
      }
//
   if (klcount==select) {	    // histogram L-K values (smallest L)
      L1=out[0];
      K1=out[1];
      delta=(L1-K1)/6;
      if (flag7!=0) {		    // rounding flag
	 j=(L1-K1)-delta*6;
	 if (j<0)
	    j=-j;
	 if (delta<0) {
	    if (j>3)
	       delta=delta-1;
	    if (j==3) {
	       if (parity1==1)
		  delta=delta-1;
	       parity1=parity1^1;
	       }
	    }
	 if (delta>0) {
	    if (j>3)
	       delta=delta+1;
	    if (j==3) {
	       if (parity1==1)
		  delta=delta+1;
	       parity1=parity1^1;
	       }
	    }
	 }
      if (((delta+offset)<0)||((delta+offset)>199)) {
	 printf("error: array not big enough (10) \n");
	 goto zskip;
	 }
      else
	 histo10[delta+offset]=histo10[delta+offset]+1;
//
      if (klcount>1) {		    // histogram L-K values (next smallest L)
	 L1=out[2];
	 K1=out[3];
	 delta=(L1-K1)/6;
	 if (flag7!=0) {	    // rounding flag
	    j=(L1-K1)-delta*6;
	    if (j<0)
	       j=-j;
	    if (delta<0) {
	       if (j>3)
		  delta=delta-1;
	       if (j==3) {
		  if (parity2==1)
		     delta=delta-1;
		  parity2=parity2^1;
		  }
	       }
	    if (delta>0) {
	       if (j>3)
		  delta=delta+1;
	       if (j==3) {
		  if (parity2==1)
		     delta=delta+1;
		  parity2=parity2^1;
		  }
	       }
	    }
	 if (((delta+offset)<0)||((delta+offset)>199)) {
	    printf("error: array not big enough (11) \n");
	    goto zskip;
	    }
	 else
	    histo11[delta+offset]=histo11[delta+offset]+1;
	 }
//
      L1=out[2*(klcount-1)];	    // histogram L-K values (largest L)
      K1=out[2*(klcount-1)+1];
      delta=(L1-K1)/6;
      if (flag7!=0) {		    // rounding flag
	 j=(L1-K1)-delta*6;
	 if (j<0)
	    j=-j;
	 if (delta<0) {
	    if (j>3)
	       delta=delta-1;
	    if (j==3) {
	       if (parity3==1)
		  delta=delta-1;
	       parity3=parity3^1;
	       }
	    }
	 if (delta>0) {
	    if (j>3)
	       delta=delta+1;
	    if (j==3) {
	       if (parity3==1)
		  delta=delta+1;
	       parity3=parity3^1;
	       }
	    }
	 }
      if (((delta+offset)<0)||((delta+offset)>199)) {
	 printf("error: array not big enough (12) \n");
	 goto zskip;
	 }
      else
	 histo12[delta+offset]=histo12[delta+offset]+1;
      }
//
   if (flag10!=0)
      start=swap;
   else
      start=0;
   mindif1=0x7fffffff;
   for (k=start; k<klcount; k++) {   // find minimum |L-K| value
      L1=out[2*k];
      K1=out[2*k+1];
      delta=L1-K1;
      if (delta<0)
	 delta=-delta;
      if (delta<mindif1) {
	 mindif1=delta;
	 L4=L1;
	 K4=K1;
	 }
      }
//
   kount=0;			    // find multiples of (L,K)
   flag1=0;			    // clear multiple flag
   for (k=0; k<klcount; k++)
      flag[k]=1;
   for (k=0; k<klcount; k++) {
      if (flag[k]==0)
	 continue;
      L1=out[2*k];
      if (L1==0)
	 continue;
      K1=out[2*k+1];
      for (i=(k+1); i<klcount; i++) {
	 L2=out[2*i];
	 K2=out[2*i+1];
	 if ((L2==(L2/L1)*L1)&&(K2==(K2/K1)*K1)) {
	    if ((L2/L1)==(K2/K1)) {
	       kount=kount+1;
	       flag[i]=0;
	       if (tflag!=0)	    // set multiple flag
		  flag1=1;
	       }
	    }
	 }
      }
//
   cnt=0;			    // find arithmetic progressions
   L1=out[0];
   if (L1==0) {
      if (swap==0) {
	 cnt=1;
	 flag1=1;		    // set multiple flag
	 flag2=1;		    // set associated cycle flag
	 mindif=0;		    // minimum |L-K| in arithmetic progressions
	 }
      goto askip;
      }
   K1=out[1];
   if (swap!=0) {		    // swap (L,K) values
      L1=out[2*swap];
      K1=out[2*swap+1];
      }
   if (L1>K1)			    // set comparison flags
      flag3=1;
   else {
      if (L1<K1)
	 flag3=0;
      else
	 flag3=2;
      }
   first=0;
   mindif=0x7fffffff;
   for (k=(1+swap); k<(klcount-1); k++) {
      L2=out[2*k];
      K2=out[2*k+1];
      if ((((L2==(L2/L1)*L1)&&(K2==(K2/K1)*K1)))&&((L2/L1)==(K2/K1)))
	 continue;
      for (i=(k+1); i<(klcount); i++) {
	 L3=out[2*i];
	 K3=out[2*i+1];
	 dL=L3-L2;
	 dK=K3-K2;
	 if ((((dL==(dL/L1)*L1)&&(dK==(dK/K1)*K1)))&&((dL/L1)==(dK/K1))) {
	    if ((dL/L1)<=limit) {
	       cnt=cnt+1;		    // increment count
	       flag1=1; 		    // set multiple flag
	       flag2=1; 		    // set associated cycle flag
	       if (first==0) {
		  if ((flag3==1)&&(L2>K2)) {
		     printf("L1>K1, L2>K2: c=%d \n",c);
		     if (owrite!=0)
			fprintf(Outfp,"L1>K1, L2>K2: c=%d \n",c);
		     mismatch=mismatch+1;
		     }
		  if ((flag3==0)&&(L2<K2)) {
		     printf("L1<K1, L2<K2: c=%d \n",c);
		     if (owrite!=0)
			fprintf(Outfp,"L1<K1, L2<K2: c=%d \n",c);
		     mismatch=mismatch+1;
		     }
		  first=1;
		  }
	       delta=L2-K2;    // minimum |L-K| in arithmetic progressions
	       if (delta<0)
		  delta=-delta;
	       if (delta<mindif) {
		  mindif=delta;
		  L5=L2;
		  K5=K2;
		  }
	       delta=L3-K3;
	       if (delta<0)
		  delta=-delta;
	       if (delta<mindif) {
		  mindif=delta;
		  L5=L3;
		  K5=K3;
		  }
	       break;
	       }
	    }
	 }
      }
   if (cnt!=0) {
      delta=L1-K1;	       // include |L1-K1| in comparisons
      if (delta<0)
	 delta=-delta;
      if (delta<mindif) {
	 mindif=delta;
	 L5=L1;
	 K5=K1;
	 }
      }
//
askip:
   if (cnt!=0) {
      if (((L5-K5+offset)<0)||((L5-K5+offset)>199)) {	// histogram smallest L-K
	 printf("error: array not big enough (14) \n");      // in arithmetic progressions
	 goto zskip;
	 }
      histo14[offset+L5-K5]=histo14[offset+L5-K5]+1;
      }
//
   if ((flag5!=0)||(cnt!=0)) {	    // histogram smallest L-K values
      if ((cnt!=0)&&(flag9!=0))
	 goto bskip;
      if (mindif1==0x7fffffff)
	 goto bskip;
      if (((L4-K4+offset4)<0)||((L4-K4+offset4)>899)) {
	 printf("error: array not big enough (15) \n");
	 printf("c=%d, L4=%d, K4=%d \n",c,L4,K4);
	 goto zskip;
	 }
      if(klcount==1)
	 histo1[offset4+L4-K4]=histo1[offset4+L4-K4]+1;
      if(klcount==2)
	 histo2[offset4+L4-K4]=histo2[offset4+L4-K4]+1;
      if(klcount==3)
	 histo3[offset4+L4-K4]=histo3[offset4+L4-K4]+1;
      if(klcount==4)
	 histo4[offset4+L4-K4]=histo4[offset4+L4-K4]+1;
      if(klcount==5)
	 histo5[offset4+L4-K4]=histo5[offset4+L4-K4]+1;
      if(klcount==6)
	 histo6[offset4+L4-K4]=histo6[offset4+L4-K4]+1;
      if(klcount==7)
	 histo7[offset4+L4-K4]=histo7[offset4+L4-K4]+1;
      if(klcount==8)
	 histo8[offset4+L4-K4]=histo8[offset4+L4-K4]+1;
      if(klcount>8)
	 histo9[offset4+L4-K4]=histo9[offset4+L4-K4]+1;
      }
//
bskip:
   if (cnt==0) {
      if (mindif1!=0x7fffffff)
	 mindif=mindif1;
      else
	 mindif=-1;
      }
   else {
      if ((mindif!=mindif1)&&(mindif1!=0x7fffffff)) {
	 if (flag11==1) {
	    delta=out[0]-out[1];
	    if (delta<0)
	       delta=-delta;
	    if (delta!=mindif1) {
	       printf("c=%d, mindif=%d, mindif1=%d \n",c,mindif,mindif1);
	       if (owrite!=0)
		  fprintf(Outfp,"c=%d, mindif=%d, mindif1=%d \n",c,mindif,mindif1);
	       }
	    }
	 else {
	    printf("c=%d, mindif=%d, mindif1=%d \n",c,mindif,mindif1);
	    if (owrite!=0)
	       fprintf(Outfp,"c=%d, mindif=%d, mindif1=%d \n",c,mindif,mindif1);
	    }
	 }
      }
   if ((flag1!=0)||(flag4!=0)) {    // output cycles
      printf(" c=%d counts=%d %d min=%d",c,kount,cnt,mindif);
      if (owrite!=0)
	 fprintf(Outfp," c=%d, counts=%d %d min=%d",c,kount,cnt,mindif);
      for (k=0; k<klcount; k++) {
	 printf(" (%d,%d)",out[2*k],out[2*k+1]);
	 if (owrite!=0)
	    fprintf(Outfp," (%d,%d)",out[2*k],out[2*k+1]);
	 }
      printf("\n");
      if (owrite!=0)
	 fprintf(Outfp,"\n");
      }
   if ((klcount==3)&&(cnt==0)&&(kount==0)) {  // count switched (L1,K1) (L2,K2)
      L1=out[0];
      K1=out[1];
      L2=out[4];
      K2=out[5];
      dL=L2-L1;
      dK=K2-K1;
      if ((dL==(out[2]*2))&&(dK==(out[3]*2)))
	 swit=swit+1;
      }
   if (klcount==3) {
      dL=out[0]-out[2];
      dK=out[1]-out[3];
      if ((dL==(out[2]-out[4]))&&(dK==(out[3]-out[5]))) {
	 if (out[2]==(out[2]/out[0])*out[0])
	    goto eskip;
	 delta=out[0]-out[1];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[2]-out[3];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[4]-out[5];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[0];
	 K6=out[1];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[2];
	    K6=out[3];
	    }
	 if (delta>delta2) {
	    L6=out[4];
	    K6=out[5];
	    }
	 if (((L6-K6+offset2)<0)||((L6-K6+offset2)>99)) {
	    printf("error: array not big enough (16), c=%d \n",c);
	    goto zskip;
	    }
	 histo15[L6-K6+offset2]=histo15[L6-K6+offset2]+1;
	 }
      }
eskip:
   if (klcount==4) {
      dL=out[0]-out[2];
      dK=out[1]-out[3];
      if ((dL==(out[2]-out[4]))&&(dK==(out[3]-out[5]))) {
	 if (out[2]==(out[2]/out[0])*out[0])
	    goto fskip;
	 delta=out[0]-out[1];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[2]-out[3];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[4]-out[5];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[0];
	 K6=out[1];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[2];
	    K6=out[3];
	    }
	 if (delta>delta2) {
	    L6=out[4];
	    K6=out[5];
	    }
	 if (((L6-K6+offset3)<0)||((L6-K6+offset3)>299)) {
	    printf("error: array not big enough (17), c=%d \n",c);
	    goto zskip;
	    }
	 histo16[L6-K6+offset3]=histo16[L6-K6+offset3]+1;
	 }
fskip:
      dL=out[0]-out[4];
      dK=out[1]-out[5];
      if ((dL==(out[4]-out[6]))&&(dK==(out[5]-out[7]))) {
	 if (out[4]==(out[4]/out[0])*out[0])
	    goto gskip;
	 delta=out[0]-out[1];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[4]-out[5];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[6]-out[7];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[0];
	 K6=out[1];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[4];
	    K6=out[5];
	    }
	 if (delta>delta2) {
	    L6=out[6];
	    K6=out[7];
	    }
	 if (((L6-K6+offset3)<0)||((L6-K6+offset3)>299)) {
	    printf("error: array not big enough (18), c=%d \n",c);
	    goto zskip;
	    }
	 histo16[L6-K6+offset3]=histo16[L6-K6+offset3]+1;
	 }
gskip:
      dL=out[2]-out[4];
      dK=out[3]-out[5];
      if ((dL==(out[4]-out[6]))&&(dK==(out[5]-out[7]))) {
	 if (out[4]==(out[4]/out[2])*out[2])
	    goto hskip;
	 delta=out[2]-out[3];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[4]-out[5];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[6]-out[7];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[2];
	 K6=out[3];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[4];
	    K6=out[5];
	    }
	 if (delta>delta2) {
	    L6=out[6];
	    K6=out[7];
	    }
	 if (((L6-K6+offset3)<0)||((L6-K6+offset3)>299)) {
	    printf("error: array not big enough (19), c=%d \n",c);
	    goto zskip;
	    }
	 histo16[L6-K6+offset3]=histo16[L6-K6+offset3]+1;
	 }
//
//  two in arithmetic progression
//
hskip:
      dL=out[2]-out[4];
      dK=out[3]-out[5];
      if ((-dL==out[0])&&(-dK==out[1])) {
	 if (out[0]==0) {
	    L6=1;
	    K6=1;
	    goto hjump;
	    }
	 if (out[2]==(out[2]/out[0])*out[0])
	    goto iskip;
	 delta=out[0]-out[1];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[2]-out[3];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[4]-out[5];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[0];
	 K6=out[1];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[2];
	    K6=out[3];
	    }
	 if (delta>delta2) {
	    L6=out[4];
	    K6=out[5];
	    }
	 if (((L6-K6+offset3)<0)||((L6-K6+offset3)>199)) {
	    printf("error: array not big enough (20), c=%d \n",c);
	    goto zskip;
	    }
hjump:	 histo17[L6-K6+offset3]=histo17[L6-K6+offset3]+1;
	 }
iskip:
      dL=out[4]-out[6];
      dK=out[5]-out[7];
      if ((-dL==out[0])&&(-dK==out[1])) {
	 if (out[0]==0)
	    goto jskip;
	 if (out[4]==(out[4]/out[0])*out[0])
	    goto jskip;
	 delta=out[0]-out[1];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[4]-out[5];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[6]-out[7];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[0];
	 K6=out[1];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[4];
	    K6=out[5];
	    }
	 if (delta>delta2) {
	    L6=out[6];
	    K6=out[7];
	    }
	 if (((L6-K6+offset3)<0)||((L6-K6+offset3)>199)) {
	    printf("error: array not big enough (21), c=%d \n",c);
	    goto zskip;
	    }
	 histo17[L6-K6+offset3]=histo17[L6-K6+offset3]+1;
	 }
jskip:
      dL=out[4]-out[6];
      dK=out[5]-out[7];
      if ((-dL==out[2])&&(-dK==out[3])) {
	 if (out[4]==(out[4]/out[2])*out[2])
	    goto kskip;
	 delta=out[2]-out[3];
	 if (delta<0)
	    delta=-delta;
	 delta1=out[4]-out[5];
	 if (delta1<0)
	    delta1=-delta1;
	 delta2=out[6]-out[7];
	 if (delta2<0)
	    delta2=-delta2;
	 L6=out[2];
	 K6=out[3];
	 if (delta>delta1) {
	    delta=delta1;
	    L6=out[4];
	    K6=out[5];
	    }
	 if (delta>delta2) {
	    L6=out[6];
	    K6=out[7];
	    }
	 if (((L6-K6+offset3)<0)||((L6-K6+offset3)>199)) {
	    printf("error: array not big enough (22), c=%d \n",c);
	    goto zskip;
	    }
	 histo17[L6-K6+offset3]=histo17[L6-K6+offset3]+1;
	 }
      }
kskip:
   if (flag2==1)			      // increment associated cycle count
      iters=iters+1;
   if (tflag!=0) {			      // flag for not counting multiples
      ncyc[klcount-kount-cnt]=ncyc[klcount-kount-cnt]+1;
      if (((klcount-kount-cnt)>6)&&(flag1==0)&&(flag4==0)) {
	 printf("large x value: c=%d, x=%d \n",klcount-kount-cnt-1);
	 fprintf(Outfp,"large x value: c=%d, x=%d \n",klcount-kount-cnt-1);
	 for (k=0; k<klcount; k++) {
	    printf(" (%d,%d)",out[2*k],out[2*k+1]);
	    fprintf(Outfp," (%d,%d)",out[2*k],out[2*k+1]);
	    }
	 printf("\n");
	 fprintf(Outfp,"\n");
	 }
      }
   else {
      ncyc[klcount-cnt]=ncyc[klcount-cnt]+1;  // don't count associated cycles
      if (((klcount-cnt)>6)&&(flag1==0)&&(flag4==0)) {
	 printf("large x value: c=%d, x=%d \n",c,klcount-cnt-1);
	 fprintf(Outfp,"large x value: c=%d, x=%d \n",c,klcount-cnt-1);
	 for (k=0; k<klcount; k++) {
	    printf(" (%d,%d)",out[2*k],out[2*k+1]);
	    fprintf(Outfp," (%d,%d)",out[2*k],out[2*k+1]);
	    }
	 printf("\n");
	 fprintf(Outfp,"\n");
	 }
      }
   if (c==1)
      klcount=klcount-3;
   if (c==11)
      klcount=klcount-1;
   if (c==49)
      klcount=klcount-1;
   if (c==179)
      klcount=klcount-1;
   if (c==601)
      klcount=klcount-1;
   if (c==1931)
      klcount=klcount-1;
   if (c==6049)
      klcount=klcount-1;
   if (c==18659)
      klcount=klcount-1;
   if (c==57001)
      klcount=klcount-1;
   if (c==173051)
      klcount=klcount-1;
   if (c==85)
      klcount=klcount-1;
   if (c==145)
      klcount=klcount-1;
   if (c==791)
      klcount=klcount-1;
   if (c==2167)
      klcount=klcount-1;
   if (c==8497)
      klcount=klcount-1;
   if (c==16133)
      klcount=klcount-1;
   if (c==29267)
      klcount=klcount-1;
   if (c==53095)
      klcount=klcount-1;
   if (c==66469)
      klcount=klcount-1;
   if (c==78313)
      klcount=klcount-1;
   if (c==117599)
      klcount=klcount-1;
   if (c==175559)
      klcount=klcount-1;
   if (c<70000)
      total=total+klcount;
   else {
      if (c<100000)
	 total1=total1+klcount;
      else {
	 if (c<150000)
	    total2=total2+klcount;
	 else
	    total3=total3+klcount;
	 }
      }
   c=c+2;
   if (c==(c/3)*3)
      c=c+2;
   }
//
printf("\n");                       // output histogram
fprintf(Outfp,"\n");
printf("actual frequency distribution (after adjusting cycle counts for associated cycles) \n");
fprintf(Outfp,"actual frequency distribution (after adjusting cycle counts for associated cycles) \n");
kount=0;
for (i=1; i<20; i++)
   kount=kount+ncyc[i];
index=0;
for (i=1; i<20; i++) {
   printf("i=%d, n=%d, f=%e \n",i-1,ncyc[i],(double)ncyc[i]/(double)kount);
   fprintf(Outfp,"i=%d, n=%d, f=%e \n",i-1,ncyc[i],(double)ncyc[i]/(double)kount);
   index=index+(i-1)*ncyc[i];
   }
lambda=(double)index/(double)kount;
printf("lambda=%e \n",lambda);
fprintf(Outfp,"lambda=%e \n",lambda);
printf("\n");
fprintf(Outfp,"\n");                // expected distribution
printf("expected frequency distribution \n");
fprintf(Outfp,"expected frequency distribution \n");
ratio=exp(-lambda);
for (i=1; i<20; i++) {
   printf("i=%d, n=%d, f=%e \n",i-1,(int)(ratio*(double)kount+0.5),ratio);
   fprintf(Outfp,"i=%d, n=%d, f=%e \n",i-1,(int)(ratio*(double)kount+0.5),ratio);
   ratio=ratio*lambda;
   if (i!=1)
      ratio=ratio/(double)i;
   }
printf("\n");
fprintf(Outfp,"\n");
printf("number of c values having associated cycles=%d\n",iters);
fprintf(Outfp,"number of c values having associated cycles=%d\n",iters);
printf("number of c values where L1>K1 and L2>K2 or L1<K1 and L2<K2=%d \n",mismatch);
fprintf(Outfp,"number of c values where L1>K1 and L2>K2 or L1<K1 and L2<K2=%d \n",mismatch);
printf("number of c values where |L-K|<=bound: %d \n",smallcnt);
fprintf(Outfp,"number of c values where |L-K|<=bound: %d \n",smallcnt);
printf("number of c values with 3 cycles and switched (L1,K1) and (L2,K2)=%d \n",swit);
fprintf(Outfp,"number of c values with 3 cycles and switched (L1,K1) and (L2,K2)=%d \n",swit);
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=1 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=1 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo1[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo1[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo1[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo1[i]);
   j=j+histo1[i];
   sum=sum+(i-offset4)*histo1[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo1[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=2 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=2 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo2[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo2[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo2[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo2[i]);
   j=j+histo2[i];
   sum=sum+(i-offset4)*histo2[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo2[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L=K|), cycle count=3 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=3 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo3[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo3[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo3[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo3[i]);
   j=j+histo3[i];
   sum=sum+(i-offset4)*histo3[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo3[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=4 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=4 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo4[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo4[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo4[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo4[i]);
   j=j+histo4[i];
   sum=sum+(i-offset4)*histo4[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo4[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=5 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=5 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo5[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo5[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo5[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo5[i]);
   j=j+histo5[i];
   sum=sum+(i-offset4)*histo5[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo5[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=6 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=6 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo6[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo6[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo6[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo6[i]);
   j=j+histo6[i];
   sum=sum+(i-offset4)*histo6[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo6[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=7 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=7 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo7[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo7[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo7[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo7[i]);
   j=j+histo7[i];
   sum=sum+(i-offset4)*histo7[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo7[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count=8 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count=8 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo8[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo8[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo8[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo8[i]);
   j=j+histo8[i];
   sum=sum+(i-offset4)*histo8[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo8[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values (for smallest |L-K|), cycle count>8 \n");
fprintf(Outfp,"histogram of L-K values (for smallest |L-K|), cycle count>8 \n");
index=899;
for (i=899; i>0; i--) {
   if (histo9[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo9[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo9[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo9[i]);
   j=j+histo9[i];
   sum=sum+(i-offset4)*histo9[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo9[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values, minimum L value, cycle count=%d \n",select);
fprintf(Outfp,"histogram of L-K values, minimum L value, cycle count=%d \n",select);
index=199;
for (i=199; i>0; i--) {
   if (histo10[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo10[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset,histo10[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset,histo10[i]);
   j=j+histo10[i];
   sum=sum+(i-offset)*histo10[i];
   sum2=sum2+(i-offset)*(i-offset)*histo10[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
if (select!=1) {
   printf("\n");
   fprintf(Outfp,"\n");
   printf("histogram of L-K values, next smallest L value, cycle count=%d \n",select);
   fprintf(Outfp,"histogram of L-K values, next smallest L value, cycle count=%d \n",select);
   index=199;
   for (i=199; i>0; i--) {
      if (histo11[i]==0)
	 index=index-1;
      else
	 break;
      }
   first=0;
   j=0;
   sum=0;
   sum2=0;
   for (i=0; i<=index; i++) {
      if ((histo11[i]==0)&&(first==0))
	 continue;
      first=1;
      printf("i=%d, %d \n",i-offset,histo11[i]);
      fprintf(Outfp,"i=%d, %d \n",i-offset,histo11[i]);
      j=j+histo11[i];
      sum=sum+(i-offset)*histo11[i];
      sum2=sum2+(i-offset)*(i-offset)*histo11[i];
      }
   mean=(double)sum/(double)j;
   var=(double)sum2*(double)j-(double)sum*(double)sum;
   var=var/(double)j/(double)(j-1);
   printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
   fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
   }
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values, maximum L value, cycle count=%d \n",select);
fprintf(Outfp,"histogram of L-K values, maximum L value, cycle count=%d \n",select);
index=199;
for (i=199; i>0; i--) {
   if (histo12[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo12[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset,histo12[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset,histo12[i]);
   j=j+histo12[i];
   sum=sum+(i-offset)*histo12[i];
   sum2=sum2+(i-offset)*(i-offset)*histo12[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values, minimum |L-K| in arithmetic progressions \n");
fprintf(Outfp,"histogram of L-K values, minimum |L-K| in arithmetic progressions \n");
index=199;
for (i=199; i>0; i--) {
   if (histo14[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo14[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset,histo14[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset,histo14[i]);
   j=j+histo14[i];
   sum=sum+(i-offset)*histo14[i];
   sum2=sum2+(i-offset)*(i-offset)*histo14[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values \n");
fprintf(Outfp,"histogram of L-K values\n");
index=899;
for (i=899; i>0; i--) {
   if (histo13[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo13[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset4,histo13[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset4,histo13[i]);
   j=j+histo13[i];
   sum=sum+(i-offset4)*histo13[i];
   sum2=sum2+(i-offset4)*(i-offset4)*histo13[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values, cycle count=3, arithmetic progression \n");
fprintf(Outfp,"histogram of L-K values, cycle count=3, arithmetic progression \n");
index=99;
for (i=99; i>0; i--) {
   if (histo15[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo15[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset2,histo15[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset2,histo15[i]);
   j=j+histo15[i];
   sum=sum+(i-offset2)*histo15[i];
   sum2=sum2+(i-offset2)*(i-offset2)*histo15[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values, cycle count=4, arithmetic progression \n");
fprintf(Outfp,"histogram of L-K values, cycle count=4, arithmetic progression \n");
index=299;
for (i=299; i>0; i--) {
   if (histo16[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo16[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset3,histo16[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset3,histo16[i]);
   j=j+histo16[i];
   sum=sum+(i-offset3)*histo16[i];
   sum2=sum2+(i-offset3)*(i-offset3)*histo16[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of L-K values, cycle count=4, arithmetic progression \n");
fprintf(Outfp,"histogram of L-K values, cycle count=4, arithmetic progression \n");
index=199;
for (i=199; i>0; i--) {
   if (histo17[i]==0)
      index=index-1;
   else
      break;
   }
first=0;
j=0;
sum=0;
sum2=0;
for (i=0; i<=index; i++) {
   if ((histo17[i]==0)&&(first==0))
      continue;
   first=1;
   printf("i=%d, %d \n",i-offset3,histo17[i]);
   fprintf(Outfp,"i=%d, %d \n",i-offset3,histo17[i]);
   j=j+histo17[i];
   sum=sum+(i-offset3)*histo17[i];
   sum2=sum2+(i-offset3)*(i-offset3)*histo17[i];
   }
mean=(double)sum/(double)j;
var=(double)sum2*(double)j-(double)sum*(double)sum;
var=var/(double)j/(double)(j-1);
printf("n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
fprintf(Outfp,"n=%d, mean=%e, std=%e \n",j,mean,sqrt(var));
//
//
// histogram of cycle counts (unmodified)
//
printf("\n");
fprintf(Outfp,"\n");
printf("histogram of cycle counts (unmodified) \n");
fprintf(Outfp,"histogram of cycle counts (unmodified) \n");
sum=0;
for (i=0; i<20; i++)
   sum=sum+ncyco[i];
for (i=0; i<20; i++) {
   printf("i=%d, n=%d, r=%e \n",i,ncyco[i],(double)ncyco[i]/(double)sum);
   fprintf(Outfp,"i=%d, n=%d, r%e \n",i,ncyco[i],(double)ncyco[i]/(double)sum);
   }
zskip:
fclose(Outfp);
return(0);
}