/* This document was written by Matt Tracy and Bart Snapp */
/* It should be pointed out that this code was merely for experimentation purposes.  */
/* Hence it may be a bit rough around the edges.  */
/* To use the code, start up pari/gp and run:  */

/* read(collatz) */



/* =================================================== */

polyMod(fg,n) = {    /* reduces the coefficients of a fraction of polys fg mod n */
	   local(f=numerator(fg),g=denominator(fg),v=[],X=[],x=variable(f));
	   for(i = 0, poldegree(f),
      	   	 v = concat(v,Mod(polcoeff(f,i),n)));
	   for(i = 0, poldegree(f),
      	   	 X = concat(X,x^i));
           X = mattranspose(X);
	   f=(v * X);
	   v = [];
	   X = [];
	   for(i = 0, poldegree(g),
      	   	 v = concat(v,Mod(polcoeff(g,i),n)));
	   for(i = 0, poldegree(g),
      	   	 X = concat(X,x^i));
           X = mattranspose(X);
	   g=(v * X);
	   return(lift(f/g));
}
addhelp(polyMod,"polyMod(f,n): computes the image of a polynomial mod a prime n.")




/* =================================================== */

polyGen(X,n) = {  /* currently only works for prime n */
	                 local(v = seqOpC(X),f=1);
             		 for(i=1, length(v),
                         	    if(v[i]==0,f=x*f,f = polyMod((f-1)/(x+1),n)));
             return(f)
}
/* this code will try to generate a polynomial if it can    */
/* HOWEVER it might output a fraction of polynomials        */

addhelp(polyGen,"polyGen(m,n): computes a polynomial with coefficients mod n, whose collatz sequence is similar to the sequence of the integer m.")

/* =================================================== */

opC(x)=if(x %2 == 0, 0, 1) /* 0 means divide by 2, 1 means mult by 3 etc */
addhelp(opC,"opC(x): displays the pattern of ups and downs of the collatz function.")

/* =================================================== */

C(x)=if(x %2 ==0, x/2, 3*x+1) /* Basic mapping */
addhelp(C,"C(x): is the collatz function.")

polyC(f,n) = { /* Basic polynomial mapping */
           local(x=variable(f),k = n - polcoeff(f,0));
	   if(polyMod(f,n) % x == 0, polyMod(f,n)/x, polyMod(f*(x+1)+k,n))} 
addhelp(polyC,"polyC(f,n): the collatz function on polynomials mod n.")


/* =================================================== */


T(x)=if(x %2 ==0, x/2, (3*x+1)/2) /* Basic mapping with the division of 2 */
addhelp(T,"T(x): the abbreviated collatz function.")


/* =================================================== */

polyT(f,n) = { /* Basic polynomial mapping with the division of x */
           local(x=variable(f),k = n - polcoeff(f,0));
	   if(polyMod(f,n) % x == 0, polyMod(f,n)/x, polyMod(f*(x+1)+k,n)/x)} 
addhelp(polyT,"polyT(f,n): the abbreviated collatz function on polynomials mod n.")

/* =================================================== */

seqC(x)={
                local(v=[]);
                while(x!=1,
                                v=concat(v,x);
                                x=C(x);
				for(i = 1, #v, if(v[i] == x,return(concat(loop,v))));
                		if(x == 1 || x == -1,return(concat(noloop,v))))
                                }        /* for instance, we will input 5, and output[noloop, 5, 16, 8, 4, 2]*/
addhelp(seqC,"seqC(x): returns the sequence given by the collatz function.")

/* =================================================== */

polySeqC(f,n)={	
		local(v=[]);
                while(poldegree(f) != 0,
				v = concat(v,f);
                                f = polyC(f,n);
                                for(i = 1, #v, if(v[i] == f,return(concat(loop,v))));
                		if(poldegree(f) == 0,return(concat(noloop,v))))
} 
					/*This will return the sequence generated by a polynomial.*/
					/*If it returns to 1, the program returns [noloop,seq]*/
					/*If it does not, the program returns [loop,seq]*/
addhelp(polySeqC,"polySeqC(x,n): returns the sequence given by the collatz function on polynomials mod n.")

/* =================================================== */


seqT(x)={
                local(v=[]);
                while(x!=1,
                                v=concat(v,x);
                                x=T(x);
				for(i = 1, #v, if(v[i] == x,return(concat(loop,v))));
                		if(x == 1 || x == -1,return(concat(noloop,v))))
                                }        /* for instance, we will input 5, and output[noloop, 5, 8, 4, 2]*/
addhelp(polySeqC,"seqT(x,n): returns the sequence given by the abbreviated collatz function.")

polySeqT(f,n)={	
		local(origf = f,v=[]);
                while(poldegree(f) != 0,
				v = concat(v,f);
                                f = polyT(f,n);
                                for(i = 1, #v, if(v[i] == f,return(concat(loop,v))));
                if(poldegree(f) == 0,return(concat(noloop,v))))
}
					/*This will return the sequence generated by a polynomial.*/
					/*If it returns to 1, the program returns [noloop,seq]*/
					/*If it does not, the program returns [loop,seq]*/
addhelp(polySeqC,"polySeqC(x,n): returns the sequence given by the abbreviated collatz function on polynomials mod n.")

seqOpC(x)={
                local(v=[]);
                while(x!=1,
                                v=concat(opC(x),v);
                                x=C(x);
                                );
                return(v)}	  /* this code outputs a vector consisting of zeros and ones  */



/* =================================================== */




lengthC(x) = #seqC(x)   /* computes the number of times it takes to reach down to 1*/


polyLengthC(f,n)={
		local(v=polySeqC(f,n));
		[v[1],#polySeqC(f,n)-1]
		}
					 /*This will return the length of a sequence generated by a polynomial.*/
					/*If it returns to 1, the program returns [noloop,count]*/
					/*If it does not, the program returns [loop,count]*/


lengthT(x) = #seqC(x)   /* computes the number of times it takes to reach down to 1*/


polyLengthT(f,n)={
		local(v=polySeqT(f,n));
		[v[1],#polySeqT(f,n)-1]
		}
					 /*This will return the length of a sequence generated by a polynomial.*/
					/*If it returns to 1, the program returns [noloop,count]*/
					/*If it does not, the program returns [loop,count]*/



/* =================================================== */

maxC(x)=vecmax(seqC(x)) /* this computes the max number of the sequence*/

maxT(x)=vecmax(seqT(x)) /* this computes the max number of the sequence*/

/* =================================================== */

/* The following two functions list i and the length of i, going all
the way up to some number x. It is not clear if they are important,
but we will use them for our next function. */

listLengthsC(x)={
                local(m=matrix(0,0));
                for(n=1,x,
                                m=concat(m,[n;lengthC(n)]);
                                );
                return(mattranspose(m))} /* Here we "flip" the final matrix on its side using "mattranspose" */


/* this won't work - we'll have to think about this one*/
polyListLengthsC(f,n)={
                local(m=matrix(0,0));
                for(a=1,n,
                                m=concat(m,[a;polyLengthC(f,n)]);
                                );
                return(mattranspose(m))} /* Here we "flip" the final matrix on its side using "mattranspose" */

listLengthsT(x)={
                local(m=matrix(0,0));
                for(n=1,x,
                                m=concat(m,[n;lengthT(n)]);
                                );
                return(mattranspose(m))} /* Here we "flip" the final matrix on its side using "mattranspose" */

polyListLengthsT(f,n)={
                local(m=matrix(0,0));
                for(a=1,n,
                                m=concat(m,[a;polyLengthT(f,n)]);
                                );
                return(mattranspose(m))} /* Here we "flip" the final matrix on its side using "mattranspose" */

/* =================================================== */

maxLengthC(x)={
                local(m = listLengthsC(x),n=1);
                while(m[n,2] != vecmax(m[,2]),
                             n=n+1;
                             );
                return(m[n,]);

        }

polyMaxLengthC(x)={
                local(m = listLengthsC(x),n=1);
                while(m[n,2] != vecmax(m[,2]),
                             n=n+1;
                             );
                return(m[n,]);

        }


maxLengthT(x)={
                local(m = listLengthsT(x),n=1);
                while(m[n,2] != vecmax(m[,2]),
                             n=n+1;
                             );
                return(m[n,]);

        }


/* =================================================== */

graphLengthC(n)=ploth(x=1,n,lengthC(floor(x)))




/* Here we copied code from Ariel Pacetti */


\\ ---------------  GP code  ---------------------------------------
\\
\\ Time-stamp: <Mon, Jan 13, 2003 - 14:17:26 - apacetti>
\\
\\ Description: Some general routines, including some vectors
\\ necessary routines.
\\
\\
\\ Original Author:         Ariel Pacetti
\\			    apacetti@math.utexas.edu
\\			    University of Texas at Austin
\\
\\ Created:             Wed Oct  4  2000
\\
\\ http://www.ma.utexas.edu/users/villegas/cnt/cnt-no-frames.html
\\======================================================================
\\ Given a number N, and a base b, gives the vector representation of
\\ N in base b.

numToBase(N,b)=
	{local(n,ans);
	ans=[];
	while(N>=b,
		ans=concat(N%b,ans);
		N=(N-N%b)/b);
	ans=concat(N,ans)}



\\ end of copied code
\\======================================================================




polyListGen(d,n) = {
		 local(v=[]);
		 for(i = n+1,n^(d+1)-1,
		 v = concat(v,Pol(numToBase(i,n))));
		 v
}


polyLoopCheck(d,n) = {
		   local(V = polyListGen(d,n),v=[]);
		   for(i = 1,#V,
		   if(polySeqT(V[i],n)[1] == loop, v = concat(v,V[i])));
}		   v


polyNoLoopCheck(d,n) = {
		   local(V = polyListGen(d,n),v=[]);
		   for(i = 1,#V,
		   if(polySeqT(V[i],n)[1] == noloop, v = concat(v,V[i])));
}		   v



polcoefficients(f) = {
		  local(v = []);		  
		  for(i = 0, poldegree(f),
      		  v = concat(polcoeff(f,i),v));
      		  v
}

baseToNum(v,b) = {
	       local(n=0,d = #v);
	       for(i = 1, d, 
	       n = n+b^(i-1)*v[d-i+1]);
	       n
}

polyLoopEval(d,n,b) = {
	   local(V = polyLoopCheck(d,n),m = #V,v = []);
	   for(i = 1,m,
	   v = concat(v,baseToNum(polcoefficients(V[i]),b)));
	   v
}

polyNoLoopEval(d,n,b) = {
	   local(V = polyNoLoopCheck(d,n),m = #V,v = []);
	   for(i = 1,m,
	   v = concat(v,baseToNum(polcoefficients(V[i]),b)));
	   v
}