/*****************************************************************************/
/*									     */
/*  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**29.   */
/*									     */
/*****************************************************************************/
#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,
0x02AAAAAA,
0x00000000,
0x00AAAAAB,
0x00000000,
0x00000000,
0x0012F685,
0x00000000,
0x000BDA13,
0x00000000,
0x00000000,
0x00069599,
0x00000000,
0x0000DB21,
0x00000000,
0x00000000,
0x0001EAED,
0x00000000,
0x00000000,
0x0000363D,
0x00000000,
0x00006E1A,
0x00000000,
0x00000000,
0x0000314F,
0x00000000,
0x00001267,
0x00000000,
0x00000000,
0x00001A65,
0x00000000,
0x00000000,
0x0000009A,
0x00000000,
0x000007A5,
0x00000000,
0x00000000,
0x00000231,
0x00000000,
0x000001C2,
0x00000000,
0x00000000,
0x0000018E,
0x00000000,
0x0000002F,
0x00000000,
0x00000000,
0x000000AF,
0x00000000,
0x00000000,
0x00000019,
0x00000000,
0x0000002D,
0x00000000,
0x00000000,
0x00000018,
0x00000000,
0x00000009,
0x00000000,
0x00000000,
0x00000011,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000004,
0x00000000,
0x00000003,
0x00000000,
0x00000000,
0x00000003,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000001,
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("out4a29.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]=0x10000000;
g=0x10000000;
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!=0))		 // 7
      flag=1;
   if ((i==4)&&(k!=0))		 // 9
      flag=1;
   if ((i==5)&&(k!=0))		 // 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!=0))		 // 22
      flag=1;
   if ((i==10)&&(k!=-1))	 // 25
      flag=1;
   if ((i==11)&&(k!=2)) 	 // 27
      flag=1;
   if ((i==12)&&(k!=1)) 	 // 30
      flag=1;
   if ((i==13)&&(k!=-3))	  // 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);
}