/*****************************************************************************/
/*									     */
/*  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**33.   */
/*									     */
/*****************************************************************************/
#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 ho[128];
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 h[4*127]={
0x00000000,
0x00000000,
0x0AAAAAAA,
0x00000000,
0x02AAAAAB,
0x00000000,
0x00000000,
0x004BDA13,
0x00000000,
0x002F684C,
0x00000000,
0x00000000,
0x001A5663,
0x00000000,
0x00036C82,
0x00000000,
0x00000000,
0x0007ABB7,
0x00000000,
0x00000000,
0x0000D8F0,
0x00000000,
0x0001B867,
0x00000000,
0x00000000,
0x0000C53A,
0x00000000,
0x000049A6,
0x00000000,
0x00000000,
0x00006999,
0x00000000,
0x00000000,
0x00000260,
0x00000000,
0x00001E9A,
0x00000000,
0x00000000,
0x000008B7,
0x00000000,
0x000006CB,
0x00000000,
0x00000000,
0x00000669,
0x00000000,
0x000000A4,
0x00000000,
0x00000000,
0x0000027E,
0x00000000,
0x00000000,
0x0000006B,
0x00000000,
0x000000AF,
0x00000000,
0x00000000,
0x00000065,
0x00000000,
0x00000023,
0x00000000,
0x00000000,
0x0000003D,
0x00000000,
0x00000000,
0x00000002,
0x00000000,
0x00000015,
0x00000000,
0x00000000,
0x0000000B,
0x00000000,
0x00000006,
0x00000000,
0x00000000,
0x00000003,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000001,
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,
0x0AAAAAAB,
0x00000000,
0x02AAAAAB,
0x00000000,
0x00000000,
0x004BDA13,
0x00000000,
0x002F684C,
0x00000000,
0x00000000,
0x001A5664,
0x00000000,
0x00036C83,
0x00000000,
0x00000000,
0x0007ABB9,
0x00000000,
0x00000000,
0x0000D8F8,
0x00000000,
0x0001B868,
0x00000000,
0x00000000,
0x0000C537,
0x00000000,
0x000049A4,
0x00000000,
0x00000000,
0x00006991,
0x00000000,
0x00000000,
0x00000262,
0x00000000,
0x00001E9A,
0x00000000,
0x00000000,
0x000008B7,
0x00000000,
0x0000075F,
0x00000000,
0x00000000,
0x00000668,
0x00000000,
0x0000009A,
0x00000000,
0x00000000,
0x00000298,
0x00000000,
0x00000000,
0x0000006D,
0x00000000,
0x000000AE,
0x00000000,
0x00000000,
0x00000069,
0x00000000,
0x00000028,
0x00000000,
0x00000000,
0x00000037,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000011,
0x00000000,
0x00000000,
0x0000000B,
0x00000000,
0x00000009,
0x00000000,
0x00000000,
0x00000002,
0x00000000,
0x00000001,
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,
0x00000001,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000000,

0x00000000,
0x00000000,
0x0AAAAAAB,
0x00000000,
0x02AAAAAB,
0x00000000,
0x00000000,
0x004BDA13,
0x00000000,
0x002F684C,
0x00000000,
0x00000000,
0x001A5662,
0x00000000,
0x00036C83,
0x00000000,
0x00000000,
0x0007ABB9,
0x00000000,
0x00000000,
0x0000D8EE,
0x00000000,
0x0001B868,
0x00000000,
0x00000000,
0x0000C538,
0x00000000,
0x000049A4,
0x00000000,
0x00000000,
0x00006998,
0x00000000,
0x00000000,
0x0000025B,
0x00000000,
0x00001E9D,
0x00000000,
0x00000000,
0x000008CF,
0x00000000,
0x00000731,
0x00000000,
0x00000000,
0x0000065C,
0x00000000,
0x000000AA,
0x00000000,
0x00000000,
0x0000028A,
0x00000000,
0x00000000,
0x00000063,
0x00000000,
0x000000CF,
0x00000000,
0x00000000,
0x00000057,
0x00000000,
0x00000026,
0x00000000,
0x00000000,
0x00000043,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000019,
0x00000000,
0x00000000,
0x0000000B,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000003,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000001,
0x00000000,
0x00000000,
0x00000000,
0x00000000,
0x00000001,
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,
0x0AAAAAAA,
0x00000000,
0x02AAAAAA,
0x00000000,
0x00000000,
0x004BDA13,
0x00000000,
0x002F684B,
0x00000000,
0x00000000,
0x001A5663,
0x00000000,
0x00036C83,
0x00000000,
0x00000000,
0x0007ABB8,
0x00000000,
0x00000000,
0x0000D8F0,
0x00000000,
0x0001B866,
0x00000000,
0x00000000,
0x0000C538,
0x00000000,
0x000049A2,
0x00000000,
0x00000000,
0x00006997,
0x00000000,
0x00000000,
0x0000025E,
0x00000000,
0x00001E9E,
0x00000000,
0x00000000,
0x000008C4,
0x00000000,
0x000006C4,
0x00000000,
0x00000000,
0x00000663,
0x00000000,
0x000000B1,
0x00000000,
0x00000000,
0x00000298,
0x00000000,
0x00000000,
0x00000067,
0x00000000,
0x000000B1,
0x00000000,
0x00000000,
0x00000065,
0x00000000,
0x00000028,
0x00000000,
0x00000000,
0x0000003A,
0x00000000,
0x00000000,
0x00000005,
0x00000000,
0x0000000B,
0x00000000,
0x00000000,
0x00000002,
0x00000000,
0x00000005,
0x00000000,
0x00000000,
0x00000003,
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,
0x00000001,
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("out4a33.dat","w");
S[0]=0;
S[1]=0;
//
// sum up histogram bins
//
for (i=0; i<127; i++) {
   ho[i]=0;
   for (j=0; j<4; j++)
      ho[i]=ho[i]+h[127*j+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]=1;
G[1]=0x00000000;
g=0x80000000;
for (i=1; i<14; i++) {
   r=(double)g;
   r=r*2.0;
   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!=0))		 // 17
      flag=1;
   if ((i==8)&&(k!=-1)) 	 // 20
      flag=1;
   if ((i==9)&&(k!=2))		 // 22
      flag=1;
   if ((i==10)&&(k!=6)) 	 // 25
      flag=1;
   if ((i==11)&&(k!=0)) 	 // 27
      flag=1;
   if ((i==12)&&(k!=4)) 	 // 30
      flag=1;
   if ((i==13)&&(k!=-7))	  // 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);
}