with(numtheory):

######################################################################
## RecurrenceTree: Save this file as  RecurrenceTree.txt to use it,  #
## open Maple and type 						     #
##	read `Dir/RecurrenceTree.txt`;				     #
## where Dir is the directory where you saved the RecurrenceTree.txt,#
## 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: March 28, 2009. Tested in Maple 10`);
print(`Current version: February 2, 2010. Tested in Maple 13`);
print(`This is RecurrenceTree, a Maple package that computes`);
print(`recurrence trees of non-linear recurrences of the form:`);
print(`sum(P_i(a(n-1),...,a(n-k))*a(n)^i,i=0..m)`);
print(`where P_i are polynomials, and with prescribed initial conditions.`);
print(``);
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/RecurrenceTree.txt`);
print(``);
print(`The data structure for a recurrence tree is as follows:`);
print(`A tree, T, is a list of levels of the tree. The value on the first level`);
print(`is T[1], i.e.- T[1] is not a list, it is just a number. The second`);
print(`level of the tree is T[2] and is a list of Deg values where Deg`);
print(`is the degree of a(n) in the recurrence. The third level of the`); 
print(`tree is T[3] and is a list of Deg lists of size Deg, etc.`);
print(`For example: where Deg=2`);
print(`T=[ a(1) , [a(2),a(2)] , [[a(3),a(3)],[a(3),a(3)]] , [[[a(4),a(4)],[a(4),a(4)]],[[a(4),a(4)],[a(4),a(4)]]],...]`);
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(`OneTreeIterOrder1, CalculateTreeRecurrenceOrder1`);
	print(`OneTreeIterOrder1NoReps, CalculateTreeRecurrenceOrder1NoReps`);
	print(`OneTreeIterOrderK, CalculateTreeRecurrecneOrderK`);
	print(`OneTreeIterOrderKNoReps, CalculateTreeRecurrenceOrderKNoReps`);
	print(`MakeSequences, JustSequences`);
	print(`CheckDiv,OnlyZ`);
	print(`For help on a specific procedure type Help(procedure_name)`);

elif nargs = 1 and args[1]=OneTreeIterOrder1 then
	print(`OneTreeIterOrder1(F,Deg,solvar,subvar,y):`);
	print(`Inputs - `);
	print(`F: the recurrence as a function in two variables, subvar and solvar`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: takes the place of a(n-1) in the recurrence F`);
	print(`Deg: the degree of solvar in F`);
	print(`y: the nth level of the recurrence tree`);
	print(`Output - `);
	print(`The n+1st level of the recurrence tree. The children of a specific`);
	print(`value, y_value in y, are the solutions of subs(subvar=y_value,F) when`);
	print(`solving for solvar WITH multiplicity. If there is no solution for`);
	print(`a specific value in y then all the children of that value are FAIL.`);	
	print(``);
	print(`For example try: OneTreeIterOrder1(Y^2*X^2+(1-4*Y)*X+(1+Y),2,X,Y,[[1, 3/4], [2, 1]])`);

elif nargs = 1 and args[1]=CalculateTreeRecurrenceOrder1 then
	print(`CalculateTreeRecurrenceOrder1(Init,F,Deg,solvar,subvar,N):`);
	print(`Inputs - `);
	print(`Init: the initial condition for the recurrence, a(1)=Init`);
	print(`F: the recurrence as a function in two variables, subvar and solvar`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: takes the place of a(n-1) in the recurrence F`);
	print(`Deg: the degree of solvar in F`);
	print(`N: the number of levels you want to compute`);
	print(`Output - `);
	print(`The recurrence tree up to level N`);
	print(``);
	print(`For example try: CalculateTreeRecurrenceOrder1(1,Y^2*X^2+(1-4*Y)*X+(1+Y),2,X,Y,5)`);

elif nargs = 1 and args[1]=OneTreeIterOrder1NoReps then
	print(`OneTreeIterOrder1NoReps(F,Deg,solvar,subvar,y,Already): `);
	print(`Inputs - `);
	print(`F: the recurrence as a function in two variables, subvar and solvar`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: takes the place of a(n-1) in the recurrence F`);
	print(`Deg: the degree of solvar in F`);
	print(`y: the nth level of the recurrence tree`);
	print(`Already: Numbers that have already been seen in the recurrence tree.`);
	print(`Output - (2 outputs)`);
	print(`1. The n+1st level of the recurrence tree. The children of a specific`);
	print(`value, y_value in y, are the solutions of subs(subvar=y_value,F) when`);
	print(`solving for solvar WITH multiplicity. If there is no solution for`);
	print(`a specific value in y then all the children of that value are FAIL.`);
	print(`If a solution has already been seen (is in the set Already) then "false"`);
	print(`is put in its place. The children of "false" are "false".`);
	print(`2. The new set for Already (Already union {all values seen in the next level)`);
	print(``);
	print(`For example try: OneTreeIterOrder1NoReps(Y^2*X^2+(1-4*Y)*X+(1+Y),2,X,Y,[1],{1})`);

elif nargs = 1 and args[1]=CalculateTreeRecurrenceOrder1NoReps then
	print(`CalculateTreeRecurrenceOrder1NoReps(Init,F,Deg,solvar,subvar,N):`);
	print(`Inputs - `);
	print(`Init: the initial condition for the recurrence, a(1)=Init`);
	print(`F: the recurrence as a function in two variables, subvar and solvar`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: takes the place of a(n-1) in the recurrence F`);
	print(`Deg: the degree of solvar in F`);
	print(`N: the number of levels you want to compute`);
	print(`Output - `);
	print(`The recurrence tree up to level N`);
	print(`Each number in the recurrence tree occurs only on one level. If it is`);
	print(`repeated then "false" is displayed instead.`);
	print(``);
	print(`For example try: CalculateTreeRecurrenceOrder1NoReps(1,Y^2*X^2+(1+5*Y)*X+(8+Y),2,X,Y,5)`);

elif nargs = 1 and args[1]=OneTreeIterOrderK then
	print(`OneTreeIterOrderK(F,n,deg,ord,solvar,subvar,TK):`);
	print(`Inputs - `);
	print(`F: the recurrence as a function in variables solvar and {subvar[n-1],...,subvar[n-ord]}`);
	print(`n: a symbol used in the set of subvars`);
	print(`Deg: the degree of solvar in F`);
	print(`ord: the order of the recurrence, how many previous terms in the tree are required`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: the variable from which the set of subvars is created`);
	print(`TK: the previous K=ord levels of the recurrence tree`);
	print(`Output - `);
	print(`The next level of the recurrence tree.`);
	print(``);
	print(`For example try: `);
	print(`F:= X^2*s[n-3]^2+s[n-1]^3*s[n-3]+s[n-1]^2*s[n-2]^2-(4*s[n-1]*s[n-2]*s[n-3]-s[n-2]^3)*X`);
	print(`OneTreeIterOrderK(F,n,2,3,X,s,[1,[1,1],[[1,1],[1,1]]],2)`);

elif nargs = 1 and args[1]=CalculateTreeRecurrenceOrderK then
	print(`CalculateTreeRecurrenceOrderK(F,n,deg,ord,solvar,subvar,Init,N):`);
	print(`Inputs - `);
	print(`F: the recurrence as a function in variables solvar and {subvar[n-1],...,subvar[n-ord]}`);
	print(`n: a symbol used in the set of subvars`);
	print(`Deg: the degree of solvar in F`);
	print(`ord: the order of the recurrence, how many terms back in the sequence are required`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: the variable from which the set of subvars is created`);
	print(`Init: The initial K=ord levels of the recurrence tree`);
	print(`N: the number of levels beyond the initial levels you want to compute`);
	print(`Output - `);
	print(`The recurrence tree up to level N+ord.`);
	print(``);
	print(`For example try:`);
	print(`F := X^2*s[n-3]^2+s[n-1]^3*s[n-3]+s[n-1]^2*s[n-2]^2-(4*s[n-1]*s[n-2]*s[n-3]-s[n-2]^3)*X`);
	print(`CalculateTreeRecurrenceOrderK(F,n,2,3,X,s,[1,[1,1],[[1,1],[1,1]]],4)`);

elif nargs = 1 and args[1]=OneTreeIterOrderKNoReps then
	print(`OneTreeIterOrderKNoReps(F,n,deg,ord,solvar,subvar,TK,Already):`);
	print(`Inputs - `);
	print(`F: the recurrence as a function in variables solvar and {subvar[n-1],...,subvar[n-ord]}`);
	print(`n: a symbol used in the set of subvars`);
	print(`Deg: the degree of solvar in F`);
	print(`ord: the order of the recurrence, how many terms back in the sequence are required`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: the variable from which the set of subvars is created`);
	print(`TK: the previous K=ord levels of the recurrence tree`);
	print(`Already: values that have already been seen in the recurrence tree`);
	print(`Output - (2 outputs)`);
	print(`1. The next level of the recurrence tree. If there are no solutions for a given`);
	print(`set of previous values the "answer" will be FAIL. If the solution has already`);
	print(`been seen then the answer will be `);
	print(`2. The new set for Already (Already union {all values seen in the next level)`);
	print(``);
	print(`For example try: `);
	print(`F:= X^2*s[n-3]^2+s[n-1]^3*s[n-3]+s[n-1]^2*s[n-2]^2-(4*s[n-1]*s[n-2]*s[n-3]-s[n-2]^3)*X`);
	print(`OneTreeIterOrderKNoReps(F,n,2,3,X,s,[1,[1,1],[[1,1],[1,1]]],{1})`);

elif nargs = 1 and args[1]=CalculateTreeRecurrenceOrderKNoReps then
	print(`CalculateTreeRecurrenceOrderKNoReps(F,n,deg,ord,solvar,subvar,Init,N,Already):`);
	print(`Inputs - `);
	print(`F: the recurrence as a function in variables solvar and {subvar[n-1],...,subvar[n-ord]}`);
	print(`n: a symbol used in the set of subvars`);
	print(`Deg: the degree of solvar in F`);
	print(`ord: the order of the recurrence, how many terms back in the sequence are required`);
	print(`solvar: takes the place of a(n) in the recurrence F`);
	print(`subvar: the variable from which the set of subvars is created`);
	print(`Init: The initial K=ord levels of the recurrence tree`);
	print(`N: the number of levels beyond the initial levels you want to compute`);
	print(`Already: the values that have already been seen in the tree`);
	print(`Output - `);
	print(`The recurrence tree up to level N+ord where each value is only seen in one level`);
	print(``);
	print(`For example try:`);
	print(`F := X^2*s[n-3]^2+s[n-1]^3*s[n-3]+s[n-1]^2*s[n-2]^2-(4*s[n-1]*s[n-2]*s[n-3]-s[n-2]^3)*X`);
	print(`CalculateTreeRecurrenceOrderKNoReps(F,n,2,3,X,s,[1,[1,1],[[1,1],[1,1]]],4,{1})`);

elif nargs = 1 and args[1]=MakeSequences then
	print(`MakeSequences(T,d):`);
	print(`Inputs - `);
	print(`T: A recurrence tree in the form specified in the initial info for this package`);
	print(`(all procedures output recurrence trees in the correct format)`);
	print(`d: the length of the sequences that will be in the output`);
	print(`Outputs - `);
	print(`All possible sequences of length d in T followed by their "path" in the recurrence tree.`);
	print(`For example the path to the 11 in the following tree is [1,2,2]`);
	print(`[1, [2,3] , [[4,5],[6,7]] , [[[8,9],[10,11]],[[12,13],[14,15]]] ]`);
	print(`and the sequence is [1,2,5,11]`);
	print(`Any sequence with a "FAIL" or a "false" are omitted`);
	print(``);
	print(`For example try:`);
	print(`MakeSequences([1, [2,3] , [[4,5],[6,7]] , [[[8,9],[10,11]],[[12,13],[14,15]]] ],4)`);

elif nargs = 1 and args[1]=JustSequences then
	print(`JustSequences(T,d)`);
	print(`Inputs - `);
	print(`T: A recurrence tree in the form specified in the initial info for this package`);
	print(`(all procedures output recurrence trees in the correct format)`);
	print(`d: the length of the sequences to be in the output`);
	print(`Outputs - `);
	print(`All possible sequences of length d in T without specifying the path.`);
	print(`Any sequence with a "FAIL" or a "false" are omitted`);
	print(``);
	print(`For example try:`);
	print(`JustSequences([1, [2,3] , [[4,5],[6,7]] , [[[8,9],[10,11]],[[12,13],[14,15]]] ],4)`);

elif nargs=1 and args[1]=CheckDiv then
	print(`CheckDiv(S,L):`);
	print(`Inputs - `);
	print(`S: a set (or list) of lists which are sequences`);
	print(`Note any 1s in any of the sequences will be disregarded`);
	print(`L: a set of numbers to check for divisibility`);
	print(`Note any 1s in L will be removed`);
	print(`Output - `);
	print(`The set of sequences for which at least one value does not have`);
	print(`a divisor from the set L`);
	print(``);
	print(`For example try: CheckDiv([[2,4,6],[2,3,5],[3,6,9]],{2,3})`);

elif nargs=1 and args[1]=OnlyZ then
	print(`OnlyZ(S):`);
	print(`Inputs - `);
	print(`S: a set (or list) of lists which are sequences`);
	print(`Output - `);
	print(`The sequences in S which consist only of integers`); 
	print(``);
	print(`For example try: OnlyZ({[2,3,4],[1/2,4,3],[4,7,8],[3,1/4,10]})`);

else
	print(`No such option`);
fi;

end:

############################### Start Procedures ######################################

#======================================================================================
#OneTreeIterOrderK(F,n,deg,ord,solvar,subvars,TK)

OneTreeIterOrderK := proc(F,n,deg,ord,solvar,subvar,TK) option remember;
local i,TKp,A,a,sub,Next,eq,ans;

if nops(TK)<>ord then
	print(`Tree wrong size, not enough or too many levels`);
	RETURN(FAIL);
fi;

if nops(TK[1])= 1 then
	A := Addresses(ord,deg); #All possible sequences of integers from 1 to deg of length ord
	Next := [];
	for a in A do
		sub := {subvar[n-ord]=TK[1],seq(subvar[n-(ord-i)] = TK[i+1][op(1..i,a)],i=1..ord-1)};
		eq := subs(sub,F);
		if subs(sub,s[n-1])=FAIL then
			Next := [op(Next),[FAIL$deg]];
		else
			ans := [solve(eq,solvar)];
			if nops(ans)=0 then
				Next := [op(Next),[FAIL$deg]];
			else
				Next := [op(Next),ans];
			fi;
		fi;
	od;
	#print(Next);
	RETURN(Paren(Next,deg)); #Paren makes sure that the answer is in the correct format adding more brackets where necessary

else
	RETURN([seq(OneTreeIterOrderK(F,n,deg,ord,solvar,subvar,[seq(TK[k][j] ,k=1..ord)]),j=1..nops(TK[1]))]);
fi;

end:

#======================================================================================
CalculateTreeRecurrenceOrderK := proc(F,n,deg,ord,solvar,subvar,Init,N) option remember;
local Ret,LastIter,i,OneIter;

Ret := [op(Init)];
LastIter := Init;
for i from 1 to N do
	OneIter := OneTreeIterOrderK(F,n,deg,ord,solvar,subvar,LastIter);
	Ret := [op(Ret),OneIter];
	LastIter := [op(2..nops(LastIter),LastIter),OneIter];
od;

RETURN(Ret);

end:

#======================================================================================
#OneTreeIterOrderKNoReps(F,n,deg,ord,solvar,subvar,TK,Already)

OneTreeIterOrderKNoReps := proc(F,n,deg,ord,solvar,subvar,TK,Already) option remember;
local i,TKp,A,a,sub,Next,eq,ans,New,P,j,S,Saw,Now;

if nops(TK)<>ord then
	print(`Tree wrong size`);
	RETURN(FAIL);
fi;

if type(TK[1],constant) then
	New := Already;
	A := Addresses(ord,deg);
	Next := [];
	for a in A do
		sub := {subvar[n-ord]=TK[1],seq(subvar[n-(ord-i)] = TK[i+1][op(1..i,a)],i=1..ord-1)};
		eq := subs(sub,F);
		if subs(sub,s[n-1])=FAIL then
			Next := [op(Next),[FAIL$deg]];
		elif subs(sub,s[n-1])=false then
			Next := [op(Next),[false$deg]];
		else
			ans := [solve(eq,solvar)];
			if nops(ans)=0 then
				Next := [op(Next),[FAIL$deg]];
			else
				P := [];
				for j from 1 to nops(ans) do
					if ans[j] in New then
						P := [op(P),false];
					else
						New := New union {ans[j]};
						P := [op(P),ans[j]];
					fi;
				od;
				Next := [op(Next),P];
			fi;
		fi;
	od;
	RETURN(Paren(Next,deg),New);

else
	S := [];
	Saw := {};
	for i from 1 to nops(TK[1]) do
		Now := OneTreeIterOrderKNoReps(F,n,deg,ord,solvar,subvar,[seq(TK[k][i] ,k=1..ord)],Already);
		S := [op(S),Now[1]];
		Saw := Saw union Now[2];
	od;
	
	RETURN(S,Saw);

fi;

end:

#======================================================================================
CalculateTreeRecurrenceOrderKNoReps := proc(F,n,deg,ord,solvar,subvar,Init,N,Already) option remember;
local Ret,LastIter,i,OneIter,New;

Ret := [op(Init)];
LastIter := Init;
New := Already;
for i from 1 to N do
	OneIter := OneTreeIterOrderKNoReps(F,n,deg,ord,solvar,subvar,LastIter,New);
	Ret := [op(Ret),OneIter[1]];
	LastIter := [op(2..nops(LastIter),LastIter),OneIter[1]];
	New := New union OneIter[2];
od;

RETURN(Ret);

end:



#======================================================================================
#OneTreeIterOrder1(F,Deg,solvar,subvar,y): Inputs a function F in two variables, solvar and 
#subvar, an integer Deg which is the degree of solvar in F, and a list of lists y of values
#for subvar. Outputs a list of lists of total length Deg*nops(y) where output[i]=answers for y[i] 
#If there is no solution then all those slots have FAIL, if there is only one solution 
#then all the slots are identical.

OneTreeIterOrder1 := proc(F,Deg,solvar,subvar,y) option remember;
local eq,R,i,ans;

if type(op(1,y),constant) then
	eq := F;
	R := [];
	
	for i from 1 to nops(y) do
		if y[i] = FAIL then
			R := [op(R),[FAIL$Deg]];
		elif subs(subvar=y[i],eq)=0 then
			R := [op(R),[FAIL$Deg]];
		else
			ans := [solve(subs(subvar=y[i],eq)=0,solvar)];
			if ans=[solvar] then
				R := [op(R),[FAIL$Deg]];
			elif nops(ans)=0 then
				R := [op(R),[FAIL$Deg]];
			elif nops(ans)=Deg then
				R := [op(R),[op(ans)]];
			elif nops(ans)=1 then
				R := [op(R),[op(ans)$Deg]];
			fi;
		fi;	
	od;
	RETURN(R);
fi;

RETURN([seq(OneTreeIterOrder1(F,Deg,solvar,subvar,y[i],X),i=1..nops(y))]);

end:

#======================================================================================
CalculateTreeRecurrenceOrder1 := proc(Init,F,Deg,solvar,subvar,N) option remember;
local OneIter,i,Ret,LastIter;

Ret := [Init];
LastIter := [Init];
for i from 1 to N do
	OneIter := OneTreeIterOrder1(F,Deg,solvar,subvar,LastIter);
	Ret := [op(Ret),op(OneIter)];
	LastIter := OneIter;
od;

RETURN(Ret);

end:

#======================================================================================
#OneTreeIterOrder1NoReps(F,Deg,solvar,subvar,y,X,Already): Inputs a function F in two variables, solvar and 
#subvar, an integer Deg which is the degree of solvar in F, and a list of lists y of values
#for subvar. Outputs a list of lists of total length Deg*nops(y) where output[i]=answers for y[i] 
#If there is no solution then all those slots have FAIL, if there is only one solution 
#then all the slots are identical. I have not yet taken care of the case where the number 
#of solutions is <Deg but >1. 

OneTreeIterOrder1NoReps := proc(F,Deg,solvar,subvar,y,Already) option remember;
local eq,R,i,ans,Saw,j,P,S,Now;

if type(op(1,y),constant) then
	eq := F;
	R := [];
	Saw := Already;
	for i from 1 to nops(y) do
		if y[i] = FAIL then
			R := [op(R),[FAIL$Deg]];
		elif y[i] = false then
			R := [op(R),[false$Deg]];
		else
			if subs(subvar=y[i],eq)=0 then
				R := [op(R),[FAIL$Deg]];
			else
				ans := [solve(subs(subvar=y[i],eq)=0,solvar)];
				if nops(ans)=0 then
					R := [op(R),[FAIL$Deg]];
				elif nops(ans)=Deg then
					P := [];
					for j from 1 to Deg do
						if ans[j] in Already then
							P := [op(P), false];
						else
							P := [op(P), ans[j]];
							Saw := Saw union {ans[j]};
						fi;
					od;
					R := [op(R),P];
				elif nops(ans)=1 then
					if ans[1] in Already then
						R := [op(R),[false$Deg]];
					else
						R := [op(R),[op(ans)$Deg]];
						Saw := Saw union {ans[1]};
					fi;
				fi;
			fi;
		fi;	
	od;
	RETURN(R,Saw);
fi;

S := [];
Saw := {};
for i from 1 to nops(y) do
	Now := OneTreeIterOrder1NoReps(F,Deg,solvar,subvar,y[i],Already);
	S := [op(S),Now[1]];
	Saw := Saw union Now[2];

od;

RETURN(S,Saw);

end:

#======================================================================================
CalculateTreeRecurrenceOrder1NoReps := proc(Init,F,Deg,solvar,subvar,N) option remember;
local OneIter,i,Ret,LastIter,Already;

Ret := [Init];
LastIter := [Init];
Already := {Init};
for i from 1 to N do
	OneIter := OneTreeIterOrder1NoReps(F,Deg,solvar,subvar,LastIter,Already);
	Ret := [op(Ret),op(OneIter[1])];
	LastIter := OneIter[1];
	Already := Already union OneIter[2];
od;

RETURN(Ret);

end:

#======================================================================================
#MakeSequences(T,d): Inputs a recurrence tree and outputs all possible sequences of 
#length d starting at the top node

MakeSequences := proc(T,d)
local SeqsInds,Bad,Pair,Ret,s,s1,i;

if d>nops(T) then
	print(`Depth too big`);
	RETURN(FAIL);
fi;

if d=1 then
	RETURN({[[T[1]],[]]});
elif d=2 then
	RETURN({seq([[T[1],T[2][i]],[i]] ,i=1..nops(T[2]))});
fi;

Bad := {false,FAIL};

SeqsInds := MakeSequences(T,d-1);
Ret := {};
for s in SeqsInds do
	for i from 1 to nops([op(T[d][op(s[2])])]) do
		if not (T[d][op(s[2]),i] in Bad) then
			s1 := [[op(s[1]),T[d][op(s[2]),i]],[op(s[2]),i]];
			Ret := Ret union {s1};
		fi;
	od;
od;

RETURN(Ret);

end:

#======================================================================================
JustSequences := proc(T,d)
local S;

S := MakeSequences(T,d);

RETURN({seq(s[1], s in S)});

end:

#======================================================================================
#Paren(L,deg): Parenthesizes (with []) the list L of length a multiple of deg into 
#groups of length deg

Paren := proc(L,deg)
local n,div,P,i;

n := nops(L);
div := n/deg;
P := [];

for i from 1 to div do
	P := [op(P),[seq(L[j],j=(i-1)*deg+1..i*deg)]];
od;

RETURN(P);

end:

#======================================================================================
#Addresses(ord,deg): all lists of integers from 1 to deg of length ord-1

Addresses := proc(ord,deg) option remember;
local A,Ao,a;

if ord=2 then
	RETURN([seq([i],i=1..deg)]);
fi;

A := Addresses(ord-1,deg);
Ao := [];
for a in A do
	Ao := [op(Ao),seq([op(a),i],i=1..deg)];
od;

RETURN(Ao);

end:

#======================================================================================
CheckDiv := proc(S,L)
local s,i,j,F,div,E,flag;

F := {};
for s in S do
	#print(s);
	flag := true;
	for i from 1 to nops(s) do
		if not s[i]=1 then
			div := divisors(s[i]);
			E := div intersect {op(L)};
			if E={} then
				flag:=false;
			fi;
		fi;
	od;
	if flag=false then
		F := F union {s};
	fi;
od;

RETURN(F);

end:

#======================================================================================
OnlyZ := proc(S)
local s,Ret;

Ret := {};

for s in S do
	if not false in {seq(type(s[i],integer),i=1..nops(s))} then
		Ret := Ret union {s};
	fi;
od;

RETURN(Ret);

end: