######################################################################
## GenerateRationalRecurrenceTrees: Save this file as  		     #
## 	GenerateRationalRecurrenceTrees.txt to use it,  	     #
##	open Maple and type 					     #
##	read `Dir/GenerateRationalRecurrenceTrees.txt`;		     #
## where Dir is the directory where you saved the file,		     #
## and then hit <enter>. Follow the instructions given once you hit  # 
## <enter>.         						     #
##                                                                   #
## Written by Emilie Hogan, Rutgers University ,                     #
##  eahogan at math dot rutgers dot edu.                             # 
######################################################################

print(`Started: November, 2009. Tested in Maple 10`);
print(`Current version: April 5, 2009`);
print(``);
print(`IMPORTANT!!!! You must have installed and read in (per the instructions) the file RecurrenceTree.txt`);
print(`which can also be found on my website at http://math.rutgers.edu/~eahogan/maple/RecurrenceTree.txt`);
print(``);
print(`This is GenerateRationalRecurrenceTrees, a Maple package that`);
print(`finds recurrences that are first order degree 2 in a(n) and `);
print(`that generate rational recurrence trees.`);
print(`This package accompanies the paper`);
print(` "Non-Linear Recurrences that Generate Rational Numbers"`);
print(`which can be found on my website`);
print(`The most current version of this package can be found at`);
print(`http://math.rutgers.edu/~eahogan/maple/GenerateRationalRecurrenceTrees.txt`);
print(``);
print(`Please report bugs and questions to Emilie Hogan at`);
print(`eahogan at math dot rutgers dot edu`);
print(`For a list of procedures in this package type "Help()"`);

Help := proc()
if args=NULL then
	print(`The MAIN procedures are:`);
	print(`FindRatRecTrees, FindRatRecTreesMoreSpec, FindRatRecTreesMostGen`);
	print(` FindRatRecTreesSmarter `);
	print(`WARNING: These all take a long time to run. FindRatRecTreesSmarter is the`);
	print(`fastest, but still slow for large values in the argument "co"`);
	print(`For help on a specific procedure type Help(procedure_name)`);
elif nargs=1 and args[1]=FindRatRecTrees then
	print(`FindRatRecTrees(d,co,x,y,N):`);
	print(`The general form of a recurrence that this procedure returns is`);
	print(`P1(y)*x^2+P2(y)*x+P3(y) where x is taking the place of a(n) and`);
	print(`y is taking the place of a(n-1)`);
	print(`Inputs - `);
	print(`d: a list of length 3 where d[i] is the degree of Pi in y`);
	print(`co: a list of length 3 where co[i] is the maximum absolute value`);
	print(`of a coefficient in the polynomial Pi`);
	print(`x,y: x is the variable to be a(n) and y is the variable to be a(n-1)`);
	print(`N: the number of iterations of the recurrence to check for rationality`);
	print(`Output - `);
	print(`All recurrences of the form P1(y)*x^2+P2(y)*x+P3(y) that generate a rational`);
	print(`recurrence trees.`);
	print(`For example try: FindRatRecTrees([2,1,1],[1,2,2],x,y,5)`);

elif nargs=1 and args[1]=FindRatRecTreesMoreSpec then
	print(`FindRatRecTreesMoreSpec(P1,P2,P3,vara,varb,varc,co,x,y,N):`);
	print(`Inputs - `);
	print(`P1: the polynomial in y which is the coefficient of x^2`);
	print(`P2: the polynomial in y which is the coefficient of x`);
	print(`P3: the polynomial in y which is the "constant" (w.r.t x)`);
	print(`vara: the variables (coefficients of y^powers) in P1`);
	print(`varb: the variables (coefficients of y^powers) in P2`);
	print(`varc: the variables (coefficients of y^powers) in P3`);
	print(`co: co[i] the largest absolute value to test for anything in (i=1) vara, (i=2) varb, (i=3) varc`);
	print(`x,y: indeterminants`);
	print(`N: the number of iterations of the recurrence to check for rationality.`);
	print(`Output - `);
	print(`All recurrences of the form P1(y)*x^2+P2(y)*x+P3(y) that generate a rational`);
	print(`recurrence tree.`);
	print(``);
	print(`For example try: FindRatRecTreesMoreSpec(a[1]+y^2,b[1]+b[2]*y,c[1]+c[2]*y,[a[1]],[b[1],b[2]],[c[1],c[2]],[1,2,2],x,y,4)`);

elif nargs=1 and args[1]=FindRatRecTreesMostGen then
	print(`FindRatRecTreesMostGen(P,vars,cos,x,y,N,seed):`);
	print(`Inputs - `);
	print(`P: a sequence of polynomials, P[i] will be coefficient for x^(i-1)`);
	print(`vars: vars[i] is the set of variables appearing in P[i] as coefficients of y^powers`);
	print(`cos: cos[i] is the max absolute value that a varaiable in vars[i] can have`);
	print(`x,y: indeterminants`);
	print(`N: the number of iterations of the recurrence to check for rationality`);
	print(`seed: the value for a(1) to use`);
	print(`Output - `);
	print(`All recurrences of the form sum(P[i]*x^(i-1), i=1..nops(P)) that generate a`);
	print(`rational recurrence tree.`);
	print(``);
	print(`For example try: FindRatRecTreesMostGen([a[1]+y^2,b[1]+b[2]*y,c[1]+c[2]*y],[[a[1]],[b[1],b[2]],[c[1],c[2]]],[1,2,2],x,y,4,1`);
	
elif nargs=1 and args[1]=FindRatRecTreesSmarter then
	print(`FindRatRecTreesSmarter(P1,P2,P3,vara,varb,varc,co,x,y,seed):`);
	print(`This procedure uses Proposition ??? to find all recurrences of the form`);
	print(`P1(y)*x^2+P2(y)*x+P3(y) that generate a rational recurrence tree.`);
	print(`No recurrence trees are calculated so this procedure is faster than FindRatRecTreesMoreSpec`);
	print(`Inputs - `);
	print(`P1: the polynomial in y which is the coefficient of x^2`);
	print(`P2: the polynomial in y which is the coefficient of x`);
	print(`P3: the polynomial in y which is the "constant" (w.r.t x)`);
	print(`vara: the variables (coefficients of y^powers) in P1`);
	print(`varb: the variables (coefficients of y^powers) in P2`);
	print(`varc: the variables (coefficients of y^powers) in P3`);
	print(`co: co[i] the largest absolute value to test for anything in (i=1) vara, (i=2) varb, (i=3) varc`);
	print(`x,y: indeterminants`);
	print(`seed: the value for a(1) to use`);
	print(`Outputs - (2 outputs) `);
	print(`1. All recurrences of the form P1(y)*x^2+P2(y)*x+P3(y) that generate a rational`);
	print(`recurrence tree that is nontrivial.`);
	print(`2. All recurrences of the form P1(y)*x^2+P2(y)*x+P3(y) that generate a rational`);
	print(`recurrence tree that is all the same number repeated over and over.`);
	print(`For example try: FindRatRecTreesSmarter(a[1]+y^2,b[1]+b[2]*y,c[1]+c[2]*y,[a[1]],[b[1],b[2]],[c[1],c[2]],[1,2,2],x,y,1)`);

fi;
end:

#===============================================================
#FindRatRecTrees(d,co,x,y,N):
#d: a list of length 3 where d[i] is the degree of the polynomial in
#   y which is the coefficient of x^(3-i)
#co: a list of length 3 where co[i] is the maximum absolute value
#    of a coefficient in the polynomial in y coefficient of x^(3-i)
#x,y: x is the variable to solve for, y is the variable that we "know"
#N: the number of iterations to check
#
#FindRatRectRees(d,co,x,y,N): returns all polynomials in x and y with
#degree(x)=2 which generate a recurrence tree of rational numbers.

FindRatRecTrees := proc(d,co,x,y,N)
local rec,P1,P2,P3,vara,varb,varc,a,b,c,posa,posb,posc,Ret,A,B,C,F,Tree,TF,L;

P1 := sum(a[i]*y^i,i=0..d[1]);
P2 := sum(b[i]*y^i,i=0..d[2]);
P3 := sum(c[i]*y^i,i=0..d[3]);

rec := P1*x^2+P2*x+P3;

vara := {seq(a[j], j=0..d[1])};
varb := {seq(b[j],j=0..d[2])};
varc := {seq(c[j],j=0..d[3])};
posa := GenerateCos(nops(vara),co[1]);
posb := GenerateCos(nops(varb),co[2]);
posc := GenerateCos(nops(varc),co[3]);

Ret := {};
for A in posa do
for B in posb do
for C in posc do
	F:=subs({seq(a[i]=A[i+1], i=0..d[1]),seq(b[i]=B[i+1],i=0..d[2]),seq(c[i]=C[i+1],i=0..d[3])},rec);
	#print(F);
	Tree := CalculateTreeRecurrenceOrder1(1,F,2,x,y,N);
	#print(Tree);
	L := Listify(Tree);
	TF := {seq(evalb(L[k]=false) or evalb(type(L[k],rational)), k=1..nops(L))};
	if TF={true} then
		Ret := Ret union {F};
	fi;
od;
od;
od;

RETURN(Ret);

end:



#===============================================================
#P1: the polynomial in y which is the coefficient of x^2
#P2: the polynomial in y which is the coefficient of x
#P3: the polynomial in y which is the "constant" (no x in it)
#vara: the variables in P1
#varb: the variables in P2
#varc: the variables in P3
#co: the largest absolute value to test for anything in vara,varb,varc
#x,y: indeterminants
#N: the number of iterations of the tree to check if the tree is rational.

FindRatRecTreesMoreSpec := proc(P1,P2,P3,vara,varb,varc,co,x,y,N)
local rec,posa,posb,posc,Ret,A,B,C,F,Tree,TF,L;

rec := P1*x^2+P2*x+P3;

posa := GenerateCos(nops(vara),co[1]);
posb := GenerateCos(nops(varb),co[2]);
posc := GenerateCos(nops(varc),co[3]);

Ret := {};
for A in posa do
for B in posb do
for C in posc do
	#print(A,B,C);
	if {op(B)}={0} or {op(C)}={0} then
		#print(`bad`);
		#do nothing
	else
		F:=subs({seq(vara[i]=A[i], i=1..nops(vara)),seq(varb[i]=B[i],i=1..nops(varb)),seq(varc[i]=C[i],i=1..nops(varc))},rec);
		#print(F);
		Tree := CalculateTreeRecurrenceOrder1(1,F,2,x,y,N);
		#print(Tree);
		L := Listify(Tree);
		TF := {seq(evalb(L[k]=false) or evalb(type(L[k],rational)), k=1..nops(L))};
		#print(TF);
		if TF={true} then
			#print(F);
			#print(Tree);
			#print(TF);
			Ret := Ret union {F};
		fi;
	fi;
od;
od;
od;

RETURN(Ret);

end:

#===============================================================
FindRatRecTreesSmarter := proc(P1,P2,P3,vara,varb,varc,co,x,y,seed)
local rec,posa,posb,posc,Ret,A,B,C,F,first,genSol,try1,try2,aTry1,aTry2,RetG,RetRep;

rec := P1*x^2+P2*x+P3;

posa := GenerateCos(nops(vara),co[1]);
posb := GenerateCos(nops(varb),co[2]);
posc := GenerateCos(nops(varc),co[3]);

RetG := {};
RetRep := {};
for A in posa do
for B in posb do
for C in posc do
	if {op(B)}={0} or {op(C)}={0} then
		#print(`bad`);
		#do nothing
	else
		F:=subs({seq(vara[i]=A[i], i=1..nops(vara)),seq(varb[i]=B[i],i=1..nops(varb)),seq(varc[i]=C[i],i=1..nops(varc))},rec);
		#print(F);
		#print(F,`First`,first);
		first := [solve(subs(y=seed,F)=0,x)];
		if (not nops(first)=0) and {seq(type(first[i],rational),i=1..nops(first))}={true} then
			genSol := [solve(F=0,x)];
			#print(genSol);
			try1 := subs(y=genSol[1],F);
			try2 := subs(y=genSol[2],F);
			#aTry1 := [solve(try1=0,x)];
			#aTry2 := [solve(try2=0,x)];
			if (subs(x=y,try1)=0 and subs(x=y,try2)=0) or {op(first)}={seed} then
				if {op(first)}={seed} then
					RetRep := RetRep union {F};
				else
					RetG := RetG union {F};
				fi; 
			else 
				#print(`F=`,F);
				#print(`GENsol=`,genSol);
				#print(CalculateTreeRecurrenceGen(1,F,2,x,y,3));
				#print(`ATry1=`,aTry1);
				#print(`ATry2=`,aTry2);
			fi;
		fi;
	fi;
od;
od;
od;

RETURN(RetG,RetRep);

end:

#===============================================================
FindRatRecTreesMostGen := proc(P,vars,cos,x,y,N,seed)
local rec,posa,posb,posc,Ret,A,B,C,F,Tree,TF,pos,i,coeffs,coef,L;

if not(evalb(nops(P)=nops(vars))) or not(evalb(nops(vars)=nops(cos))) then
	print(`P,vars,cos need to be all the same size`);
fi;

rec := sum(P[i]*x^(i-1),i=1..nops(P));

pos := [seq([op(GenerateCos(nops(vars[i]),cos[i]))], i=1..nops(vars))];

coeffs := OneFromEach(pos);

Ret := {};
for coef in coeffs do
	if true in {seq(evalb({op(coef[i])}={0}), i=1..nops(coef))} then
		#print(`bad`);
		#do nothing
	else
		F:=subs({seq(seq(vars[j][k]=coef[j][k] ,k=1..nops(coef[j])) ,j=1..nops(coef))},rec);
		#print(F);
		Tree := CalculateTreeRecurrenceOrder1(seed,F,nops(P)-1,x,y,N);
		#print(evalf(Tree));
		L := Listify(Tree);
		TF := {seq(evalb(L[k]=false) or evalb(type(L[k],rational)), k=1..nops(L))};
		#print(TF);
		if TF={true} then
			#print(F);
			#print(Tree);
			#print(TF);
			Ret := Ret union {F};
		fi;
	fi;

od;

RETURN(Ret);

end:

#===============================================================
#n: the length of each of the lists in the output
#co: the max absolute value of each of the entries of the output lists

GenerateCos := proc(n,co)
local CosOneLess,Ret,coe;

if n=0 then
	RETURN([]);
fi;
if n=1 then
	RETURN({seq([i], i=-co..co)});
fi;

CosOneLess := GenerateCos(n-1,co);

Ret := {};
for coe in CosOneLess do
	Ret := Ret union {seq([op(coe),i],i=-co..co)};
od;

RETURN(Ret);

end:

#===============================================================
OneFromEach := proc(sets)
local prev,Ret,p;

if nops(sets)=1 then
	RETURN([seq([s],s in sets[1])]);
fi;

prev := OneFromEach([seq(sets[i], i=1..nops(sets)-1)]);

Ret := [];
for p in prev do
	Ret := [op(Ret),seq([op(p),s],s in sets[nops(sets)])];
od;

RETURN(Ret);

end:

#===============================================================
Listify := proc(Tree)
local L,level;

L := [];
for level in Tree do
	if type(level,constant) then
		L := [op(L),level];
	else
		L := [op(L),op(Listify(level))];
	fi;
od;

RETURN({op(L)});

end: