
<< ct.m 

<< Singular.m

f = Together[2y(1-x)x(x+1)(x+2)(t+x)(x y-y-t^4)/(1-x(2-x+(x+1)(x+2)(t+x)(x y-y-t^4)^2))];
(* f = t/(1-x-y-y^2-t) *)

L = ctRat6[f, t, x, y][[2, 1]];

Lf = Together[ApplyOreOperator[L, f]];

{r, ab} = HermiteReduce[Lf, y];
Print["hermite passed: ", Short[Together[Lf - (D[r, y] + ab)]]];

betaMP = Factor[Denominator[ab]];  (* now y represents a root of betaMP, viz. a singularity of f *)
betaMP = Last[betaMP]; 
dy = -D[betaMP, x]/D[betaMP, y];

(* y1 == t^4/(x-1) - 1/sqrt(x(x+1)(x+2)(x+t)) *)
(* y2 == t^4/(x-1) + 1/sqrt(x(x+1)(x+2)(x+t)) == 2t^4/(x-1) - y1 *)
flip = {y -> u - y} /. First[Solve[Thread[CoefficientList[(betaMP/.y->u-y)-betaMP, y]==0],{u}]];
Print["flip passed: ", Short[Together[betaMP - (betaMP /. flip)]]];

alphaMP = Factor[Resultant[Numerator[ab] - z D[Denominator[ab], y], betaMP, y]];
alphaMP = Last[alphaMP]; 

(* choose a root alpha of alphaMP in K(x,y) *)
alpha = Factor[u + v y /. First[Solve[Thread[SingularMinAssGTZ[Numerator[Together[
  CoefficientList[PolynomialRemainder[alphaMP /. z -> u + v y, betaMP, y], y]]], {u, v}][[1]]==0], {u, v}]]]

Print["alpha passed: ", Short[Together[PolynomialRemainder[alphaMP /. z -> alpha, betaMP, y]]]];

gamma = Last[ctAlgNew[alpha, betaMP, t, x, y]];
Print["gamma passed: ", Short[Together[PolynomialRemainder[D[gamma, x] + D[gamma,y]dy - alpha, betaMP, y]]]];

g0 = ExpandDenominator[Together[gamma/(z - y) + (gamma/(z - y) /. flip)]];
g0 = Together[PolynomialRemainder[Numerator[g0], betaMP, y]/PolynomialRemainder[Denominator[g0], betaMP, y]];
g0 = g0 /. z -> y;

h0 = ExpandDenominator[Together[(r/.y->z) + gamma*dy/(z - y) + (gamma*dy/(z - y) /. flip)]];
h0 = Together[PolynomialRemainder[Numerator[h0], betaMP, y]/PolynomialRemainder[Denominator[h0], betaMP, y]];
h0 = h0 /. z -> y;

check = Together[Lf - (D[g0, x] + D[h0, y])];

Print["check passed: ", Short[check]];
