# $Source: /u/maple/research/lib/DEtools/diffop/src/RCS/rightdivision,v $ # $Notify: hoeij@sci.kun.nl $ # Procedures for users # rightdivision, leftdivision, GCRD, LCLM, adjoint, mult, l_p # Procedures for use only in the rest of the diffop package: # none # Procedures for use only in this file: # lift_rightdivision, lift_adjoint # (only for l_p): xplusa, lift_xplusa, lift_xplusa2, kx_lin_ser, # lift_linser, L_p # macro's for g_ext: macro( infosubs=`DEtools/diffop/infosubs`, modulus=`DEtools/diffop/modulus`, x=`DEtools/diffop/x`, DF=`DEtools/diffop/DF`, xDF=`DEtools/diffop/xDF` ): macro( g_conversion1=`algcurves/g_conversion1`, g_conversion2=`algcurves/g_conversion2`, rootof=`algcurves/rootof`, degree_ext=`algcurves/degree_ext`, g_evala=`algcurves/g_evala`, g_evala_rem=`algcurves/g_evala_rem`, g_zero_of=`algcurves/g_zero_of`, v_ext_m=`algcurves/v_ext_m`, ext_to_coeffs=`algcurves/e_to_coeff`, g_ext_r=`algcurves/g_ext_r`, g_gcdex=`algcurves/gcdex`, truncate=`algcurves/truncate`, g_factors=`algcurves/g_factors` ); macro( new_laurent2=`DEtools/diffop/new_laurent2`, differentiate=`DEtools/diffop/differentiate`, g_ext=`DEtools/diffop/g_ext`, g_factors=`DEtools/diffop/g_factors`, g_normal=`DEtools/diffop/g_normal`, set_laurents=`DEtools/diffop/set_laurents`, description_laurent=`DEtools/diffop/description_laurent`, modm=`DEtools/diffop/modm`, g_solve=`DEtools/diffop/g_solve`, g_quotient=`DEtools/diffop/g_quotient`, g_rem=`DEtools/diffop/g_rem`, g_quo=`DEtools/diffop/g_quo`, truncate_x=`DEtools/diffop/truncate_x`, g_expand=`DEtools/diffop/g_expand` ): # macro's for eval_laurent: macro( accuracy_laurent=`DEtools/diffop/accuracy_laurent`, value_laurent=`DEtools/diffop/value_laurent`, LL=_LL ): macro( lowerbound_val=`DEtools/diffop/lowerbound_val`, max0_val=`DEtools/diffop/max0_val`, new_laurent=`DEtools/diffop/new_laurent`, nmterms_laurent=`DEtools/diffop/nmterms_laurent`, nterms_expression=`DEtools/diffop/nterms_expression`, nthterm_laurent=`DEtools/diffop/nthterm_laurent`, ramification_laur=`DEtools/diffop/ramification_laur`, lift_ramification_laur=`DEtools/diffop/lift_ramification_laur`, eval_laurent=`DEtools/diffop/eval_laurent`, subsvalueslaurents=`DEtools/diffop/subsvalueslaurents`, series_val=`DEtools/diffop/pseries_val`, valuation_laurent=`DEtools/diffop/valuation_laurent`, nt=`DEtools/diffop/nt` ): # macro's for factor_op and substitute macro( Newtonpolygon=`DEtools/diffop/Newtonpolygon`, coefs_operator=`DEtools/diffop/coefs_operator`, factor_newton=`DEtools/diffop/factor_newton`, factor_newton2=`DEtools/diffop/factor_newton2`, factor_op=`DEtools/diffop/factor_op`, factor_riccati=`DEtools/diffop/factor_riccati`, lift_newton=`DEtools/diffop/lift_newton`, lift_ramification=`DEtools/diffop/lift_ramification`, lift_substitute=`DEtools/diffop/lift_substitute`, nm_mult=`DEtools/diffop/nm_mult`, nm_block=`DEtools/diffop/nm_block`, op_with_slope=`DEtools/diffop/op_with_slope`, ram_laur=`DEtools/diffop/ram_laur`, ramification_of=`DEtools/diffop/ramification_of`, rem_lift=`DEtools/diffop/rem_lift`, skipped_factors=`DEtools/diffop/skipped_factors`, substitute=`DEtools/diffop/substitute` ): macro( rightdivision=`DEtools/diffop/rightdivision`, leftdivision=`DEtools/diffop/leftdivision`, GCRD=`DEtools/diffop/GCRD`, LCLM=`DEtools/diffop/LCLM`, adjoint=`DEtools/diffop/adjoint`, lift_adjoint=`DEtools/diffop/lift_adjoint` ): macro( lift_rightdivision=`DEtools/diffop/lift_rightdivision`, mult=readlib(`DEtools/diffop/mult`), solve=readlib(`solve/linear`) ): # result: [a,b] such that a*right+b=f # right must be monic # If the 4th argument slope is specified we get a faster division, but no # remainder. We can use this 4'th argument only if this slope is the # only slope. rightdivision:=proc(f,right,ext,slope) local a,b,t; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if nargs=4 then # `DEtools/diffop/readlibs`(factor_op): # division of local differential operators a:=op_with_slope(degree(f,xDF)-degree(right,xDF),slope,0, [lift_rightdivision,f,right,ext]); # rem_lift[a]:=0; I used this variable for computation of the remainder # but that introduced a bug. I didn't need the remainder anyway for # this case so now the remainder is not implemented. RETURN([a,`not implemented`]) elif nargs=2 then a:=collect(right/lcoeff(right,DF),DF,g_normal); t:=g_ext([args]); a:=rightdivision(op(subs(g_conversion1,[f,a])),t); RETURN(subs(g_conversion2,[mult(a[1],1/lcoeff(right,DF)),a[2]])) fi; a:=0; b:=f; if member(DF,indets([args])) then # global factor if lcoeff(right,DF)<>1 then ERROR(right,`must be monic`) fi; if right=1 then RETURN([f,0]) fi; userinfo(5,'diffop',`right division`,infosubs(f,right)); while degree(b,DF)>=degree(right,DF) do a:=a+lcoeff(b,DF)*DF^(degree(b,DF)-degree(right,DF)); b:=collect(b-lcoeff(b,DF)*mult(seq(DF ,t=1..degree(b,DF)-degree(right,DF)),right,ext),DF,g_normal) od; userinfo(5,'diffop',`done right division`) else if lcoeff(right,xDF)<>1 then ERROR(right,`must be monic`) fi; if right=1 then RETURN([f,0]) fi; while degree(b,xDF)>=degree(right,xDF) do a:=a+lcoeff(b,xDF)*xDF^(degree(b,xDF)-degree(right,xDF)); b:=eval_laurent(expand(b-lcoeff(b,xDF)*mult(seq(xDF ,t=1..degree(b,xDF)-degree(right,xDF)),right,ext)),ext,xDF) od fi; [a,b] end: lift_rightdivision:=proc(llaur,order_l,slope,f,right,ext,left,acc) global value_laurent, accuracy_laurent; local lau,l,r,lr,n_lifts,n_known,mult_args,i,l_extra,le,fe,r_low; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; n_lifts:=acc-accuracy_laurent[llaur]; if n_lifts<=0 then RETURN() fi; n_known:=accuracy_laurent[coeff(left,xDF,0)]+degree(left,xDF)*slope; if n_known=0 then l:=0 else l:=expand(subsvalueslaurents(left)) fi; r:=nm_block(right,0,n_known+n_lifts,slope,ceil); r_low:=coefs_operator(r,slope,-numer(slope)*degree(right,xDF),0); mult_args:=ceil(n_known-slope*degree(f,xDF)), ceil(n_known+n_lifts)-1,ext; # lr:=rem_lift[left]-nm_block(f,n_known,n_known+n_lifts,slope,ceil) # +nm_mult(l,r,mult_args[2]+1-n_lifts,mult_args[2..3]); lr:=-nm_block(f,n_known,n_known+n_lifts,slope,ceil) +nm_mult(l,r,mult_args); le:=-numer(slope)*degree(left,xDF); fe:=-numer(slope)*degree(f,xDF); for i from n_known*denom(slope) to (n_known+n_lifts)*denom(slope)-1 do l_extra:=coefs_operator(g_quo(-coefs_operator(lr ,slope,i+fe,0),r_low,x,ext),slope,i+le,1); l:=l+l_extra; lr:=lr+nm_mult(l_extra,r,mult_args) od; # rem_lift[left]:=lr; l:=collect(l,xDF); for i from 0 to degree(left,xDF)-1 do lau:=coeff(left,xDF,i); value_laurent[lau]:=coeff(l,xDF,i); accuracy_laurent[lau]:=accuracy_laurent[lau]+n_lifts od; NULL end: leftdivision:=proc(f,L,ext) local a,b; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if nargs=2 then a:=mult(L,1/lcoeff(L,DF)); b:=g_ext([args]); a:=procname(op(subs(g_conversion1,[f,a])),b); RETURN(subs(g_conversion2,[a[1]/lcoeff(L,DF),a[2]])) fi; userinfo(5,'diffop',`left division`,infosubs(f,L)); a:=0; b:=f; while degree(b,DF)>=degree(L,DF) do a:=a+lcoeff(b,DF)*DF^(degree(b,DF)-degree(L,DF)); b:=collect(b-mult(L,lcoeff(b,DF),ext)*DF^ (degree(b,DF)-degree(L,DF)),DF,g_normal) od; userinfo(5,'diffop',`done left division`); [a,b] end: GCRD:=proc(f,g) global DF; local a; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if f=0 or g=0 then f+g else a:=collect(g/lcoeff(g,DF),DF,g_normal); procname(a,readlib(rightdivision)(f,a)[2]) fi end: # Input: sequence of operators # or: # sequence of operators,`and conjugates`,ext # ext: optional, if not present then ext=[] # Output: the LCLM of the operators and their conjugates over Q(ext) # RootOf syntax LCLM:=proc(L) local ext,j,i,ind,sol,ansatz,a; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; ext:=g_ext([args]); if has([args],`and conjugates`) then if not type(args[nargs],list) then # ground field is Q RETURN(procname(args,[])) fi; if nargs>3 then RETURN(procname(seq(procname(args[j],args[nargs-1..nargs]),j=1..nargs-2))) fi; j:=subs(g_conversion1,args[3]); ansatz:=convert([seq(a[i]*DF^i,i=0..degree(L,DF) *degree_ext(ext,j))],`+`); sol:={seq(coeffs(i,DF),i=[ext_to_coeffs(g_expand(subs(g_conversion1, numer(readlib(rightdivision)(ansatz,L)[2])),ext),j)])} else if type(args[nargs],list) then # ext should not be given RETURN(procname(args[1..nargs-1])) elif nargs=1 then RETURN(collect(L/lcoeff(L,DF),DF,g_normal)) fi; ansatz:=convert([seq(a[i]*DF^i,i=0..convert([seq(degree(j,DF),j=args)] ,`+`))],`+`); sol:={seq(coeffs(g_expand(numer(readlib(rightdivision)( ansatz,subs(g_conversion1,i))[2]),ext),DF),i=args)} fi; ind:=indets(ansatz) minus {DF}; sol:=solve(subs(g_conversion2,sol),ind); i:=degree(ansatz,DF); while nops(indets(subs(sol,ansatz)) intersect ind)>1 do sol:=solve({a[i]=0,op(sol)},ind); i:=i-1 od; sol:=solve({a[i]=1,op(sol)},ind); collect(subs(sol,ansatz),DF,g_normal) end: # Compute the adjoint, both local and global case. adjoint:=proc(f) global x,xDF,DF; local i,ext,v,D,g; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; userinfo(5,'diffop',`computing the adjoint of`,infosubs(f)); if nargs=1 then ext:=g_ext(f) else ext:=args[2] fi; if has(f,xDF) and has(f,set_laurents) then v:=lowerbound_val(f); g:= convert([xDF^degree(f,xDF),seq(new_laurent(v,0, [lift_adjoint,i,f,ext])*xDF^i,i=0..degree(f,xDF)-1)],`+`) else if has(f,xDF) then D:=xDF else D:=DF fi; v:=convert([seq(mult(D^i*(-1)^(i+degree(f,D)),coeff(f,D,i) ,ext),i=0..degree(f,D))],`+`); if type(v,polynom(anything,x)) then g:= g_expand(v,ext) else g:= collect(v,D,g_normal) fi fi; procname(g,args[2..nargs]):=f; g end: # L is the n'th coefficient in the adjoint of f. # ext seems to be unused. This is to make sure that ext appears in the # descriptions of the Laurent series. lift_adjoint:=proc(L,n,f,ext,acc) global x,xDF; local a,res,i; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; a:=acc,accuracy_laurent[L]; res:=0; for i from degree(f,xDF)-1 by -1 to n do res:=expand((-1)^(i+degree(f,xDF))*nmterms_laurent( coeff(f,xDF,i),a)+x*diff(res,x)*(i+1)/(i+1-n)) od; value_laurent[L]+res end: #savelib('rightdivision','lift_rightdivision','`DEtools/diffop/rightdivision.m`'): #savelib('leftdivision','`DEtools/diffop/leftdivision.m`'): #savelib('GCRD','`DEtools/diffop/GCRD.m`'): #savelib('LCLM','`DEtools/diffop/LCLM.m`'): #savelib('adjoint','lift_adjoint','`DEtools/diffop/adjoint.m`'): macro( mult=`DEtools/diffop/mult` ): # Multiplication works for both syntaxes xDF and DF, but these syntaxes # must not be mixed. # Multiplication of local operators is not optimized. Could be done faster # using a liftable result, lifted by nm_mult. mult:=proc() global x,xDF,DF; local a,b,i,Dib,result,D,xx,ext; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if not type(args[nargs],list) then RETURN(mult(args,g_ext([args]))) fi; if has([args],RootOf) then RETURN(subs(g_conversion2,mult(op(subs(g_conversion1,[args]))))) fi; if nargs=2 then RETURN(args[1]) fi; ext:=args[nargs]; if nargs>3 then RETURN(mult(args[1],mult(args[2..nargs]),ext)) fi; if member(xDF,indets([args])) then xx:=x; D:=xDF else xx:=1; D:=DF fi; if type(args[1],polynom(anything,x)) then a:=g_expand(args[1],ext) else a:=collect(args[1],D,g_normal) fi; b:=args[2]; result:=0; for i from 0 to degree(a,D) do if i=0 then Dib:=b else Dib:=xx*differentiate(Dib)+Dib*D; if has(Dib,set_laurents) then Dib:=eval_laurent(expand(Dib),ext,D) elif type(Dib,polynom(anything,x)) then Dib:=g_expand(Dib,ext) else Dib:=collect(Dib,D,g_normal) fi fi; result:=result+coeff(a,D,i)*Dib od; if has(result,set_laurents) then result:=eval_laurent(expand(result),ext,D) elif type(result,polynom(anything,x)) or D=xDF then result:=g_expand(result,ext) else result:=collect(result,D,g_normal) fi; result end: #savelib('mult','`DEtools/diffop/mult.m`'): macro( substitute=readlib(`DEtools/diffop/substitute`), l_p=`DEtools/diffop/l_p`, L_p=`DEtools/diffop/L_p`, xplusa=`DEtools/diffop/xplusa`, lift_xplusa=`DEtools/diffop/lift_xplusa`, lift_xplusa2=`DEtools/diffop/lift_xplusa2`, kx_lin_ser=`DEtools/diffop/kx_lin_ser`, lift_linser=`DEtools/diffop/lift_linser` ): ############################################ # Localization of differential operators #------------ ############################################ # Input: f in k(x)[DF], p in P^1(k_bar). In RootOf or rootof notation. # Output: l_p(f) in k(p)((x))[xDF], i.e. f localized in p and monic. l_p:=proc() global `DEtools/diffop/modulus`; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if type(`DEtools/diffop/modulus`,posint) then # turn modular computations off (the fact that this happens # indicates that the computation was interrupted. # The following code has already been loaded because that file # is the only spot where the modulus is set to something nonzero. # I placed this in a seperate procedure in order to have an # options remember in the rest of the l_p code, in L_p. `DEtools/diffop/compute_modp`() fi; L_p(args) end: L_p:=proc(f,p) global DF,x,xDF; local g,ext; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if type(p,equation) then ERROR(`wrong input`) fi; ext:=g_ext([args]); if args[nargs]<>ext then # this is to make the options remember in l_p effective in # more cases: RETURN(procname(op(subs(g_conversion1,[args])),ext)) fi; g:=g_expand(subs(g_conversion1,numer(f)),ext); if nargs>2 and args[3]=`global` and p<>infinity then g:=subs(g_conversion2,g_expand(subs(x=x+p,g),ext)); RETURN(collect(g/lcoeff(g,DF),DF,g_normal)) fi; if p=infinity then g:=substitute(-x*xDF,expand(numer(subs({DF=xDF,x=1/x},g))),ext); if nargs>2 and args[3]=`global` then RETURN(collect(subs(g_conversion2,g/lcoeff(g,DF)),DF,g_normal)) fi; if not has(subs(g_conversion2,lcoeff(g,xDF)),RootOf) then g:=collect(g/lcoeff(g,xDF),xDF,g_normal) fi; # if lcoeff does contain RootOf's then such a division could # be quite costly so in that case we first apply eval_laurent: g:=eval_laurent(g,ext,xDF) elif p=0 then g:=collect(g,DF,eval_laurent); g:=expand(numer(substitute(1/x*xDF,subs(DF=xDF,g),ext))); g:=kx_lin_ser(g,ext) else g:=xplusa(g,p,ext); g:=expand(numer(substitute(1/x*xDF,subs(DF=xDF,g),ext))); g:=kx_lin_ser(g,ext) fi; eval_laurent(collect(g/lcoeff(g,xDF),xDF),ext,xDF) end: # Input: P in k[x] # a in k_bar # Output: a series representing subs(x=x+a,P) xplusa:=proc(P,a,ext) global DF,x; local i; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if has(P,DF) then convert([seq(procname(coeff(P,DF,i),a,ext)*DF^i,i=0..degree(P,DF))],`+`) elif member(a,{op(ext)} minus {op(g_ext(P))}) then new_laurent(0,0,[lift_xplusa2,P,a, g_expand(subs(_Z=x,subs(g_conversion1,op(subs(g_conversion2,a)))),ext) ,ext]) else new_laurent(0,0,[lift_xplusa,P,a,ext]) fi end: lift_xplusa:=proc(L,P,a,ext,n) global x; local m,i,j; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; m:=accuracy_laurent[L]; g_expand(value_laurent[L]+convert([ seq(seq(coeff(P,x,i)*binomial(i,j)*x^j*g_evala_rem(a^(i-j)) ,j=m..min(n-1,i)),i=m..degree(P,x))],`+`),ext,x) end: lift_xplusa2:=proc(L,P,a,pa,ext,n) global description_laurent,x; local p,v,m,i; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; m:=accuracy_laurent[L]; v:=0; p:=P; for i from m to n-1 do # if not divide(p,pa) then v:=v+subs(x=a,rem(p,pa,x))*x^i; # fi; p:=collect(diff(p,x),x)/(i+1) od; description_laurent[L]:=[lift_xplusa2,L,p,a,pa,ext]; g_expand(value_laurent[L]+v,ext,x) end: # k[x] linear combination of laurent series kx_lin_ser:=proc(a,ext) global xDF; local v; options remember, `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; if has(a,xDF) then RETURN(convert([seq(procname(coeff(a,xDF,v),ext)*xDF^v,v =0..degree(a,xDF))],`+`)) fi; if member(a,set_laurents) then RETURN(a) fi; new_laurent(0,nt(a,0), [lift_linser,a,[op(indets(a) intersect set_laurents)],ext]) end: lift_linser:=proc(L,a,v,ext,n) global x; local i,s,m,c; option `Copyright (c) 1997 Waterloo Maple Inc. All rights reserved.`; s:=0; m:=accuracy_laurent[L]; for i in v do c:=coeff(a,i,1); s:=s+nmterms_laurent(i,n-ldegree(c,x),m-degree(c,x))*c od; s:=g_expand(s,ext,x); value_laurent[L]+convert([seq(coeff(s,x,i)*x^i,i=m..n-1)],`+`) end: #savelib('L_p','l_p','xplusa','lift_xplusa','lift_xplusa2','kx_lin_ser',\ 'lift_linser','`DEtools/diffop/l_p.m`'):