### INDICIAL EQUATION valpol:=proc( P :: depends(polynom(anything, t)), Q :: depends(polynom(anything, t)), t :: name, co ) local i,q,p,cont,valco; p:=P;q:=primpart(Q,t,'cont'); if normal(P)=0 then return infinity fi; if type(p,polynom(rational)) and type(q,polynom(rational)) then for i while divide(p,q,'p') do od else for i while rem(p,q,t,'p')=0 do od fi; if nargs=4 then valco:=rem(p,q,t)/cont^(i-1); # ASSERT(expand(rem(normal(P/Q^(i-1)),Q,t)-valco)=0); co:=valco fi; i-1 end: intsols := proc(pol, s :: name) local fpol; fpol := factor(pol); if fpol = 0 then return {s}; elif type(fpol, `*`) then fpol := select(p -> indets(p) = {s}, fpol); elif indets(fpol) <> {s} then fpol := 1; fi; return map(subs, {isolve(fpol)}, s); end proc: oldindicialeq0 := proc( L :: depends(polynom(anything, [Dvar, tvar])), pol :: depends(Or(identical(infinity), polynom(anything, tvar))), var, { Dvar :: name := NULL, tvar :: name := NULL }) local q, r, v, lc, shift, retshift, retlc; if L = 0 then return 1, -infinity; fi; r := rem(L, Dvar, Dvar, 'q'); v := valpol(r, pol, tvar); lc, shift := thisproc(q, pol, var+1, _options); retshift := max(shift+1, -v); retlc := 0; if -v = retshift then retlc := quo(r, pol^v, tvar); fi; if shift+1 = retshift then retlc := retlc + lc*(-var)*diff(pol, tvar); fi; retlc := rem(retlc, pol, tvar); return retlc, retshift; end proc: indicialeq0 := proc( L :: depends(polynom(anything, [Dvar, tvar])), pol :: depends(Or(identical(infinity), polynom(anything, tvar))), var, { Dvar :: name := NULL, tvar :: name := NULL }) local retshift, lvals, co, i, j, res, dq, pow, q, cont; if L = 0 then return 1, -infinity; fi; lvals:=[seq(i-valpol(coeff(L,Dvar,i),pol,tvar,co[i]),i=0..degree(L,Dvar))]; retshift:=max(op(lvals)); q:=primpart(pol,tvar,'cont'); res:=0;dq:=diff(q,tvar);pow:=1; for i from 0 to degree(L,Dvar) do if i>=1 then pow:=rem(pow*dq,q,tvar) fi; if lvals[i+1]=retshift then res:=res+cont^i*rem(pow*co[i],q,tvar)*mul(-var-j,j=0..i-1) fi od; # ASSERT(map(expand,[oldindicialeq0(args)])=map(expand,[res,retshift])); res,retshift end proc: indicialeq0inf := proc( L :: depends(polynom(anything, [Dvar, tvar])), var, { Dvar :: name := NULL, tvar :: name := NULL }) local q, r, v, lc, shift, retshift, retlc; if L = 0 then return 1, -infinity; fi; r := rem(L, Dvar, Dvar, 'q'); v := degree(r, tvar); lc, shift := thisproc(q, var-1, _options); retshift := max(shift-1, v); retlc := 0; if v = retshift then retlc := lcoeff(r, tvar); fi; if shift-1 = retshift then retlc := retlc + lc*var; fi; return retlc, retshift; end proc: # input: # L, diff. op. in Dvar and tvar # pol, a polynomial in tvar, or infinity (think pol=1/tvar) # var, a variable name # output: # - a set of integers # - the corresponding formal leading coeff Q(s, tvar), # - the shift # They satisfy # L(1/pol^s) = Q/pol^(s+shift) # and the set of integers are the values of s for which Q is not invertible # modulo pol. indicialeq1 := proc( L :: depends(polynom(anything, [Dvar, tvar])), pol :: depends(Or(identical(infinity), polynom(anything, tvar))), var :: name, { Dvar :: name := NULL, tvar :: name := NULL, dofactor :: boolean := false }) option remember; local lc, shift, indeq, integers; if pol = infinity then lc, shift := indicialeq0inf(L, var, _options); else lc, shift := indicialeq0(L, pol, var, _options); fi; if dofactor or pol = infinity then # we can assume that pol is irreducible indeq := foldl(gcd, coeffs(collect(lc, tvar), tvar)); else indeq := resultant(lc, pol, tvar); fi; integers := intsols(indeq, var); userinfo(3,'redct',sprintf("integral exponents at %.300a: % a",pol,integers)); userinfo(3,'redct',sprintf("shift there: % a",shift)); return integers, normal(lc), shift end proc: # Returns pol as a product of L-irreducible polynomials. splitpol := proc( pol :: depends(polynom(anything, tvar)), L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, dofactor :: boolean := true }) option remember; local sing, Q, Qev, _, p1, p2, sval; # La factorisation semble beaucoup plus rapide que splitpot... if dofactor then return factor(pol); fi; if type(pol, `*`) then return map(splitpol, pol, L, _options) elif type(pol,`^`) then return splitpol(op(1,pol),L,_options)^op(2,pol) elif degree(pol, tvar) <= 0 then return pol; fi; sing, Q, _ := indicialeq1(L, pol, indvar_, _options); for sval in sing do Qev := normal(subs(indvar_ = sval, Q)); if Qev <> 0 then ASSERT(degree(Qev, tvar) > 0); p1 := gcd(pol, numer(Qev)); p2 := normal(pol/p1); ASSERT(degree(p1, tvar) < degree(pol, tvar)); ASSERT(degree(p2, tvar) < degree(pol, tvar)); return thisproc(p1, L, _options)*thisproc(p2, L, _options); fi; end do; return pol; end proc: # Returns a right multiple of L such that the shift associated to any factor of # singp is null. zeroshift := proc( L :: depends(polynom(anything, [Dvar, tvar])), singp, { Dvar :: name := NULL, tvar :: name := NULL }) option remember; local s, f, Lden, Lnum, facts, shiftpol, shift, splsing, _; Lden := denom(normal(L)); Lnum := numer(normal(L)); # singp := Lden * lcoeff(collect(Lnum, Dvar), Dvar); splsing := splitpol(singp, Lnum, _options); _, facts := op(sqrfree(splsing, tvar)); shiftpol := 1; for f in facts do _, _, shift := indicialeq1(Lnum, f[1], s, _options); shiftpol := shiftpol*f[1]^(shift+valpol(Lden, f[1], tvar)); end do; return DETools['mult'](L, shiftpol, [Dvar,tvar]), shiftpol; end proc: # Returns L(A) applydop := proc( A, L :: depends(polynom(anything, [Dvar,tvar])), { Dvar :: name := NULL, tvar :: name := NULL }) local Lc, fn, ret, k; Lc := PolynomialTools['CoefficientVector'](L, Dvar); fn := A; ret := 0; for k from 1 to LinearAlgebra['Dimension'](Lc) do ret := ret + Lc[k]*fn; fn := diff(fn, tvar); od; return normal(ret); end: # If A = B + C*CERT, checks that B + L(C) = 0 checkcert := proc( A, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert :: name := 0 }) local lhs,rhs, res; lhs := subs(cert=0, A); rhs := coeff(A, cert); res := lhs + applydop(rhs, L, _options); return evalb(normal(res) = 0); end proc: rembezout := proc(P, Q, t) option remember; local A, B; gcdex(P, Q, t, A, B); return A, B; end proc: # return A such that # A*P*Q + B*R = S # for some B customgcdex := proc(P, Q, R, S, t) local U, V, ret; U, V := rembezout(P, R, t); ASSERT(normal(U*P+V*R-1)=0); ret := powmod(Q, -1, R, t); ret := rem(U*rem(rem(S, R, t)*ret, R, t), R, t); ASSERT(rem(ret*P*Q - S, R, t) = 0); return ret; end proc: # Input: # - A, a rational function # - denfact, a polynomial # - L, a differential operator # Output # U, V, C # where A = U + V + L(C), # U is the reduced part, V is the part to be reduced further and C is the certificate. localreduction_base := proc( A, denfact, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert :: anything := 0, absolute :: boolean := false }) local sing, lcgen, shift, den, exponent, lcev, redpart, redA, num, R, U, reductor, ccert, g, invden; sing, lcgen, shift := indicialeq1(L, denfact, indvar_, _options); if indvar_ in sing then error "Split %1 or shift the operator %2.", denfact, L; fi; den := denom(normal(A)); exponent := valpol(den, denfact, tvar); if exponent = 0 then return A, 0, 0; fi; lcev := normal(subs(indvar_ = exponent - shift, lcgen)); redA := normal(A*denfact^exponent); num := numer(redA); den := denom(redA); invden:=powmod(den,-1,denfact,tvar); ASSERT(degree(gcd(den, denfact), tvar) = 0); if exponent - shift in sing then if lcev <> 0 then error "Please split the denominator."; fi; if absolute then # Compute U such that A = diff(U/P,t\$(exponent-1)) + O(P^(-exponent + 1)) U := rem( num*invden*powmod(-diff(denfact, tvar), -exponent+1, denfact, tvar)/(exponent-1)!, denfact, tvar ); if exponent = 1 then redpart := U/denfact; else redpart := diff(U/denfact, tvar\$(exponent-1)); fi; else U := rem(invden*rem(num, denfact, tvar), denfact, tvar); redpart := U/denfact^exponent; fi; ASSERT(valpol(denom(normal(A - redpart)), denfact, tvar) < exponent); return normal(redpart), normal(A-redpart), 0; else # gcdex(den*lcev, denfact, num, tvar, 'R', 'U'); g:=gcdex(lcev, denfact, tvar, 'R'); ASSERT(g=1); # R:=rem(R*invden*num,denfact,tvar); R:=rem(R*rem(invden*rem(num,denfact,tvar),denfact,tvar),denfact,tvar); ASSERT(degree(R, tvar) < degree(denfact, tvar)); ccert := R/denfact^(exponent - shift); reductor := applydop(ccert, L, _options); ASSERT(valpol(denom(normal(A - reductor)), denfact, tvar) < exponent); return 0, normal(A-reductor), ccert; end if; end proc: localreduction := proc( A, denfact, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert :: Or(identical(0), name) := 0 }) local red, C, V, U, Cpart; red := 0; C := 0; V := normal(A); while V <> 0 do U, V, Cpart := localreduction_base(V, denfact, L, _options); red := red + U; C := C + Cpart; end do; return normal(red), C; end proc: # Performs the localreduction on the factors of the denominator. hermitered0 := proc( A, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert := 0, dofactor :: boolean := true }) local splpol, sqfden, _, red, C, f, Cpart; splpol := splitpol(denom(A), L, _options); _, sqfden := op(sqrfree(splpol)); # NB: sqrfree preserves factorization as much as possible sqfden := select(has, sqfden, tvar); red := A; C := 0; for f in sqfden do red, Cpart := localreduction(red, f[1], L, _options); C := C + Cpart; od; return red, C; end proc: # Hermite reduction at infinity. hermitered1 := proc( A, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert := 0 }) local numq, numr, An, sing, lcgen, shift, red, ccert, deg, lcnumq, lcev, reductor; An := normal(A); numr := rem(numer(An), denom(An), tvar, 'numq'); sing, lcgen, shift := indicialeq1(L, infinity, indvar_, _options); red := 0; ccert := 0; while numq <> 0 do deg := degree(numq, tvar); lcnumq := lcoeff(numq, tvar); if not deg - shift in sing and deg >= shift then lcev := subs(indvar_ = deg - shift, lcgen); reductor := applydop(lcnumq/lcev*tvar^(deg - shift), L, _options); ccert := ccert + lcnumq/lcev*tvar^(deg - shift); numq := normal(numq - reductor); else numq := normal(numq - lcnumq*tvar^deg); red := red + lcnumq*tvar^deg fi; ASSERT(degree(numq, tvar) < deg); end do; return numr/denom(An) + red, ccert; end proc: hermitered := proc( A, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert := 0 }) local red, C0, C1; red, C0 := hermitered0(A, L, _options); red, C1 := hermitered1(red, L, _options); return red + cert*(C0+C1); end proc: exceptionalbasis := proc( L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert :: Or(identical(0), name) := 0, liouvilleform :: boolean := false }) option remember; local lc, sqflc, basis, f, sing, Q, shift, Qev, g, cden, pbasis, new, num, den, co1, co2, newbasis, _, s, k, R, st; st:=time(); userinfo(2,'redct',"entering exceptionalbasis at time",time()); basis := {}; if not liouvilleform then lc := splitpol(lcoeff(normal(numer(L)), Dvar), L, _options); _, sqflc := op(sqrfree(lc)); sqflc := select(has, sqflc, tvar); for f in sqflc do sing, Q, shift := indicialeq1(L, f[1], indvar_, _options); for s in select(`>`, sing, 0) do Qev := subs(indvar_ = s, Q); g := normal(f[1]/gcd(f[1], Qev)); # Qev*h = 0 mod f iff h is divided by g. ASSERT(degree(g, tvar) < degree(f[1], tvar)); # Otherwise s is not really a singularity basis := basis union {seq(tvar^k*g/f[1]^s, k=0..degree(f[1], tvar)-degree(g, tvar)-1)}; end do; for s in {seq(s, s=1..-shift)} do Qev := subs(indvar_ = s, Q); g := normal(f[1]/gcd(f[1], Qev)); # Qev*h = 0 mod f iff h is divided by g. ASSERT(type(applydop(R, L, _options), polynom(anything, tvar))); basis := basis union {seq(tvar^k/f[1]^s, k=0..degree(f[1], tvar)-1)}; end do; end do; fi; sing, Q, shift := indicialeq1(L, infinity, indvar_, _options); basis := basis union {seq(tvar^s, s in select(`>=`, sing, 0))}; for R in basis do new:=normal(hermitered(applydop(R, L, _options), L, _options)); if new=0 then newbasis[R]:=NULL else num:=primpart(numer(new),tvar,'co1'); den:=primpart(denom(new),tvar,'co2'); newbasis[R]:=num/den-co2/co1*cert*R fi od; basis:={seq(newbasis[R],R=basis)}; # basis:={seq(hermitered(applydop(R, L, _options), L, _options)-cert*R,R=basis)}; cden := lcm(seq(denom( `if`(type(cert, name), normal(subs(cert=0, R)), R) ), R in basis)); pbasis := [seq( `if`(type(cert, name), collect(cden*R, cert, normal), normal(cden*R)) , R in basis)]; pbasis:=echelonpols(pbasis, _options); if pbasis<>[] then userinfo(2,'redct',"dimension of the exceptional set",nops(pbasis)); userinfo(2,'redct',"degrees of its elements",(map(degree,pbasis,tvar))); userinfo(2,'redct',"sizes of its elements",(map(length,pbasis))); fi; userinfo(2,'redct',"exceptional basis computed in ",time()-st,"sec."); return cden, pbasis end proc: echelonpols := proc( pols :: list, { tvar :: name := NULL, cert :: Or(identical(0), name) := 0 }) local npols, imax, echelon, rpol, spol, co; if nops(pols) <= 1 then return pols; fi; if type(cert, name) then npols := map(normal, subs(cert=0, pols)); else npols := map(normal, pols); fi; ASSERT(type(npols, list(polynom(anything, tvar)))); imax := max[index](map(degree, npols, tvar)); echelon := echelonpols(subsop(imax=NULL, pols), _options); rpol := redmodpols(pols[imax], echelon, _options); if type(cert,'name') then spol:=coeff(rpol,cert,0) else spol:=rpol fi; spol:=primpart(spol,tvar,'co'); if type(cert,'name') then spol:=spol+coeff(rpol,cert,1)/co fi; return remove(`=`, [spol, op(echelon)], 0); end proc: redmodpols := proc( p, pols :: list, { tvar :: name := NULL, cert :: Or(identical(0), name) := 0 }) local spol, L, degs, deg, j, ccert, certL, lc; L:=pols; spol := p; if type(cert, name) then ccert := coeff(spol,cert); spol:=subs(cert=0,spol); certL:=map(coeff,L,cert); L:=subs(cert=0,L) else ccert:=0 fi; degs:=map(degree,L,tvar); while L<>[] do j:=max['index'](degs); deg:=degs[j]; lc:=coeff(spol,tvar,deg); if lc<>0 then spol:=collect(spol-lc*L[j]/coeff(L[j],tvar,deg),tvar,normal); if type(cert,name) then ccert:=ccert-lc/coeff(L[j],tvar,deg)*certL[j] fi; fi; degs:=subsop(j=NULL,degs); L:=subsop(j=NULL,L); if type(cert,name) then certL:=subsop(j=NULL,certL) fi end do; return spol + cert*ccert; end proc: redmodfracs := proc( A, cden, basis, { tvar :: name := NULL, cert :: Or(identical(0), name) := 0 }) local numq, B, numr; B := normal(A*cden); numr := rem(numer(B), denom(B), tvar, 'numq'); numq := redmodpols(numq, basis, _options); return numq/cden + numr/cden/denom(B); end proc: reduction := proc( A, L :: depends(polynom(anything, [Dvar, tvar])), { Dvar :: name := NULL, tvar :: name := NULL, cert := 0 }) local cden, basis, red; cden, basis := exceptionalbasis(L, _options, _rest); red := hermitered(A, L, _options, _rest); red := redmodfracs(red, cden, basis, _options, _rest); return normal(red); end proc: