(* ::Package:: *)

(* ::Text:: *)
(*StirlingTrees.m is a package written by *)
(*	Luis A. Medina (lmedina@math.rutgers.edu, medinarivera@gmail.com)*)
(*	Laine Noble (noble@math.ohio-state.edu) *)
(*	*)
(*This packages accompanies the article :*)
(*"The p-adic valuation of Stirling Numbers"*)
(*	by Ana Berrizbeitia, Luis A. Medina, Alexander C. Moll, Victor H. Moll, and Laine Noble*)
(**)
(*Please save StilingTrees.m in a folder that is readable by Mathematica.  For example, *)
(*if you are a Windows user, you can save the file StirlingTrees.m in the following location:  *)
(*C:\ProgramData\Mathematica\Applications.*)
(**)
(*Please report bug to : lmedina at math dot rutgers dot edu.*)
(* *)
(*The most current version of this  package and article*)
(*are  available from http://www.math.rutgers.edu/~lmedina.*)
(**)
(*Created: 8/19/2009.*)


BeginPackage["StirlingTrees`"];


(* ::Text:: *)
(*Below you can find the functions that are public.*)


Mensaje::usage=" Created: August 19, 2009.

StirlingTrees.m was written by 
  Luis A. Medina
  Rutgers University and University of Puerto Rico
  lmedina at math dot rutgers dot edu

  Laine Noble
  Ohio State University
  noble at math dot ohio-state dot edu

This package accompanies the article:
'The p-adic valuation of Stirling Numbers'
	by Ana Berrizbeitia, Luis A. Medina, Alexander C. Moll, Victor H. Moll, and Laine Noble

Please report bug to: lmedina at math dot rutgers dot edu.

The most current version of this  package and article
are  available from http://www.math.rutgers.edu/~lmedina.

For a list of the procedures simply input ?Help to see them."


Help::usage="The functions included in this package are:
  (1) ModTree
  (2) pAdicTree
  (3) Test

If you want to know more information about these functions, 
simpliy input ?NameOfFunction.  For example, try ?pAdicTree."


ModTree::usage = "ModTree[k, p, max] returns a tree of StirlingNumbersS2[n,k] modulo powers of the prime p.\n
 In this case k fixed and n is variable.";


pAdicTree::usage = "pAdicTree[k, p, max] returns a tree of the p-adic valuation of StirlingNumbersS2[n,k].
 This procedure is heuristic. \n
 In this case k fixed and n is variable.";


Test::usage = "Test[w,p,k,j]";


Begin["Private`"];


(* ::Section:: *)
(*Auxiliary Functions*)


(* ::Text:: *)
(*S(n,k) is the Stirling number of the second kind.*)


S[n_,k_]:=Sum[(-1)^i (k-i)^n/(k-i)!/i!,{i,0,k-1}]


(* ::Text:: *)
(*L is the period length of S(n,k) mod p^m for k fixed and n varies among all integers bigger or equal to k.*)


NewFloor[p_,k_]:=If[IntegerQ[Log[p,k]],Log[p,k]-1,Floor[Log[p,k]]];
L[k_,p_,m_]:=L[k,p,m]=(p-1)*p^(NewFloor[p,k]+m-1);


\[Nu][n_Integer,p_]:=Module[{k=0,m=n},
	While[Mod[m,p]==0,(m=m/p;k++)];k];
\[Nu][0,p_]=Infinity;


(* ::Section:: *)
(*ModTrees*)


(* ::Text:: *)
(*ZeroCongruenceClassesList makes a list of congruence classes mod L_m where S\[Congruent] 0 mod p^m*)


ZeroCongruenceClassesList[k_,p_,1]:=Module[{PL=L[k,p,1]},
Flatten[
 Position[Table[Mod[S[j+PL,k],p],{j,0,PL-1}],0]-1
 ]
]


ZeroCongruenceClassesList[k_,p_,m_]:=ZeroCongruenceClassesList[k,p,m]=Block[{$RecursionLimit=Infinity},
 Module[{PL=L[k,p,m-1],NL=L[k,p,m],
 LL,Z=ZeroCongruenceClassesList[k,p,m-1],c=0,j,list={}},
 LL=Length[Z];
 While[(c=c+1)<=LL,
 j=Z[[c]];
 list=Join[list,j + Flatten[
  Position[
   Table[
    Mod[S[j+i PL+ NL, k], p^m],{i,0,p-1}],0]-1] PL]
 ];
 list
 ]
]


(* ::Text:: *)
(*EdgeList uses zcclist to compile a list of rules for branching from one generation to the next, labeled with the modular congruence of S*)


EdgeList[k_,p_,1]:=With[{LL=L[k,p,1]},
    If[p==2,
	  Table[{\[DoubleStruckCapitalN]->{1,j},If[Mod[S[j+(p+1)LL,k],p]==0,"", 0]},{j,0,LL-1}],
      Table[{\[DoubleStruckCapitalN]->{1,j},If[Mod[S[j+2 LL,k],p]==0,"", 0]},{j,0,LL-1}]
    ]
]


EdgeList[k_,p_,m_]:=Module[{c=0,j,M, Z=ZeroCongruenceClassesList[k,p,m-1],LL, list={}, PL = L[k,p,m-1], NL = L[k,p,m]},
 LL=Length[Z];
 While[(c=c+1)<=LL,
 j=Z[[c]];
 list = Join[list,
			Table[{{m-1, j}->{m,j+ i PL}, If[Mod[S[j + i PL + NL,k],p^m]==0,"",m-1]},{i,0,p-1}]]
 ];
list
]


(* ::Text:: *)
(*TreeList puts together several generations of edgelists*)


TreeList[k_,p_,max_]:=Table[EdgeList[k,p,i],{i,1,max}]


(* ::Text:: *)
(*ModTree plots max levels of the tree associated with a k,p pair.*)
(*THIS FUNCTION IS PUBLIC!!!!!*)


ModTree[k_Integer,p_?PrimeQ,max_Integer]/;(k>0 && max>0):=TreePlot[Flatten[TreeList[k,p,max],1],Top,\[DoubleStruckCapitalN],VertexLabeling->Tooltip];


(* ::Section:: *)
(*P-adic trees*)


(* ::Text:: *)
(*AllEqual checks whether the Stirling numbers of a particular congruence class Subscript[D, m,j] all have the same p-adic valuation, using a bubble sort with depth 2p.*)


AllEqual[k_,p_,m_,j_]:=Module[{soFar=True,a,b,i=0,LL},
 LL = L[k,p,m];
 While[(i=i+1)<=2p && soFar,
  a = \[Nu][S[j+i LL,k],p];
  b = \[Nu][S[j+(i+1) LL,k],p];
  If[a!=b, 
   soFar = False;
  ]
 ];
 If[soFar,
  {soFar,b},
  {soFar,-1}
 ]
] 


(* ::Text:: *)
(*This generates trees using dual criteria: first checks modular congruence of S, then checks constancy of p-adic valuations.*)


Clear[pZeroCongruenceClassesList]


pZeroCongruenceClassesList[k_,p_,1]:=
Select[ZeroCongruenceClassesList[k,p,1],Not[First[AllEqual[k,p,1,#]]]&]


pZeroCongruenceClassesList[k_,p_,m_]:=pZeroCongruenceClassesList[k,p,m]=Block[{$RecursionLimit=Infinity},
Module[{PL=L[k,p,m-1],NL=L[k,p,m],
 LL,Z=pZeroCongruenceClassesList[k,p,m-1],c=0,j,list={},i},
 LL=Length[Z];
 While[(c=c+1)<=LL,
  i=-1;
  j=Z[[c]];
  While[(i=i+1)<=p-1,
   If[Mod[S[j+ i PL + NL,k],p^m]==0 && Not[First[AllEqual[k,p,m,j+i PL]]],
    list=Join[list,{j+i PL}]
   ]
  ]
 ];
list
]
]


(* ::Text:: *)
(*Similar to EdgeList. *)


pEdgeList[k_,p_,1]:=With[{LL=L[k,p,1]},
	Table[{\[DoubleStruckCapitalN]->{1,j},If[Length[Union[Table[IntegerExponent[S[j+ i LL,k],p],{i,p,2p-1}]]]==1,IntegerExponent[S[j+(p+1)LL,k],p], ""]},{j,0,LL-1}]
]


pEdgeList[k_,p_,m_]:=Module[{c=0,j,M, Z=pZeroCongruenceClassesList[k,p,m-1],LL, list={}, PL = L[k,p,m-1], NL = L[k,p,m]},
 LL=Length[Z];
 While[(c=c+1)<=LL,
 j=Z[[c]];
 list = Join[list, 
			Table[{{m-1, j}->{m,j+ i PL}, If[Length[Union[Table[IntegerExponent[S[j + i PL + w NL,k],p],{w,p,2p-1}]]]==1,IntegerExponent[S[j + i PL + NL,k],p]]},{i,0,p-1}]]
 ];
list
]
(*M=v[S[(j2)+L[k,p,m],k],p]*)


(* ::Text:: *)
(*Similar to TreeList*)


pTreeList[k_,p_,max_]:=Table[pEdgeList[k,p,i],{i,1,max}]


(* ::Text:: *)
(*pAdicTree plots max levels of the tree associated with a k,p pair.*)
(*THIS FUNCTION IS PUBLIC!!!!!*)


pAdicTree[k_Integer,p_?PrimeQ,mmax_Integer]/;(k>0&&mmax>0):=TreePlot[Flatten[pTreeList[k,p,mmax],1],Top,\[DoubleStruckCapitalN],VertexLabeling->Tooltip]


(* ::Section:: *)
(*Test*)


F[r_,p_,w_]:=IntegerDigits[Mod[r^EulerPhi[p^(w+1)],p^(2w+2)],p,2w+2];
SubTest[w_,p_,k_,j_]:=Module[{list,\[Beta]},
list=Table[F[i,p,w],{i,1,p-1}];
\[Beta][i_,r_]:=Reverse[list[[r]][[1;;w+1]]][[i]];
Mod[Sum[(-1)^r Binomial[k,r] r^j Total[Table[\[Beta][i,r] p^(i-1),{i,1,w+1}]],{r,1,k}],p^(w+1)]!=0
];
Test[w_,p_,k_,j_]:=SubTest[w-1,p,k,j]


?Mensaje


End[];
EndPackage[];