/******************************************************************************/
/* */
/* COMPUTE 3N+1 OR 3N-1 SEQUENCE (tree for k<=26, more efficient algorithm) */
/* 01/18/10 (dkc) */
/* */
/* This C program generates sequence vectors of live limbs. The error in */
/* approximating z is computed. The maximum and minimum y values are */
/* computed for sequence vectors. The maximum and minimum z values are */
/* also computed for sequence vectors. For limb lengths of 4, 9, 14, 22, */
/* 27, 35, 40, ..., the maximum z value, minimum y value ratio is less than */
/* 3/2. For limb lengths of 7, 12, 17, 20, 25, 30, 33, 38, ..., the minimum */
/* z value, maximum y value ratio is greater than 3/4. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void add64(unsigned int *a, unsigned int *b);
void shift(unsigned int *a, unsigned int *b, unsigned int c);
int main () {
//
unsigned int order=50331648; // order
unsigned int tflag=0; // A, B, C, or D (0, 1, 2, 3) (k=26)
// // A or B (0 or 1) (k=25)
int c=1; // c
unsigned int l=22; // specified limb length
unsigned int a=9; // number of elements in sv
unsigned int b=9; // number of sequence vectors
//
unsigned int d,e,f,g,h,i,j,n,t,u,mask,count,flag,offset,first; // indices
unsigned int savea,saveb;
int k,m;
unsigned int X[2],Y[2];
unsigned int v[66]; // sequence vector
unsigned int maxy[66]; // maximum y value array
unsigned int miny[66]; // minimum y value array
unsigned int maxz[66]; // maximum z value array
unsigned int minz[66]; // minimum z value array
unsigned int he[66]; // sequence vector histogram
unsigned int s[131072]; // duplicate array, minimum size=(order/12)/32
unsigned int ho[512]; // histogram of lengths
unsigned int dist[200]; // histogram of z values
int num[19]={1,-1,7,5,-17,47,13,-217,295,-139,1909,1631,-3299,
13085,6487,-46075,84997,-7153,502829};
int den[19]={4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,
32768,65536,131072,262144,524288,2097152};
unsigned int len[19]={2,4,5,7,9,10,12,14,15,17,18,20,22,23,25,27,28,30,33};
unsigned int aa[19]= {0,0,0,1,1, 0, 2, 2, 0, 3, 0, 4, 4, 0, 5, 5, 0, 6, 7};
unsigned int bb[19]= {1,2,0,2,3, 0, 3, 4, 0, 4, 0, 4, 5, 0, 5, 6, 0, 6, 6};
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
//int z[1]={2}; // length=2
//int z[1]={10}; // length=4
//int z[1]={58}; // length=7
//int z[1]={206}; // length=9
//unsigned int z[2]={ // length=12, live only
// 986,842};
//unsigned int z[2]={ // length=14, live only
// 2782,3214};
//unsigned int z[5]={ // length=17, live only
// 14314,11434,12586,10138,11290};
//unsigned int z[5]={ // length=20, live only
// 54718,49534,50110,44926,61630};
unsigned int z[9]={ // length=22, live only
156794,158522,172346,193082,145130,142970,132602,131306,120938};
//unsigned int z[19]={ // length=25, live only
// 666542,532910,569774,574382,728750,611246,636590,673454,811694,619886,486254,
// 523118,527726,564590,451262,492734,584894,488126,529598};
//unsigned int z[23]={ // length=27, live only
// 1788682,1654330,2065162,1975306,2251786,1774858,1524298,1899274,1648714,
// 2500618,1925194,1541578,1326010,1634890,1419322,1759306,1529914,1543738,
// 1820218,1664266,1430986,2085898,1436602};
//unsigned int z[66]={ // length=30, live only
// 10356766,9361438,8697886,8614942,8255518,8055070,7960606,7951390,7635166,
// 7509022,7453726,7391518,7320238,7214110,7011358,6971614,6949150,6893854,
// 6716446,6679582,6656686,6654238,6529246,6520606,6473950,6451486,6384670,
// 6234334,6214318,6159022,6156574,6119710,6100702,6078238,6031582,5919406,
// 5824798,5820766,5785774,5783326,5746462,5736670,5716654,5699806,5658334,
// 5505838,5451550,5421742,5404894,5384878,5378398,5363422,5343406,5326558,
// 5089966,5083486,5063470,5048494,5046622,5031646,5011630,4768558,4751710,
// 4731694,4716718,4436782};
//unsigned int z[66]={ // length=33, live only
// 38508634,35854426,35522650,34084954,33283162,32905306,32868442,31603546,
// 31098970,30877786,30343834,29919322,28949338,28859482,28638298,27928666,
// 27781210,27689626,27679834,27179866,27145306,26958682,26868826,26601562,
// 26000218,25920154,25698970,25689178,25541722,25465690,25189210,24740506,
// 24362074,24345946,24196186,24048730,24009562,23862106,23696218,23086234,
// 22869082,22749850,22682458,22602394,22576474,22436506,22369114,21422746,
// 21396826,21316762,21256858,21249370,21189466,21109402,20137114,20069722,
// 19989658,19929754,18810010,42489946,30628954,29108314,25375834,24205978,
// 23929498,22516570};
//unsigned int sv[1*1]={ // length=2, live only
// 1};
//unsigned int sv[2*1]={ // length=4, live only
// 1, 1};
//unsigned int sv[3*1]={ // length=7, live only
// 1, 2, 1};
//unsigned int sv[4*1]={ // length=9, live only
// 1, 2, 1, 1};
//unsigned int sv[5*2]={ // length=12, live only
// 1, 2, 2, 1, 1, // 986
// 1, 2, 1, 2, 1}; // 842
//unsigned int sv[6*2]={ // length=14, live only
// 1, 2, 2, 1, 1, 1, // 3214
// 1, 2, 1, 2, 1, 1}; // 2782
//unsigned int sv[7*5]={ // length=17, live only
// 1, 2, 2, 2, 1, 1, 1, // 14314
// 1, 2, 2, 1, 2, 1, 1, // 12586
// 1, 2, 2, 1, 1, 2, 1, // 11434
// 1, 2, 1, 2, 2, 1, 1, // 11290
// 1, 2, 1, 2, 1, 2, 1}; // 10138
//unsigned int sv[8*5]={ // length=20, live only
// 1, 2, 2, 2, 2, 1, 1, 1, // 61630
// 1, 2, 2, 2, 1, 2, 1, 1, // 54718
// 1, 2, 2, 2, 1, 1, 2, 1, // 50110
// 1, 2, 2, 1, 2, 2, 1, 1, // 49534
// 1, 2, 2, 1, 2, 1, 2, 1}; // 44926
unsigned int sv[9*9]={ // length=22, live only
1, 2, 2, 2, 2, 1, 1, 1, 1, // 193082
1, 2, 2, 2, 1, 2, 1, 1, 1, // 172346
1, 2, 2, 2, 1, 1, 2, 1, 1, // 158522
1, 2, 2, 1, 2, 2, 1, 1, 1, // 156794
1, 2, 1, 2, 2, 2, 1, 1, 1, // 145130
1, 2, 2, 1, 2, 1, 2, 1, 1, // 142970
1, 2, 2, 1, 1, 2, 2, 1, 1, // 132602
1, 2, 1, 2, 2, 1, 2, 1, 1, // 131306
1, 2, 1, 2, 1, 2, 2, 1, 1}; // 120938
/* unsigned int sv[10*19]={ // length=25, live only
1, 2, 2, 2, 2, 2, 1, 1, 1, 1, // 811694
1, 2, 2, 2, 2, 1, 2, 1, 1, 1, // 728750
1, 2, 2, 2, 2, 1, 1, 2, 1, 1, // 673454
1, 2, 2, 2, 1, 2, 2, 1, 1, 1, // 666542
1, 2, 2, 2, 2, 1, 1, 1, 2, 1, // 636590
1, 2, 2, 1, 2, 2, 2, 1, 1, 1, // 619886
1, 2, 2, 2, 1, 2, 1, 2, 1, 1, // 611246
1, 2, 1, 2, 2, 2, 2, 1, 1, 1, // 584894
1, 2, 2, 2, 1, 2, 1, 1, 2, 1, // 574382
1, 2, 2, 2, 1, 1, 2, 2, 1, 1, // 569774
1, 2, 2, 1, 2, 2, 1, 2, 1, 1, // 564590
1, 2, 2, 2, 1, 1, 2, 1, 2, 1, // 532910
1, 2, 1, 2, 2, 2, 1, 2, 1, 1, // 529598
1, 2, 2, 1, 2, 2, 1, 1, 2, 1, // 527726
1, 2, 2, 1, 2, 1, 2, 2, 1, 1, // 523118
1, 2, 1, 2, 2, 2, 1, 1, 2, 1, // 492734
1, 2, 1, 2, 2, 1, 2, 2, 1, 1, // 488126
1, 2, 2, 1, 2, 1, 2, 1, 2, 1, // 486254
1, 2, 1, 2, 2, 1, 2, 1, 2, 1}; // 451262 */
/* unsigned int sv[11*23]={ // length=27, live only
1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, // 2500618
1, 2, 2, 2, 2, 1, 2, 1, 1, 1, 1, // 2251786
1, 2, 2, 2, 2, 1, 1, 2, 1, 1, 1, // 2085898
1, 2, 2, 2, 1, 2, 2, 1, 1, 1, 1, // 2065162
1, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, // 1975306
1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, // 1925194
1, 2, 2, 2, 1, 2, 1, 2, 1, 1, 1, // 1899274
1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, // 1820218
1, 2, 2, 2, 1, 2, 1, 1, 2, 1, 1, // 1788682
1, 2, 2, 2, 1, 1, 2, 2, 1, 1, 1, // 1774858
1, 2, 2, 1, 2, 2, 1, 2, 1, 1, 1, // 1759306
1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 1, // 1664266
1, 2, 1, 2, 2, 2, 1, 2, 1, 1, 1, // 1654330
1, 2, 2, 1, 2, 2, 1, 1, 2, 1, 1, // 1648714
1, 2, 2, 1, 2, 1, 2, 2, 1, 1, 1, // 1634890
1, 2, 1, 2, 2, 2, 1, 1, 2, 1, 1, // 1543738
1, 2, 2, 1, 1, 2, 2, 2, 1, 1, 1, // 1541578
1, 2, 1, 2, 2, 1, 2, 2, 1, 1, 1, // 1529914
1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, // 1524298
1, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, // 1436602
1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 1, // 1430986
1, 2, 1, 2, 2, 1, 2, 1, 2, 1, 1, // 1419322
1, 2, 1, 2, 1, 2, 2, 1, 2, 1, 1}; // 1326010 */
/* unsigned int sv[12*66]={ // length=30, live only
1, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, // 10356766
1, 2, 2, 2, 2, 2, 1, 2, 1, 1, 1, 1, // 9361438
1, 2, 2, 2, 2, 2, 1, 1, 2, 1, 1, 1, // 8697886
1, 2, 2, 2, 2, 1, 2, 2, 1, 1, 1, 1, // 8614942
1, 2, 2, 2, 2, 2, 1, 1, 1, 2, 1, 1, // 8255518
1, 2, 2, 2, 1, 2, 2, 2, 1, 1, 1, 1, // 8055070
1, 2, 2, 2, 2, 2, 1, 1, 1, 1, 2, 1, // 7960606
1, 2, 2, 2, 2, 1, 2, 1, 2, 1, 1, 1, // 7951390
1, 2, 2, 1, 2, 2, 2, 2, 1, 1, 1, 1, // 7635166
1, 2, 2, 2, 2, 1, 2, 1, 1, 2, 1, 1, // 7509022
1, 2, 2, 2, 2, 1, 1, 2, 2, 1, 1, 1, // 7453726
1, 2, 2, 2, 1, 2, 2, 1, 2, 1, 1, 1, // 7391518
1, 2, 1, 2, 2, 2, 2, 2, 1, 1, 1, 1, // 7320238
1, 2, 2, 2, 2, 1, 2, 1, 1, 1, 2, 1, // 7214110
1, 2, 2, 2, 2, 1, 1, 2, 1, 2, 1, 1, // 7011358
1, 2, 2, 1, 2, 2, 2, 1, 2, 1, 1, 1, // 6971614
1, 2, 2, 2, 1, 2, 2, 1, 1, 2, 1, 1, // 6949150
1, 2, 2, 2, 1, 2, 1, 2, 2, 1, 1, 1, // 6893854
1, 2, 2, 2, 2, 1, 1, 2, 1, 1, 2, 1, // 6716446
1, 2, 2, 2, 2, 1, 1, 1, 2, 2, 1, 1, // 6679582
1, 2, 1, 2, 2, 2, 2, 1, 2, 1, 1, 1, // 6656686
1, 2, 2, 2, 1, 2, 2, 1, 1, 1, 2, 1, // 6654238
1, 2, 2, 1, 2, 2, 2, 1, 1, 2, 1, 1, // 6529246
1, 2, 2, 2, 1, 1, 2, 2, 2, 1, 1, 1, // 6520606
1, 2, 2, 1, 2, 2, 1, 2, 2, 1, 1, 1, // 6473950
1, 2, 2, 2, 1, 2, 1, 2, 1, 2, 1, 1, // 6451486
1, 2, 2, 2, 2, 1, 1, 1, 2, 1, 2, 1, // 6384670
1, 2, 2, 1, 2, 2, 2, 1, 1, 1, 2, 1, // 6234334
1, 2, 1, 2, 2, 2, 2, 1, 1, 2, 1, 1, // 6214318
1, 2, 1, 2, 2, 2, 1, 2, 2, 1, 1, 1, // 6159022
1, 2, 2, 2, 1, 2, 1, 2, 1, 1, 2, 1, // 6156574
1, 2, 2, 2, 1, 2, 1, 1, 2, 2, 1, 1, // 6119710
1, 2, 2, 1, 2, 1, 2, 2, 2, 1, 1, 1, // 6100702
1, 2, 2, 2, 1, 1, 2, 2, 1, 2, 1, 1, // 6078238
1, 2, 2, 1, 2, 2, 1, 2, 1, 2, 1, 1, // 6031582
1, 2, 1, 2, 2, 2, 2, 1, 1, 1, 2, 1, // 5919406
1, 2, 2, 2, 1, 2, 1, 1, 2, 1, 2, 1, // 5824798
1, 2, 2, 1, 1, 2, 2, 2, 2, 1, 1, 1, // 5820766
1, 2, 1, 2, 2, 1, 2, 2, 2, 1, 1, 1, // 5785774
1, 2, 2, 2, 1, 1, 2, 2, 1, 1, 2, 1, // 5783326
1, 2, 2, 2, 1, 1, 2, 1, 2, 2, 1, 1, // 5746462
1, 2, 2, 1, 2, 2, 1, 2, 1, 1, 2, 1, // 5736670
1, 2, 1, 2, 2, 2, 1, 2, 1, 2, 1, 1, // 5716654
1, 2, 2, 1, 2, 2, 1, 1, 2, 2, 1, 1, // 5699806
1, 2, 2, 1, 2, 1, 2, 2, 1, 2, 1, 1, // 5658334
1, 2, 1, 2, 1, 2, 2, 2, 2, 1, 1, 1, // 5505838
1, 2, 2, 2, 1, 1, 2, 1, 2, 1, 2, 1, // 5451550
1, 2, 1, 2, 2, 2, 1, 2, 1, 1, 2, 1, // 5421742
1, 2, 2, 1, 2, 2, 1, 1, 2, 1, 2, 1, // 5404894
1, 2, 1, 2, 2, 2, 1, 1, 2, 2, 1, 1, // 5384878
1, 2, 2, 1, 1, 2, 2, 2, 1, 2, 1, 1, // 5378398
1, 2, 2, 1, 2, 1, 2, 2, 1, 1, 2, 1, // 5363422
1, 2, 1, 2, 2, 1, 2, 2, 1, 2, 1, 1, // 5343406
1, 2, 2, 1, 2, 1, 2, 1, 2, 2, 1, 1, // 5326558
1, 2, 1, 2, 2, 2, 1, 1, 2, 1, 2, 1, // 5089966
1, 2, 2, 1, 1, 2, 2, 2, 1, 1, 2, 1, // 5083486
1, 2, 1, 2, 1, 2, 2, 2, 1, 2, 1, 1, // 5063470
1, 2, 1, 2, 2, 1, 2, 2, 1, 1, 2, 1, // 5048494
1, 2, 2, 1, 1, 2, 2, 1, 2, 2, 1, 1, // 5046622
1, 2, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, // 5031646
1, 2, 1, 2, 2, 1, 2, 1, 2, 2, 1, 1, // 5011630
1, 2, 1, 2, 1, 2, 2, 2, 1, 1, 2, 1, // 4768558
1, 2, 2, 1, 1, 2, 2, 1, 2, 1, 2, 1, // 4751710
1, 2, 1, 2, 1, 2, 2, 1, 2, 2, 1, 1, // 4731694
1, 2, 1, 2, 2, 1, 2, 1, 2, 1, 2, 1, // 4716718
1, 2, 1, 2, 1, 2, 2, 1, 2, 1, 2, 1}; // 4436782 */
/* unsigned int sv[13*66]={ // length=33, live only
1,2,2,2,2,2,2,2,1,1,1,1,1, // 42489946
1,2,2,2,2,2,2,1,2,1,1,1,1, // 38508634
1,2,2,2,2,2,2,1,1,2,1,1,1, // 35854426
1,2,2,2,2,2,1,2,2,1,1,1,1, // 35522650
1,2,2,2,2,2,2,1,1,1,2,1,1, // 34084954
1,2,2,2,2,1,2,2,2,1,1,1,1, // 33283162
1,2,2,2,2,2,2,1,1,1,1,2,1, // 32905306
1,2,2,2,2,2,1,2,1,2,1,1,1, // 32868442
1,2,2,2,1,2,2,2,2,1,1,1,1, // 31603546
1,2,2,2,2,2,1,2,1,1,2,1,1, // 31098970
1,2,2,2,2,2,1,1,2,2,1,1,1, // 30877786
1,2,2,2,2,1,2,2,1,2,1,1,1, // 30628954
1,2,2,1,2,2,2,2,2,1,1,1,1, // 30343834
1,2,2,2,2,2,1,2,1,1,1,2,1, // 29919322
1,2,2,2,2,2,1,1,2,1,2,1,1, // 29108314
1,2,2,2,1,2,2,2,1,2,1,1,1, // 28949338
1,2,2,2,2,1,2,2,1,1,2,1,1, // 28859482
1,2,2,2,2,1,2,1,2,2,1,1,1, // 28638298
1,2,2,2,2,2,1,1,2,1,1,2,1, // 27928666
1,2,2,2,2,2,1,1,1,2,2,1,1, // 27781210
1,2,2,1,2,2,2,2,1,2,1,1,1, // 27689626
1,2,2,2,2,1,2,2,1,1,1,2,1, // 27679834
1,2,2,2,1,2,2,2,1,1,2,1,1, // 27179866
1,2,2,2,2,1,1,2,2,2,1,1,1, // 27145306
1,2,2,2,1,2,2,1,2,2,1,1,1, // 26958682
1,2,2,2,2,1,2,1,2,1,2,1,1, // 26868826
1,2,2,2,2,2,1,1,1,2,1,2,1, // 26601562
1,2,2,2,1,2,2,2,1,1,1,2,1, // 26000218
1,2,2,1,2,2,2,2,1,1,2,1,1, // 25920154
1,2,2,1,2,2,2,1,2,2,1,1,1, // 25698970
1,2,2,2,2,1,2,1,2,1,1,2,1, // 25689178
1,2,2,2,2,1,2,1,1,2,2,1,1, // 25541722
1,2,2,2,1,2,1,2,2,2,1,1,1, // 25465690
1,2,2,2,2,1,1,2,2,1,2,1,1, // 25375834
1,2,2,2,1,2,2,1,2,1,2,1,1, // 25189210
1,2,2,1,2,2,2,2,1,1,1,2,1, // 24740506
1,2,2,2,2,1,2,1,1,2,1,2,1, // 24362074
1,2,2,2,1,1,2,2,2,2,1,1,1, // 24345946
1,2,2,1,2,2,1,2,2,2,1,1,1, // 24205978
1,2,2,2,2,1,1,2,2,1,1,2,1, // 24196186
1,2,2,2,2,1,1,2,1,2,2,1,1, // 24048730
1,2,2,2,1,2,2,1,2,1,1,2,1, // 24009562
1,2,2,1,2,2,2,1,2,1,2,1,1, // 23929498
1,2,2,2,1,2,2,1,1,2,2,1,1, // 23862106
1,2,2,2,1,2,1,2,2,1,2,1,1, // 23696218
1,2,2,1,2,1,2,2,2,2,1,1,1, // 23086234
1,2,2,2,2,1,1,2,1,2,1,2,1, // 22869082
1,2,2,1,2,2,2,1,2,1,1,2,1, // 22749850
1,2,2,2,1,2,2,1,1,2,1,2,1, // 22682458
1,2,2,1,2,2,2,1,1,2,2,1,1, // 22602394
1,2,2,2,1,1,2,2,2,1,2,1,1, // 22576474
1,2,2,2,1,2,1,2,2,1,1,2,1, // 22516570
1,2,2,1,2,2,1,2,2,1,2,1,1, // 22436506
1,2,2,2,1,2,1,2,1,2,2,1,1, // 22369114
1,2,2,1,2,1,2,2,2,1,2,1,1, // 21316762
1,2,2,1,2,2,1,2,2,1,1,2,1, // 21256858
1,2,2,2,1,1,2,2,1,2,2,1,1, // 21249370
1,2,2,2,1,2,1,2,1,2,1,2,1, // 21189466
1,2,2,1,2,2,2,1,1,2,1,2,1, // 21422746
1,2,2,2,1,1,2,2,2,1,1,2,1, // 21396826
1,2,2,1,2,2,1,2,1,2,2,1,1, // 21109402
1,2,2,1,2,1,2,2,2,1,1,2,1, // 20137114
1,2,2,2,1,1,2,2,1,2,1,2,1, // 20069722
1,2,2,1,2,1,2,2,1,2,2,1,1, // 19989658
1,2,2,1,2,2,1,2,1,2,1,2,1, // 19929754
1,2,2,1,2,1,2,2,1,2,1,2,1}; // 18810010 */
int r,del,savex;
double x,dx,maxr,minr,rat,ratio;
double maxmin[2*66];
unsigned int savem;
FILE *Outfp;
Outfp = fopen("out6a.dat","w");
if ((l==5)||(l==10)||(l==15)||(l==18)||(l==23)||(l==28)) {
printf("dead limbs \n");
goto zskip;
}
if (l>33) {
printf("length too big \n");
goto zskip;
}
if (order>50331648*4) {
printf("order too big \n");
goto zskip;
}
if (c<-1) {
printf("c too small \n");
goto zskip;
}
if (c>1) {
printf("c too big \n");
goto zskip;
}
if (c==1)
offset=8;
else
offset=4;
for (i=0; i<200; i++) // clear distribution of e values
dist[i]=0;
for (i=0; i<512; i++) // clear histogram of lengths
ho[i]=0;
for (i=0; i<131072; i++) // clear duplicate array
s[i]=0;
for (i=0; i<66; i++) { // clear maximum y values
maxy[i]=0;
miny[i]=0x7fffffff;
maxz[i]=0;
minz[i]=0x7fffffff;
he[i]=0; // clear histogram of sv's
}
for (i=0; i<66; i++) {
maxmin[2*i]=0.0;
maxmin[2*i+1]=1000000.0;
}
fprintf(Outfp,"ORDER=%d c=%d \n",order,c);
//
// limbs in A, B, C, and D
//
g=0; // clear counter
count=0; // clear counter
maxr=0.0;
minr=1000000.0;
for (i=order/2; i<order; i+=2) { // possible starting elements
if ((i-offset)==((i-offset)/12)*12) // check for elements of U
continue;
k=i; // back-track
if (k!=(k/3)*3) { // check for dead limb
if ((k-c)!=((k-c)/3)*3) { // check for beginning of limb
goto askip;
}
k=(k-c)/3; // previous number in sequence
if (k==(k/3)*3) // check for dead limb
goto bskip;
k=k*2; // previous number in sequence
while ((k<(int)(order/2))||((k-c)==((k-c)/3)*3)) { // check for beginning
if ((k-c)==((k-c)/3)*3) {
k=(k-c)/3; // previous number in sequence
if (k==(k/3)*3) // check for dead limb
goto bskip;
k=k*2; // previous number in sequence
}
else
k=k*2; // previous number in sequence
}
}
askip:
m=k; // save beginning of limb
n=1; // set length
while (k==(k/2)*2) { // check for even element
k=k/2; // next element of limb
n=n+1; // increment length
}
for (j=1; j<1000000; j++) { // iterate until end of limb
h=3*k+c; // next element of limb
if (h>order) { // check for end of limb
if (k<(int)(order/3))
goto bskip;
if (order<=50331648)
t=((k-(order/3))-1)/2; // index into array
if (order==50331648*2) {
if (tflag==0) {
if ((k-1)!=((k-1)/4)*4)
goto bskip;
}
if (tflag==1) {
if ((k-3)!=((k-3)/4)*4)
goto bskip;
}
if (tflag>1)
printf("error: flag=%d \n",tflag);
t=((k-(order/3))-1)/4; // index into array
}
if (order==50331648*4) {
if (tflag==0) {
if ((k-1)!=((k-1)/8)*8)
goto bskip;
}
if (tflag==1) {
if ((k-3)!=((k-3)/8)*8)
goto bskip;
}
if (tflag==2) {
if ((k-5)!=((k-5)/8)*8)
goto bskip;
}
if (tflag==3) {
if ((k-7)!=((k-7)/8)*8)
goto bskip;
}
if (tflag>3)
printf("error: flag=%d \n",tflag);
t=((k-(order/3))-1)/8; // index into array
}
u=t-(t/32)*32; // fractional portion of index
t=t/32; // integer portion of index
mask=1<<u; // set mask
if ((s[t]&mask)==0) // check if bit not set
s[t]=s[t]|mask; // set bit in array
else // already set
goto bskip;
g=g+1; // increment counter
if (n<512) // check length
ho[n]=ho[n]+1; // histogram length
else {
printf("error: histogram array not big enough \n");
goto zskip;
}
if (n!=l) // continue if not specified length
goto bskip;
if (m==(m/3)*3)
goto bskip;
r=(int)(m)-(int)(h/2); // difference of first elements
savem=m;
for (f=0; f<19; f++) {
if (n==len[f]) {
savea=aa[f];
saveb=bb[f];
savex=den[f];
x=(double)m;
x=x*(double)num[f];
dx=(double)r;
dx=dx*(double)den[f];
x=x-dx;
del=(int)x;
flag=0;
for (e=0; e<b; e++) {
if (del==c*(int)z[e]) {
flag=z[e];
goto fskip;
}
}
printf("error: no solution of Diophantine equation found \n");
goto zskip;
fskip: if (l==2) {
v[0]=1;
f=1;
goto dskip;
}
f=0;
if ((m&1)==0) {
count=0;
while ((m&1)==0) {
count=count+1;
m=m/2;
}
v[0]=count;
f=1;
}
first=1;
while (m!=k) {
if ((m&1)==0) {
m=m/2;
count=count+1;
}
else {
m=3*m+c;
if (first==0) {
v[f]=count;
f=f+1;
if (f>16) {
printf("error: sequence vector array not big enough \n");
goto zskip;
}
}
else
first=0;
count=0;
}
}
v[f]=count;
f=f+1;
dskip: if (f!=a) {
printf("error: incorrect sequence vector length \n");
goto zskip;
}
x=0.5;
for (d=0; d<b; d++) {
for (j=0; j<a; j++) {
if (v[j]!=sv[j+a*d])
goto eskip;
if (sv[j+a*d]==1)
x=x*1.5;
else
x=x*0.75;
}
//
// delta=z-(y/2)((3/4)**a)*((3/2)**b), then multiply by X/y
// (same as X((z/y)-(1/2)((3/4)**a)*((3/2)**b))
//
rat=(double)savex*((double)(h/2)/(double)savem-x);
if (dist[d]==0) {
printf("limb=(%d, %d)=>%d, z'=%e, del=%e \n",
savem,k,h/2,x*(double)savem,(double)(h/2)-x*(double)savem);
printf("del/y=%e, r=%e, index=%d \n",
(double)(h/2)/(double)savem-x,rat,d);
fprintf(Outfp,"del/y=%e, r=%e, index=%d, z=%d \n",
(double)(h/2)/(double)savem-x,rat,d,flag);
printf("\n");
}
dist[d]=dist[d]+1;
if (rat<0.0)
rat=-rat;
if (rat>maxr)
maxr=rat;
if (rat<minr)
minr=rat;
if (rat>maxmin[2*d])
maxmin[2*d]=rat;
if (rat<maxmin[2*d+1])
maxmin[2*d+1]=rat;
if (savem>maxy[d])
maxy[d]=savem;
if (savem<miny[d])
miny[d]=savem;
if ((h/2)>maxz[d])
maxz[d]=h/2;
if ((h/2)<minz[d])
minz[d]=h/2;
he[d]=he[d]+1;
//
// check if o=[(y/2)(3/4)**a]+0,1
//
if (d==0) {
if ((l==2)||(l==4))
goto bskip;
X[0]=0;
X[1]=savem/2;
for (j=0; j<savea; j++) {
Y[0]=X[0];
Y[1]=X[1];
add64(Y, X);
add64(Y, X);
}
shift(X, X, 2*savea);
m=(int)X[1];
if ((c==1)&&(m!=(m/2)*2)) {
printf("error: o is odd \n");
goto zskip;
}
if ((c==-1)&&(m==(m/2)*2)) {
printf("error: o is even \n");
goto zskip;
}
if (c==1)
m=m+1;
m=m+c;
for (j=0; j<saveb; j++) {
m=m/2;
m=m*3;
}
m=m-c;
if (m!=(int)(h/2)) {
printf("error: y=%d, o=%d, z=%d \n",savem,m,h/2);
goto zskip;
}
}
goto bskip;
eskip: x=0.5;
}
printf("error: sequence vector not found \n");
fprintf(Outfp,"error: sequence vector not found \n");
goto zskip;
}
}
goto bskip;
}
k=h; // next element
if (k==(k/8)*8) // not valid limb, start over
goto bskip; // (prevents taking even path at nodes)
n=n+1; // increment length
while (k==(k/2)*2) { // check for even element
k=k/2; // next element
n=n+1; // increment length
}
}
bskip:
n=0;
}
printf("maxr=%e, minr=%e \n",maxr,minr);
//
// compute upper bound
//
if (l>4) {
printf("\n");
x=2.0;
j=4;
for (i=0; i<a; i++) {
if (sv[i]==1) {
x=x*1.5;
j=j<<1;
}
else
j=j<<2;
}
x=x-1.0;
x=x*(double)j;
x=x/(double)order;
printf("upper bound of X(delta/y)=%e \n",x);
}
//
// DISTRIBUTION OF Z VALUES
//
printf("\n");
fprintf(Outfp,"\n");
printf("distribution of e values \n");
for (i=0; i<b; i++) {
printf("i=%d, d=%d, max=%e, min=%e \n",i,dist[i],maxmin[2*i],maxmin[2*i+1]);
fprintf(Outfp,"i=%d, d=%d, max=%e, min=%e \n",i,dist[i],maxmin[2*i],maxmin[2*i+1]);
}
printf("\n");
printf("maximum y values \n");
for (i=0; i<b; i++)
printf(" i=%d, max=%d, n=%d, ratio=%e \n",i,maxy[i],he[i],(double)maxy[i]/(double)order);
printf("\n");
printf("minimum y values \n");
for (i=0; i<b; i++)
printf(" i=%d, min=%d, n=%d, ratio=%e \n",i,miny[i],he[i],(double)miny[i]/(double)order);
printf("\n");
printf("maximum z values \n");
for (i=0; i<b; i++)
printf(" i=%d, max=%d, n=%d, ratio=%e \n",i,maxz[i],he[i],(double)maxz[i]/(double)order);
printf("\n");
printf("minimum z values \n");
for (i=0; i<b; i++)
printf(" i=%d, min=%d, n=%d, ratio=%e \n",i,minz[i],he[i],(double)minz[i]/(double)order);
if ((l==7)||(l==12)||(l==17)||(l==20)||(l==25)||(l==30)||(l==33))
flag=1;
else
flag=0;
printf("\n");
printf("maximum z/minimum y values \n");
for (i=0; i<b; i++) {
ratio=(double)maxz[i]/(double)miny[i];
printf("maximum z/minimum y=%e \n",ratio);
if ((flag==0)&&(ratio>1.5)) {
printf("error: ratio>1.5 \n");
goto zskip;
}
}
printf("\n");
printf("minimum z/maximum y values \n");
for (i=0; i<b; i++) {
ratio=(double)minz[i]/(double)maxy[i];
printf("minimum z/maximum y=%e \n",ratio);
if ((flag==1)&&(ratio<.75)) {
printf("error: ratio<.75 \n");
goto zskip;
}
}
//
// check for valid lengths
//
for (i=0; i<40; i++) {
if (ho[y[i]]!=0) {
printf("error: n=%d \n",y[i]);
goto zskip;
}
}
//
// HISTOGRAM OF LENGTHS
//
if (order<=50331648) {
if (g!=(order/12)) {
printf("error: incorrect number of entries in histogram \n");
printf("g=%d, order/12=%d \n",g,order/12);
goto zskip;
}
}
if (order==(50331648*2)) {
if (g!=(order/24)) {
printf("error: incorrect number of entries in histogram \n");
printf("g=%d, order/12=%d \n",g,order/12);
goto zskip;
}
}
if (order==(50331648*4)) {
if (g!=(order/48)) {
printf("error: incorrect number of entries in histogram \n");
printf("g=%d, order/12=%d \n",g,order/12);
goto zskip;
}
}
fprintf(Outfp,"\n");
printf("\n");
fprintf(Outfp,"HISTOGRAM: number of entries=%d \n",g);
printf("HISTOGRAM: number of entries=%d \n",g);
for (h=0; h<64; h++) {
fprintf(Outfp," %d %d %d %d %d %d %d %d \n",ho[8*h],ho[8*h+1],ho[8*h+2],
ho[8*h+3],ho[8*h+4],ho[8*h+5],ho[8*h+6],ho[8*h+7]);
}
for (h=0; h<16; h++) {
printf(" %d %d %d %d %d %d %d %d \n",ho[8*h],ho[8*h+1],ho[8*h+2],
ho[8*h+3],ho[8*h+4],ho[8*h+5],ho[8*h+6],ho[8*h+7]);
}
zskip:
fclose(Outfp);
return(0);
}