# GeneralizedAbel v2.4 -- Maple 2026, exact arithmetic over Q(x).
# v2.4: lazy pole descriptors, specialization before symbolic products,
#       and jet refinement in a pre-solved Wronskian space.
# v2: section-7 linear screening, local jets, exact monomial-rank pruning.
# v2.1: verbose=true displays the active branch and current computation stage.
# v2.2: the combined entry computes a polynomial inverse integrating factor
#       directly for the n=2 stage, following Full algorithm.txt's ansatz.
# v2.3: direct Wronskian linear reduction and unique-candidate verification.
# Load with read "C:/.../GeneralizedAbel_Maple.txt":
# IMPORTANT: ReduceHigher is complete for n >= 3.
# ReduceRiccatiRegular retains the earlier transformation-search interface.
# Reduce's n=2 stage now searches polynomial inverse integrating factors.
# Reduce is NOT a complete decision procedure for unrestricted n >= 2.
# A failed combined search returns UNKNOWN, never an unjustified NO.

GeneralizedAbel := module()
option package;
export Reduce, ReduceHigher, ReduceRiccatiRegular, Verify, Pullback,
       RiccatiInverseFactor, InverseFactor;
local CheckInput, Equations, RationalPoints, Fit, MakeResult,
      HigherCore, RiccatiCore, MakeQuadratic, GoodPoint, JetFactor,
      LinearJetSpace, MonomialInconsistent, ExactScreen, GroebnerCalls,
      PoleDivisors, InverseFactorCore, WronskianSpace;
GroebnerCalls:=0;

CheckInput := proc(M,N,x,y,u)
    if nops({x,y,u}) <> 3 or not type([x,y,u],list(name)) then
        error "x, y, u must be distinct unassigned names"
    end if;
    if N=0 or not type(M,polynom(rational,[x,y]))
           or not type(N,polynom(rational,[x,y])) then
        error "M,N must be polynomials over Q, with N nonzero"
    end if;
    NULL
end proc;

Equations := proc(p,y)
    local q;
    q := expand(numer(normal(p)));
    if q=0 then return {} end if;
    {coeffs(q,y)} minus {0}
end proc;

# All K-rational points of a ZERO-DIMENSIONAL algebraic set, K=Q(x).
# This deliberately does not accept an arbitrary solve() parameterization.
# A missing eliminant is an error, not a proof of nonexistence.
RationalPoints := proc(eqs,vs::list,x,verbose)
    local es,G,z,rest,g,h,fac,roots,r,ans,sol,eqn,t0;
    es := {seq(numer(normal(eqn)),eqn in eqs)} minus {0};
    if nops(vs)=0 then
        if es={} then return [{}] else return [] end if
    end if;
    if es={} then error "unexpected positive-dimensional candidate set" end if;
    if verbose then
        printf("[GA][points] Groebner start: variables=%d, equations=%d, eliminate=%a\n",
               nops(vs),nops(es),vs[-1])
    end if;
    t0:=time();
    GroebnerCalls:=GroebnerCalls+1;
    G := Groebner:-Basis([op(es)],plex(op(vs)));
    if verbose then printf("[GA][points] Groebner done: basis=%d, CPU=%.3f s\n",
                           nops(G),time()-t0) end if;
    if member(1,G) then
        if verbose then printf("[GA][points] inconsistent system; no rational points\n") end if;
        return []
    end if;
    z := vs[-1]; rest := vs[1..-2]; g := 0;
    for h in G do
        if degree(h,z)>0 and not has(h,{op(rest)}) then
            if g=0 or degree(h,z)<degree(g,z) then g:=h end if
        end if
    end do;
    if g=0 then error "unexpected positive-dimensional candidate set" end if;
    fac := factors(numer(normal(g)))[2]; roots := {};
    for h in fac do
        if degree(h[1],z)=1 then
            r := normal(-coeff(h[1],z,0)/coeff(h[1],z,1));
            if type(r,ratpoly(rational,x)) then roots:=roots union {r} end if
        end if
    end do;
    if verbose then printf("[GA][points] rational roots for %a: %d\n",z,nops(roots)) end if;
    ans := [];
    for r in roots do
        for sol in thisproc(subs(z=r,G),rest,x,verbose) do
            ans := [op(ans),sol union {z=r}]
        end do
    end do;
    if verbose then printf("[GA][points] complete: variables=%d, points=%d\n",
                           nops(vs),nops(ans)) end if;
    ans
end proc;

# Degree-aware subset generation. Return [degree,factor-descriptors] only;
# no symbolic denominator is constructed before specialization screening.
PoleDivisors := proc(eligible,lower,y)
    local weights,suffix,count,ii,ans,walk;
    count:=nops(eligible);
    weights:=[seq(degree(eligible[ii][1],y)*eligible[ii][3],ii=1..count)];
    suffix:=Array(1..count+1,0);
    for ii from count by -1 to 1 do suffix[ii]:=suffix[ii+1]+weights[ii] end do;
    ans:=[];
    walk:=proc(pos,deg,chosen)
        local h;
        if deg+suffix[pos]<lower then return NULL end if;
        if pos>count then
            ans:=[op(ans),[deg,chosen]];
            return NULL
        end if;
        thisproc(pos+1,deg,chosen);
        thisproc(pos+1,deg+weights[pos],[op(chosen),eligible[pos]]);
        NULL
    end proc;
    walk(1,0,[]);
    sort(ans,(ll,rr)->evalb(ll[1]<rr[1]))
end proc;

# Choose an exact, regular specialization. No floating point decisions.
GoodPoint := proc(M,N,fl,x,y)
    local q,v,M0,N0,rad,p,h,ok;
    q:=0;
    do
        if q=0 then v:=0 elif irem(q,2)=1 then v:=iquo(q+1,2)
        else v:=-iquo(q,2) end if;
        q:=q+1;
        try
            M0:=subs(x=v,M); N0:=subs(x=v,N); rad:=1; ok:=true;
            if degree(M0,y)<>degree(M,y) or degree(N0,y)<>degree(N,y) then next end if;
            for h in fl do
                p:=normal(subs(x=v,h[1]));
                if not type(p,polynom(rational,y)) or degree(p,y)<>degree(h[1],y) then
                    ok:=false; break
                end if;
                rad:=expand(rad*p)
            end do;
            if ok and degree(gcd(rad,diff(rad,y)),y)=0
                  and degree(gcd(M0,N0),y)=0 then return v end if
        catch:
            # A pole of an input coefficient is a bad point, not a rejection.
        end try
    end do
end proc;

# J is the truncated local primitive: (p^(k+1)*J)' = N/M mod p^(k+e).
JetFactor := proc(M,Np,p,k,e,y)
    option remember;
    local jj,vs,J,ss,ii,eq;
    vs:=[seq(jj[ii],ii=0..e*degree(p,y)-1)];
    J:=add(jj[ii]*y^ii,ii=0..e*degree(p,y)-1);
    eq:=Equations(rem(expand(M*((k+1)*diff(p,y)*J+p*diff(J,y))-Np),p^e,y),y);
    ss:=[solve(eq,{op(vs)})];
    if nops(ss)<>1 then error "local linear jet system did not have a unique solution" end if;
    J:=normal(subs(ss[1],J));
    if has(J,{op(vs)}) then error "undetermined local jet coefficient" end if;
    J
end proc;

# Necessary homogeneous space W=T*C, T=N/B^(n-2), over Q(x).
# Use polynomial coefficients before elimination, avoiding remainders by
# a large monic critical factor with rational-function coefficients.
WronskianSpace := proc(N,B,n,a,y)
    local b,d,BB,T,ids,zz,qq,A,C,vs,cs,es,ss,pars,ii;
    b:=degree(B,y); d:=a+(n-1)*b-1-degree(N,y);
    if d<0 then return [0,[]] end if;
    BB:=numer(normal(B)); T:=normal(N/B^(n-2));
    if degree(denom(T),y)>0 then return FAIL end if;
    T:=numer(T);
    ids:=remove(ii->evalb(ii=b),[seq(ii,ii=0..a)]);
    vs:=[seq(zz[ii],ii in ids)]; cs:=[seq(qq[ii],ii=0..d)];
    A:=add(zz[ii]*y^ii,ii in ids); C:=add(qq[ii]*y^ii,ii=0..d);
    es:={coeffs(expand(BB*diff(A,y)-A*diff(BB,y)-T*C),y)} minus {0};
    ss:=[solve(es,{op(vs),op(cs)})];
    if nops(ss)<>1 or not type(ss[1],set) then return FAIL end if;
    A:=collect(subs(ss[1],A),y,normal);
    pars:=select(ii->has(A,ii),[op(vs),op(cs)]);
    [A,pars]
end proc;

# Homogeneous linear space for A, with the translation coefficient A_b=0.
# Works over Q for screening and over Q(x) for exact parameter reduction.
LinearJetSpace := proc(M,N,fl,eps,B,n,a,einf,y,seed:=FAIL)
    local b,ii,ids,zz,vs,A,W,es,h,p,k,e,C,Np,J,z,Mr,Nr,Ar,Cr,
          eqs,mat,ns,pp,pars,Ard,rr,cc;
    b:=degree(B,y);
    ids:=remove(ii->evalb(ii=b),[seq(ii,ii=0..a)]);
    if seed=FAIL then
        vs:=[seq(zz[ii],ii in ids)];
        A:=add(zz[ii]*y^ii,ii in ids)
    else
        A:=seed[1]; vs:=seed[2];
        if vs=[] then return [0,[]] end if
    end if;
    es:={};
    if seed=FAIL then
    W:=expand(B*diff(A,y)-A*diff(B,y));
    for h in fl do
        if rem(B,h[1],y)<>0 then
            es:=es union Equations(rem(W,h[1]^h[2],y),y)
        end if
    end do;
    end if;
    for h in eps do
        p:=h[1]; k:=h[2]; e:=h[3];
        if e>1 then
            C:=normal(B/p^e); Np:=normal(N/p^k);
            J:=JetFactor(M,Np,p,k,e,y);
            es:=es union Equations(rem(numer(normal(
                   (n-1)*J*(C*diff(A,y)-A*diff(C,y))+diff(J,y)*C*A)),p^(e-1),y),y)
        end if
    end do;
    if einf>1 then
        Mr:=-expand(z^degree(M,y)*subs(y=1/z,M));
        Nr:=expand(z^degree(N,y)*subs(y=1/z,N));
        Ar:=expand(z^a*subs(y=1/z,A)); Cr:=expand(z^b*subs(y=1/z,B));
        k:=(n-1)*einf-1;
        J:=JetFactor(Mr,Nr,z,k,einf,z);
        es:=es union Equations(rem(numer(normal(
               (n-1)*J*(Cr*diff(Ar,z)-Ar*diff(Cr,z))+diff(J,z)*Cr*Ar)),z^(einf-1),z),z)
    end if;
    eqs:=[op(es)];
    mat:=Matrix(nops(eqs),nops(vs),(rr,cc)->normal(coeff(eqs[rr],vs[cc],1)));
    ns:=[op(LinearAlgebra:-NullSpace(mat))];
    pars:=[seq(pp[ii],ii=1..nops(ns))];
    Ard:=subs({seq(vs[rr]=add(pp[ii]*ns[ii][rr],ii=1..nops(ns)),rr=1..nops(vs))},A);
    [collect(Ard,y,normal),pars]
end proc;

# If even independent values for all occurring monomials cannot satisfy the
# equations, there is no algebraic solution. This is only a rejection test.
MonomialInconsistent := proc(eqs,vs)
    local e,cs,ms,rows,mons,cols,h,mat,rr,j,idx,eq;
    rows:=[]; mons:={};
    for eq in eqs do
        e:=expand(numer(normal(eq)));
        if e=0 then next end if;
        cs:=[coeffs(e,{op(vs)},'ms')]; ms:=[ms];
        rows:=[op(rows),[cs,ms]];
        mons:=mons union {op(ms)}
    end do;
    if rows=[] then return [false,0,0,0] end if;
    mons:=[op(mons minus {1})]; cols:=table();
    for j to nops(mons) do cols[mons[j]]:=j end do;
    cols[1]:=nops(mons)+1;
    mat:=Matrix(nops(rows),nops(mons)+1,0);
    for rr to nops(rows) do
        for j to nops(rows[rr][1]) do
            mat[rr,cols[rows[rr][2][j]]]:=rows[rr][1][j]
        end do
    end do;
    if mons=[] then return [true,0,0,1] end if;
    rr:=LinearAlgebra:-Rank(LinearAlgebra:-SubMatrix(mat,1..nops(rows),1..nops(mons)));
    j:=LinearAlgebra:-Rank(mat);
    [evalb(j>rr),nops(mons),rr,j]
end proc;

# Work over the algebraic closure at a good rational x value. We set f_n=1
# by an algebraic scaling of r, not by assuming f_n was originally 1.
# All principal parts are regular there; the independent-monomial test is
# therefore a sound necessary condition for a generic rational solution.
ExactScreen := proc(M,N,fl,eps,B,n,a,einf,y)
    local model,A,vs,es,h,p,k,e,C,Np,pol,check;
    # At the rational screening point, first solve the nonpole constraints.
    # Empty eps and zero infinity order postpone all jet computations.
    model:=LinearJetSpace(M,N,fl,[],B,n,a,0,y);
    model:=LinearJetSpace(M,N,fl,eps,B,n,a,einf,y,model);
    A:=model[1]; vs:=model[2];
    if vs=[] then return ["REJECT_LINEAR",0,0,0,0] end if;
    # Keep screening bounded in large, highly nongeneric surviving spaces.
    if binomial(nops(vs)+n-2,n-1)>500 then return ["PASS",nops(vs),-1,-1,-1] end if;
    es:={};
    for h in eps do
        p:=h[1]; k:=h[2]; e:=h[3];
        C:=normal(B/p^e); Np:=normal(N/p^k);
        pol:=expand(Np*A^(n-1)+e*M*diff(p,y)*C^(n-1));
        es:=es union Equations(rem(pol,p,y),y)
    end do;
    if einf>0 then
        es:=es union {coeff(A,y,a)^(n-1)-einf*lcoeff(M,y)/lcoeff(N,y)}
    end if;
    check:=MonomialInconsistent(es,vs);
    if check[1] then return ["REJECT_LIFTED",nops(vs),check[2],check[3],check[4]] end if;
    ["PASS",nops(vs),check[2],check[3],check[4]]
end proc;

Verify := proc(M,N,r,P,x,y,u)
    local res;
    CheckInput(M,N,x,y,u);
    if not type(r,ratpoly(rational,[x,y])) or normal(diff(r,y))=0 then
        return table(["status"="INVALID_TRANSFORMATION"])
    end if;
    if not type(P,polynom(ratpoly(rational,x),u)) or degree(P,u)<2 then
        return table(["status"="INVALID_TARGET"])
    end if;
    res:=normal(diff(r,x)+(M/N)*diff(r,y)-subs(u=r,P));
    table(["status"=ifelse(res=0,"VERIFIED","FAILED"),
           "transformation"=normal(r),"polynomial"=collect(P,u,normal),
           "n"=degree(P,u),"residual"=res])
end proc;

Pullback := proc(r,P,x,y,u)
    if normal(diff(r,y))=0 then error "r_y must be nonzero" end if;
    normal((subs(u=r,P)-diff(r,x))/diff(r,y))
end proc;

# Exact membership test: D(r) in span_K(1,r,...,r^n).
Fit := proc(M,N,r,n,x,y,u)
    local A,B,W,L,ff,vs,es,ss,P,j;
    A:=numer(normal(r)); B:=denom(normal(r));
    W:=expand(B*diff(A,y)-A*diff(B,y));
    if W=0 then return FAIL end if;
    vs:=[seq(ff[j],j=0..n)];
    L:=expand((N*(B*diff(A,x)-A*diff(B,x))+M*W)*B^(n-2)
              -N*add(ff[j]*A^j*B^(n-j),j=0..n));
    es:=Equations(L,y);
    ss:=[solve(es,{op(vs)})];
    if nops(ss)<>1 then return FAIL end if;
    P:=normal(subs(ss[1],add(ff[j]*u^j,j=0..n)));
    if not type(P,polynom(ratpoly(rational,x),u)) then return FAIL end if;
    if normal(diff(r,x)+M/N*diff(r,y)-subs(u=r,P))<>0 then
        error "internal membership verification failed"
    end if;
    collect(P,u,normal)
end proc;

MakeResult := proc(M,N,r,P,x,y,u,branch)
    local ans,j;
    ans:=Verify(M,N,r,P,x,y,u);
    if ans["status"]<>"VERIFIED" then error "internal certificate failure" end if;
    ans["status"]:="FOUND";
    ans["branch"]:=branch;
    ans["coefficients"]:=[seq(normal(coeff(P,u,j)),j=0..degree(P,u))];
    ans["equation"]:=diff(u(x),x)=subs(u=u(x),P);
    ans["normalization"]:="D = d/dx + (M/N) d/dy; x is unchanged";
    ans
end proc;

HigherCore := proc(M,N,x,y,u,verbose,pruning)
    local R,MM,NN,m,j,fl,h,p,k,s,nmax,n,eligible,bs,bb,B,b,
          choices,aa,ainf,ii,A,av,c,sat,vars,H,es,ep,L,t,sols,sol,
          r,P,ans,systems,tested,resB,cc,unknownc,pol,W,required,nonpoles,Cp,Np,
          point,M0,N0,fl0,eps0,B0,screen,model,apars,pivot,lc,stat,lower,
          started,bi,ci,branchid,selected_ids,lastdegree;
    started:=time();
    if verbose then printf("[GA][higher] start: normalize M/N\n") end if;
    forget(JetFactor);
    R:=normal(M/N); MM:=numer(R); NN:=denom(R);
    GroebnerCalls:=0;
    stat:=table(["degree_candidates"=0,"pruned_linear"=0,"pruned_lifted"=0,
                 "pruned_symbolic_linear"=0,"screen_passed"=0,"groebner_calls"=0,
                 "wronskian_unique"=0,"fixed_candidates"=0,
                 "specialized_denominators"=0,"symbolic_denominators"=0,
                 "raw_denominator_descriptors"=0,
                 "enabled"=pruning,"branches"=[]]);
    if MM=0 then
        if verbose then printf("[GA][higher] zero slope: no degree >=3 target\n") end if;
        return table(["status"="NO_N_GE_3","systems"=0,"candidates"=0,"screening"=eval(stat)])
    end if;
    m:=degree(MM,y)-degree(NN,y); j:=degree(NN,y);
    if verbose then printf("[GA][higher] factor denominator: deg_y(N)=%d, deg_y(M)-deg_y(N)=%d\n",j,m) end if;
    fl:=[];
    for h in factors(NN)[2] do
        if degree(h[1],y)>0 then
            p:=normal(h[1]/lcoeff(h[1],y));
            fl:=[op(fl),[p,h[2]]]
        end if
    end do;
    if verbose then
        for ii to nops(fl) do
            printf("[GA][higher] pole p%d: deg_y=%d, ord_N=%d, polynomial=%a\n",
                   ii,degree(fl[ii][1],y),fl[ii][2],fl[ii][1])
        end do
    end if;
    if verbose then printf("[GA][higher] denominator factors=%d; possible target degrees=3..%d\n",
                           nops(fl),max(2,m,op([seq(h[2]+2,h in fl)]))) end if;
    if pruning then
        if verbose then printf("[GA][higher] find regular x-specialization for exact screening\n") end if;
        point:=GoodPoint(MM,NN,fl,x,y);
        M0:=subs(x=point,MM); N0:=subs(x=point,NN); fl0:=subs(x=point,fl);
        stat["specialization_point"]:=point;
        if verbose then printf("[GA][higher] screening point x=%a\n",point) end if
    end if;
    nmax:=max(2,m,op([seq(h[2]+2,h in fl)]));
    systems:=0; tested:=0;
    for n from 3 to nmax do
        eligible:=[];
        for h in fl do
            if irem(h[2]+1,n-1)=0 then
                eligible:=[op(eligible),[h[1],h[2],iquo(h[2]+1,n-1)]]
            end if
        end do;
        if verbose then printf("[GA][higher] degree n=%d/%d: eligible pole factors=%d\n",
                               n,nmax,nops(eligible)) end if;
        # degree(C)=a+(n-1)*b-1-j >= 0, together with the infinity case.
        if pruning then
            if m>1 then
                lower:=max(m-1,ceil((m+j)/n));
                if irem(m-1,n-1)=0 then
                    lower:=min(lower,max(0,ceil((j+1-iquo(m-1,n-1))/n)))
                end if
            else lower:=ceil((j+2)/n)
            end if;
            if verbose then printf("[GA][higher] n=%d: enumerate B with degree lower bound %d\n",n,lower) end if;
            bs:=PoleDivisors(eligible,lower,y)
        else
            bs:=[[1,[]]];
            for h in eligible do
                bs:=[op(bs),seq([expand(bb[1]*h[1]^h[3]),[op(bb[2]),h]],bb in bs)]
            end do;
            bs:=sort(bs,(ll,rr)->evalb(degree(ll[1],y)<degree(rr[1],y)))
        end if;
        stat["raw_denominator_descriptors"]:=stat["raw_denominator_descriptors"]+nops(bs);
        if verbose then printf("[GA][higher] n=%d: %d pole descriptors; screening precedes symbolic expansion\n",n,nops(bs)) end if;
        lastdegree:=-1;
        for bi to nops(bs) do
            bb:=bs[bi];
            if pruning then
                b:=bb[1];
                if verbose and b<>lastdegree then
                    printf("[GA][higher] degree frontier b=%d: %d descriptors\n",b,nops(select(h->evalb(h[1]=b),bs)))
                end if;
                lastdegree:=b;
                B0:=expand(mul(subs(x=point,h[1])^h[3],h in bb[2]));
                eps0:=subs(x=point,bb[2]);
                stat["specialized_denominators"]:=stat["specialized_denominators"]+1
            else B:=bb[1]; b:=degree(B,y)
            end if;
            if verbose then
                selected_ids:=[];
                for ii to nops(fl) do
                    for h in bb[2] do
                        if fl[ii][1]=h[1] then
                            selected_ids:=[op(selected_ids),[ii,h[3]]]
                        end if
                    end do
                end do;
                printf("[GA][higher] n=%d B=%d/%d: deg_y(B)=%d, pole IDs and B powers=%a\n",
                       n,bi,nops(bs),b,selected_ids)
            end if;
            if m>1 then
                if b>=m-1 then choices:=[[b-m+1,0]] else choices:=[] end if
            else
                choices:=[seq([ii,0],ii=0..b-1)]
            end if;
            if m>1 and irem(m-1,n-1)=0 then
                ainf:=iquo(m-1,n-1);
                choices:=[op(choices),[b+ainf,ainf]]
            end if;
            nonpoles:=[]; required:=0;
            for h in fl do
                if not member(h[1],[seq(ep[1],ep in bb[2])]) then
                    nonpoles:=[op(nonpoles),h];
                    required:=required+h[2]*degree(h[1],y)
                end if
            end do;
            required:=required+add((h[3]-1)*degree(h[1],y),h in bb[2]);
            if verbose then printf("[GA][higher] n=%d B=%d/%d: %d numerator-degree choices, critical degree bound=%d\n",
                                   n,bi,nops(bs),nops(choices),required) end if;
            for ci to nops(choices) do
                aa:=choices[ci];
                if aa[1]+b-1<required then
                    if verbose then printf("[GA][branch] n=%d B=%d/%d A=%d/%d (a=%d, b=%d): reject degree bound\n",
                                           n,bi,nops(bs),ci,nops(choices),aa[1],b) end if;
                    next
                end if;
                stat["degree_candidates"]:=stat["degree_candidates"]+1;
                branchid:=stat["degree_candidates"];
                if verbose then printf("[GA][branch %d] n=%d B=%d/%d A=%d/%d: a=%d, b=%d, infinity order=%d; screen start\n",
                                       branchid,n,bi,nops(bs),ci,nops(choices),aa[1],b,aa[2]) end if;
                if pruning then
                    screen:=ExactScreen(M0,N0,fl0,eps0,B0,n,aa[1],aa[2],y);
                    stat["branches"]:=[op(stat["branches"]),[n,aa[1],b,op(screen)]];
                    if verbose then
                        printf("[GA][branch %d] screen=%s; linear parameters=%d, lifted monomials=%d\n",
                               branchid,screen[1],screen[2],screen[3])
                    end if;
                    if screen[1]="REJECT_LINEAR" then
                        stat["pruned_linear"]:=stat["pruned_linear"]+1; next
                    elif screen[1]="REJECT_LIFTED" then
                        stat["pruned_lifted"]:=stat["pruned_lifted"]+1; next
                    end if;
                    stat["screen_passed"]:=stat["screen_passed"]+1;
                    B:=expand(mul(h[1]^h[3],h in bb[2]));
                    stat["symbolic_denominators"]:=stat["symbolic_denominators"]+1;
                    if verbose then printf("[GA][branch %d] exact Wronskian linear reduction over Q(x)\n",branchid) end if;
                    model:=WronskianSpace(NN,B,n,aa[1],y);
                    if model=FAIL or nops(model[2])>1 then
                        if verbose then printf("[GA][branch %d] remaining freedom: refine with exact linear jets\n",branchid) end if;
                        model:=LinearJetSpace(MM,NN,fl,bb[2],B,n,aa[1],aa[2],y,model)
                    elif nops(model[2])=1 then
                        stat["wronskian_unique"]:=stat["wronskian_unique"]+1;
                        if verbose then printf("[GA][branch %d] Wronskian space has one scaling parameter\n",branchid) end if
                    end if;
                    A:=model[1]; apars:=model[2]; lc:=normal(coeff(A,y,aa[1]));
                    if lc=0 then
                        if verbose then printf("[GA][branch %d] reject symbolic leading coefficient\n",branchid) end if;
                        stat["pruned_symbolic_linear"]:=stat["pruned_symbolic_linear"]+1; next
                    end if;
                    for pivot in apars do if coeff(lc,pivot,1)<>0 then break end if end do;
                    A:=normal(subs(pivot=normal((1-subs(pivot=0,lc))/coeff(lc,pivot,1)),A));
                    apars:=remove(ii->evalb(ii=pivot),apars)
                else
                if verbose then printf("[GA][branch %d] pruning disabled; build full numerator ansatz\n",branchid) end if;
                # Affine normalization: A monic and coefficient of y^b zero.
                av:=[seq(ii,ii=0..aa[1]-1)];
                av:=remove(ii->evalb(ii=b),av);
                A:=y^aa[1]+add(t[ii]*y^ii,ii in av);
                apars:=[seq(t[ii],ii in av)]
                end if;
                # A unique normalized A needs no auxiliary top-coefficient
                # system. Full membership and identity verification decide it.
                if pruning and apars=[] then
                    stat["fixed_candidates"]:=stat["fixed_candidates"]+1;
                    tested:=tested+1; r:=normal(A/B);
                    if verbose then printf("[GA][branch %d] unique normalized candidate: direct target fit and verification\n",branchid) end if;
                    P:=Fit(MM,NN,r,n,x,y,u);
                    if P<>FAIL and degree(P,u)=n then
                        ans:=MakeResult(M,N,r,P,x,y,u,"n>=3: exhaustive pole search");
                        ans["systems"]:=systems; ans["candidates"]:=tested;
                        stat["groebner_calls"]:=GroebnerCalls;
                        ans["screening"]:=eval(stat);
                        if verbose then printf("[GA][higher] FOUND n=%d; unique-candidate path, CPU=%.3f s\n",n,time()-started) end if;
                        return ans
                    end if;
                    if verbose then printf("[GA][branch %d] unique candidate rejected by full target identity\n",branchid) end if;
                    next
                end if;
                if verbose then printf("[GA][branch %d] build exact polynomial conditions\n",branchid) end if;
                if aa[2]>0 then
                    # Leading term at infinity fixes the top target coefficient.
                    cc:=normal(aa[2]*lcoeff(MM,y)/lcoeff(NN,y));
                    unknownc:=[]
                else
                    cc:=c; unknownc:=[c]
                end if;
                vars:=[op(unknownc),op(apars)];
                W:=expand(B*diff(A,y)-A*diff(B,y));
                H:=expand(MM*W*B^(n-2)
                          -cc*NN*A^n);
                es:={};
                for h in nonpoles do
                    es:=es union Equations(rem(W,h[1]^h[2],y),y)
                end do;
                for ep in bb[2] do
                    Cp:=normal(B/ep[1]^ep[3]);
                    Np:=normal(NN/ep[1]^ep[2]);
                    # Divide the leading local equation by A, which is a unit.
                    # Its nonzero RHS itself enforces c != 0 and gcd(A,B)=1.
                    pol:=numer(normal(cc*Np*A^(n-1)
                         +ep[3]*MM*diff(ep[1],y)*Cp^(n-1)));
                    es:=es union Equations(rem(pol,ep[1],y),y);
                    if ep[3]>1 then
                        pol:=numer(normal(H/ep[1]^ep[2]));
                        es:=es union Equations(rem(pol,ep[1]^ep[3],y),y)
                    end if
                end do;
                if aa[2]>0 then
                    L:=j+n*b;
                    es:=es union {seq(numer(normal(coeff(H,y,ii))),
                          ii=L+(n-1)*aa[2]+1..L+n*aa[2])}
                end if;
                systems:=systems+1;
                if verbose then printf("[GA][branch %d] rational points: system=%d, equations=%d, variables=%d\n",
                                       branchid,systems,nops(es),nops(vars)) end if;
                sols:=RationalPoints(es,vars,x,verbose);
                if verbose then printf("[GA][branch %d] rational points found=%d; verify candidates\n",
                                       branchid,nops(sols)) end if;
                for sol in sols do
                    r:=normal(subs(sol,A)/B); tested:=tested+1;
                    if verbose then printf("[GA][branch %d] fit target for candidate %d\n",branchid,tested) end if;
                    P:=Fit(MM,NN,r,n,x,y,u);
                    if P<>FAIL and degree(P,u)=n then
                        if verbose then printf("[GA][higher] FOUND n=%d; systems=%d, candidates=%d, CPU=%.3f s\n",
                                               n,systems,tested,time()-started) end if;
                        ans:=MakeResult(M,N,r,P,x,y,u,"n>=3: exhaustive pole search");
                        ans["systems"]:=systems; ans["candidates"]:=tested;
                        stat["groebner_calls"]:=GroebnerCalls;
                        ans["screening"]:=eval(stat);
                        return ans
                    end if
                end do
            end do
        end do
    end do;
    stat["groebner_calls"]:=GroebnerCalls;
    if verbose then printf("[GA][higher] NO_N_GE_3; branches=%d, systems=%d, candidates=%d, Groebner=%d, CPU=%.3f s\n",
                           stat["degree_candidates"],systems,tested,GroebnerCalls,time()-started) end if;
    table(["status"="NO_N_GE_3","systems"=systems,"candidates"=tested,"screening"=eval(stat),
           "meaning"="No rational nonconstant transformation to a polynomial of actual degree >=3"])
end proc;

ReduceHigher := proc(M,N,x,y,u,{seconds::nonnegint:=0,verbose::truefalse:=false,pruning::truefalse:=true})
    CheckInput(M,N,x,y,u);
    try
        if seconds=0 then return HigherCore(M,N,x,y,u,verbose,pruning)
        else return timelimit(seconds,HigherCore(M,N,x,y,u,verbose,pruning)) end if
    catch "time expired":
        if verbose then printf("[GA][higher] time limit reached after %d s; result is UNKNOWN_TIMEOUT\n",seconds) end if;
        return table(["status"="UNKNOWN_TIMEOUT","scope"="n>=3"])
    end try
end proc;

# Turn a target of degree <=2 into actual degree 2 by a K-Mobius map.
MakeQuadratic := proc(r,P,x,y,u)
    local s,v,Q;
    if degree(P,u)=2 then return [r,P] end if;
    for s in [0,1,x,2*x] do
        v:=normal(1/(r+s));
        Q:=normal(-u^2*(subs(u=1/u-s,P)+diff(s,x)));
        if degree(Q,u)=2 then return [v,collect(Q,u,normal)] end if
    end do;
    error "internal quadratic normalization failure"
end proc;

# Complete only when W=B*A_y-A*B_y has no canceled finite factor.
# Equivalently, after reduced input, W is a K* multiple of N in K[y].
RiccatiCore := proc(M,N,x,y,u,verbose)
    local R,MM,NN,N0,j,d,b,A,B,a,v,vars,es,sols,sol,W,rr,P,
          sat,rB,ans,ii,q,systems,tested,started;
    started:=time();
    if verbose then printf("[GA][Riccati] start: no finite ramification cancellation\n") end if;
    R:=normal(M/N); MM:=numer(R); NN:=denom(R);
    N0:=normal(NN/lcoeff(NN,y)); j:=degree(NN,y);
    systems:=0; tested:=0;
    if verbose then printf("[GA][Riccati] normalized deg_y(N)=%d; map degrees=%d..%d\n",
                           j,ceil((j+2)/2),j+1) end if;
    for d from ceil((j+2)/2) to j+1 do
        b:=j+1-d;
        if b<0 or b>=d then next end if;
        A:=y^d+add(a[ii]*y^ii,ii in remove(ii->evalb(ii=b),[seq(ii,ii=0..d-1)]));
        B:=y^b+add(v[ii]*y^ii,ii=0..b-1);
        vars:=[sat,seq(a[ii],ii in remove(ii->evalb(ii=b),[seq(ii,ii=0..d-1)])),
                   seq(v[ii],ii=0..b-1)];
        W:=expand(B*diff(A,y)-A*diff(B,y));
        es:=Equations(W-(d-b)*N0,y);
        rB:=resultant(A,B,y);
        es:=es union {sat*rB-1};
        systems:=systems+1;
        if verbose then printf("[GA][Riccati] branch %d: map degree=%d, deg_y(B)=%d, equations=%d, variables=%d; rational points start\n",
                               systems,d,b,nops(es),nops(vars)) end if;
        sols:=RationalPoints(es,vars,x,verbose);
        if verbose then printf("[GA][Riccati] branch %d: rational points found=%d; verify candidates\n",
                               systems,nops(sols)) end if;
        for sol in sols do
            rr:=normal(subs(sol,A/B)); tested:=tested+1;
            if verbose then printf("[GA][Riccati] fit target for candidate %d\n",tested) end if;
            P:=Fit(MM,NN,rr,2,x,y,u);
            if P<>FAIL then
                q:=MakeQuadratic(rr,P,x,y,u);
                if verbose then printf("[GA][Riccati] FOUND; systems=%d, candidates=%d, CPU=%.3f s\n",
                                       systems,tested,time()-started) end if;
                ans:=MakeResult(M,N,q[1],q[2],x,y,u,"Riccati: no finite ramification cancellation");
                ans["systems"]:=systems; ans["candidates"]:=tested;
                return ans
            end if
        end do
    end do;
    if verbose then printf("[GA][Riccati] NO_REGULAR_RICCATI; systems=%d, candidates=%d, CPU=%.3f s\n",
                           systems,tested,time()-started) end if;
    table(["status"="NO_REGULAR_RICCATI","systems"=systems,"candidates"=tested,
           "meaning"="Only the no-finite-ramification-cancellation Riccati class is excluded"])
end proc;

ReduceRiccatiRegular := proc(M,N,x,y,u,{seconds::nonnegint:=0,verbose::truefalse:=false})
    CheckInput(M,N,x,y,u);
    try
        if seconds=0 then return RiccatiCore(M,N,x,y,u,verbose)
        else return timelimit(seconds,RiccatiCore(M,N,x,y,u,verbose)) end if
    catch "time expired":
        if verbose then printf("[GA][Riccati] time limit reached after %d s; result is UNKNOWN_TIMEOUT\n",seconds) end if;
        return table(["status"="UNKNOWN_TIMEOUT","scope"="regular Riccati"])
    end try
end proc;

# Direct polynomial inverse integrating factor of N*d_x+M*d_y.
# This is the degree bound and coefficient-ODE template in Full algorithm.txt.
# A failure of dsolve or this bounded ansatz is not a proof of nonexistence.
InverseFactorCore := proc(M,N,x,y,verbose)
    local bound,V,vs,eqs,sol,res,z,divergence,started;
    started:=time();
    bound:=max(degree(M,y),degree(N,y)+1);
    divergence:=diff(N,x)+diff(M,y);
    if verbose then printf("[GA][IIF] start: polynomial ansatz deg_y(V)<=%d for original field [N,M]\n",bound) end if;
    if divergence=0 then
        if verbose then printf("[GA][IIF] divergence is zero; V=1, residual=0\n") end if;
        return table(["status"="FOUND_IIF","result_type"="inverse_integrating_factor",
                      "inverse_factor"=1,"residual"=0,"degree_bound"=bound,
                      "vector_field"=[N,M],"method"="direct polynomial ansatz"])
    end if;
    V:=add(vcoef[z](x)*y^z,z=0..bound);
    vs:={seq(vcoef[z](x),z=0..bound)};
    eqs:={coeffs(expand(V*divergence-M*diff(V,y)-N*diff(V,x)),y)} minus {0};
    if verbose then printf("[GA][IIF] coefficient equations=%d, unknown functions=%d; dsolve start\n",
                           nops(eqs),nops(vs)) end if;
    try
        sol:=dsolve(eqs,vs)
    catch "time expired":
        error "time expired"
    catch:
        if verbose then printf("[GA][IIF] dsolve could not resolve the coefficient system\n") end if;
        return table(["status"="UNKNOWN_IIF","degree_bound"=bound,
                      "reason"="coefficient ODE solver did not finish with a usable solution"])
    end try;
    if not (type(sol,set) or type(sol,list)) then
        if verbose then printf("[GA][IIF] dsolve returned no usable coefficient solution\n") end if;
        return table(["status"="UNKNOWN_IIF","degree_bound"=bound,
                      "reason"="no usable solution returned for the bounded polynomial ansatz"])
    end if;
    V:=subs(sol,V);
    if normal(V)=0 then
        if verbose then printf("[GA][IIF] only the zero polynomial was returned\n") end if;
        return table(["status"="UNKNOWN_IIF","degree_bound"=bound,
                      "reason"="only the zero polynomial was returned within the chosen bound"])
    end if;
    if verbose then printf("[GA][IIF] candidate obtained; verify inverse-factor identity\n") end if;
    res:=simplify(normal(N*diff(V,x)+M*diff(V,y)-divergence*V));
    if res<>0 then
        if verbose then printf("[GA][IIF] residual could not be certified as zero\n") end if;
        return table(["status"="UNKNOWN_IIF","degree_bound"=bound,
                      "reason"="candidate did not pass symbolic residual verification"])
    end if;
    if verbose then printf("[GA][IIF] FOUND_IIF; residual=0, CPU=%.3f s\n",time()-started) end if;
    table(["status"="FOUND_IIF","result_type"="inverse_integrating_factor",
           "inverse_factor"=V,"residual"=0,"degree_bound"=bound,
           "vector_field"=[N,M],"method"="direct polynomial ansatz",
           "condition"="Choose integration constants so that V is nonzero"])
end proc;

InverseFactor := proc(M,N,x,y,{seconds::nonnegint:=0,verbose::truefalse:=false})
    if nops({x,y})<>2 or not type([x,y],list(name)) then
        error "x and y must be distinct unassigned names"
    end if;
    if N=0 or not type(M,polynom(rational,[x,y]))
           or not type(N,polynom(rational,[x,y])) then
        error "M,N must be polynomials over Q, with N nonzero"
    end if;
    try
        if seconds=0 then return InverseFactorCore(M,N,x,y,verbose)
        else return timelimit(seconds,InverseFactorCore(M,N,x,y,verbose)) end if
    catch "time expired":
        if verbose then printf("[GA][IIF] time limit reached after %d s; result is UNKNOWN_TIMEOUT\n",seconds) end if;
        return table(["status"="UNKNOWN_TIMEOUT","scope"="polynomial inverse integrating factor"])
    end try
end proc;

Reduce := proc(M,N,x,y,u,{seconds::nonnegint:=0,verbose::truefalse:=false,pruning::truefalse:=true})
    local h,q,R,show,cap,do_prune;
    CheckInput(M,N,x,y,u);
    show:=verbose; cap:=seconds; do_prune:=pruning;
    R:=normal(M/N);
    if verbose then printf("[GA][Reduce] start: check polynomial target\n") end if;
    # Preserve the immediate transformation for a polynomial target of degree >=3.
    # A quadratic polynomial input follows the direct IIF path below.
    if type(R,polynom(ratpoly(rational,x),y)) and degree(R,y)>=3 then
        if verbose then printf("[GA][Reduce] polynomial input: identity transformation, degree=%d\n",degree(R,y)) end if;
        return MakeResult(M,N,y,subs(y=u,R),x,y,u,"polynomial input")
    end if;
    if verbose then printf("[GA][Reduce] enter complete n>=3 search\n") end if;
    # Call the positional cores here: Maple evaluates a nested keyword name
    # against the caller's keyword local, which can silently restore defaults.
    try
        if cap=0 then h:=HigherCore(M,N,x,y,u,show,do_prune)
        else h:=timelimit(cap,HigherCore(M,N,x,y,u,show,do_prune)) end if
    catch "time expired":
        if show then printf("[GA][higher] time limit reached after %d s; result is UNKNOWN_TIMEOUT\n",cap) end if;
        h:=table(["status"="UNKNOWN_TIMEOUT","scope"="n>=3"])
    end try;
    if h["status"]="FOUND" then return h end if;
    if verbose then printf("[GA][Reduce] higher-degree status=%s; enter direct polynomial IIF search\n",h["status"]) end if;
    try
        if cap=0 then q:=InverseFactorCore(M,N,x,y,show)
        else q:=timelimit(cap,InverseFactorCore(M,N,x,y,show)) end if
    catch "time expired":
        if show then printf("[GA][IIF] time limit reached after %d s; result is UNKNOWN_TIMEOUT\n",cap) end if;
        q:=table(["status"="UNKNOWN_TIMEOUT","scope"="polynomial inverse integrating factor"])
    end try;
    if q["status"]="FOUND_IIF" then
        q["higher_status"]:=h["status"];
        return q
    end if;
    if verbose then printf("[GA][Reduce] final status=UNKNOWN; higher=%s, IIF=%s\n",
                           h["status"],q["status"]) end if;
    table(["status"="UNKNOWN","higher_status"=h["status"],"iif_status"=q["status"],
           "reason"="No polynomial IIF was certified within the template bound; this is not a proof against all n=2 reductions"])
end proc;

# Inverse integrating factor for the ORIGINAL polynomial field N*d_x+M*d_y.
# w,E are unassigned function names. No claim that the IIF is rational.
RiccatiInverseFactor := proc(M,N,r,P,x,y,u,w,E)
    local a,b,c,V,we,ee,res,vf;
    vf:=Verify(M,N,r,P,x,y,u);
    if vf["status"]<>"VERIFIED" or vf["n"]<>2 then
        error "a verified quadratic reduction is required"
    end if;
    a:=coeff(P,u,2); b:=coeff(P,u,1); c:=coeff(P,u,0);
    we:=diff(w(x),x$2)=(b+diff(a,x)/a)*diff(w(x),x)-a*c*w(x);
    ee:=diff(E(x),x)=-b*E(x);
    V:=normal(N/diff(r,y))*E(x)*(w(x)*r+diff(w(x),x)/a)^2;
    res:=normal(subs({we,ee},N*diff(V,x)+M*diff(V,y)
                           -(diff(N,x)+diff(M,y))*V));
    if res<>0 then error "inverse integrating factor verification failed" end if;
    table(["status"="VERIFIED","inverse_factor"=V,"w_ODE"=we,"E_ODE"=ee,
           "E_expression"=exp(-Int(b,x)),"residual"=res,
           "condition"="Choose a nonzero solution w and a nonzero solution E",
           "vector_field"=[N,M]])
end proc;

end module:

# Minimal usage:
# R := GeneralizedAbel:-Reduce(y^6+x,2*y,x,y,u):
# For a quadratic outcome: R["inverse_factor"]; R["residual"];
# For n>=3: R["transformation"]; R["equation"]; R["residual"];
# I := GeneralizedAbel:-InverseFactor(M,N,x,y,verbose=true):
# H := GeneralizedAbel:-ReduceHigher(M,N,x,y,u):
# H := GeneralizedAbel:-ReduceHigher(M,N,x,y,u,verbose=true):
# verbose=true prints live stages and branches; default false is silent.
# C := GeneralizedAbel:-Verify(M,N,known_r,known_P,x,y,u):
# IFAC := GeneralizedAbel:-RiccatiInverseFactor(M,N,R["transformation"],
#          R["polynomial"],x,y,u,w,E):  # for an independently given quadratic map