/******************************************************************************/
/* */
/* COMPUTE 3N+1 OR 3N-1 SEQUENCE */
/* 02/13/10 (dkc) */
/* */
/* This C program finds long limbs in S. Sequence vectors are generated. */
/* */
/******************************************************************************/
#include <stdio.h>
#include <math.h>
void div256_32(unsigned int a0, unsigned int a1, unsigned int a2,
unsigned int a3, unsigned int a4, unsigned int a5,
unsigned int a6, unsigned int a7, unsigned int *quotient,
unsigned int d7);
void add256(unsigned int *a, unsigned int *b);
void sub256(unsigned int *c, unsigned int *d);
void shift256(unsigned int *a, unsigned int *b, unsigned int shift);
void lshift256(unsigned int *a, unsigned int *b, unsigned int shift);
void copy256(unsigned int *a, unsigned int *b);
void set256(unsigned int *a, unsigned int b);
int main () {
//
int c=-1; // c
//
unsigned int y[178]={0,3,6,8,11, // I (invalid lengths)
13,16,19,21,24, // I
26,29,32,34,37, // I
39,42,44,47,50, // J
52,55,57,60,63, // J
65,68,70,73,75, // K
78,81,83,86,88, // K
91,94,96,99,101, // K
104,106,109,112,114, // L
117,119,122,125,127, // L
130,132,135,138,140, // L
143,145,148,150,153, // M
156,158,161,163,166, // M
171,174,176,179, // M
184,187,189,192,194, // I'
197,200,202,205,207, // I'
210,212,215,218,220, // J'
223,225,228,231,233, // J'
236,238,241,243,246, // K'
249,251,254,256,259, // K'
262,264,267,269,272, // K'
275,277,280,282,285, // K'
287,290,293,295,298, // L'
300,303,306,308,311, // L'
313,316,318,321,324, // M'
326,329,331,334,337, // M'
339,342,344,347,349, // N
352,355,357,360,362, // N
365,368,370,373,375, // N
378,380,383,386,388, // O
391,393,396,399,401, // O
404,406,409,412,414, // O
417,419,422,424,427, // P
430,432,435,437,440, // P
443,445,448,450,453, // P
458,461,463,466}; // I
unsigned int g,h,i,count,scount,temp; // indices, etc.
unsigned int O[8],P[8],Q[8],R[8],U[8],V[8],X[8],W[8],Y[8],Z[8];
unsigned int first,sum,save,len,sv[2048];
FILE *Outfp;
Outfp = fopen("out11.dat","w");
if (c==1)
set256(Z, 1);
else {
Z[0]=0xffffffff;
Z[1]=0xffffffff;
Z[2]=0xffffffff;
Z[3]=0xffffffff;
Z[4]=0xffffffff;
Z[5]=0xffffffff;
Z[6]=0xffffffff;
Z[7]=0xffffffff;
}
//
// begin loop
//
for (h=2; h<256; h++) {
set256(Y, 1);
for (i=0; i<h; i++) { // for (i=0; i<h; i++)
copy256(Y, X);
add256(Y, X);
add256(X, Y); // n=n*3;
if (Y[0]>0x40000000) {
printf("h=%d, order too big \n",h);
goto zskip;
}
}
//
copy256(Z, V);
sub256(Y, V);
shift256(V, Y, 1); // n=(n-c)/2;
if ((Y[7]&1)==0)
continue;
//
count=2*(h+1);
g=h; // power of 3 count
lshift256(Y, U, 2);
set256(O, 3);
temp=0;
copy256(U, W);
sub256(O, W);
while ((W[0]&0x80000000)!=0) { // while (order<(4*n))
add256(O, O); // order=order*2;
temp=temp+1;
copy256(U, W);
sub256(O, W);
}
//
div256_32(O[0],O[1],O[2],O[3],O[4],O[5],O[6],O[7],W, 3); // order/3
copy256(W, U);
sub256(Y, U);
while ((U[0]&0x80000000)!=0) { // while (n<(order/3)) {
g=g+1;
set256(U, 1);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0) // if (n==1)
goto askip; // goto askip;
if (c==-1) {
set256(U, 5);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 7);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 17);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 25);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 37);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 55);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 41);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 61);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
set256(U, 91);
sub256(Y, U);
if ((U[0]|U[1]|U[2]|U[3]|U[4]|U[5]|U[6]|U[7])==0)
goto askip;
}
copy256(Y, U);
add256(Y, Y); // 2*n
add256(U, Y); // 3*n
add256(Z, Y); // n=3*n+c
count=count+1;
if ((Y[7]&0x7)==0)
goto askip; // goto askip;
while ((Y[7]&1)==0) {
shift256(Y, Y, 1); // n=n/2;
count=count+1;
} // }
copy256(W, U);
sub256(Y, U);
} // }
//
g=g+1;
copy256(Y, P);
set256(Y, 1);
lshift256(Y, Y, h);
copy256(Z, W);
sub256(Y, W);
copy256(W, Y);
add256(Y, Y); // n=n*2;
count=count+1;
shift256(O, W, 1); // order/2
copy256(Y, V);
sub256(W, V);
while ((V[0]&0x80000000)==0) { // while (n<(order/2)) {
copy256(Z, U);
sub256(Y, U); // n-c
div256_32(U[0],U[1],U[2],U[3],U[4],U[5],U[6],U[7],Q,3); // (n-c)/3
copy256(Q, R);
add256(Q, R); // ((n-c)/3)*2
add256(Q, R); // ((n-c)/3)*3
sub256(U, R);
if ((R[0]|R[1]|R[2]|R[3]|R[4]|R[5]|R[6]|R[7])==0) {
g=g+1;
copy256(Q, Y);
add256(Y, Y); // n=n*2;
count=count+2;
} // }
else {
add256(Y, Y); // n=n*2;
count=count+1;
}
copy256(Y, V);
sub256(W, V);
}
printf("order=(%#10x %#10x %#10x %#10x \n",O[0],O[1],O[2],O[3]);
printf(" %#10x %#10x %#10x %#10x), length=%d, h=%d, g=%d \n",O[4],O[5],O[6],O[7],count,h,g);
printf("limb= (%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
printf(" %#10x %#10x %#10x %#10x, \n",Y[4],Y[5],Y[6],Y[7]);
printf(" %#10x %#10x %#10x %#10x \n",P[0],P[1],P[2],P[3]);
printf(" %#10x %#10x %#10x %#10x) \n",P[4],P[5],P[6],P[7]);
fprintf(Outfp,"order=(%#10x %#10x %#10x %#10x \n",O[0],O[1],O[2],O[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x), length=%d, h=%d, g=%d \n",O[4],O[5],O[6],O[7],count,h,g);
fprintf(Outfp,"limb= (%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x, \n",Y[4],Y[5],Y[6],Y[7]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",P[0],P[1],P[2],P[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x) \n",P[4],P[5],P[6],P[7]);
for (i=0; i<178; i++) {
if (count==y[i]) {
printf("invalid length=%d \n",count);
fprintf(Outfp,"invalid length=%d \n",count);
goto zskip;
}
}
scount=count+1;
//
// COMPUTE SEQUENCE VECTOR
//
i=0;
if ((Y[7]&1)==0) {
count=1;
while ((Y[7]&1)==0) {
count=count+1;
shift256(Y, Y, 1);
}
sv[0]=count;
i=1;
}
first=1;
copy256(P, W);
sub256(Y, W);
while ((W[0]|W[1]|W[2]|W[3]|W[4]|W[5]|W[6]|W[7])!=0) {
if ((Y[7]&1)==0) {
shift256(Y, Y, 1);
count=count+1;
}
else {
copy256(Y, W);
add256(Y, W);
add256(W, Y);
add256(Z, Y);
if (first==0) {
sv[i]=count;
i=i+1;
if (i>2047) {
printf("error: sequence vector array not big enough \n");
goto zskip;
}
}
else
first=0;
count=0;
}
copy256(P, W);
sub256(Y, W);
}
sv[i]=count;
i=i+1;
//
// COMPUTE X, Y, AND Z
//
len=i;
if (len!=g) {
printf("error: len=%d, g=%d \n",len,g);
goto zskip;
}
sum=0;
for (i=0; i<len; i++)
sum=sum+sv[i];
save=sum;
if ((sum+len)!=scount) {
printf("error: length=%d, sum=%d \n",scount-1,sum+len);
goto zskip;
}
set256(X, 1);
for (i=0; i<sum; i++)
add256(X, X);
set256(Y, 1);
for (i=0; i<len; i++) {
copy256(Y, W);
add256(Y, W);
add256(W, Y);
}
set256(U, 0);
for (h=0; h<len; h++) {
set256(W, 1);
for (i=0; i<h; i++) {
copy256(W, V);
add256(W,V);
add256(V,W);
}
sum=0;
for (i=0; i<(len-h); i++)
sum=sum+sv[i];
lshift256(W, W, sum);
add256(W, U);
}
shift256(U, U, 1);
printf(" X=%#10x %#10x %#10x %#10x \n",X[0],X[1],X[2],X[3]);
printf(" %#10x %#10x %#10x %#10x, exp=%d \n",X[4],X[5],X[6],X[7],save);
printf(" Y=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
printf(" %#10x %#10x %#10x %#10x, exp=%d \n",Y[4],Y[5],Y[6],Y[7],len);
printf(" Z=%#10x %#10x %#10x %#10x \n",U[0],U[1],U[2],U[3]);
printf(" %#10x %#10x %#10x %#10x \n",U[4],U[5],U[6],U[7]);
fprintf(Outfp," X=%#10x %#10x %#10x %#10x \n",X[0],X[1],X[2],X[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x,exp=%d \n",X[4],X[5],X[6],X[7],save);
fprintf(Outfp," Y=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x, exp=%d \n",Y[4],Y[5],Y[6],Y[7],len);
fprintf(Outfp," Z=%#10x %#10x %#10x %#10x \n",U[0],U[1],U[2],U[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",U[4],U[5],U[6],U[7]);
sub256(X, Y);
printf(" numer=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
printf(" %#10x %#10x %#10x %#10x \n",Y[4],Y[5],Y[6],Y[7]);
fprintf(Outfp," numer=%#10x %#10x %#10x %#10x \n",Y[0],Y[1],Y[2],Y[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",Y[4],Y[5],Y[6],Y[7]);
temp=temp-1;
shift256(U, U, temp); // Z/(order/6)
div256_32(U[0],U[1],U[2],U[3],U[4],U[5],U[6],U[7],W,3); // Z/(order/2)
if ((Y[0]&0x80000000)!=0) {
set256(X, 0);
sub256(X, Y);
}
sub256(Y, W);
printf(" delta=%#10x %#10x %#10x %#10x \n",W[0],W[1],W[2],W[3]);
printf(" %#10x %#10x %#10x %#10x \n",W[4],W[5],W[6],W[7]);
fprintf(Outfp," delta=%#10x %#10x %#10x %#10x \n",W[0],W[1],W[2],W[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",W[4],W[5],W[6],W[7]);
printf("\n");
fprintf(Outfp,"\n");
askip:
count=0;
}
zskip:
fclose(Outfp);
return(0);
}