AdditiveDecomposition.m

(* ::Package:: *)

BeginPackage["AdditiveDecomposition`"]

(* Welcome Box *)
If[$Notebooks,
   CellPrint[Cell[
     # <> "\[LongRightArrow] Type  ?AdditiveDecomposition for help",
     "Text", FontColor -> RGBColor[0, 0, 0], CellFrame -> 0.5,
     Background -> RGBColor[0.796887, 0.789075, 0.871107]]],
   Print[# <> "--> Type  ?AdditiveDecomposition for help"]]&["AdditiveDecomposition.m is written by Hao Du and Elaine Wong, Austrian Academy of Sciences, RICAM, Version 0.2 (July 8, 2020)\n"];

(* Package Usage Notes *)
AdditiveDecomposition::usage="Welcome to AdditiveDecomposition.m. This package accompanies the paper **An Additive Decomposition in Logarithmic Towers and Beyond** by Hao Du, Jing Guo, Ziming Li and Elaine Wong.

The main objective of this package is an implementation of the algorithm for decomposing a function in an S-primitive tower into its integrable part and a remainder that is minimal in some sense, 
both of which are in the same field. If the tower is also logarithmic, then we are able to determine an extension field such that the function can be decomposed further. A function in an S-primitive
tower is integrable in the tower if and only if the remainder is equal to zero.

The following commands serve the above objectives: AddDecompInField and WellGenLogTower.

Common subtasks in the above algorithm are Hermite reduction and computing the matryoshka decomposition. These are also provided in this implementation via HermiteReduceGeneral and
ProperDecomposition, respectively. Please refer to the example notebook for usage examples.

Some other functions that might be useful: ExtendedEuclidean, HeadTerm, MonomialIndicator, AddDecompLogTower, ToGenNames.
";


(* ::Subsection:: *)
(*Usage Notes Public*)


Euclidean::usage="Euclidean[a, b, var] computes the monic gcd of a and b with respect to var.";
ExtendedEuclidean::usage="ExtendedEuclidean[a, b, c, var] computes {r,s} such that r a + s b = c and the degree of r with respect to var is lower than that of b.";
DiffPoly::usage="DiffPoly[poly, var, gen, der] outputs the derivative of poly (a polynomial with respect to the last generator of a primitive tower) with respect to var.";
DiffRational::usage="DiffRational[ratfunc, var, gen, der] outputs the derivative of ratfunc (a rational function in the primitive tower represented by the generator and derivative lists gen and der, respectively) with respect to var.";
HermiteReduceGeneral::usage="HermiteReduceGeneral[f, var, gen, der] outputs {g, h, p}, such that f = g' + h + p, where h is simple and p is a polynomial with respect to the last generator in gen.";
ProperDecomposition::usage="ProperDecomposition[f, {t_1, ..., t_n}] outputs a list of functions {f_0, f_1, ..., f_n} such that f = f_0 + f_1 + ... + f_n, f_i in K(t_1, ..., t_i)[t_{i+1}, ...t_n] with t_i-proper coefficients for all i with i > 0 and f_0 in K[t_1, ..., t_n].";
AssMat::usage="AssMat[gen, der] outputs the associated matrix for the tower generated by the generator list gen characterized by their derivatives der.";
CLinearlyIndependentQ::usage="CLinearlyIndependentQ[f, var, gen, list] ouputs True if f is linearly independent with respect to the list (over a field independent of var and gen); otherwise, it outputs {False, clist}, where clist is the coefficient list of f with respect to the given list.";
CLinearBasis::usage="CLinearBasis[list, var, gen] outputs a linear basis of list (over a field independent of var and gen) and the list of coefficient lists of the functions in the list with respect to this basis.";
SigIndex::usage="SigIndex[f, gen] outputs the position of the last nozero projection with respect to the proper decomposition of f.";
GenSort::usage="GenSort[gen, der] outputs a list of new derivatives of the generators whose significant vector is larger and the permutation order of der that would give this new list.";
IndVec::usage="IndVec[posintlist, a positive integer k] outputs an ordered list of zeroes and ones indicating the elements of {1, ..., k} that are present in posintlist (marked by a one if present, and zero otherwise).";
SigComponentList::usage="SigComponentList[var, gen, der] outputs the significant component list of the primitive tower represented by the generator and derivative lists gen and der, respectively.";
GenReady::usage="GenReady[var, gen, der] takes generators from a primitive tower and outputs a list of new derivatives for generators in a well-sorted primitive tower and a transfer matrix between old and new generators.";
WellGenLogTower::usage="WellGenLogTower[var, gen, der] takes input from a primitive tower and outputs {new generators, the transfer matrix between old and new generators, new derivatives, associated matrix} of a well-generated tower.";
Leading::usage="LeadingTerm[poly, gen] outputs the leading coefficient and leading monomial of poly (a multivariable polynomial) with respect to the purely lexicographic ordering of generators.";
HeadTerm::usage="HeadTerm[f, gen] outputs the highest term among all of the leading terms of each projection of f with respect to the proper decomposition in the form of {coefficient, monomial}.";
MonomialIndicator::usage="MonomialIndicator[monomial, gen] outputs the index of the lowest generator (with respect to the generator list, gen) appearing in the monomial. It outputs the length of gen if the monomial is 1.";
AddDecompInField::usage="AddDecompInField[f, var, gen, der] takes a function in an S-primitive tower and outputs the pair {integrable part of f, minimal remainder (0, if f is integrable)} from the same tower.";
AddDecompWellGen::usage="AddDecompWellGen[f, var, gen, der, associated matrix] takes a function in a well-generated S-primitive tower and outputs the pair {integrable part of f, its minimal remainder (0, if f is integrable)}.";
RationalFunctionQ::usage="RationalFunctionQ[f, var, gen] outputs True if f is a rational function with respect to the variable and generators.";
SimpleQ::usage="RationalFunctionQ[f, var, gen, der] outputs True if f is simple with respect to the variable and generators.";
UserInputCheck::usage="UserInputCheck[f, var, gen, der] outputs True if f satisfies all input requirements for AddDecompLogTower.";
AddDecompLogTower::usage="AddDecompLogTower[f, var, gen, der] takes a function in a well-generated S-primitive tower and outputs the pair {integrable part of f, its minimal remainder (0 if f is integrable)}
 and all of the information needed to transform the old generators to the new ones in the form of a triple {new generator names, transfer matrix from old to new generators, new derivatives},";
CheckAddDecompLogTower::usage="CheckAddDecompLogTower[f, var, gen, results] takes the results from AddDecompLogTower[f, var, gen, der] and outputs True if the decomposition is correct.";
GenSubstitute::usage="GenSubstitute[var, gen, der] outputs a list of substitutions from the generators to its corresponding function of x.";
ToGenNames::usage="ToGenNames[f, var, gen] takes a function in variable var and a list of generators, renames the generators to t1,t2,etc..., and outputs the same function with the new names substituted, along with the generator and derivative lists with the new names as a list {fnew, genlist, derlist}. Set f=1 if you only want to convert the generator list with variable var.";
(* Usage Notes for Christoph's PQR *)
MyPolynomialQuotientRemainder::usage = "MyPolynomialQuotientRemainder[p, q, x] gives a list of the quotient and remainder of p and q, treated as polynomials in x.";


(* ::Subsection::Closed:: *)
(*Begin Private*)


Begin["`Private`"]


(* ::Subsection::Closed:: *)
(*Euclidean*)


(* Input: Two rational functions a and b which are polynomials in the variable var. *)
(* Output: The monic gcd of the numerators of a and b with respect to var. *)
Euclidean[a_,b_,var_Symbol]:=Module[{g},
(* Observe that we are taking the GCD of their numerators (which are polynomials). *)
g=PolynomialGCD@@({a,b}//Together//Numerator);
(* The main difference in this function compared to PolynomialGCD is that we make the gcd monic with respect to the desired variable. *)
g=Cancel[g/Coefficient[g,var,Exponent[g,var]]]
];


(* ::Subsection::Closed:: *)
(*ExtendedEuclidean*)


(* Input: Three polynomials a, b and c in var. *)
(* Output: Two polynomials r and s such that r a + s b = c and degree of r in var is lower than that of b. *)
ExtendedEuclidean[a_,b_,c_,var_Symbol]:=Module[{a1=Together[a],b1=Together[b],c1=Together[c],g,r,s,h,r1,rem,q},
{g,{r,s}}=PolynomialExtendedGCD[a1,b1,var];
{h,rem}=MyPolynomialQuotientRemainder[c1,g,var];
If[rem=!=0,Throw["Error in ExtendedEuclidean: c is not in the ideal generated by a and b."],
{r,s}=Cancel[h {r,s}];
(* In this case we take care to keep track of the main variable to give the correct output in terms of the main variable. *)
If[!(r===0)&&Exponent[b1,var]<= Exponent[r,var],
(* An ad hoc trick: we only compute the quotient and remainder of the numerators in order to save time (because we know that it is a polynomial in terms of var). *)
{q,r1}=MyPolynomialQuotientRemainder[Numerator[r],Numerator[b1],var];
(* We remember to return the denominator here! *)
s=Together[q a1 Denominator[b]/Denominator[r]+s];
{r1/Denominator[r],s},
{r,s}]]];


(* ::Subsection::Closed:: *)
(*DiffPoly & DiffRational*)


(* A general derivation operator. Generators and their derivatives must be known and the tower must be primitive! *)
(* Input: p \in C(var)(t_1, ..., t_{n-1})[t_n] (primitive tower), 
          gen = {t_1, ..., t_n},
          der = {t_1', ..., t_n'};
Output:  p';
Note that ' = d/d(var)
*)
DiffPoly[p_,var_Symbol,gen_List,der_List] := Module[{glast,dlast,q,e,i,coe,r=0},
If[Length[gen]===0,D[p, var], 
glast=gen[[-1]];
dlast=der[[-1]];
q=Collect[p,glast];
(* The max function prevents the value of e to be negative infinity and allows us to use the Do loop. *)
e=Max[Exponent[q,glast],-1];
(* This part requires the assumption that we are in a primitive tower every time we run through this function! We address this in our user input check function. *)
Do[
coe=Together[Coefficient[q,glast,i]];
(* Quotient rule; we use recursion here. *)
r=r+DiffRational[coe,var,Most[gen],Most[der]]*glast^i ,{i,0,e}]; 
r=r+Cancel[D[q,glast]*dlast]
]
]


(* This function is used in AddDecompWellGen and CheckAddDecompLogTower *)
DiffRational[f_,var_Symbol,gen_List,der_List] := Module[{f1=Together[f],num,den},
num=Numerator[f1];
den=Denominator[f1];
(* Quotient rule *)
Together[(DiffPoly[num,var,gen,der]*den-num*DiffPoly[den,var,gen,der])/den^2]
]


(* ::Subsection::Closed:: *)
(*HermiteReduceGeneral*)


(* This will decompose any function into its integrable part, gen[last]-simple part, and polynomial part. Generators and their derivatives must be known and the tower must be primitive! *)
(* Input: f in (C(var)(gen), d/d(var)), where gen={t_1,...t_n} is a list of generators and der is a list of their derivatives *)
(* Output: {g,h,p} such that f = g'+h+p, where g and h are in C(var)(gen), h is t_n-simple, and p is in C(var)(t_1,...t_{n-1})[t_n] *)
(* See Bronstein's book, page 44. *)
(* We use Mack's version here...it may be beneficial to try to implement the quadratic version to see if the speed can be improved. *)
(* Another option to improve speed: program Extended Euclidean. *)
HermiteReduceGeneral[f_,var_Symbol,gen_List,der_List] := Module[{v,f1=Together[f],den,p1,p2,a1,a2,d1,d2,e1,e2,e3,s,t,int=0},
If[Length[gen]===0,v=var,v=Last[gen]];
den=Denominator[f1];
(* f=f1=p1+(a1/den) *)
{p1,a1}=Together[MyPolynomialQuotientRemainder[Numerator[f1],den,v]];
(* d1=gcd(den,den') *)
d1=Euclidean[den,DiffPoly[den,var,gen,der],v];
(* d2=den/d1=squarefree part of den *)
d2=Together[den/d1];
(* If the last generator does not appear in d1,then the denominator of the remainder is already squarefree with respect to the last generator and we don't need to use this loop. *)
While[Exponent[d1,v] > 0, 
(* e1=gcd(d1,d1') *)
e1=Euclidean[d1,DiffPoly[d1,var,gen,der],v];
(* e2=d1/e1=squarefree part of d1 *)
e2=Together[d1/e1];
(* e3 = -d2*d1'/d1 *)
e3=Together[-d2*DiffPoly[d1,var,gen,der]/d1];
(* s e3 + t e2 = a1, which implies that 
a1/den = (s/d1)'+(t-s'd2/e2)/(d2e1) *)
{s,t}=Catch[ExtendedEuclidean[e3, e2, a1, v]];
(* We separate the integrable part and define the new a1 and d1 (because d2 is already squarefree). *)
int=int+Cancel[s/d1];
a1=Together[t-DiffPoly[s,var,gen,der]*d2/e2];
d1=e1
];
{p2,a1}=MyPolynomialQuotientRemainder[a1,d2,v];
p1=p1+p2;
If[Length[gen]===0,int = int+Integrate[p1,var];
{int,Cancel[a1/d2],0},
{int,Cancel[a1/d2],p1}
]
]


(* ::Subsection::Closed:: *)
(*ProperDecomposition*)


(* Input: gen = {t_1, ..., t_n}, a list of generators (no need to be primitive);
          f, an element in K(t_1, ..., t_n) where K is a differential field (could be "anything"); *)
(* Output: {f_0, f_1, ..., f_n} such that f_i \in K(t_1, ..., t_i)[t_{i+1}, ...t_n] with t_i-proper coefficients and f = Sum_{0\[LessEqual]i\[LessEqual]n}f_i (this decomposition is unique) *)
ProperDecomposition[f_,gen_List] := Module[{glast,num,den,q,r,p,i,j,coe,numcoe,p1},
If[Length[gen]===0,{Together[f]},
glast=Last[gen];
{num,den}={Numerator[#],Denominator[#]}&[Together[f]];
{q,r}=MyPolynomialQuotientRemainder[num,den,glast];
(* Here, we decompose the function according to the main generator and then decompose its coefficients (in the previous field) recursively (Matryoshka viewpoint). *)
If[q===0,p=ConstantArray[0,Length[gen]],
 coe=Table[ProperDecomposition[Coefficient[q,glast,j],Drop[gen,-1]],{j,0,Exponent[q,glast]}]; 
(* Old code: be careful about the index shift. We start at 0, but Mathematica starts at 1. *)
(* p=Table[Sum[coe[[j+1,i]] glast^j,{j,0,Exponent[q,glast]}],{i,Length[gen]}]; *)
(* We transpose so that we can use the dot product instead. *)
numcoe=Exponent[q,glast]+1;
p=Table[Take[Transpose[coe][[i]],numcoe]. glast^(Range[numcoe]-1),{i,Length[gen]}];
];
p=Append[p,Cancel[r/den]]
]
]


(* ::Subsection::Closed:: *)
(*AssMat*)


(* Input: generator and derivative lists *)
(* Output: The associated matrix for the tower *)
AssMat[gen_List, der_List]:=Module[{},
Transpose[ProperDecomposition[#,Drop[gen,-1]]&/@der]//MatrixForm
]


(* ::Subsection::Closed:: *)
(*CLinearlyIndependentQ*)


(* Input: f is in C(var)(gen) and row is a list of already C-linearly independent functions *)
(* Output: True if the f is C-linearly independent with respect to the elements in the row; 
		{False, coefficient list (clist)} otherwise. *)
CLinearlyIndependentQ[f_,var_Symbol,gen_List,row_List]:=Module[{list=Together[row],denlcm,clhs,clist,sol,c},
(* Extreme cases: If there is nothing in the row, then linear independence will depend on if the function that we check is zero or not.
1. Zero function against an empty row = coefficients should be zero and we return false;
2. Non-Zero function against an empty row = then this is linearly independent without the need to solve a system; returns true
*)
If[row==={},If[f===0,{False,{}},{True}],
(* If the row is not empty, then there is actually something to check. *)
(* We first compute the lcm of all the denominators in row including f. *)
(* It will only return the coefficients if the system has a solution. *)
denlcm=PolynomialLCM@@(Denominator/@Append[list,Together[f]]);
clist=Array[c,Length[list]];
(* Find all coefficients of (Sum c[j] row[j] - f) *)
clhs=CoefficientList[Together[(clist.list -f)denlcm],Prepend[gen,var]];
(* Set these coefficients = 0 and solve the system for c[j] for all j *)
If[sol=Solve[Flatten[clhs]==0,clist];sol==={},{True},clist=clist/.sol[[1]];{False,clist}] 
]
]


(* ::Subsection::Closed:: *)
(*CLinearBasis*)


(* Input: row, a list of rational functions in C(var)(gen);
         var, a variable;
         gen, a list of generators; *)
(* Output:  {basis, res} where basis is a list of C-linear basis of row; and res is a list coefficient lists, with res[[i]] begin a list of coefficients of row[[i]] w.r.t. basis
*)
CLinearBasis[row_List,var_Symbol,gen_List]:=Module[{list=Together[row],len=Length[row],basis = {},res={}, check,i},
Do[check=CLinearlyIndependentQ[list[[i]],var,gen,basis];
(* We check each element against all previous elements in the list to determine the basis for that list; if linearly independent, then we add it to the current basis list and coefficients of the form {0,...,0,1}; otherwise, we simply use the coefficients that were derived from CLinearlyIndependentQ. *)
If[check[[1]],
res=Append[res,Append[ConstantArray[0,Length[basis]],1]];basis = Append[basis, list[[i]]],
res=Append[res,check[[2]]]
],{i,len}];
(* Unify the length of each list in res, which is equal to the length of basis *)
res = Table[PadRight[res[[i]],Max[Length[basis],1]],{i,len}];
{basis,res}
]


(* ::Subsection::Closed:: *)
(*SigIndex*)


(* Determine the highest generator present in a function given ANY ordered list of generators (the position is returned). *)
(* Input: A function f in K(gen), and ordered list of generators, where K can be any field. *)
(* Output: Highest position attained from the list of generators (and 0, if f is in K). *)
SigIndex[f_,gen_List]:=Max[Position[Rest[ProperDecomposition[f,gen]],a_/;a=!=0,{1}]]


(* ::Subsection::Closed:: *)
(*GenSort*)


(* Input: A list of generators and a corresponding list of their derivatives. *)
(* Output: The list of derivatives sorted according to the given generator order (with substitutions), and the permutation corresponding to the new order of generators (with respect to the old ordering). *)
GenSort[gen_List,der_List]:=Module[{newder=Together[der],len=Length[gen],transmat,posorder,perm},
posorder=First/@SortBy[Table[{i,SigIndex[newder[[i]],gen]},{i,len}],Last];
perm = posorder;
While[posorder=!=Range[len],
newder=(newder/.Thread[gen->gen[[InversePermutation[posorder]]]])[[posorder]];
posorder=First/@SortBy[Table[{i,SigIndex[newder[[i]],gen]},{i,len}],Last];
perm = PermutationProduct[posorder,perm]
];
Return[{newder,perm}]
]


(* ::Subsection::Closed:: *)
(*IndVec*)


IndVec[list_List,len_Integer]:=If[MemberQ[list,#],1,0]&/@ Range[len]


(* ::Subsection::Closed:: *)
(*SigComponentList*)


(* Input: main variable, generator list, derivative list *)
(* Output: the significant component list of the generators (which form the "treads" of the associated matrix) *)
SigComponentList[var_,gen_List,der_List]:=Module[{len=Length[gen],num=Numerator[Together[der]],den=Denominator[Together[der]],genvar=Prepend[gen,var],i},
Table[Together[PolynomialRemainder[num[[i]],den[[i]],genvar[[SigIndex[der[[i]],genvar]]]]/den[[i]]],
{i,len}]
]


(* ::Subsection::Closed:: *)
(*GenReady*)


(* C-Linearly Independent Computation Script (all rounds except for the first round!): *)
(* Identify the positions of the basis elements in sclist, which gives us basispos. *)
(* In basis[[2]], we collect the coefficients of elements in treads with respect to the basis. *)
(* Some of them will look like elementary basis vectors such as (1,0,0,0) or (0,0,1,0,0,0). *)
(* We now determine which of them (i.e. positions) are actually basis vectors simply by matching them to the real elementary basis vectors, which comes from the identity matrix. *)
(* If there are more than one, we just take the first one because the others can be written as a C-linear combination of the first one anyways. *)
(* The output of basispos will be sorted because the inputs were sorted. *)
(* transmat collects all of the information to perform the correct column operations in order to remove (subtracting the C-linear combinations that we computed in basis[[2]])the dependent sclist elements (i.e. make them zero) and also makes sure that the columns that DON'T have dependent sclist elements is kept the same (adding IndVec matrix). *)
(* Now we substitute to get the new derivatives after this transformation to the newgen. *)
(* newder = "oldder" * transmat *)
(* newgen = "oldgen" * transmat *)
(* BUT, we want to use the SAME names of the generators, so if we want to replace the "oldgen" to "newgen", we have to multiply "newgen" with "transmat inverse" on the right. *)
(* In the substitution below, the left gen is old gen, and the right gen is newgen but we keep the same names. *)
(* Sorting will now give us a new upper triangular matrix (which look like a staircase). *)
(* We are moving all columns where the sclist element had turned to zero to the left in such a way, that the significant index list of the sclist decrease in value. We are guaranteed this because the left-most column with a certain signficant index is moved to a place that pushes the element in the sclist with the greater significant index one column over to the right, so in that position, the significant index must decrease. *)
(* Important! We are guaranteed such a decrease in every loop and at some point, the algorithm must stop because we have a minimal element (all zeros) according to this ordering. *)


(* Basis Computation Script (every round!): *)
(* We compute the proper parts of each derivative with respect to its significant generator, and x is included in that list of generators. This list of proper parts will form a collection of projections that can be grouped by their significant indices, and these will form the significant component list of the associated matrix. *)
(* Next, we compute the basis list for all of the sclist at once. *)
(* Note that certain groups of elements in sclist will always be C-linearly independent, because the proper decomposition gives us a direct sum. *)
(* basis[[1]] is a subset of sclist. *)
(* basis[[2]] is the list of coefficients for each element in the sclist (there are len of them) with respect to basis[[1]], which has lenbasis elements. Thus, basis[[2]] is a len by lenbasis matrix. *)
(* If there are C-linearly dependent elements in the sclist, then the number of basis elements will be less than the total number of generators, and this forces us to continue the loop. *)


(* Input: a list of generators and their derivatives using variables we like *)
(* Output: a list of new derivatives in which the elements of its significant component list are sorted and C-linearly independent; a transfer matrix from old to new generators; does not guarantee simple *)
GenReady[var_,gen_List,der_List]:=Module[{genvar=Prepend[gen,var],newder=der,len=Length[gen],lenbasis,sclist,transmat,transmatall,basis,basispos,sort,di,si,rem,i,firstround=True},
While[!(len===lenbasis),
If[firstround,
(* The first round sorts the original input. *)
sort = GenSort[gen,der];
newder=sort[[1]];
(* The first round only stores the sort, because that is all we have done so far. *)
transmatall=Transpose[IdentityMatrix[len][[sort[[2]]]]];
firstround=False
,(* else: all subsequent rounds... *)
(* This is the script that produces C-Linearly independent treads. *)
basispos = Position[basis[[2]],#][[1,1]]&/@IdentityMatrix[lenbasis];
transmat = Transpose[Together[IdentityMatrix[len]-basis[[2]].IdentityMatrix[len][[basispos]]+
DiagonalMatrix[IndVec[basispos,len]]]];
newder=Together[(newder.transmat)/.Table[gen[[i]]->(gen.Inverse[transmat])[[i]],{i,len}]];
(* Again, we apply the "sort" and "basis computation" script. *)
sort=GenSort[gen,newder];
newder = sort[[1]];
(* now we store the column operations and then the sort operation, in that order *)
transmatall=transmatall.transmat.Transpose[IdentityMatrix[len][[sort[[2]]]]]; 
];
(* This is the "basis computation" script that must be done for all rounds. *)
sclist=SigComponentList[var,gen,newder];
basis = CLinearBasis[sclist,var,gen];
lenbasis = Length[basis[[1]]]
];
{newder,transmatall}
]


(* ::Subsection::Closed:: *)
(*WellGenLogTower*)


(*
Input: der = {t_1', ..., t_n'}, a list of derivatives and der[i] is in C(var)(gen[1,2,..i-1]);
       var, a variable and the derivation is d/d(var);
       gen = {t_1, ..., t_n}, a list of logarithmic generators names; 
Output: newgen, the name of new generators;
        transmat, the transfer matrix from gen to newgen such that 
                    gen = newgen * transmat;
        newder, the derivatives of newgen;
        wmat, a floating staircase matrix such that each row is C-linear independent gives rise to a well-generated log tower
*)
WellGenLogTower[var_Symbol,gen_List,der_List] := Module[{n=Length[der],genready,mat,data,len, transmat,m,newgen,wmat,basis,i,j,rownum,colnum,newder},
If[Length[gen]===0,{{},{{}},{},{{}}},
genready=GenReady[var,gen,der];
(* Original n by n matrix (pi_i(t_j^\prime)) with 0 \[LessEqual] i \[LessEqual] n-1 and 1 \[LessEqual] j \[LessEqual] n. *)
mat=Transpose[ProperDecomposition[#,Drop[gen,-1]]&/@genready[[1]]];
(* data[[i,1]] is the basis of i-th row of mat and data[[i,2]] is the C-linear relationship. 
Moreover, delete the zero rows because we don't want to create a new constant. *)
data =  DeleteCases[Table[CLinearBasis[mat[[i]],var,Take[gen,i-1]],{i,n}],{{},_}];
len=Length[data];
(* Define the new generators of the well generated tower and compute the transfer matrix between gen and newgen by collecting the j-th element of all rows from the second part of data to be the j-th column of the transfer matrix. *)
transmat= Transpose[Table[Flatten[Table[data[[i,2,j]],{i,len}]],{j,n}]];
m = Length[transmat];
(* Ascribe a new name to the generators. *)
newgen = Table[Symbol["u"<>ToString[j]],{j,m}];
(* wmat is an m by m matrix associated to the well generated tower, where the derivatives of generators are the elements in data[[1,i,1]](basis of each row) for all i. This is what we will consider to be our "floating" matrix in the end. *)
(* We initialize with zeros. *)
wmat=Table[0,{i,m},{j,m}];
(* We now convert the generators to the new basis. *)
basis=Table[data[[i,1]]/.Table[gen[[j]]-> (newgen.transmat)[[j]],{j,n}],{i,len}];
(* Note that there should be no zeros in basis, so we can just flatten. *)
newder = Flatten[basis];
colnum=1;
 (* Next, we construct our floating matrix from the information in basis by making sure the i-th row has only u_{i-1}-simple elements and that every column only has one non-zero entry. *)
Do[
(* The row number is the significant index of the functions in the i-th part of basis. *)
rownum=SigIndex[basis[[i,1]],newgen]+1;
(* Replace part of i-th row of wmat by changing zero into the exact elements of the i-th part. *)
wmat[[rownum,colnum;;colnum+Length[basis[[i]]]-1]]=basis[[i]];
(* The next column number is determined by the total length of the previous rows up to i. *)
colnum=colnum+Length[basis[[i]]],
{i,1,len}
];
{newgen,transmat.Inverse[genready[[2]]],newder,wmat}]
]


(* ::Subsection::Closed:: *)
(*Leading & HeadTerm*)


(* Input: A polynomial in K[gen] where gen is a list of generators. *)
(* Output: Leading term with respect to the plex ordering of the generators in the form of {coefficient, monomial}. *)
(* We offer two versions of the code. *) 
(* Version 1: Determining the leading term via functional programming. *)
Leading[f_,gen_List]:=
If[f===0,{0,0},Fold[With[{e=Exponent[#1[[1]],#2]},{Coefficient[#1[[1]],#2,e],#1[[2]]*#2^e}]&,{f,1},Reverse[gen]]]
(* Version 2: Determining the leading term in a more straightforward way. *)
LeadingAlt[f_,gen_List]:=Module[{a=Expand[f],e,M=1,i},
If[a===0,{0,0},
For[i=Length[gen],i>0,i--,e=Exponent[a,gen[[i]]];
M=M*gen[[i]]^e;
a=Coefficient[a,gen[[i]],e]];
{a,M}]
]


(* Input: A (rational) function in K(gen) where K is any field, and a list of generators. *)
(* Output: The highest term among all of the leading terms of each projection in the form {coefficient, monomial}. *)
HeadTerm[f_,gen_List]:=Module[{pd,len=Length[gen],sum=0,c=0,M,i,h},
pd=ProperDecomposition[f,gen];
Do[
h[i]=Leading[pd[[i]],gen[[i;;len]]];
sum=sum+h[i][[2]],
{i,1,len+1}
];
M=Leading[sum,gen][[2]];
Do[
Which[h[i][[2]]===M,c=c+h[i][[1]]], (*CK: Which \[Rule] If *)
{i,1,len+1}
];
{c,M}
]


(* ::Subsection::Closed:: *)
(*MonomialIndicator*)


(* Input: a monomial and a list of generator list (in any order) *)
(* Output: an integer that identifies the position of the lowest generator in the monomial based on the order *)
(* Exceptions: If generator list is empty, it returns 1. If constant with respect to generator list, then returns length of list, which is equivalent to saying the lowest generator is the last one. *)
MonomialIndicator[m_, gen_List] := Module[{s=1},
If[gen==={},1,
If[MatchQ[Exponent[m,gen],{(0)..}],
Length[gen], 
While[Exponent[m, gen[[s]]] == 0, s = s+1];s
]
]
]


(* ::Subsection::Closed:: *)
(*AddDecompInField*)


(* Input: Function in a logarithmic tower, main variable, list of generators, and the derivatives of generators *)
(* Output: Integrable part and the remainder (minimal) in the same field. *)
(* Compare first part with WellGenLogTower *)
(* See detailed documentation in the function AddDecompWellGen *)

AddDecompInField[f_, var_Symbol, gen_List, der_List] :=
 Module[{len = Length[gen], hc, hm, s, e, pd, hr, g = 0, hp = 0, pos,
   check, cs = 0, i, j, int, r, low, res},
  If[Together[f] === 0, {0, 0},
   If[len === 0,
    hr = HermiteReduceGeneral[f, var, gen, der]; {hr[[1]], hr[[2]]},
    (* else if *)
    {hc, hm} = HeadTerm[f, gen];
    s = MonomialIndicator[hm, gen];
    e = Exponent[hm, gen[[s]]];
    pd = ProperDecomposition[hc, gen];
    Do[
     hr=HermiteReduceGeneral[pd[[i + 1]], var, gen[[1 ;; i]],der[[1 ;; i]]];
     (* Throw the integrable part together. *)
     check=CLinearlyIndependentQ[hr[[2]],var,gen[[1;;s]], der[[1;;s]]];
     If[check[[1]]||(check[[2]]==={}), g = g + hr[[1]]; hp=hp+hr[[2]],
     (* else if *)
     cs=cs+check[[2, -1]];
     g=g+hr[[1]]+Sum[check[[2,j]]*gen[[j]], {j,1,s-1}]
     ], {i, 0, s}
    ];
    low = Together[f-hc*hm-g*DiffRational[hm, var, gen, der]-(cs/(e + 1))*gen[[s]]^(e+1)*DiffRational[Cancel[hm/gen[[s]]^e], var, gen, der]];
    res = AddDecompInField[low, var, gen, der];
    int = g*hm + (cs/(e + 1))*gen[[s]]*hm + res[[1]];
    r = hp*hm+res[[2]];
    {int, Together[r]}
    ]]]


(* ::Subsection::Closed:: *)
(*AddDecompWellGen*)


(* Input: Function in a well-generated S-primitive tower, main variable, list of generators, the derivatives of generators, and the associated matrix. *)
(* Output: Integrable part and the remainder (minimal) in a well-generated tower. *)
AddDecompWellGen[f_,var_Symbol,gen_List, der_List,assmat_List]:=Module[{len = Length[gen],hc,hm,s,e,pd,hr,g=0,hp=0,pos,check,cs=0,i,j,int,r,low,res},
If[Together[f]===0,{0,0},
If[len===0,hr=HermiteReduceGeneral[f,var,gen,der];{hr[[1]],hr[[2]]},
(* leading coefficient and leading monomial of f according to our NEW ordering *)
{hc,hm}= HeadTerm[f,gen];
(* The following is information extracted from our highest head monomial. *)
(* In particular, we take out the s-th generator and its power, hm = gen_s^e * N. *)
s = MonomialIndicator[hm,gen];
e = Exponent[hm,gen[[s]]];
(* From this point on, we extract information about our highest head coefficient. *)
(* First, we decompose our highest head coefficient. *)
pd = ProperDecomposition[hc,gen];
Do[
(* Here, the i-th projection of hc is hr[[1]]'+hr[[2]], because pd[[i+1]] is gen[[i]]-proper *)
(* Which means there is no polynomial part. *)
hr=HermiteReduceGeneral[pd[[i+1]],var,gen[[1;;i]],der[[1;;i]]];
(* Throw the integrable part together. *)
g=g+hr[[1]];
(* We are checking to see if CLinearlyIndependentQ is necessary. If we are already at the end of the pd list (or the hermitian part is already 0), then we don't need to check any C-Linear independence; we just add the hermitian part. *)
If[i===len||hr[[2]]===0,hp=hp+hr[[2]],
(* Identify all non-zero positions of the gen_i row. *)
pos=Complement[Range[s],Flatten[Position[assmat[[i+1,1;;s]],0]]];
(* If there is a zero row, then we move on. *)
(* Otherwise, we check the C-linear independence of hr[[2]] with respect to the non-zero elements of the gen-i row that occur before the indicator column. *)
(* If it IS linearly independent, then we just throw all the hp together. *)
If[pos==={},hp=hp+hr[[2]],
check=CLinearlyIndependentQ[hr[[2]],var,gen[[1;;i]],assmat[[i+1,pos[[1]];;pos[[-1]]]]];
If[check[[1]],hp=hp+hr[[2]],
(* If it is NOT linearly independent, then we should add the C-linear combination of those generators using the coefficients from CLinearlyIndependentQ up to the s-1 position to the integrable part g. *)
(* We also store the coefficients of the s position to cs. *)
If[pos[[-1]]===s,
g=g+Sum[check[[2]][[j-pos[[1]]+1]]gen[[j]],{j,pos[[1]],s-1}];
cs=cs+check[[2]][[s-pos[[1]]+1]],
g=g+Sum[check[[2]][[j-pos[[1]]+1]]gen[[j]],{j,pos[[1]],pos[[-1]]}]
]]]],
{i,0, len}
];
(* After the for loop, we have the following information: *)
(* hc = g' + hp + cs * gen_s' *)
(* From here, we have three different kinds of terms of hc * hm (our highest term): *)
(* Term 1 comes from g'*hm = (g*hm)'-g*(hm)' (the second is lower order) *)
(* Term 2 comes from hp*hm *) 
(* Term 3 comes from (cs*gen_s')*hm=(cs/(e+1)*gen_s*hm)'- cs*gen_s^(e+1)/(e+1)*(N)' (the second is lower order) *)
(* We view f = hc * hm + (f - hc * hm) *)
(* Now we substitute the first hc * hm, and regroup the integrable part together and the lower order terms together *)
(* Integrable parts: g * hm + cs/e+1 * gen_s * hm + integrable parts from the recursion *)
(* Lower Order Terms: f-hc*hm-g*(hm)'- cs*gen_s^(e+1)/(e+1)*(N)'; Note that we do the recursion on these terms!!! *)
low = Together[f-hc*hm-g*DiffRational[hm,var,gen,der]-(cs/(e+1))*gen[[s]]^(e+1)* DiffRational[Cancel[hm/gen[[s]]^e],var,gen,der]];
res = AddDecompWellGen[low,var,gen, der,assmat];
int = g*hm + (cs/(e+1))*gen[[s]]*hm + res[[1]];
(* We keep all the Term 2 from the recursions as part of the remainder *)
r = hp*hm + res[[2]];
{int,Together[r]}
]]]


(* ::Subsection::Closed:: *)
(*UserInputCheck*)


RationalFunctionQ[f_,var_Symbol,gen_List]:=Module[{g=Together[f],genvar=Prepend[gen,var]},
PolynomialQ[Numerator[g],genvar]&&PolynomialQ[Denominator[g],genvar]
]


(* Input: a function to test, main variable (usually x), a generator list and their corresponding derivatives. *)
(* Output: True if the function is simple with respect to the generators; False otherwise *)
(* Simple Definition: i-th projection is t_i-simple for all i *)
SimpleQ[f_,var_Symbol,gen_List,der_List]:=Module[{pd,len=Length[gen],i=1,p,den,v,check=True},
pd=ProperDecomposition[f,Prepend[gen,var]];
(* Since we are including the var in the proper decomposition, if the first projection is non-zero, then the x-projection of f is not x-proper. *)
If[pd[[1]]=!=0,False,
(* Here, we "start" at the second decomposition, which is the x-projection of f. *)
Do[
p=Together[pd[[i+1]]];
den=Denominator[p];
If[i===1,v=var,v=gen[[i-1]]];
Which[Not[Intersection[Variables[p],gen[[i;;len]]]==={}&&Euclidean[den,DiffPoly[den,var,gen[[1;;i-1]],der[[1;;i-1]]],v]===1],check=False;Break[]]
,{i,len+1}
];
check
]
]


(* Sanitizing Inputs *)
(* We check just enough conditions on the derivatives for which our additive decomposition function can work. *)
(* The simplified conditions include:
1. t_i' belongs to the previous field, i.e. C(var)(t_1,...,t_{i-1}) (primitive tower condition);
  2. t_i' is simple (i.e. the j-th projection of t_i' is t_j-simple for all j=0,...,i-1);
  3. der_List are C-linearly independent (using CLinearlyIndependentQ)
*)
(* Input: function and var with generator and derivative list;
Output: Some kind of error OR True
*)
(* Another way to specify inputs: gen:{(_Symbol)...} *)
Clear[UserInputCheck];
UserInputCheck[f_,var_Symbol,gen_List,der_List]:=Module[{len=Length[gen]},
(* First we check the function itself, and the generator and derivative lists. *)
Which[
Not[gen==={}||Union[Head/@gen]==={Symbol}],
Throw["Error: The generators are not all variables in the form of symbols."],
Not[RationalFunctionQ[f,var,gen]],Throw["Error: f is not a rational function of the variable and generators."],
Not[len===Length[der]],Throw["Error: The number of generators does not match the number of derivatives."]
];
(* Next, we check the derivatives to see it satisfies the minimum tower requirements for our function to produce a result. *)
Do[Which[
Not[RationalFunctionQ[der[[i]],var,gen[[1;;i-1]]]&&Intersection[gen[[i;;len]],Variables[der[[i]]]]==={}],
Throw["Error: The input is not a primitive tower."],
Not[SimpleQ[der[[i]],var,gen[[1;;i-1]],der[[1;;i-1]]]],
Throw["Error: One of the derivatives is not simple."],
Not[CLinearlyIndependentQ[der[[i]],var,gen[[1;;i-1]],der[[1;;i-1]]][[1]]],
Throw["Error: The derivatives are C-linearly dependent."]],
{i,len}
];
True
]


(* ::Subsection::Closed:: *)
(*AddDecompLogTower*)


(* Input: Function in a well-sorted S-primitive tower OR a logarithmic tower, main variable, list of generators, the derivatives of generators. *)
(* Output: 
   1. integrable part and the remainder (minimal) in a well-generated tower(both in terms of the new generators);
	2. new generator list; transfer matrix; derivatives of the new generators
*)
AddDecompLogTower[f_,var_Symbol,gen_List,der_List]:=Module[{len=Length[gen],hr,genready,newgen,transmat,newder,assmat,newf=f,adwg},
If[len===0,  hr=HermiteReduceGeneral[f,var,gen,der];{{hr[[1]],hr[[2]]},{{},{{}},{}}},
Catch[
UserInputCheck[f,var,gen,der];
{newgen,transmat,newder,assmat}=WellGenLogTower[var,gen,der];
newf = newf/.Thread[gen->(newgen.transmat)];
adwg=AddDecompWellGen[newf,var,newgen,newder,assmat];
{adwg,{newgen,transmat,newder}} 
]
]
]


(* ::Subsubsection:: *)
(*CheckAddDecompLogTower*)


(* Check the output of AddDecompLogTower; Only use this function if there is no error! *)
CheckAddDecompLogTower[f_,var_Symbol,gen_List,results_List]:=Module[{check,ad, newgen},
ad=results[[1]];
newgen=results[[2]];
check=Together[DiffRational[ad[[1]],var,newgen[[1]],newgen[[3]]]+ad[[2]]-f/.Thread[gen -> newgen[[1]].newgen[[2]]]];
Return[check===0];
]


(* ::Subsection::Closed:: *)
(*GenSubstitute*)


(* Input: Main variable, generator list, and derivative list *)
(* Output: A list of substitutions from the generator name to its corresponding function of x. *)
GenSubstitute[var_,gen_,der_]:=Module[{len=Length[gen],derinitial=der,dersub={}},
If[len===0,{},
Do[
dersub=Append[dersub,Integrate[derinitial[[i]],var]];
derinitial[[i+1;;len]]=derinitial[[i+1;;len]]//.{gen[[i]]->dersub[[i]]},
{i,len-1}
];
dersub=Append[dersub,Integrate[derinitial[[len]],var]];
Thread[gen->dersub]
]
]


(* ::Subsection::Closed:: *)
(*ToGenNames*)


(* This helps the user to automate converting their function to be used with AddDecompInField. *)
(* Input: A function, the main variable, and its corresponding list of proposed generators (user choice). *)
(* Output: The same function with the newly named generators, the generator and derivative lists with the new names. *)
ToGenNames[fold_,var_Symbol,genloglist_List]:=Module[{derloglist,gentlist,dertlist,subst,fnewt},
derloglist=D[genloglist,var]//Together;
gentlist=Table[Symbol["t"<>ToString[i]],{i,Length[genloglist]}];
subst=Thread[genloglist->gentlist];
fnewt=fold//.subst;
dertlist=derloglist//.subst;
{fnewt,gentlist,dertlist}
]


(* ::Subsection::Closed:: *)
(*Christoph's PQR*)


(* This section represents all of Christoph's PQR, which we use in our code. *)

(* NormalizeVector takes a list of rational functions and removes the content.
 * i.e., the output is a list of polynomials whose common gcd is 1 which is
 * a rational function multiple of the input. *)
 
Options[NormalizeVector] = {Modulus -> 0,
                            Extension -> None,
                            Extended -> False,
                            Denominator -> True,
                            CoefficientNormal -> Expand};
                            
NormalizeVector[vec1_List, opts:((_Rule|_RuleDelayed)...)] :=
Module[{vec = vec1, test, len, t, lcm, gcd, ftl, modulus, ext, extended, denominator,
        coeffNormal, factor, norm, myTogether},
Catch[
  If[(test = Complement[First /@ {opts}, First /@ Options[NormalizeVector]]) =!= {},
    Message[NormalizeVector::optx, test]; Throw[$Failed]];
  {modulus, ext, extended, denominator, coeffNormal} = {Modulus, Extension, Extended, Denominator, CoefficientNormal}
    /. {opts} /. Options[NormalizeVector];

  If[MatchQ[vec, {(0)...}], Return[If[extended === True, {1, vec}, vec]]];
  Switch[(len = Length[vec]),
    0, Return[If[extended === True, {1, {}}, {}]],
    1, Return[If[extended === True, {First[vec], {1}}, {1}]]];

  myTogether = If[MatchQ[ext, AlgebraicNumber[_,_]], MyTogether, Together];
  norm = False; factor = 1;

  (* If Denominator -> True then we multiply by the common denominator. *)
  If[denominator === True,
    norm = True;
    vec = myTogether[vec, Modulus -> modulus, Extension -> ext];
    lcm = Function[v, Fold[PolynomialLCM[#1, #2, Modulus -> modulus, Extension -> ext]&, First[v], Rest[v]]][
            SortBy[Denominator[vec], ByteCount]];
    If[lcm =!= 1,
      factor = 1 / lcm;
      vec = (Numerator[#] * myTogether[lcm / Denominator[#], Modulus -> modulus, Extension -> ext])& /@ vec;
    ];
  ];

  If[ext === None, (* if no extension, we can use FactorTermsList *)
    If[$NormalizeMethod === "ftl",
      If[Length[vec] < 10 ||
         Function[v, Fold[PolynomialGCD[#1, #2, Modulus -> modulus]&, First[v], Rest[v]]][SortBy[vec, ByteCount]] =!= 1,

        (* a special case in which FactorTermsList does not produce the desired result *)
        If[MatchQ[vec, {_, (0)..}],
          norm = True;
          factor *= First[vec];
          vec = {1};
        ,
          ftl = FactorTermsList[vec . t^Range[0, Length[vec]-1], t, Modulus -> modulus];
          If[Not[MatchQ[Most[ftl], {(1)..}]],
            norm = True;
            factor *= Times @@ Most[ftl];
            vec = CoefficientList[Last[ftl], t];
          ];
        ];
      ];
    ,
      (* alternatively, use PolynomialGCD and Together:*)
      gcd = Function[v, Fold[PolynomialGCD[#1, #2, Modulus -> modulus]&, First[v], Rest[v]]][SortBy[vec, ByteCount]];
      If[gcd =!= 1,
        norm = True;
        factor *= gcd;
        vec = Together[# / gcd, Modulus -> modulus]& /@ vec;
      ];
    ];

  , (* Compute the polynomial gcd in some algebraic extension. *)
    gcd = Function[v, Fold[PolynomialGCD[#1, #2, Modulus -> modulus, Extension -> ext]&, First[v], Rest[v]]][
            SortBy[vec, ByteCount]];
    If[gcd =!= 1,
      norm = True;
      factor *= gcd;
      vec = myTogether[# / gcd, Modulus -> modulus, Extension -> ext]& /@ vec;
    ];

    (* Unfortunately, PolynomialGCD does not give common integer factors when using an algebraic extension.
       Hence we perform a second step in order to cancel out these. *)
    gcd = Function[v, Fold[PolynomialGCD[#1, #2, Modulus -> modulus]&, First[v], Rest[v]]][SortBy[vec, ByteCount]];
    If[gcd =!= 1,
      norm = True;
      factor *= gcd;
      vec /= gcd;
    ];
  ];

  If[norm, factor = coeffNormal[factor]; vec = coeffNormal /@ vec];
  vec = PadRight[vec, len];
  Return[If[extended === True, {factor, vec}, vec]];
]];


(* MyPolynomialDivision takes two polynomials in K[a1, ..., ak, x] and computes
 *   r\in K[a1, ..., ak, x],  c\in K(a1, ..., ak),  and q\in K(a1, ..., ak)[x]
 *   such that c*p1 = q*p2 + r.
 * The polynomials must be given as CoefficientLists w.r.t. x. *)
 
Options[MyPolynomialDivision] = {Modulus -> 0};

MyPolynomialDivision[p1_List, p2_List, opts:((_Rule|_RuleDelayed)...)] :=
Module[{r, q, c, len1, len2, f1, f2, gcd, mod},

  mod = Modulus /. {opts} /. Options[MyPolynomialDivision];

  r = p1; q = {}; c = 1;
  {len1, len2} = Length /@ {r, p2};
  While[len1 >= len2,
    {f1, f2} = NormalizeVector[{Last[r], Last[p2]}, Modulus -> mod];
    r = Together[f2 * Most[r] - Join[Table[0, {len1-len2}], f1 * Most[p2]], Modulus -> mod];
    len1--;
    q = Prepend[f2 * q, f1];
    c *= f2;

    (* If the leading coefficient turns out to be zero. *)
    While[len1 >= len2 && Last[r] === 0, r = Most[r]; len1--; PrependTo[q, 0]];

    If[MatchQ[r, {(0)..}], Return[{c, q, {}}]];

    (* Remove content *)
    {gcd, r} = NormalizeVector[r, Denominator -> False, Modulus -> mod, Extended -> True];
    If[gcd =!= 1, {q, c} = Together[{q, c} / gcd, Modulus -> mod]];
  ];

  While[len1 > 0 && Last[r] === 0, r = Most[r]; len1--];

  Return[{c, q, r}];
];

(* MyPolynomialQuotientRemainder performs polynomial division with remainder.
 * It is optimized for input with parameters, i.e., p1, p2 in K(a1, ..., ak)[x].
 * However, if the leading coefficient is a constant then MMA is faster. *)

Options[MyPolynomialQuotientRemainder] = {Modulus -> 0};

MyPolynomialQuotientRemainder[p1_, p2_, x_, opts:((_Rule|_RuleDelayed)...)] :=
Module[{cl2 = CoefficientList[p2, x], c, q, r, mod},

  mod = Modulus /. {opts} /. Options[MyPolynomialQuotientRemainder];
  (* In MMA lower than Version 6.0, PolynomialQuotientRemainder was not yet available. *)
  If[Variables[Last[cl2]] =!= {} || $VersionNumber < 6,
    {c, q, r} = MyPolynomialDivision[Together[CoefficientList[p1, x]], Together[cl2], Modulus -> mod];
    Return[{Together[(q / c), Modulus -> mod] . Table[x^(i-1), {i, Length[q]}],
            Together[(r / c), Modulus -> mod] . Table[x^(i-1), {i, Length[r]}]}];
  ,
    Return[Collect[PolynomialQuotientRemainder[p1, p2, x, Modulus -> mod], x, Together[#, Modulus -> mod]&]];
  ];
];


(* ::Subsection::Closed:: *)
(*End Credits*)


End[]


EndPackage[]