/*  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);
}