/******************************************************************************/
/* */
/* 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);
}