(* ::Package:: *)

(* Pacakge: EulerianTrees
   Author:  Luis A. Medina
   Institution: University of Puerto Rico, Rio Piedras
   email: luis.medina17@upr.edu
*)


BeginPackage["EulerianTrees`"];


(* ::Text:: *)
(*Message*)


Print["This program was written by 

Luis A. Medina 
University of Puerto Rico, Rio Piedras
San Juan, PR

Please report bugs to: luis.medina17@upr.edu

The most current version of this  package is available 
from http://emmy.uprrp.edu/lmedina/

Below you can see the list of functions included in this
package.  For information about a particular function, 
simply click on it or type ?NameOfFunction, e.g. ?LinRec.
Type ?EulerianTrees`* if you want to see again the list of 
functions included in the package."]


(* ::Text:: *)
(*Functions visible to the user*)


EulerianTree::usage="EulerianTree[p,k,M]: Given a prime 'p' and a positive integer 'k', returns
the tree up to level 'M' for the the p-adic valuation of the Eulerian number A(n,k) 
when n runs through the positive integers bigger than k."


EulerianTreeNode::usage="EulerianTreeNode[node,p,k,M]: Given a prime 'p', a positive integer 'k', 
and a vertex 'node' in the Eulerian tree for p and k, it returns the branch of 'node' 
from its position down M levels. Note, 'node' should have the form {elm level}.";


Begin["Private`"];


(* ::Text:: *)
(*Explicit definition of the Eulerian numbers*)


A[n_,k_]:=Sum[(-1)^j Binomial[n+1,j](k-j+1)^n,{j,0,k}]


(* ::Text:: *)
(*A number that has the same p-adic valuation as A(n, k).*)


Ap[n_,k_,p_?PrimeQ]:=Module[{list},
list=Select[Table[j,{j,0,k}],Mod[k-#+1,p]!=0&];
Sum[(-1)^j Binomial[n+1,j](k-j+1)^n,{j,list}]
]


(* ::Text:: *)
(*Period of the Eulerian numbers modulo (p^m).*)


LP[p_,k_,m_]:=(p-1)p^(Floor[Log[p,k]]+m)


(* ::Text:: *)
(*Level one of the tree.*)


level[p_,k_,1]:=Thread[Naturals->Table[{i,1},{i,0,(p-1)p^(Floor[Log[p,k]]+1)-1}]]


(* ::Text:: *)
(*Level M of the tree.*)


level[p_,k_,M_]:=Module[{PreviousLevel=level[p,k,M-1],count=0,L,el,list={},nodestocheck,elementstocheck,tocheck},
nodestocheck=Quiet[Select[PreviousLevel,Last[Last[#]]==M-1&]];
elementstocheck=First/@Last/@nodestocheck;
L=Length[nodestocheck];
While[(count=count+1)<=L,
	el=elementstocheck[[count]];
	tocheck=Table[IntegerExponent[Ap[el+LP[p,k,M-1]m,k,p],p],{m,1,p+1}];
		If[Length[Union[tocheck]]>1,
			list=Append[list,el],
			PreviousLevel=PreviousLevel/.{(nodestocheck[[count]])->{nodestocheck[[count]],ToString[IntegerExponent[Ap[el+LP[p,k,M-1],k,p],p]]}}
		]
];
L=Length[list];
Join[PreviousLevel,Flatten[Table[Table[{list[[i]],M-1}->{list[[i]]+r LP[p,k,M-1],M},{r,0,p-1}],{i,1,L}],1]]
]


(* ::Text:: *)
(*Level one of the node.*)


nodelevel[node_List,p_,k_,1]:=Module[{previouslevel,currentlevel,nlevel,initialvertex,newnode},
newnode=node;
currentlevel=Last[node];
initialvertex=First[node];
Table[newnode->{initialvertex+r LP[p,k,currentlevel],currentlevel+1},{r,0,p-1}]
]


(* ::Text:: *)
(*M levels down the node.*)


nodelevel[node_List,p_,k_,M_]:=Module[{PreviousLevel,count=0,L,el,list={},nodestocheck,elementstocheck,tocheck,currentlevel},
currentlevel=Last[node];
PreviousLevel=nodelevel[node,p,k,M-1];
nodestocheck=Quiet[Select[PreviousLevel,Last[Last[#]]==currentlevel+M-1&]];
elementstocheck=First/@Last/@nodestocheck;
L=Length[nodestocheck];
While[(count=count+1)<=L,
	el=elementstocheck[[count]];
	tocheck=Table[IntegerExponent[Ap[el+LP[p,k,currentlevel+M-1]m,k,p],p],{m,1,p+1}];
	If[Length[Union[tocheck]]>1,
		list=Append[list,el],
		PreviousLevel=PreviousLevel/.{(nodestocheck[[count]])->{nodestocheck[[count]],ToString[IntegerExponent[Ap[el+LP[p,k,currentlevel+M-1],k,p],p]]}}
	]
];
L=Length[list];
Join[PreviousLevel,Flatten[Table[Table[{list[[i]],currentlevel+M-1}->{list[[i]]+r LP[p,k,M-1+currentlevel],currentlevel+M},{r,0,p-1}],{i,1,L}],1]]
]


(* ::Text:: *)
(*Eulerian p-adic Tree (guess).*)


EulerianTree[p_?PrimeQ,k_Integer,M_Integer]/;(k>=1&&M>=1):=TreePlot[level[p,k,M],Automatic,Naturals]


(* ::Text:: *)
(*Branch of the tree down a certain node.*)


EulerianTreeNode[node_List,p_?PrimeQ,k_Integer,M_Integer]/;(k>=1&&M>=1&&(Length[node]==2)&&IntegerQ[First[node]]&&IntegerQ[Last[node]]):=Module[{bool,currentlevel,el},
{el,currentlevel}=node;
bool=Table[IntegerExponent[Ap[el+LP[p,k,currentlevel]r,k,p],p],{r,1,p+1}];
If[Length[Union[bool]]>1,
	TreePlot[nodelevel[node,p,k,M],Automatic,node],
	Print["The node ", node, " is either a constant node in the tree or it belongs to a constant node in a previous level."]
]
]


?EulerianTrees`*


End[];
EndPackage[];


(* ::Input:: *)
(**)
