﻿ /*******************************
```/*****************************************************************************/
/*									     */
/*  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**30.   */
/*									     */
/*****************************************************************************/
#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]={
0x00000000,
0x00000000,
0x05555555,
0x00000000,
0x01555555,
0x00000000,
0x00000000,
0x0025ED09,
0x00000000,
0x0017B426,
0x00000000,
0x00000000,
0x000D2B32,
0x00000000,
0x0001B642,
0x00000000,
0x00000000,
0x0003D5DC,
0x00000000,
0x00000000,
0x00006C78,
0x00000000,
0x0000DC33,
0x00000000,
0x00000000,
0x000062A1,
0x00000000,
0x000024D1,
0x00000000,
0x00000000,
0x000034CC,
0x00000000,
0x00000000,
0x0000012C,
0x00000000,
0x00000F50,
0x00000000,
0x00000000,
0x00000463,
0x00000000,
0x00000388,
0x00000000,
0x00000000,
0x0000032E,
0x00000000,
0x0000005A,
0x00000000,
0x00000000,
0x0000014A,
0x00000000,
0x00000000,
0x0000002F,
0x00000000,
0x0000004E,
0x00000000,
0x00000000,
0x0000003A,
0x00000000,
0x0000000D,
0x00000000,
0x00000000,
0x00000020,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000007,
0x00000000,
0x00000000,
0x00000003,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000003,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000};

int k;
unsigned int S[2],T[2];
double r;
FILE *Outfp;
Outfp = fopen("out4a30.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];
}
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]=0x20000000;
g=0x20000000;
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!=0))		 // 2 (k is rounded)
flag=1;
if ((i==2)&&(k!=0))		 // 4
flag=1;
if ((i==3)&&(k!=0))		 // 7
flag=1;
if ((i==4)&&(k!=0))		 // 9
flag=1;
if ((i==5)&&(k!=0))		 // 12
flag=1;
if ((i==6)&&(k!=-1)) 	 // 14
flag=1;
if ((i==7)&&(k!=0))		 // 17
flag=1;
if ((i==8)&&(k!=1))		 // 20
flag=1;
if ((i==9)&&(k!=1))		// 22
flag=1;
if ((i==10)&&(k!=-4))	 // 25
flag=1;
if ((i==11)&&(k!=1)) 	 // 27
flag=1;
if ((i==12)&&(k!=0)) 	 // 30
flag=1;
if ((i==13)&&(k!=2)) 	 // 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);
}
```