
(* Manuel Kauers, 2025-07-10 *)

(* << HolonomicFunctions.m *)
<< Asymptotics.m

SmallestRoot[p_, x_] := Module[ {roots, abs},
   roots = x /. Solve[p==0, x] // RootReduce;
   abs = Abs /@ roots;
   roots = Cases[Transpose[{abs, roots}], {Min[abs], _}];
   If[ Length[roots] != 1, Throw["no unique dominant finite singularity"] ];
   Return[ToNumberField[roots[[1,2]]]];   
];

DominantTerm[p_, x_, y_] := Module[ {alpha, ser, exponents},
   alpha = SmallestRoot[Discriminant[p, y], x];
   ser = Series[y /. First[Solve[0==p /. x -> alpha*(1-x), y]], {x, 0, 1}];
   exponents = Range[ser[[-3]], ser[[-2]]]/ser[[-1]];
   Do[
     If[ exponents[[i]] < 0 || !IntegerQ[exponents[[i]]],
       Return[{ser[[3,i]]/Gamma[-exponents[[i]]], 1/alpha, -1-exponents[[i]]}]
     ]
   , {i, 1, Length[exponents]}]
];

AlgebraicAsymptotics[p_, x_, y_, r_Integer] := Module[ {init, rec, expansion},
   init = DominantTerm[p, x, y];
   rec = ApplyOreOperator[First[DFiniteDE2RE[Annihilator[Evaluate[y /. First[Solve[p==0, y]]], {Der[x]}], x, n]], f[n]];
   expansion = FurtherTerms[rec, f[n], {0, init[[2]], {}, init[[3]], {1}}, r];
   init[[1]] expansion[[2]]^n n^expansion[[4]] Sum[expansion[[5, i]]n^(1-i), {i, 1, Length[expansion[[5]]]}]
]