(*************************************************************************)
(*************************************************************************)
(**** : JFractions.m:                                 ********************)
(**** : Author:  Maxie D. Schmidt (maxieds@gmail.com) ********************)
(**** : Created: 09.19.2018                           ********************) 
(*************************************************************************)
(*************************************************************************)

(**** : Clear out any old definitions and stored data on reload: ****) 
ClearAll["JFractions`*"];
BeginPackage["JFractions`"];

(***************************************************************************)
(**** : Package configuration for the core functionality offered here : ****) 
(***************************************************************************)

ReturnedParameters::usage = "ABSequenceIndex := 1 || CSequenceIndex := 2 || " <> 
     "LambdaParameterIndex := 3 || MhParameterIndex := 4";
DisplaySequenceParameters::usage = "Display the parameters computed from a call to ComputeSequenceParameters[...] " <> 
     "in a nicely formatted notebook box summarizing these results.";
     
PhzLocal::usage = "Internal package function.";
QhzLocal::usage = "Internal package function.";
ConvLocal::usage = "Internal package function.";
LambdaLocal::usage = "Internal package function.";
MhLocal::usage = "Internal package function.";

ComputeSequenceParameters::usage = "ComputeSequenceParameters[seqFunc, maxN] || ComputeSequenceParametersByList[seqList]\n" <> 
     "Explicitly compute the sequence parameters for the J-fraction generating the input sequence terms.";
VerifyParameters::usage = "VerifyParameters[seqTerms] || VerifyParameters[seqFunc, maxN]: \n" <> 
     "Computationally obtain a certificate that the approximate generating function involving the previously " <> 
     "computed parameters generated by a call to ComputeSequenceParameters[...] does in fact accurately " <> 
     "generate the input sequence. This method is not needed to establish the validity of this method, but " <> 
     "rather is a testing tool used to develop and spot logic errors introduced when developing the package source.";
HthConvergentPhz::usage = "HthConvergentPhz[seqFunc, h, z]: \n " <> 
     "Generate the hth element of the numerator sequence of rational J-fraction approximations which generate the input " <> 
     "sequence. Observe that the computed sequence is polynomial, of degree at most (or typcially) h in z.";
HthConvergentQhz::usage = "HthConvergentQhz[seqFunc, h, z]: \n " <> 
     "Generate the hth element of the denominator sequence of rational J-fraction approximations which generate the input " <> 
     "sequence. Observe that the computed sequence is polynomial, of degree at most (or typcially) h in z. " <> 
     "These denominator polynomials satisfy the same second-order difference equation as the numerator sequence, but " <> 
     "predictably, and inexplicably, tend to be significantly easier to recognize in closed-form than the first case.";
HthConvergent::usage = "HthConvergent[seqFunc, h, z]: \n " <> 
     "Computes the hth rational convergent to the infinite J-fraction that formally generates the input sequence. " <> 
     "Note that wrapping this call by Series[..., {z, 0, 2 * h - 1}] allows one to compute the coefficients we have " <> 
     "input to the function. Also, for small h, say h := 2, it is possible to use Mathematica to solve for " <> 
     "exact algebraic formulas, often providing pairity in many cases, for the associated congruences satisfied by the " <> 
     "input sequence modulo each M_h, or divisors of this sequence-specific parameter. For example, " <> 
     "modulo M_2 (which is espcially useful if 2 | M_2), the input sequence is congruent to the precise formulaation " <> 
     "given by SeriesCoefficient[HthConvergent[seqFunc, 2, z], {z, 0, n}, Assumptions->None].";
EnumerateCongruenceSeries::usage = "EnumerateCongruence[seqFunc, h]:\n" <> 
     "Apply's Flajolet's observation that the approximate coefficients generated by the hth convergent functions " <> 
     "yield exact new congruences for the input sequence modulo M_h, which may icidentally be polynomial in some " <> 
     "other parameters which define the sequence. It is typically the case that strictly numerical input integer " <> 
     "sequences will lead to congruences which only hold when a rational valued M_h assumes integer values. " <> 
     "This limitation is undesirable and is a limitation of our methods. Nonetheless, in some cases the resulting " <> 
     "moduli for which these congruences hold *are* integer-valued, e.g., for generalized factorial product sequences, " <> 
     "or are nearly integer enough in form so that we can still obtain combinatorially useful new results from " <> 
     "the J-fraction methods which hold for wide subsets of the natural numbers. Properties of the resulting " <> 
     "congruences obtained here which are indeed accurate moduli some rational bases are intentionally unexplored. " <> 
     "Related congruences of the GF-of-GFs construction outlined in the examples from the " <> 
     "end user documentation can be obtained using PolynomialMod[gfCf, q] in the secondary series variable. " <>  
     "Wrapping calls to this method by Series[..., {z, 0, 2 * h - 1}], or by " <> 
     "Collect[Series[..., {z, 0, 2 * h - 1}], {z, z^_}, PolynomialMod[#, SomeFuncOrInteger]&]; may suggest " <> 
     "further intuition to the congruence properties of the convergent function sequences over h. "; 
     
OEISDatabaseLoad::usage = "Parses the internal format structure of a HTML page fetched live from the " <> 
     "OEIS database online listings. This function should be combined with a call to the other " <> 
     "OEIS*[...] package functions defined below.";
OEISSequenceLoad::usage = "Returns a list of integers which correspond to the input AXXXXXX sequence. " <> 
     "Use this method to obtain live sequence data for parsing if you are unsure of what the other OEIS " <> 
     "functions provided here are doing or how they work.";
OEISPrintf::usage = "A pretty printing printf-like function which respects the internal format " <> 
     "specifiers outlined at https://oeis.org/eishelp1.html. Example: " <> 
     "OEISPrintf[\"#%I (%N) \[RightTeeArrow] {%S,%T,%U}\", OEISDatabaseLoad[\"A290529\"].";

ABSequenceIndex=1;
CSequenceIndex=2;
LambdaParameterIndex=3;
MhParameterIndex=4;

ReturnedParameters[] := {
     {ToExpression["ab_n", TeXForm], ABSequenceIndex}, 
     {ToExpression["c_n", TeXForm], CSequenceIndex}, 
     {ToExpression["\\lambda_h", TeXForm], LambdaParameterIndex}, 
     {ToExpression["M_h", TeXForm], MhParameterIndex}
};
     
DisplaySequenceParameters[seqParams_] := Module[{paramSpec, itemizedList}, 
     paramSpec = ReturnedParameters[];
     itemizedList = Map[{"\[RightTriangle] ", #1[[1]], "\[RightTeeArrow] ", TraditionalForm[seqParams[[ #1[[2]] ]]]}&, 
                        paramSpec];
     PrintNotebookNotification[LightBlue][BoxData[ToBoxes[TableForm[itemizedList]]]];
     Return[];
] 

PhzLocal[h_, z_] := PhzLocal[h, a] = If[h <= 1, 
          KroneckerDelta[h == 1, True], 
          Simplify[(1 - c[h] * z) * PhzLocal[h - 1, z] - ab[h] * (z^2) * PhzLocal[h - 2, z]]
     ]
     
QhzLocal[h_, z_] := QhzLocal[h, z] = If[h <= 1, 
          KroneckerDelta[h == 0, True] + KroneckerDelta[h == 1, True] * (1 - c[1] * z), 
          Simplify[(1 - c[h] * z) * QhzLocal[h - 1, z] - ab[h] * (z^2) * QhzLocal[h - 2, z]] 
     ]
     
ConvLocal[h_, z_] := ConvLocal[h, z] = Simplify[PhzLocal[h, z] / QhzLocal[h, z]]

LambdaLocal[h_] := Product[ab[i], {i, 2, h}] 

MhLocal[h_] := LambdaLocal[h] 
     
ComputeSequenceParameters[seqFunc_, maxN_] := Module[{z, upperConvFunc, jFracEqns, seqVars, seqSols, abSeq, cSeq, lambdaParams, MhParams, convH, numVars}, 
     
     convH = If[OddQ[maxN], (maxN - 1) / 2, maxN / 2];
     numVars = If[EvenQ[maxN], maxN / 2, (maxN - 1) / 2];
     upperConvFunc = ConvLocal[convH, z];
     jFracEqns = Table[SeriesCoefficient[upperConvFunc, {z, 0, n}] == FunctionExpand[ seqFunc[n] ], {n, 1, 2 * numVars - 2}];
     seqVars = Flatten[Table[{c[n - 1], ab[n]}, {n, 2, numVars}]];
     seqSols = First[Solve[jFracEqns, seqVars]];
     abSeq = Table[ab[n] -> (ab[n] /. seqSols), {n, 2, numVars}];
     cSeq = Table[c[n] -> (c[n] /. seqSols), {n, 1, numVars - 1}];
     lambdaParams = Map[(LambdaLocal[#1] /. seqSols)&, Range[2, convH]];
     MhParams = Map[(MhLocal[#1] /. seqSols)&, Range[2, convH]];
     Return[{abSeq, cSeq, lambdaParams, MhParams}];
] 
     
VerifyParameters[seqFunc_, maxN_] := VerifyParametersByList[Table[seqFunc[n], {n, 1, maxN}]]

HthConvergentPhz[seqFunc_, h_, z_] := HthConvergentPhz[seqFunc, h, z] = 
     Module[{seqParams}, 
     
     seqParams = ComputeSequenceParameters[seqFunc, 2 * h - 1];
     seqSubsts = Union[seqParams[[ABSequenceIndex]], seqParams[[CSequenceIndex]]];
     Return[PhzLocal[h, z] /. seqSubsts];
]
     
HthConvergentQhz[seqFunc_, h_, z_] := HthConvergentQhz[seqFunc, h, z] = 
     Module[{seqParams}, 
     
     seqParams = ComputeSequenceParameters[seqFunc, 2 * h - 1];
     seqSubsts = Union[seqParams[[ABSequenceIndex]], seqParams[[CSequenceIndex]]];
     Return[QhzLocal[h, z] /. seqSubsts];
]
     
HthConvergent[seqFunc_, h_, z_] := HthConvergent[seqFunc, h, z] = 
     Module[{seqParams}, 
     
     seqParams = ComputeSequenceParameters[seqFunc, 2 * h - 1];
     seqSubsts = Union[seqParams[[ABSequenceIndex]], seqParams[[CSequenceIndex]]];
     Return[ConvLocal[h, z] /. seqSubsts];
]
     
EnumerateCongruenceSeries[seqFunc_, h_, z_, printSummary_:True] := 
     Module[{hthConvFunc, MhModulus, extendedSeqTerms, listItems}, 
     
     hthConvFunc = HthConvergent[seqFunc, h, z];
     seqParams = ComputeSequenceParameters[seqFunc, 2 * h - 1];
     seqSubsts = Union[seqParams[[ABSequenceIndex]], seqParams[[CSequenceIndex]]];
     MhModulus = MhLocal[h - 1] /. seqSubsts;
     If[!printSummary, 
          Return[{hthConvFunc, MhModulus}];
     ];
     extendedSeqTerms = Table[SeriesCoefficient[hthConvFunc, {z, 0, n}], {n, 0, 4 * h}];
     listItems = {"\[RightTriangle] ", Simplify[hthConvFunc], Mh -> MhModulus, SeqTerms -> extendedSeqTerms};
     PrintNotebookNotification[Cyan][BoxData[ToBoxes[TableForm[listItems]]]];
     Return[{hthConvFunc, MhModulus}];
]

(*******************************************************************************)
(**** : Routines for importing and formatting sequence data from the OEIS : ****) 
(*******************************************************************************)

OEISDatabaseLoad[seqANumber_] := Module[{entryHTML, percentMarkedLines}, 
     entryHTML = Import["https://oeis.org/" <> seqANumber <> "/internal", "Source"];
     percentMarkedLines = StringCases[entryHTML, RegularExpression["%([A-Z]+) (.*?)\\n"] -> {"$1", "$2"}];
     Return[percentMarkedLines];
] 

OEISSequenceLoad[seqANumber_] := 
     Module[{fullDBEntry, seqData, splicedSeqData, seqString, intSeqTerms}, 
     
     fullDBEntry = OEISDatabaseLoad[seqANumber]; 
     seqData = Select[fullDBEntry, MemberQ[{"S", "T", "U"}, First[#]]&]; 
     seqData = Map[#1[[2]]&, seqData]; 
     splicedSeqData = Flatten[Map[{ #1, "\[RawComma]" }&, Drop[seqData, -1]]];
     seqData = Append[splicedSeqData, Last[seqData]]; 
     seqString = StringReplace[StringJoin @@ seqData, " " -> ""]; 
     seqString = "{" <> seqString <> "}";
     intSeqTerms = Flatten@ToExpression@StringSplit[seqString, "\[RawComma]"];
     Return[intSeqTerms];
]

OEISPrintf[fmtStr_, fullSeqEntry_] := Module[{}, 
     formatSpecs = StringCases[fmtStr, RegularExpression["%([A-Z]+)"] -> ("$1" -> Select[fullSeqEntry, First[#] == "$1"&])];
     Print[StringReplace[fmtStr, formatSpecs]];
     Return[];
]

ExtractSequenceFuncFromList[lst_] := If[#1 <= 0 || #1 > Length[lst], If[#1 < 0, 0, 1] * s[#1], lst[[#1]]]&

(*******************************************************************************)
(**** : Sequence transformations (Largely for future imporovements) :       ****) 
(*******************************************************************************)

SequenceTransformationsSpec = {
     {"AERATE", 0&}, 
     {"ALLDIFFS", 0&}, 
     {"ANDCONV", 0&}, 
     {"BINOMIAL", Sum[Binomial[#2, k] * #1[[k]], {k, 0, #2}]&}, 
     {"BINOMIALi", 0&}, 
     {"BIN1", 0&}, 
     {"BISECT(a, 0)", 0&}, 
     {"BISECT(a, 1)", 0&}, 
     {"BOUS", 0&}, 
     {"BOUS2", 0&}, 
     {"BOUS2i", 0&}, 
     {"CHAR", 0&}, 
     {"COMP", 0&}, 
     {"COMPl", 0&}, 
     {"COMPOSE", 0&}, 
     (* ... : See the complete listing at https://oeis.org/transforms.txt ... *)
     {"rev", SeriesCoefficient[InverseSeries[1 + Sum[#1[n] * Power[z, n], {n, 1, #2}], z], {z, 0, #2}]&}, 
     {"inv", SeriesCoefficient[1 / (1 + Sum[#1[n] * Power[z, n], {n, 1, #2}]), {z, 0, #2}]&}, 
     {"exp", SeriesCoefficient[Exp[Sum[#1[n] * Power[z, n], {n, 1, #2}]], {z, 0, #2}]&}, 
     {"log", SeriesCoefficient[Log[1 + Sum[#1[n] * Power[z, n], {n, 1, #2}]], {z, 0, #2}]&}, 
     {"int", #1[#2 - 1] / #2&}, 
     {"der", (n + 1) * #1[#2]&}, 
     {"logder", SeriesCoefficient[D[Log[1 + Sum[#1[n] * Power[z, n], {n, 1, #2}]], z], {z, 0, #2}]&}, 
     {"bgf", Mod[#1[#2], 2, 0]&}
};
     

(*************************************************************************)
(**** : Text processing, display and general utility functions      : ****) 
(*************************************************************************)

GetDatestamp[] := ToString[StringForm["`1`/`2`/`3` : `4`:`5`:`6`", ##]] @@ 
     Map[DateValue, {"Month", "Day", "Year", "Hour", "Minute", "Second"}]

PrintNotebookNotificationDetailed[textColor_, backgroundColor_, borderColor_] := 
Function[cellText, 
     If[$Notebooks,
     CellPrint[Cell[cellText, "Print", 
                    FontColor -> textColor, 
                    FontSize -> 14, 
                    CellFrame -> 2.5, 
                    CellFrameColor -> borderColor, 
                    Background -> backgroundColor, 
                    CellFrameMargins -> 14, 
                    ShowCellBracket -> False, 
                    PageWidth -> WindowWidth
                   ]
              ],
     Print[cellText]; 
    ]
];

PrintNotebookNotification[baseColor_:Green] := 
     PrintNotebookNotificationDetailed[Darker[baseColor, 0.8], Lighter[baseColor, 0.5], Black];

PrintError[errorText_] := 
     PrintNotebookNotification[Red][StringJoin["\[ScriptX]\n", " !!! ERROR : >>> @ ", 
                                               GetDatestamp[], "\n", errorText, "\n\[ScriptX]"]]

(*************************************************************************)
(**** : Perform pretty printing of the package details and          : ****) 
(**** : revsision information if the package is loaded in a notebook: ****)
(*************************************************************************)

(** Local package loaded announcements: **)
packageAnnouncement = ToString["JFractions`"] <> " Package (Rev. 2018.09.20-v1) Loaded!\n\n";
authorAnnouncement = "Author: Maxie D. Schmidt (maxieds@gmail.com)\n"; 
usageAnnouncement = "Usage: See https://github.com/maxieds/INTEGERS2018-JFractions for complete documentation.\n\n";
descAnnouncement = "This package is a combinatorially insired tribute to combinatorial number theory and experiemental mathemetics " <> 
     "created to accompany the authors talk at the INTEGERS 2018 conference. It is my aim that the computational data " <> 
     "generated by this package will yield new and interesting results when applied to sequence classes which I have " <> 
     "not touched on in my own research related to J-fraction expansions. Enjoy the package, and please " <> 
     "contibute your most interesting finds with this package to the working list of examples maintained " <> 
     "on the package source repository page!";
PrintNotebookNotification[][packageAnnouncement <> authorAnnouncement <> usageAnnouncement <> descAnnouncement]

EndPackage[]

(*************************************************************************)
(*************************************************************************)
(*************************************************************************)
(*************************************************************************)
