﻿ /*  FACTOR COLLATZ NUMBERS
```/*  FACTOR COLLATZ NUMBERS						      */
/*  06/09/13 (dkc)							      */
/*									      */
/*  This C program determines the prime factors of 2^(K+L)-3^K. 	      */
/*									      */
/******************************************************************************/
#include <math.h>
#include <stdio.h>
#include "table4.h"
void div128_64(unsigned int a0, unsigned int a1, unsigned int a2, unsigned int a3,
unsigned int *quotient, unsigned int d0, unsigned int d1);
void lshiftn(unsigned int *a, unsigned int *b, unsigned int shift, unsigned int n);
void addn(unsigned int *a, unsigned int *b, unsigned int n);
void subn(unsigned int *c, unsigned int *d, unsigned int n);
void copyn(unsigned int *a, unsigned int *b, unsigned int n);
void setn(unsigned int *a, unsigned int b, unsigned int n);
unsigned int orn(unsigned int *a, unsigned int n);
int main () {
unsigned int KLlim=10;	   // set to various values
unsigned int maxKL=10;
unsigned int L,K,count,histo[20];
unsigned int h,i,k,l,m;
unsigned int A[4],B[4],C[4],Q[4],Z[4],M[4],o;
FILE *Outfp;
Outfp = fopen("test23b.dat","w");
if (KLlim>maxKL) {
printf("limit too big; Collatz numbers not completely factored \n");
goto zskip;
}
for (i=0; i<20; i++)
histo[i]=0;
m=4;
setn(M, 0, m);
setn(Z, 0, m);
for (K=1; K<=KLlim; K++) {
for (L=1; L<=KLlim; L++) {
setn(A, 1, m);
for (h=0; h<K; h++) {	    // 3**K
copyn(A, B, m);
addn(A, A, m);
addn(B, A, m);
}
setn(B, 1, m);
lshiftn(B, C, K+L, m);	    // 2**(K+L)
subn(C, A, m);		    // 2**(K+L)-3**K
if ((A[0]&0x80000000)!=0)     // check if negative
subn(Z, A, m);
copyn(A, C, m);		    // find maximum
subn(Z, C, m);
addn(M, C, m);
if ((C[0]&0x80000000)!=0)
copyn(A, M, m);
count=0;
printf("L=%d, K=%d, %#10x %#10x %#10x %#10x \n",L,K,A[0],A[1],A[2],A[3]);
fprintf(Outfp,"L=%d, K=%d, %#10x %#10x %#10x %#10x \n",L,K,A[0],A[1],A[2],A[3]);
for (k=1; k<17893; k++) {     // divide by primes in LUT
l=table[k];
if (l>65536)		    // related to maximum
break;
setn(C, l, m); 	    // check if prime is larger
subn(A, C, m);
if ((C[0]&0x80000000)!=0)
break;
div128_64(A[0], A[1], A[2], A[3], Q, 0, l);
setn(B, 0, m);
for (i=0; i<l; i++)
addn(Q, B, m);
subn(A, B, m);
o=orn(B, m);		    // check if prime is factor
if (o==0) {
printf(" %d",l);
fprintf(Outfp," %d",l);
copyn(Q, A, m);
count=count+1;
loop:	    div128_64(A[0], A[1], A[2], A[3], Q, 0, l);
setn(B, 0, m);
for (i=0; i<l; i++)
addn(Q, B, m);
subn(A, B, m);
o=orn(B, m);
if (o==0) {
printf(" %d",l);
fprintf(Outfp," %d",l);
copyn(Q, A, m);
count=count+1;
goto loop;
}
}
}			    //	check for 1
setn(C, 1, m);
subn(A, C, m);
o=orn(C, m);
if (o!=0)
count=count+1;
histo[count]=histo[count]+1;
printf("\n");
fprintf(Outfp,"\n");
printf(" %#10x %#10x %#10x %#10x \n",A[0],A[1],A[2],A[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",A[0],A[1],A[2],A[3]);
printf("\n");
fprintf(Outfp,"\n");
}
}
for (i=0; i<10; i++) {
printf("i=%d, n=%d \n",i,histo[i]);
fprintf(Outfp,"i=%d, n=%d \n",i,histo[i]);
}
printf("\n");
fprintf(Outfp,"\n");
printf("maximum \n");
fprintf(Outfp,"maximum \n");
printf(" %#10x %#10x %#10x %#10x \n",M[0],M[1],M[2],M[3]);
fprintf(Outfp," %#10x %#10x %#10x %#10x \n",M[0],M[1],M[2],M[3]);
zskip:
fclose(Outfp);
return(0);
}
```