/******************************************************************************/
/* */
/* COMPUTE 3N+1 OR 3N-1 SEQUENCE */
/* 01/27/10 (dkc) */
/* */
/* This C program generates limbs where the first element of the limb is */
/* even, less than the order, and greater than the order over 2. The last */
/* element of the limb is odd and greater than the order/3. Only live limbs */
/* are generated. The limbs must have at least one element that is divisible*/
/* by 8. */
/* */
/* Note: The order must be less than 3*2**31. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void add64(unsigned int *a, unsigned int *b);
void sub64(unsigned int *c, unsigned int *d);
void sub128(unsigned int *c, unsigned int *d);
void add128(unsigned int *c, unsigned int *d);
void shift(unsigned int *a, unsigned int *b, unsigned int shift);
void mul64_32(unsigned int a0, unsigned int a2, unsigned int m0,
unsigned int *product);
unsigned int subr2(unsigned int A, unsigned int B, unsigned int *P);
int main () {
//
int c=-1; // c
unsigned int lshift=20; // order shift
unsigned int l=17; // limb length
unsigned int p=10; // number of elements in e array
unsigned int a=7; // number of elements in sv
unsigned int b=10; // number of sequence vectors
//
unsigned int d,e,f,g,h,i,j,n,iters,count,first,savee; // indices
unsigned int flag,order;
unsigned int v[2*4096]; // sub-sequence arrays
unsigned int t[200];
unsigned int ho[512]; // histogram of limb lengths
unsigned int he[512];
unsigned int maxy[200*2]; // maximum y value array
int del;
unsigned int T[2],U[2],V[2],W[2],X[2],Y[2],Z[2];
unsigned int A[3],B[4],C[4];
unsigned int dist[200]; // distribution of e values
int num[15]={1,-1,7,5,-17,47,13,-217,295,-139,1909,1631,-3299,
13085,6487};
int den[15]={4,8,16,32,64,128,256,512,1024,2048,4096,8192,16384,
32768,65536};
unsigned int len[15]={2,4,5,7,9,10,12,14,15,17,18,20,22,23,25};
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]={1202}; // length=12, live only
// int z[1]={3862}; // length=14, live only
// int z[2]={ // length=15, live only
// 7724,5780};
int z[10]={ // length=17, live only
24196,22738,18850,18364,16906,16258,14530,14476,13378,13018};
// int z[10]={ // length=18, live only
// 48392,45476,37700,36728,33812,32516,29060,28952,26756,26036};
// int z[45]={ // length=20, live only
// 149272,140524,136150,117196,114280,112822,105532,101644,101158,99700,
// 97270,95326,91276,90952,86902,85606,84364,84148,82204,79990,79756,79774,
// 77830,76372,75400,75382,75238,73780,71998,69406,68326,67462,66868,66652,
// 63718,62494,62278,62260,60820,60550,57886,56446,55942,55366,50758};
// int z[63]={ // length=22, live only
// 456008,429764,416642,359780,351032,346658,324788,313124,311666,307292,
// 300002,294170,282020,281048,268898,265010,261284,260636,254804,248162,
// 247514,247460,241682,234392,234338,237308,233906,229532,225644,224186,
// 216410,213170,212522,210578,208796,208148,203288,199346,195674,195026,
// 194972,194540,190652,189842,181850,181418,178988,177530,177044,176018,
// 174290,173804,165866,163922,160682,160466,159980,159548,150098,147884,
// 146858,146426,134762};
// int z[98]={ // length=23
// 912016,859528,833284,719560,702064,693316,649576,626248,623332,614584,
// 600004,588340,570844,564040,562096,537796,530020,522568,521272,509608,
// 496324,495028,494920,483364,477532,476488,474616,468784,468676,467812,
// 459064,451288,450244,448372,432820,430876,426340,425044,421156,417592,
// 416296,415324,407548,406576,398692,391348,390052,389944,389080,381304,
// 380260,379684,373852,371512,368668,363700,362836,357976,355060,354088,
// 352036,348580,347608,346204,345340,345268,337564,333604,331732,327844,
// 327772,327196,321364,320932,319960,319096,314236,303868,302500,301528,
// 300196,299548,296092,295768,293716,292852,281764,281116,276220,275356,
// 275284,269524,268444,257788,252028,250012,247708,229276};
//unsigned int sv[5*1]={ // length=12, live only
// 1, 3, 1, 1, 1}; // 1202
//unsigned int sv[6*1]={ // length=14, live only
// 1, 3, 1, 1, 1, 1}; // 3862
//unsigned int sv[6*2]={ // length=15, live only
// 2, 3, 1, 1, 1, 1, // 7724
// 2, 1, 3, 1, 1, 1};// 5780
unsigned int sv[7*10]={ // length=17, live only
2, 3, 1, 1, 1, 1, 1, // 24196
1, 4, 1, 1, 1, 1, 1, // 22738
1, 3, 2, 1, 1, 1, 1, // 18850
2, 1, 3, 1, 1, 1, 1, // 18364
1, 2, 3, 1, 1, 1, 1, // 16906
1, 3, 1, 2, 1, 1, 1, // 16258
1, 3, 1, 1, 2, 1, 1, // 14530
2, 1, 1, 3, 1, 1, 1, // 14476
1, 3, 1, 1, 1, 2, 1, // 13378
1, 2, 1, 3, 1, 1, 1};// 13018
/* unsigned int sv[7*10]={ // length=18, live only
3, 3, 1, 1, 1, 1, 1, // 48392
2, 4, 1, 1, 1, 1, 1, // 45476
2, 3, 2, 1, 1, 1, 1, // 37700
3, 1, 3, 1, 1, 1, 1, // 36728
2, 2, 3, 1, 1, 1, 1, // 33812
2, 3, 1, 2, 1, 1, 1, // 32516
2, 3, 1, 1, 2, 1, 1, // 29060
3, 1, 1, 3, 1, 1, 1, // 28952
2, 3, 1, 1, 1, 2, 1, // 26756
2, 2, 1, 3, 1, 1, 1};// 26036 */
/* unsigned int sv[8*45]={ // length=20, live only
3, 3, 1, 1, 1, 1, 1, 1, // 149272
2, 4, 1, 1, 1, 1, 1, 1, // 140524
1, 5, 1, 1, 1, 1, 1, 1, // 136150
2, 3, 2, 1, 1, 1, 1, 1, // 117196
3, 1, 3, 1, 1, 1, 1, 1, // 114280
1, 4, 2, 1, 1, 1, 1, 1, // 112822
2, 2, 3, 1, 1, 1, 1, 1, // 105532
2, 3, 1, 2, 1, 1, 1, 1, // 101644
1, 3, 3, 1, 1, 1, 1, 1, // 101158
2, 1, 4, 1, 1, 1, 1, 1, // 99700
1, 4, 1, 2, 1, 1, 1, 1, // 97270
1, 2, 4, 1, 1, 1, 1, 1, // 95326
2, 3, 1, 1, 2, 1, 1, 1, // 91276
3, 1, 1, 3, 1, 1, 1, 1, // 90952
1, 4, 1, 1, 2, 1, 1, 1, // 86902
1, 3, 2, 2, 1, 1, 1, 1, // 85606
2, 3, 1, 1, 1, 2, 1, 1, // 84364
2, 1, 3, 2, 1, 1, 1, 1, // 84148
2, 2, 1, 3, 1, 1, 1, 1, // 82204
1, 4, 1, 1, 1, 2, 1, 1, // 79990
1, 2, 3, 2, 1, 1, 1, 1, // 79774
2, 3, 1, 1, 1, 1, 2, 1, // 79756
1, 3, 1, 3, 1, 1, 1, 1, // 77830
2, 1, 2, 3, 1, 1, 1, 1, // 76372
3, 1, 1, 1, 3, 1, 1, 1, // 75400
1, 4, 1, 1, 1, 1, 2, 1, // 75382
1, 3, 2, 1, 2, 1, 1, 1, // 75238
2, 1, 3, 1, 2, 1, 1, 1, // 73780
1, 2, 2, 3, 1, 1, 1, 1, // 71998
1, 2, 3, 1, 2, 1, 1, 1, // 69406
1, 3, 1, 2, 2, 1, 1, 1, // 67462
1, 3, 2, 1, 1, 2, 1, 1, // 68326
2, 1, 3, 1, 1, 2, 1, 1, // 66868
2, 2, 1, 1, 3, 1, 1, 1, // 66652
1, 3, 2, 1, 1, 1, 2, 1, // 63718
1, 2, 3, 1, 1, 2, 1, 1, // 62494
1, 3, 1, 1, 3, 1, 1, 1, // 62278
2, 1, 3, 1, 1, 1, 2, 1, // 62260
2, 1, 2, 1, 3, 1, 1, 1, // 60820
1, 3, 1, 2, 1, 2, 1, 1, // 60550
1, 2, 3, 1, 1, 1, 2, 1, // 57886
1, 2, 2, 1, 3, 1, 1, 1, // 56446
1, 3, 1, 2, 1, 1, 2, 1, // 55942
1, 3, 1, 1, 2, 2, 1, 1, // 55366
1, 3, 1, 1, 2, 1, 2, 1}; // 50758 */
/* unsigned int sv[9*63]={ // length=22, live only
3, 3, 1, 1, 1, 1, 1, 1, 1, // 456008
2, 4, 1, 1, 1, 1, 1, 1, 1, // 429764
1, 5, 1, 1, 1, 1, 1, 1, 1, // 416642
2, 3, 2, 1, 1, 1, 1, 1, 1, // 359780
3, 1, 3, 1, 1, 1, 1, 1, 1, // 351032
1, 4, 2, 1, 1, 1, 1, 1, 1, // 346658
2, 2, 3, 1, 1, 1, 1, 1, 1, // 324788
2, 3, 1, 2, 1, 1, 1, 1, 1, // 313124
1, 3, 3, 1, 1, 1, 1, 1, 1, // 311666
2, 1, 4, 1, 1, 1, 1, 1, 1, // 307292
1, 4, 1, 2, 1, 1, 1, 1, 1, // 300002
1, 2, 4, 1, 1, 1, 1, 1, 1, // 294170
2, 3, 1, 1, 2, 1, 1, 1, 1, // 282020
3, 1, 1, 3, 1, 1, 1, 1, 1, // 281048
1, 4, 1, 1, 2, 1, 1, 1, 1, // 268898
1, 3, 2, 2, 1, 1, 1, 1, 1, // 265010
2, 3, 1, 1, 1, 2, 1, 1, 1, // 261284
2, 1, 3, 2, 1, 1, 1, 1, 1, // 260636
2, 2, 1, 3, 1, 1, 1, 1, 1, // 254804
1, 4, 1, 1, 1, 2, 1, 1, 1, // 248162
1, 2, 3, 2, 1, 1, 1, 1, 1, // 247514
2, 3, 1, 1, 1, 1, 2, 1, 1, // 247460
1, 3, 1, 3, 1, 1, 1, 1, 1, // 241682
2, 1, 2, 3, 1, 1, 1, 1, 1, // 237308
3, 1, 1, 1, 3, 1, 1, 1, 1, // 234392
1, 4, 1, 1, 1, 1, 2, 1, 1, // 234338
1, 3, 2, 1, 2, 1, 1, 1, 1, // 233906
2, 1, 3, 1, 2, 1, 1, 1, 1, // 229532
2, 1, 1, 4, 1, 1, 1, 1, 1, // 225644
1, 2, 2, 3, 1, 1, 1, 1, 1, // 224186
1, 2, 3, 1, 2, 1, 1, 1, 1, // 216410
1, 3, 2, 1, 1, 2, 1, 1, 1, // 213170
1, 2, 1, 4, 1, 1, 1, 1, 1, // 212522
1, 3, 1, 2, 2, 1, 1, 1, 1, // 210578
2, 1, 3, 1, 1, 2, 1, 1, 1, // 208796
2, 2, 1, 1, 3, 1, 1, 1, 1, // 208148
3, 1, 1, 1, 1, 3, 1, 1, 1, // 203288
1, 3, 2, 1, 1, 1, 2, 1, 1, // 199346
1, 2, 3, 1, 1, 2, 1, 1, 1, // 195674
1, 3, 1, 1, 3, 1, 1, 1, 1, // 195026
2, 1, 3, 1, 1, 1, 2, 1, 1, // 194972
2, 1, 1, 3, 2, 1, 1, 1, 1, // 194540
2, 1, 2, 1, 3, 1, 1, 1, 1, // 190652
1, 3, 1, 2, 1, 2, 1, 1, 1, // 189842
1, 2, 3, 1, 1, 1, 2, 1, 1, // 181850
1, 2, 1, 3, 2, 1, 1, 1, 1, // 181418
2, 1, 1, 2, 3, 1, 1, 1, 1, // 178988
1, 2, 2, 1, 3, 1, 1, 1, 1, // 177530
2, 2, 1, 1, 1, 3, 1, 1, 1, // 177044
1, 3, 1, 2, 1, 1, 2, 1, 1, // 176018
1, 3, 1, 1, 2, 2, 1, 1, 1, // 174290
2, 1, 1, 3, 1, 2, 1, 1, 1, // 173804
1, 2, 1, 2, 3, 1, 1, 1, 1, // 165866
1, 3, 1, 1, 1, 3, 1, 1, 1, // 163922
1, 2, 1, 3, 1, 2, 1, 1, 1, // 160682
1, 3, 1, 1, 2, 1, 2, 1, 1, // 160466
2, 1, 1, 3, 1, 1, 2, 1, 1, // 159980
2, 1, 2, 1, 1, 3, 1, 1, 1, // 159548
1, 3, 1, 1, 1, 2, 2, 1, 1, // 150098
2, 1, 1, 2, 1, 3, 1, 1, 1, // 147884
1, 2, 1, 3, 1, 1, 2, 1, 1, // 146858
1, 2, 2, 1, 1, 3, 1, 1, 1, // 146426
1, 2, 1, 2, 1, 3, 1, 1, 1};// 134762 */
/* unsigned int sv[9*98]={ // length=23, live only
4, 3, 1, 1, 1, 1, 1, 1, 1, // 912016
3, 4, 1, 1, 1, 1, 1, 1, 1, // 859528
2, 5, 1, 1, 1, 1, 1, 1, 1, // 833284
3, 3, 2, 1, 1, 1, 1, 1, 1, // 719560
4, 1, 3, 1, 1, 1, 1, 1, 1, // 702064
2, 4, 2, 1, 1, 1, 1, 1, 1, // 693316
3, 2, 3, 1, 1, 1, 1, 1, 1, // 649576
3, 3, 1, 2, 1, 1, 1, 1, 1, // 626248
2, 3, 3, 1, 1, 1, 1, 1, 1, // 623332
3, 1, 4, 1, 1, 1, 1, 1, 1, // 614584
2, 4, 1, 2, 1, 1, 1, 1, 1, // 600004
2, 2, 4, 1, 1, 1, 1, 1, 1, // 588340
2, 1, 5, 1, 1, 1, 1, 1, 1, // 570844
3, 3, 1, 1, 2, 1, 1, 1, 1, // 564040
4, 1, 1, 3, 1, 1, 1, 1, 1, // 562096
2, 4, 1, 1, 2, 1, 1, 1, 1, // 537796
2, 3, 2, 2, 1, 1, 1, 1, 1, // 530020
3, 3, 1, 1, 1, 2, 1, 1, 1, // 522568
3, 1, 3, 2, 1, 1, 1, 1, 1, // 521272
3, 2, 1, 3, 1, 1, 1, 1, 1, // 509608
2, 4, 1, 1, 1, 2, 1, 1, 1, // 496324
2, 2, 3, 2, 1, 1, 1, 1, 1, // 495028
3, 3, 1, 1, 1, 1, 2, 1, 1, // 494920
2, 3, 1, 3, 1, 1, 1, 1, 1, // 483364
2, 1, 4, 2, 1, 1, 1, 1, 1, // 477532
3, 3, 1, 1, 1, 1, 1, 2, 1, // 476488
3, 1, 2, 3, 1, 1, 1, 1, 1, // 474616
4, 1, 1, 1, 3, 1, 1, 1, 1, // 468784
2, 4, 1, 1, 1, 1, 2, 1, 1, // 468676
2, 3, 2, 1, 2, 1, 1, 1, 1, // 467812
3, 1, 3, 1, 2, 1, 1, 1, 1, // 459064
3, 1, 1, 4, 1, 1, 1, 1, 1, // 451288
2, 4, 1, 1, 1, 1, 1, 2, 1, // 450244
2, 2, 2, 3, 1, 1, 1, 1, 1, // 448372
2, 2, 3, 1, 2, 1, 1, 1, 1, // 432820
2, 1, 3, 3, 1, 1, 1, 1, 1, // 430876
2, 3, 2, 1, 1, 2, 1, 1, 1, // 426340
2, 2, 1, 4, 1, 1, 1, 1, 1, // 425044
2, 3, 1, 2, 2, 1, 1, 1, 1, // 421156
3, 1, 3, 1, 1, 2, 1, 1, 1, // 417592
3, 2, 1, 1, 3, 1, 1, 1, 1, // 416296
2, 1, 4, 1, 2, 1, 1, 1, 1, // 415324
2, 1, 2, 4, 1, 1, 1, 1, 1, // 407548
4, 1, 1, 1, 1, 3, 1, 1, 1, // 406576
2, 3, 2, 1, 1, 1, 2, 1, 1, // 398692
2, 2, 3, 1, 1, 2, 1, 1, 1, // 391348
2, 3, 1, 1, 3, 1, 1, 1, 1, // 390052
3, 1, 3, 1, 1, 1, 2, 1, 1, // 389944
3, 1, 1, 3, 2, 1, 1, 1, 1, // 389080
3, 1, 2, 1, 3, 1, 1, 1, 1, // 381304
2, 3, 2, 1, 1, 1, 1, 2, 1, // 380260
2, 3, 1, 2, 1, 2, 1, 1, 1, // 379684
2, 1, 4, 1, 1, 2, 1, 1, 1, // 373852
3, 1, 3, 1, 1, 1, 1, 2, 1, // 371512
2, 1, 3, 2, 2, 1, 1, 1, 1, // 368668
2, 2, 3, 1, 1, 1, 2, 1, 1, // 363700
2, 2, 1, 3, 2, 1, 1, 1, 1, // 362836
3, 1, 1, 2, 3, 1, 1, 1, 1, // 357976
2, 2, 2, 1, 3, 1, 1, 1, 1, // 355060
3, 2, 1, 1, 1, 3, 1, 1, 1, // 354088
2, 3, 1, 2, 1, 1, 2, 1, 1, // 352036
2, 3, 1, 1, 2, 2, 1, 1, 1, // 348580
3, 1, 1, 3, 1, 2, 1, 1, 1, // 347608
2, 1, 4, 1, 1, 1, 2, 1, 1, // 346204
2, 1, 2, 3, 2, 1, 1, 1, 1, // 345340
2, 2, 3, 1, 1, 1, 1, 2, 1, // 345268
2, 1, 3, 1, 3, 1, 1, 1, 1, // 337564
2, 3, 1, 2, 1, 1, 1, 2, 1, // 333604
2, 2, 1, 2, 3, 1, 1, 1, 1, // 331732
2, 3, 1, 1, 1, 3, 1, 1, 1, // 327844
2, 1, 4, 1, 1, 1, 1, 2, 1, // 327772
2, 1, 3, 2, 1, 2, 1, 1, 1, // 327196
2, 2, 1, 3, 1, 2, 1, 1, 1, // 321364
2, 3, 1, 1, 2, 1, 2, 1, 1, // 320932
3, 1, 1, 3, 1, 1, 2, 1, 1, // 319960
3, 1, 2, 1, 1, 3, 1, 1, 1, // 319096
2, 1, 2, 2, 3, 1, 1, 1, 1, // 314236
2, 1, 2, 3, 1, 2, 1, 1, 1, // 303868
2, 3, 1, 1, 2, 1, 1, 2, 1, // 302500
3, 1, 1, 3, 1, 1, 1, 2, 1, // 301528
2, 3, 1, 1, 1, 2, 2, 1, 1, // 300196
2, 1, 3, 2, 1, 1, 2, 1, 1, // 299548
2, 1, 3, 1, 2, 2, 1, 1, 1, // 296092
3, 1, 1, 2, 1, 3, 1, 1, 1, // 295768
2, 2, 1, 3, 1, 1, 2, 1, 1, // 293716
2, 2, 2, 1, 1, 3, 1, 1, 1, // 292852
2, 3, 1, 1, 1, 2, 1, 2, 1, // 281764
2, 1, 3, 2, 1, 1, 1, 2, 1, // 281116
2, 1, 2, 3, 1, 1, 2, 1, 1, // 276220
2, 1, 3, 1, 1, 3, 1, 1, 1, // 275356
2, 2, 1, 3, 1, 1, 1, 2, 1, // 275284
2, 2, 1, 2, 1, 3, 1, 1, 1, // 269524
2, 1, 3, 1, 2, 1, 2, 1, 1, // 268444
2, 1, 2, 3, 1, 1, 1, 2, 1, // 257788
2, 1, 2, 2, 1, 3, 1, 1, 1, // 252028
2, 1, 3, 1, 2, 1, 1, 2, 1, // 250012
2, 1, 3, 1, 1, 2, 2, 1, 1, // 247708
2, 1, 3, 1, 1, 2, 1, 2, 1};// 229276 */
FILE *Outfp;
Outfp = fopen("out3b.dat","w");
if (c==1) {
W[0]=0; // load 1
W[1]=1;
}
else {
W[0]=0xffffffff; // load -1
W[1]=0xffffffff;
}
for (i=0; i<512; i++) { // clear histogram of limb lengths
ho[i]=0;
he[i]=0;
}
for (i=0; i<200; i++) // clear histogram of e values
dist[i]=0;
for (i=0; i<200*2; i++) // clear maximum y values
maxy[i]=0;
order=3*(1<<lshift); // order
Z[0]=0;
Z[1]=order;
X[0]=0;
X[1]=1<<lshift; // order/3
h=0; // clear limb count
iters=0; // clear display count
//
// LIVE LIMBS (ending in an odd natural number)
//
for (i=order; i>(order/2); i-=2) {
V[0]=0; // load k
V[1]=i;
flag=subr2(V[0], V[1], T);
if (flag==1)
continue;
g=0; // clear segment count
v[0]=V[0]; // save sequence value
v[1]=V[1];
n=1; // element count
while ((V[1]&1)==0) {
shift(V, V, 1);
v[2*n]=V[0];
v[2*n+1]=V[1];
n=n+1; // increment number of terms in sequence
}
for (j=1; j<100000; j++) {
if ((V[0]==0)&&(V[1]==1))
goto mskip;
if (c==-1) {
U[0]=0;
U[1]=5;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=7;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=17;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=25;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=37;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=55;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=41;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=61;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
U[0]=0;
U[1]=91;
sub64(V, U);
if ((U[0]|U[1])==0)
goto mskip;
}
flag=subr2(V[0], V[1], T); // check if 3 divides k
if (flag==1)
goto mskip;
U[0]=X[0]; // check if k>(order/3)
U[1]=X[1];
sub64(V, U);
if ((U[0]&0x80000000)==0)
goto lskip;
U[0]=V[0];
U[1]=V[1];
add64(V, V);
add64(U, V);
add64(W, V); // next sequence value
U[0]=Z[0]; // check if k>order
U[1]=Z[1];
sub64(V, U);
if ((U[0]&0x80000000)==0)
goto mskip;
if ((V[0]&0xc0000000)!=0)
goto mskip; // check for overflow
if ((V[1]&0x7)==0)
g=g+1; // increment primary connection point count
v[2*n]=V[0]; // save sequence value
v[2*n+1]=V[1];
n=n+1; // increment element count
while ((V[1]&1)==0) {
shift(V, V, 1);
v[2*n]=V[0]; // save sequence value
v[2*n+1]=V[1];
n=n+1; // increment number of terms in sequence
if (n>4095) {
printf("error: array not big enough \n");
goto mskip;
}
}
}
goto mskip;
//
// check array
//
lskip:
if (g==0) // skip segment counts of zero
goto mskip;
h=h+1;
if (n<512)
ho[n]=ho[n]+1;
else {
printf("error: ho array not big enough \n");
goto zskip;
}
if (n<12) {
printf("error: unexpected limb length \n");
goto zskip;
}
if (n!=l) // check for specified limb length
goto mskip;
if (iters<=20) { // output a few limbs
printf("limb=(%d, %d), g=%d \n",v[1],v[2*n-1],g);
iters=iters+1;
}
U[0]=v[2*n-2];
U[1]=v[2*n-1];
T[0]=U[0];
T[1]=U[1];
add64(U, U);
add64(T, U);
add64(W, U);
shift(U, U, 1);
T[0]=v[0];
T[1]=v[1];
sub64(T, U);
for (f=0; f<15; f++) {
if (n==len[f]) {
if (num[f]>0)
mul64_32(v[0], v[1], num[f], A);
else
mul64_32(v[0], v[1], -num[f], A);
B[0]=0;
B[1]=A[0];
B[2]=A[1];
B[3]=A[2];
mul64_32(U[0], U[1], den[f], A);
C[0]=0;
C[1]=A[0];
C[2]=A[1];
C[3]=A[2];
if (num[f]>0)
sub128(B, C);
else {
add128(B, C);
B[0]=0;
B[1]=0;
B[2]=0;
B[3]=0;
sub128(B, C);
}
del=(int)C[3];
for (e=0; e<p; e++) {
if (del==c*z[e]) {
savee=z[e];
dist[e]=dist[e]+1;
goto bskip;
}
}
printf("error: No solution of Diophantine equation %d \n",del);
printf("limb=%d, %d, n=%d \n",v[1],v[2*n-1],n);
goto zskip;
bskip:
//
// COMPUTE SEQUENCE VECTOR
//
Y[0]=v[0];
Y[1]=v[1];
count=0;
while ((Y[1]&1)==0) {
count=count+1;
shift(Y, Y, 1);
}
t[0]=count;
f=1;
first=1;
U[0]=v[2*n-2];
U[1]=v[2*n-1];
sub64(Y, U);
while ((U[0]|U[1])!=0) {
if ((Y[1]&1)==0) {
shift(Y, Y, 1);
count=count+1;
}
else {
U[0]=Y[0];
U[1]=Y[1];
add64(Y, U);
add64(U, Y);
add64(W, Y);
if (first==0) {
t[f]=count;
f=f+1;
if (f>199) {
printf("error: sequence vector array not big enough \n");
goto zskip;
}
}
else
first=0;
count=0;
}
U[0]=v[2*n-2];
U[1]=v[2*n-1];
sub64(Y, U);
}
t[f]=count;
f=f+1;
if (f!=a) {
printf("error: incorrect sequence vector length=%d \n",f);
goto zskip;
}
for (d=0; d<b; d++) {
for (j=0; j<a; j++) {
if (t[j]!=sv[j+a*d])
goto eskip;
}
he[d]=he[d]+1;
T[0]=v[0];
T[1]=v[1];
U[0]=maxy[2*d];
U[1]=maxy[2*d+1];
sub64(T, U);
if ((U[0]&0x80000000)==0) {
maxy[2*d]=T[0];
maxy[2*d+1]=T[1];
}
goto mskip;
eskip: j=0;
}
printf("error: sequence vector not found \n");
goto zskip;
}
}
mskip:n=0;
}
//
// check for valid lengths
//
for (i=0; i<40; i++) {
if (ho[y[i]]!=0) {
printf("error: n=%d \n",y[i]);
goto zskip;
}
}
//
// output maximum y values
//
printf("\n");
printf("maximum y values \n");
for (i=0; i<p; i++) {
T[0]=maxy[2*i];
T[1]=maxy[2*i+1];
if (lshift>16) {
shift(T, T, lshift-16);
flag=subr2(T[0], T[1], U);
}
else {
U[0]=0;
U[1]=0;
}
printf("max=%#010x %#010x, ratio=%#010x %#010x \n",maxy[2*i],maxy[2*i+1], U[0], U[1]);
}
//
// output distribution of e values
//
printf("\n");
printf("distribution of e values \n");
for (i=0; i<p; i++)
printf("count=%d, %d \n",dist[i],he[i]);
//
// output histogram
//
fprintf(Outfp,"\n");
printf("\n");
printf("expected number of entries=%d \n",order/144);
fprintf(Outfp,"HISTOGRAM: number of entries=%d \n",h);
printf("HISTOGRAM: number of entries=%d \n",h);
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<48; 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);
}