/*****************************************************************************/
/* */
/* HISTOGRAM OF LIMB LENGTHS */
/* 11/28/09 (dkc) */
/* */
/* This C program compares the expected number of limbs of a given length */
/* to the actual number. The limbs must be alive. The order is 3*2**17. */
/* */
/*****************************************************************************/
#include <stdio.h>
#include <math.h>
void add64(unsigned int *a0, unsigned int *a1);
void mul64_32(unsigned int a0, unsigned int a2, unsigned int m0,
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);
int main () {
unsigned int P[3],Q[4],G[2];
unsigned int v[14]={1,2,4,7,9,12,14,17,20,22,25,27,30,33};
unsigned int y[40]={0,3,6,8,11, // B (invalid lengths)
13,16,19,21,24, // B
26,29,32,34,37, // B
39,42,44,47,50, // C
52,55,57,60,63, // C
65,68,70,73,75, // D
78,81,83,86,88, // D
91,94,96,99,101}; // D
//
// fractions
//
unsigned int z[14*2]={1,6, // 1 1/2*3
1,6, // 2 1/(2**1)*(3**1)
1,24, // 4 1/(2**3)*(3**1)
1,216, // 7 1/(2**3)*(3**3)
5,1728, // 9 5/(2**6)*(3**3)
25,15552, // 12 25/(2**6)*(3**5)
13,62208, // 14 13/(2**8)*(3**5)
233,497664, // 17 233/(2**11)*(3**5)
695,13436928, // 20 695/(2**11)*(3**8)
6349,1024*59049, // 22 6349/(2**10)*(3**10)
3791,80621568, // 25 3791/(2**12)*(3**9)
2831,161243136, // 27 2831/(2**13)*(3**9)
3425,256*531441, // 30 3425/(2**8)*(3**12)
7359,1594323}; // 33 7359/(2**13)*(3**13)
//
// l=33 7359/(2**13)*(3**13) or 345/(2**7)*(3**14) or 8279/(2**10)*(3**15)
// l=35 2979/(2**8)*(3**13) or 8937/(2**8)*(3**14)
// Lots of solutions since histogram value is small.
//
unsigned int g,i,j,flag;
unsigned int ho[127]={
0,0,10922,0,2731,0,0,304,
0,189,0,0,106,0,14,0,
0,30,0,0,4,0,8,0,
0,2,0,2,0,0,1,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0};
int k;
unsigned int S[2],T[2];
double r;
FILE *Outfp;
Outfp = fopen("out4a17.dat","w");
//
// sum up histogram bins
//
S[0]=0;
S[1]=0;
for (i=0; i<127; i++) {
T[0]=0;
T[1]=ho[i];
add64(T, S);
}
printf("sum=%#10x %#10x \n",S[0],S[1]);
printf("\n");
for (i=0; i<40; i++) {
if (ho[y[i]]!=0) {
printf("error: invalid length=%d \n",y[i]);
goto zskip;
}
}
G[0]=0;
G[1]=0x00010000;
g=0x00010000;
for (i=1; i<14; i++) {
r=(double)g;
r=r*(double)z[2*i];
r=r/(double)z[2*i+1];
if (i==13)
r=r/(double)8192;
j=(unsigned int)r;
j=j+1;
mul64_32(G[0], G[1], z[2*i], P);
div128_64(P[0], P[1], P[2], 0, Q, 0, z[2*i+1]);
if (i==13)
div128_64(Q[0], Q[1], Q[2], Q[3], Q, 0, 8192);
if (((j-Q[2])!=0)&&((j-Q[2])!=1)) {
printf("rounding error: j=%#10x, Q=%#10x %#10x %#10x %#10x\n",j,Q[0],Q[1],Q[2],Q[3]);
// goto zskip;
}
k=Q[2]+(Q[3]>>31)-(int)ho[v[i]];
flag=0;
if ((i==1)&&(k!=1)) // 2 (k is rounded)
flag=1;
if ((i==2)&&(k!=0)) // 4
flag=1;
if ((i==3)&&(k!=-1)) // 7
flag=1;
if ((i==4)&&(k!=1)) // 9
flag=1;
if ((i==5)&&(k!=-1)) // 12
flag=1;
if ((i==6)&&(k!=0)) // 14
flag=1;
if ((i==7)&&(k!=1)) // 17
flag=1;
if ((i==8)&&(k!=-1)) // 20
flag=1;
if ((i==9)&&(k!=-1)) // 22
flag=1;
if ((i==10)&&(k!=1)) // 25
flag=1;
if ((i==11)&&(k!=-1)) // 27
flag=1;
if ((i==12)&&(k!=1)) // 30
flag=1;
if ((i==13)&&(k!=0)) // 33
flag=1;
if (flag!=0) {
printf("error: j=%d, h=%d, v=%d, i=%d, r=%f \n",Q[2]+(Q[3]>>31),ho[v[i]],v[i],i,r);
// goto zskip;
}
else {
printf("l=%d, j=%d, h=%d, r=%d, d=%d \n",v[i],Q[2]+(Q[3]>>31),ho[v[i]],Q[3]>>31,Q[2]+(Q[3]>>31)-ho[v[i]]);
}
}
zskip:
printf("\n");
for (j=0; j<100; j++) {
if (ho[j]!=0)
printf(" %d %d \n",j,ho[j]);
}
fclose(Outfp);
return(0);
}