# Author: Mark van Hoeij # University of Nijmegen, Netherlands # e-mail: hoeij@sci.kun.nl # If you change/improve this program, or if you find a bug, please let me know # diffop package: This package provides factorizations of local (i.e. power # series coefficients) and global (i.e. rational functions coefficients) # linear differential operators. It also provides formal solutions for systems # of linear differential equations (matrix differential equations). # lprint(`Note: preliminary version, the global factorization is currently`); # lprint(`being rewritten.`); lprint(`Additional documentation and the most recent version of this program is`); lprint(`available via WWW from: http://www-math.sci.kun.nl/math/compalg/diffop/`); x:=evaln(x): # We use these 3 variables as global variables. xDF:=evaln(xDF): # xDF = x*d/dx DF:=evaln(DF): # DF = d/dx # The following are used globally. Uncomment these lines to change that. # macro( # LL=`diffop/LL`, # nt=`diffop/nt`, # rootof=`diffop/rootof` # ); infolevel[diffop]:=1: ############################# # groundfield computation #------------------------------------------------- ############################# diffop_groundfield # Note: the procedures in this section are copied from the file IntBasis. # For brief descriptions see IntBasis. # The main purpose is to have ground field independent procedures like # g_expand (groundfield dependent expand) which work at a reasonable # efficiency. The field is always given by the argument ext. macro( bug=`diffop/bug`, degree_ext=`diffop/degree_ext`, ext_to_coeffs=`diffop/ext_to_coeffs`, g_conversion1=`diffop/g_conversion1`, g_conversion2=`diffop/g_conversion2`, g_evala=`diffop/g_evala`, g_evala_rem=`diffop/g_evala_rem`, g_expand=`diffop/g_expand`, normal_coeffs=`diffop/normal_coeffs`, g_ext=`diffop/g_ext`, g_ext_r=`diffop/g_ext_r`, g_factors=`diffop/g_factors`, g_gcdex=`diffop/g_gcdex`, g_normal=`diffop/g_normal`, g_quo=`diffop/g_quo`, g_quotient=`diffop/g_quotient`, g_rem=`diffop/g_rem`, g_solve=`diffop/g_solve`, g_zero_of=`diffop/g_zero_of`, modm=`diffop/modm`, modulus=`diffop/modulus`, truncate=`diffop/truncate`, truncate_x=`diffop/truncate_x`, v_ext_m=`diffop/v_ext_m`, factors=readlib('factors') ): g_conversion1:={}: # RootOf syntax -> my own syntax g_conversion2:={}: # my syntax -> RootOf syntax modulus:=0: # For modular computation modm:=proc() if modulus=0 then args else args mod modulus fi end: # modm:=proc() args end: g_solve:=proc() local v,w,res,i,a,ext; if nargs<3 then # The "if" is to avoid that this message gets repeated in the recursion userinfo(10,diffop,`solving linear equations`) fi; if modulus=0 then subs(g_conversion1,evala(solve(op(subs(g_conversion2,[args[1 ..min(2,nargs)]]))))) else # Solving linear equations modulo a number, not very efficiently written. v:=args[1]; if nargs=3 then ext:=args[3] else ext:=g_ext(v) fi; if nargs>=2 then w:=args[2] else w:=indets(v,`name`) minus {op(ext)} fi; v:=modm(map(g_expand,v,ext)) minus {0}; if v={} then RETURN(v) fi; a:=v[1]; i:=indets(a,`name`) intersect w; if i=[] then RETURN(NULL) fi; i:=i[1]; res:=g_normal(-coeff(a,i,0)/coeff(a,i,1)); a:=g_solve(subs(i=res,v minus {a}),w minus {i},ext); if a=NULL then a else a union {i=modm(g_expand(subs(a,res),ext))} fi fi end: # Compute the quotient of two truncated power series. g_quotient:=proc(aa,bb,n,ext) local a,b,res,i,g,la,lb; a:=modm(g_expand(aa,ext,x)); b:=modm(g_expand(bb,ext,x)); if a=0 or b=0 then RETURN(0) fi; g:=g_normal(1/tcoeff(b,x)); res:=0; la:=ldegree(a,x); lb:=ldegree(b,x); for i from la-lb to n-1 do res:=res+g*x^i*coeff(a,x,i+lb); a:=a-g_expand(b*g*x^i*coeff(a,x,i+lb),ext) od; modm(g_expand(res,ext,x)) end: # Same as in the algcurves package: g_gcdex:=proc(a,b,c,x,ext) global g_conversion1,g_conversion2; local r,ss,tt,s,t,i,j; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; r:=subs(g_conversion2,[a,b,c]); for i from 1 to 3 do if indets(denom(r[i]))<>{} then # Clear the denominator: r:=[seq(evala(Expand(g_normal(j*denom(r[i])))),j=r)] fi od; if ext=[] then r:=gcdex(op(r),x,'ss','tt') else r:=subs(g_conversion1,evala(Gcdex(r[1],r[2],x,'ss','tt'))) fi; s:=subs(g_conversion1,ss); t:=subs(g_conversion1,tt); [c,s,t] end: g_rem:=proc(a,b,y,ext) g_expand(rem(g_expand(a,ext),b,y),ext,y) end: g_quo:=proc(a,b,y,ext) modm(g_expand(quo(g_expand(a,ext),modm(g_expand(b,ext)),y),ext,y)) end: # Same as in the algcurves package, but with the difference that the factors are made monic. # Factorization over the groundfield # Both RootOf and rootof syntax # ext: describes the groundfield g_factors:=proc(f,x,ext) global g_conversion1,g_conversion2; local v,w,i; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if not has(f,x) then RETURN([]) elif ext<>[] and not has(ext,RootOf) then # rootof syntax input, so rootof syntax output RETURN(subs(g_conversion1, procname(op(subs(g_conversion2,[args]))))) fi; userinfo(5,diffop,` Factorizing`,f); v:=traperror(factors(numer(f),{op(ext)})); if v=lasterror then # Now try evala(Factor()), it can handle more # cases than factors g_ext([args]); v:=subs(g_conversion2,factors(subs(g_conversion1 ,evala(Factor(numer(f),{op(ext)}))))) fi; v:=v[2]; userinfo(5,diffop,`Done factorization`); w:=NULL; for i in v do if has(i,x) then w:=w,[collect(i[1]/lcoeff(i[1],x),x,g_normal),i[2]] fi od; [w] end: # Same as in the algcurves package: # Gives the algebraic degree of a over b degree_ext:=proc(aa,bb) global g_conversion2; local a,b,v,i,all,d,var; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; a:=subs(g_conversion2,aa); b:=subs(g_conversion2,bb); v:=indets(a,RootOf) minus indets(b,RootOf); all:=[op(indets([a,b],RootOf))]; all:={seq(all[i]=var[i],i=1..nops(all))}; d:=1; for i in v do d:=d*degree(subs(all,op(i)),_Z) od; d end: # Same as in the algcurves package, with normal_tcoeff replaced by normal_coeffs # evala(Expand( )) # Give 3'rd argument x for normalizing the coefficients in x # So then ldegree(result,x) gives the right answer g_expand:=proc(a,ext) global g_conversion1,g_conversion2; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if type(a,polynom) then g_evala(expand(a),ext) elif nargs=3 and type(a,polynom(anything,ext)) then normal_coeffs(g_evala(expand(a),ext),args[3]) elif nargs=3 then normal_coeffs(procname(a,ext),args[3]) else subs(g_conversion1,evala(Expand(subs(g_conversion2,a)) ,'independent')) fi end: # Normalize all coefficients normal_coeffs:=proc(a,x) if type(a,polynom) then a else collect(a,x,g_normal) fi end: # Same as in the algcurves package: # Input: A polynomial in the rootofs in ext # Output: a evala normalization g_evala:=proc(a,ext) local i,e; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if nops(ext)=0 then RETURN(a) elif nops(ext)=1 then e:=ext[1]; expand(convert([seq(coeff(a,e,i)*g_evala_rem(e^i) ,i=0..degree(a,e))],`+`)) else e:=g_evala(a,[op(2..nops(ext),ext)]); g_evala(expand(convert([seq(coeff(e,ext[1],i)*g_evala_rem(ext[1]^i) ,i=0..degree(e,ext[1]))],`+`)),[op(2..nops(ext),ext)]) fi end: # Same as in the algcurves package: g_evala_rem:=proc() global g_conversion1,g_conversion2; options remember, `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; expand(subs(g_conversion1,evala(Expand(subs(g_conversion2,args)) ,'independent'))) end: g_normal:=proc(aa) local a; if indets(aa,RootOf)<>{} then RETURN(evala(Normal(aa))) fi; if modulus<>0 then a:=subs(g_conversion2,Normal(aa) mod modulus); if has(a,RootOf) then subs(g_conversion1,evala(Normal(a))) else a fi else a:=subs(g_conversion2,aa); if has(a,RootOf) then subs(g_conversion1,evala(Normal(a))) else normal(a) fi fi end: # Same as in the algcurves package: # Input : an irreducible polynomial kk in x, not necessarily monic # Output: a root of k # If an algebraic extension is needed it will be placed in ext. g_zero_of:=proc(k,x,ext) global g_conversion1,g_conversion2; local a,v; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; if degree(k,x)=1 then ext:=NULL; RETURN(g_normal(-coeff(k,x,0)/coeff(k,x,1))) fi; a:=RootOf(subs(g_conversion2,k),x); if not member(_Z,indets(op(subs(g_conversion1,a)))) then a:=subs(g_conversion1,a); ext:=a; RETURN(a) fi; if g_conversion1={} then g_conversion1:=NULL fi; v:=nops(g_conversion2); g_conversion1:=a=rootof.v,g_conversion1; g_conversion2:={rootof.v=a,op(g_conversion2)}; ext:=rootof.v end: # Same as in the algcurves package: # Gives the zeros of the factors, their multiplicities and algebraic extensions # The input of this procedure is the output of g_factors v_ext_m:=proc(f,x) local ext,nulp,i,res; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; res:=NULL; for i in f do nulp:=g_zero_of(i[1],x,'ext'); ext:=eval(ext); res:=res,[nulp,i[2],[ext],degree(i[1],x)] od; [res] end: # Same as in the algcurves package: # ext_to_coeffs does basically: coeffs(expression,RootOf( .. )) # We need this procedure because coeffs doesn't always work this way. ext_to_coeffs:=proc(a,ext) global g_conversion2; local i,aa; option `Copyright (c) 1996 Waterloo Maple Inc. All rights reserved.`; aa:=(indets(a) minus indets(ext)) intersect {seq(rootof.i,i=0..nops(g_conversion2))}; coeffs(a,[op(aa)]) end: # Takes the lowest coefficients truncate:=proc(aa,n,y,ext) local dummy,a; a:=collect(aa,y); a:=expand(convert([seq(y^dummy*coeff(a,y,dummy) ,dummy=ldegree(a,y)..n-1)],`+`)); modm(g_evala(a,ext)) end: truncate_x:=proc(a,f,n) local dummy; convert([seq(x^dummy*coeff(f,x,dummy),dummy=ldegree(f,x)..n-1)],`+`) end: # a very useful procedure! bug:=proc() lprint(`Bug alert: please send this example to hoeij@sci.kun.nl`); ERROR(args) end: ################################ # computation in k((x)) #---------------------------------------------- ################################ diffop_laurent # where k is the groundfield, an algebraic extension of Q # variables: # LL.1 LL.2 ... these stand for the Laurent series in x # with finite pole order. # accuracy_laurent value_laurent description_laurent: these are # tables `help/text/eval_laurent`:=TEXT(``, `FUNCTION: eval_laurent - Convert an expression to a Laurent series`,``, `CALLING SEQUENCE:`,` eval_laurent(e);`,``,`PARAMETERS:`, ` e - an expression that can be converted to a Laurent series in x`,``, `SYNOPSIS:`,`- A Laurent series in this package is a series of the form`, ``,` infinity`, ` -----`, ` \\ i`, ` ) a[i] x`, ` /`,` -----`, ` i = - N`,``, ` for some integer N and coefficients a[i]. Currently the coefficients a[i]`, ` are restricted in diffop to algebraic numbers, although in the future`, ` other coefficient fields will be supported as well.`,``, ` A Laurent series is denoted in diffop as a variable LL.m where m is some`, ` integer. Since such a series can have infinitely many coefficients the`, ` value can be given only up to some finite accuracy a (i.e. modulo x^a,`,` \ so all terms up to x^(a-1)). This is done with the procedure nt (= n terms).`, ``,`EXAMPLES:`,`> l:=eval_laurent( (x^2+x+1)/(x^5+x^17) );`,``, ` l := LL3`,``,`> nt(l,8); `,``, ` 1 1 1 7`, ` ---- + ---- + ---- - x`, ` 5 4 3`, ` x x x`,`> nt(l,10);`,``, ` 1 1 1 7 8 9`, ` ---- + ---- + ---- - x - x - x`, ` 5 4 3`, ` x x x`,``, `> l:=eval_laurent( (l^3*x^10+4)/sin(x) );`,``, ` l := LL12`,``,`> nt(l,5); `,``,` 1\ 3 37 15 2527 329 6407 733 25561 2 91381 3`, `---- + ---- + ---- + ---- + ------ + ---- + ---- + --- x + ----- x + ------ \ x`, ` 6 5 4 3 2 40 x 3024 560 86400 604800`, ` x x 6 x 2 x 360 x`,``,` 3955031 4`, ` + --------- x`,` 119750400`): `help/text/nt`:=TEXT(``, `FUNCTION: nt - n terms of an expression in Laurent series`,``, `CALLING SEQUENCE:`,` nt(e,n);`,``,`PARAMETERS:`, ` e - an expression`, ` n - an integer`,``, `SYNOPSIS:`,`- This procedure evaluates all Laurent series in the expression`, ` e up to accuracy n (i.e. mod x^n), substitutes these values in e and`, ` returns e mod x^n.`,``,`EXAMPLES: For examples see ?eval_laurent`): # The purpose of this section is computation with infinite series. # This section uses the following variables macro( accuracy_laurent=`diffop/accuracy_laurent`, description_laurent=`diffop/description_laurent`, set_laurents=`diffop/set_laurents`, value_laurent=`diffop/value_laurent` ): set_laurents:={}: # It uses the global variable x and the procedures: macro( differentiate=`diffop/differentiate`, lowerbound_val=`diffop/lowerbound_val`, max0_val=`diffop/max0_val`, new_laurent=`diffop/new_laurent`, new_laurent2=`diffop/new_laurent2`, nmterms_laurent=`diffop/nmterms_laurent`, nterms_expression=`diffop/nterms_expression`, nthterm_laurent=`diffop/nthterm_laurent`, ramification_laur=`diffop/ramification_laur`, lift_ramification_laur=`diffop/lift_ramification_laur`, subsvalueslaurents=`diffop/subsvalueslaurents`, series_val=`diffop/pseries_val`, valuation_laurent=`diffop/valuation_laurent` ): # g_ext: gives a list of the algebraic extensions. g_ext_r:=proc(a) local v,vv,i,tail; options remember; v:=indets(a,RootOf); if nops(v)=0 then RETURN([]) fi; vv:={}; for i in v do vv:=vv union indets(op(i),RootOf) od; tail:=g_ext_r(vv); v:=[op(v minus vv)]; [op(v),op(tail)] end: # Gives the algebraic extensions appearing in aa. # For laurent series, the groundfield must be given in their descriptions. g_ext:=proc(aa) global g_conversion1, g_conversion2; local v,i,result,ii,vv; options remember; v:=aa; if indets(aa) intersect set_laurents<>{} then for i in indets(aa) intersect set_laurents do v:=subs(i=description_laurent[i],v) od; fi; v:=g_ext_r(subs(g_conversion2,v)); vv:=subs(g_conversion1,v); result:=NULL; for i from 0 to nops(g_conversion2) do if member(rootof.i,vv) then result:=rootof.i,result fi od; for i from nops(v) by -1 to 1 do if not member(subs(g_conversion1,v[i]) ,{seq(rootof.ii,ii=0..nops(g_conversion2))}) then if g_conversion1={} then g_conversion1:=NULL fi; vv:=nops(g_conversion2); g_conversion1:=v[i]=rootof.vv,g_conversion1; g_conversion2:={rootof.vv=v[i],op(g_conversion2)}; result:=subs(g_conversion1,v[i]),result fi od; [result] end: new_laurent:=proc() global set_laurents,accuracy_laurent,value_laurent,description_laurent; local i; i:=nops(set_laurents)+1; set_laurents:=set_laurents union {LL.i}; if nargs=3 then accuracy_laurent[LL.i]:=args[1]; value_laurent[LL.i]:=args[2]; description_laurent[LL.i]:= [args[3][1..min(1,nops(args[3]))],LL.i ,args[3][2..nops(args[3])]] else accuracy_laurent[LL.i]:=NULL; value_laurent[LL.i]:=NULL; description_laurent[LL.i]:=NULL fi; LL.i end: subsvalueslaurents:=proc(a) local i,res; res:=a; for i in indets(a) intersect set_laurents do res:=subs(i=value_laurent[i],res) od; res end: # Computes a laurent series modulo x^n # Input: laurent series, upperbound n, (optional) lower bound m nmterms_laurent:=proc() global accuracy_laurent,value_laurent; local g,procedure,dummy; if accuracy_laurent[args[1]]0. ramification_laur:=proc(L,r,ext) options remember; if r=x then RETURN(L) fi; new_laurent(max0_val(L)*degree(r,x),0,[lift_ramification_laur,L,r,ext]) end: lift_ramification_laur:=proc(laur,L,r,ext,ac) modm(g_expand(subs(x=r,nmterms_laurent(L,ceil(ac/degree(r,x)) ,ceil(accuracy_laurent[laur]/degree(r,x)))),ext))+ value_laurent[laur] end: # Converts an expression to a Laurent series. # Arguments: f, optional ext, # optional var (then f is a polynomial in var) eval_laurent:=proc() local v,f,ext,i; options remember; f:=args[1]; if nargs=1 then ext:=g_ext(f) else ext:=args[2] fi; if nargs=3 then # Consider the expression as a polynomial over Laurent series v:=args[3]; if lcoeff(f,v)=1 then # Don't apply eval_laurent on the constant 1. RETURN(expand(convert([seq(eval_laurent(coeff(f,v,i) ,ext)*v^i,i=0..degree(f,v)-1),v^degree(f,v)],`+`))) else RETURN(expand(convert([seq(eval_laurent(coeff(f,v,i) ,ext)*v^i,i=0..degree(f,v))],`+`))) fi fi; if has(f,ExpInt) or has(f,log(x)) then for i in indets(f,anything) do if (type(i,function) and op(0,i)=ExpInt) or i=log(x) then RETURN(expand(convert([seq(i^v*eval_laurent(coeff(f,i,v),ext) ,v=ldegree(f,i)..degree(f,i))],`+`))) fi od fi; if member(f,set_laurents) then RETURN(f) fi; v:=indets(f) intersect set_laurents; if v={} then if type(f,ratpoly(rational,[x,op(ext)])) then f:=[g_expand(numer(f),ext,x),g_expand(denom(f),ext,x)]; if type(f[2],polynom(rational,x)) and ldegree(f[2],x)=degree(f[2],x) then f:=expand(f[1]/f[2]); RETURN(new_laurent(0,truncate_x(f,f,0),[truncate_x,f])) fi; RETURN(new_laurent2([seq(new_laurent(0,truncate_x(i,i,0), [truncate_x,i]),i=f)],`/`,ext)) elif type(f,procedure) then RETURN(new_laurent2(f,procedure,ext)) else RETURN(new_laurent2(evala(Normal(subs(g_conversion2,f))),maple,ext)) fi fi; if type(f,polynom(rational,v)) and degree(f,v)<=2 and type(f,`+`) and coeffs(f,v)=1 then RETURN(new_laurent2([f,v],`deg 2 polynom`,ext)) fi; v:=[op(v)]; if not type(f,polynom(anything,v)) then v:=g_normal(f); RETURN(new_laurent2(map(eval_laurent,[numer(v),denom(v)]) ,`/`,ext)) fi; v:=v[1]; if f=x^degree(f,x)*v then RETURN(new_laurent2([degree(f,x),v],`*x^n`,ext)) fi; f:=collect(f,v); if f=v^degree(f,v) then f:=eval_laurent(f/v); RETURN(new_laurent2([f,v],`*`,ext)) fi; new_laurent2([eval_laurent(coeff(f,v,0)),new_laurent2( [eval_laurent(expand((f-coeff(f,v,0))/v)),v],`*`,ext)],`+`,ext) end: # For use with eval_laurent. new_laurent2:=proc() local l; options remember; if member(`deg 2 polynom`,[args]) then l:=lowerbound_val([args]); l:=min(l,2*l); new_laurent(l,0,[nterms_expression,args]) elif nargs>1 and args[2]=`/` then new_laurent(max0_val(args[1][1])-valuation_laurent(args[1][2],infinity) ,0,[nterms_expression,args]) elif nargs>1 and args[2]=`*` then new_laurent(max0_val(args[1][1])+max0_val(args[1][2]) ,0,[nterms_expression,args]) else new_laurent(-10,nterms_expression(0,args,-10),[nterms_expression,args]) fi end: # For use with eval_laurent # This procedure computes an expression mod x^n. nterms_expression:=proc(lau,aa,what,ext,n) local a,s,d,v1,v2; if what=maple then s:=series(aa,x=0,max(n+series_val(denom(aa)), series_val(denom(aa))+2)); if order(s)<>infinity then d:=0; while order(s)<>infinity and order(s)aa do d:=d+1; s:=series(aa,x=0,max(n+series_val(denom(aa)) ,series_val(denom(aa))+2)+d) od fi; s:=evala(Expand(convert(s,polynom))); normal_coeffs(subs(g_conversion1,convert( [seq(coeff(s,x,d)*x^d,d=ldegree(s,x)..n-1)],`+`)),x) elif what=`+` then if lau<>0 and value_laurent[lau]<>0 then convert(map(nmterms_laurent,aa,n,accuracy_laurent[lau]),`+`) +value_laurent[lau] else normal_coeffs(convert(map(nmterms_laurent,aa,n),`+`),x) fi elif what=`*` then v2:=max0_val(aa[2]); a:=nmterms_laurent(aa[1],n-v2); modm(g_expand(value_laurent[lau]+convert([seq(coeff(a,x,d)*x^d *nmterms_laurent(aa[2],n-d,accuracy_laurent[lau]-d) ,d=ldegree(a,x)..degree(a,x))],`+`),ext)) elif what=`*x^n` then expand(x^aa[1]*nmterms_laurent(aa[2],n-aa[1])) elif what=`/` then v1:=max0_val(aa[1]); v2:=valuation_laurent(aa[2],infinity); if lau=0 or value_laurent[lau]=0 then g_quotient(nmterms_laurent(aa[1],v2+n),nmterms_laurent( aa[2],2*v2+n-v1),n,ext) else a:=nmterms_laurent(aa[2],2*v2+n-v1); value_laurent[lau]+g_quotient(nmterms_laurent(aa[1],n+v2 ,accuracy_laurent[lau]+v2)- convert([seq(convert([seq(x^s*coeff(value_laurent[lau],x,s) ,s=accuracy_laurent[lau]-d+v2..n-d-1+v2)],`+`)*x^d *coeff(a,x,d),d=ldegree(a,x)..degree(a,x))],`+`),a,n,ext) fi elif what=differentiate then diff(nmterms_laurent(aa,n+1),x) elif what=`deg 2 polynom` then # Now aa[1] is a sum normal_coeffs( convert(map(nterms_expression,[op(aa[1])],accuracy_laurent[lau], `deg 2 1 term`,ext,n),`+`),x) elif what=`deg 2 1 term` then if type(lau,`*`) then v2:=max0_val(op(2,lau)); a:=nmterms_laurent(op(1,lau),n-v2); modm(g_expand(convert([seq(coeff(a,x,d)*x^d *nmterms_laurent(op(2,lau),n-d,aa-d) ,d=ldegree(a,x)..degree(a,x))],`+`),ext)) else nmterms_laurent(lau,n,aa) fi elif what=procedure then convert([value_laurent[lau],seq(aa(d)*x^d ,d=accuracy_laurent[lau]..n)],`+`) else bug() fi end: # Gives the derivative of an expression differentiate:=proc(ff) local f,dummy,v; v:=indets(ff) intersect set_laurents; if v={} then if has(subs(g_conversion2,subs(x=dummy,ff)),x) then RETURN(subs(g_conversion1,diff(subs(g_conversion2,ff),x))) else RETURN(diff(ff,x)) fi fi; v:=[op(v)]; v:=v[1]; f:=collect(ff,v); if type(f,polynom(anything,v)) then expand(convert([seq(differentiate(coeff(f,v,dummy))*v^dummy +dummy*v^(dummy-1)*coeff(f,v,dummy)*new_laurent2(v ,differentiate,g_ext(v)),dummy=0..degree(f,v))],`+`)) else f:=g_normal(f); g_normal((denom(f)*differentiate(numer(f))-numer(f) *differentiate(denom(f)))/denom(f)^2) fi end: # valuation of a. It looks at most as far as bound valuation_laurent:=proc(a,bound) global value_laurent; local i; options remember; i:=value_laurent[a]; if i<>0 then if g_normal(tcoeff(i,x))=0 then value_laurent[a]:=collect(i,x,g_normal); RETURN(procname(args)) fi; RETURN(min(bound,ldegree(i,x))) fi; i:=accuracy_laurent[a]-1; while (bound=infinity or i=bound then bound else procname(args) fi end: # This procedure computes the valuation of an expression. series_val:=proc(a) local d,n; options remember; if a=0 then ERROR(`division by zero`) fi; d:=1; n:=0; while n=0 do d:=d+1; n:=evala(Expand(convert(series(a,x=0,d),polynom))) od; ldegree(n,x) end: # Returns a lower bound for the valuation of the Laurent series in f. lowerbound_val:=proc(f) local i; min(seq( # op([accuracy_laurent[i],ldegree(value_laurent[i],x)]), max0_val(i), i=indets(f) intersect set_laurents)) end: # A lower bound (not much too low) for the valuation max0_val:=proc(a) local v; v:=value_laurent[a]; if v=0 then v:=accuracy_laurent[a]; if v<0 then valuation_laurent(a,0) else v fi else ldegree(v,x) fi end: # nt (n terms), compute f mod x^n # Only for expressions linear in Laurent series, like the standard form # for local differential operators. nt:=proc(f,n) local a,i; if has(f,log(x)) then RETURN(subs(a=log(x),nt(subs(log(x)=a,f),n))) fi; if has(f,ExpInt) then for i in indets(f,anything) do if type(i,function) and op(0,i)=ExpInt then RETURN(subs(a=i,nt(subs(i=a,f),n))) fi od fi; if type(f,list) then RETURN(map(nt,f,n)) fi; map(nmterms_laurent,indets(f) intersect set_laurents,n); truncate(subsvalueslaurents(f),n,x,[]) end: ############################## # Computation in k((x))[xDF] #------------------------------------------------ ############################## diffop_local # Everything from here until the help section deals with # the local factorization. The syntax for local operators is: # f=xDF^n + LL.? * xDF^(n-1) + .. # So f is monic, and has elements of set_laurents as coefficients. # Furthermore these operators must use my own syntax for algebraic numbers, # and must use xDF (differentiation followed by a multiplication by x) # as a syntax for operators, instead of DF. macro( Newtonpolygon=`diffop/Newtonpolygon`, coefs_operator=`diffop/coefs_operator`, factor_newton=`diffop/factor_newton`, factor_newton2=`diffop/factor_newton2`, factor_op=`diffop/factor_op`, factor_riccati=`diffop/factor_riccati`, faster_riccati_split=`diffop/faster_riccati_split`, indeterminate=`diffop/indeterminate`, indeterminates_op=`diffop/indeterminates_op`, lift_newton=`diffop/lift_newton`, lift_ramification=`diffop/lift_ramification`, lift_rightfactor=`diffop/lift_rightfactor`, lift_rsplit=`diffop/lift_rsplit`, lift_rsplit2=`diffop/lift_rsplit2`, lift_substitute=`diffop/lift_substitute`, make_rightfactor=`diffop/make_rightfactor`, nm_mult=`diffop/nm_mult`, nm_block=`diffop/nm_block`, nm_block2=`diffop/nm_block2`, op_with_slope=`diffop/op_with_slope`, ram_laur=`diffop/ram_laur`, ramification_of=`diffop/ramification_of`, rem_lift=`diffop/rem_lift`, skipped_factors=`diffop/skipped_factors`, substitute=`diffop/substitute` ): # Gives terms of a multiplication in k[x,xDF] from x^low to x^high # So the accuracy is high+1 nm_mult:=proc(l,r,low,high,ext) local j,i; modm(g_expand(convert([seq(coeff(r,x,j)*x^j*subs(xDF=xDF+j ,convert([seq(x^i*coeff(l,x,i),i=max(ldegree(l,x),low-j)..min(degree( l,x),high-j))],`+`)),j=max(ldegree(r,x),low-degree(l,x))..min(high -ldegree(l,x),degree(r,x)))],`+`),ext)) end: # 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() local a,b,i,Dib,result,D,xx,ext; options remember; 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: # ff is f with xDF + a for xDF substituted. This procedure lifts ff. # Note: the resulting Laurents may have a higher degree then their # accuracy. The terms higher than the accuracy are not yet correct. lift_substitute:=proc(l,d,slope,a,f,ext,ff,acc) global value_laurent,accuracy_laurent; local i,res,n_lifts; n_lifts:=acc-accuracy_laurent[l]; if n_lifts<=0 then RETURN() fi; i:=accuracy_laurent[coeff(ff,xDF,0)]+d*slope; res:=modm(substitute(a,nm_block(f,i,i+n_lifts,slope,ceil),ext)); for i from 0 to degree(f,xDF)-1 do value_laurent[coeff(ff,xDF,i)]:=coeff(res,xDF,i) +value_laurent[coeff(ff,xDF,i)]; accuracy_laurent[coeff(ff,xDF,i)]:= accuracy_laurent[coeff(ff,xDF,i)]+n_lifts od; NULL end: # This procedure substitutes a for xDF in an operator f in internal form. # a must be xDF + something of valuation at least -slope # arguments: a,f,slope,shift,ext # Other call syntax: a,f,ext substitute:=proc() local i,ide,res,f,D; if type(args[2],list) then f:=args[2]; if nops(f)=4 and f[4]=`alg over k((x))` then RETURN([substitute(ramification_of( args[1],f[3],args[nargs]),f[1],args[3..nargs]),f[2..4]]) else RETURN([seq(substitute(args[1],i,args[3..nargs]),i=f)]) fi elif nargs=5 then op_with_slope(degree(args[2],xDF),args[3..4], [lift_substitute,args[1..2],args[5]]) elif nargs=3 then if has(args[2],DF) then D:=DF else D:=xDF fi; res:=coeff(args[2],D,0); for i from 1 to degree(args[2],D) do if i=1 then ide:=args[1] else ide:=mult(args[1],ide,args[3]) fi; res:=res+coeff(args[2],D,i)*ide od; if type(res,polynom(anything,x)) or D=xDF then g_expand(res,args[3]) else collect(res,D,g_normal) fi elif nargs=2 or nargs=4 then substitute(args,g_ext([args])) fi end: # Returns an operator having a given slope. op_with_slope:=proc(order,slope,shift,d) global description_laurent; local i,result; result:=xDF^args[1]; for i from 0 to args[1]-1 do result:=result+new_laurent(ceil((i-args[1])* slope+shift),0,[])*xDF^i # Note: shift <= 0. od; if d=[] then RETURN(result) elif member(d[1],{lift_newton,lift_rsplit}) then result:=[result,op_with_slope(d[2],slope,0,[])] fi; for i in indets(result) minus {xDF} do description_laurent[i]:= [d[1],i,args[1..2],d[2..nops(d)],result] od; result end: # Gives all terms from x^n to x^m of f if slope=0. nm_block:=proc(f,n,m,slope,round_off) local i,res; if round_off(n)<=0 and 0)`, ``, `PARAMETERS:`, ` f - a differential operator`, ``, `SYNOPSIS:`, ``, ` Various conversions can be performed with this command, see the examples. For`, ` local operators a different syntax is used than for global operators. Local`, ` operators are represented as non-commutative polynomials in xDF with Laurent`, ` series (always with finite pole order) as coefficients. Conversion of a global`, ` operator to a local operator in a point p can be done with the argument x=p.`, ` A truncation of this local operator is given using the argument ``truncated``.`, ` The result will be represented as a truncated power series in x, instead of a`, ` series in x-p.`, ``, ` The name xDF stands for x times DF, i.e. a differentiation followed by a`, ` multiplication by x.`, ``, `EXAMPLES:`, ` f:=2*DF^2-1;`, ` 2`, ` f := 2 DF - 1`, ``, ` convert(f,diffop,monic); # makes f monic`, ` 2`, ` DF - 1/2`, ``, ` convert(f,diffop,axiom); # first use: y:=operator y in Axiom`, ``, ` - y(x) + 2 D(y(x), x, 2) = 0`, ``, ` convert(f,diffop,maple);`, ``, ` / 2 \`, ` | d |`, ` 2 |----- y(x)| - y(x) = 0`, ` | 2 |`, ` \ dx /`, ``, ` convert(f,diffop,reduce); # Apply lprint after these latter conversions to`, ` # get machine readable expressions`, ``, ` - Y + 2 DF(Y, X, 2) = 0`, ``, ` convert(",diffop);`, ` 2`, ` 2 DF - 1`, ``, ` convert(f,diffop,x=0);`, ` # Converts f to a local operator at x=0, and also makes the result monic`, ` 2`, ` LL4 + LL5 xDF + xDF`, ``, ` convert(",diffop,truncated,10); # gives f modulo x^10`, ``, ` 2 2`, ` - 1/2 x - xDF + xDF`, ``, `SEE ALSO: diffop` ): `convert/diffop`:=proc() global x,y; local f,i,l; f:=args[1]; if nargs=1 then if member(xDF,indets(f)) then f:=substitute(x*DF,args) elif member(X,indets(f)) then # reduce's syntax f:=subs({X=x,Y=1},eval(subs(DF =proc() if nargs=2 then xDF else xDF^args[3] fi end,f))); f:=op(1,sort(subs(xDF=DF,f),[x,DF])) elif f<>eval(subs(y(x)=1,f)) or has(convert(f,D),D) then if not has(f,y(x)) or not type(f,equation) or not op(2,f)=0 then ERROR(`wrong type of arguments`) fi; # Maple's syntax f:=convert(f,D); i:=1; l:=diff(y(x),x); while has(f,D) do f:=subs(convert(l,D)=DF^i,f); l:=diff(l,x); i:=i+1 od; f:=op(1,subs(y(x)=1,f)) fi; RETURN(collect(f,DF,g_normal)) elif nargs>=2 then if nargs=3 and args[2]=truncated then RETURN(subs(g_conversion2,nt(f,args[3]))) fi; f:=convert(f,diffop); if type(args[2],equation) and op(args[2])[1]=x then RETURN(l_p(f,op(args[2])[2])) elif args[2]=reduce then RETURN(subs(x=X,coeff(f,DF,0)*Y+sum(coeff(f,DF,i)*DF(Y,X,i) ,i=1..degree(f,DF)))=0) elif args[2]=axiom then RETURN(coeff(f,DF,0)*y(x)+sum(coeff(f,DF,i)*''D''(y(x),x,i) ,i=1..degree(f,DF))=0) elif args[2]=maple then RETURN(sum('coeff(f,DF,i)*diff(y(x),x$i)',i=1..degree(f ,DF))+coeff(f,DF,0)*y(x)=0) elif args[2]=monic then RETURN(collect(f/lcoeff(f,DF),DF,g_normal)) fi fi; ERROR(`wrong arguments`) end: # 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(f,p) local g,i,ext; if type(p,equation) then ERROR(`wrong input`) fi; ext:=g_ext([args]); 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(convert(subs(g_conversion2,g),diffop,monic)) 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) local i; options remember; 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) local m,i,j; 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; local p,v,m,i,aa; m:=accuracy_laurent[L]; v:=0; p:=P; for i from m to n-1 do v:=v+subs(x=a,rem(p,pa,x))*x^i; p:=collect(diff(p,x)/(i+1),x) 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) local v; options remember; 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) local i,s,m,c; 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: ######################################### # Factorization in k((x))[xDF] # ######################################### # Input: a monic operator f in xDF, in Laurent form # a string what: # `split over k((x))` gives a factorization over k((x)) # `all right factors` gives all safe right factors over k((x)) # `all alg factors` gives all safe right factors over alg. closure of k((x)) # `alg factor` gives one safe right factor over alg. closure k((x)) # ext: gives the constant field factor_op:=proc() local r,dummy; options remember; userinfo(6,diffop,`factorizing`,args); if nargs=1 then # shouldn't use factor_op with 1 argument RETURN(factor_op(eval_laurent(args[1],[],xDF) ,`split over k((x))`,[])) elif nargs=2 then RETURN(factor_op(args,g_ext([args]))) fi; if degree(args[1],xDF)<=1 then RETURN([args[1]]) fi; r:=factor_newton(args[1],args[2..nargs]); [seq(op(factor_riccati(dummy,args[2..nargs])),dummy=r)] end: ################################################# # Factorization using the Newton method # ################################################# # Input: an operator f, and (optional) a second argument # Output (if called with 2 arguments): a list giving 3 things: # 1) the coordinates of the extreme points of the Newton polygon # 2) the slope of a point to the next point # 3) the Newton polynomial of this slope # If called with only 1 argument then only the extreme points will be given. Newtonpolygon:=proc() global Newtonpolygon; local f,n,vals,dummy,res,m,i,powd, val_powd,npg,slope; options remember; userinfo(8,diffop,`Computing Newton polygon of`,args); f:=args[1]; n:=degree(f,xDF); vals:=[seq(valuation_laurent(coeff(f,xDF,dummy),0),dummy=0..n-1),0]; powd:=0; val_powd:=min(op(vals)); npg:=[[powd,val_powd]]; while powd(i-powd)*m do i:=i-1 od; powd:=i;val_powd:=vals[i+1]; npg:=[op(npg),[i,vals[i+1]]] od; if nargs=1 then RETURN(npg) else Newtonpolygon(args[1]):=npg fi; res:=NULL; for i from 1 to nops(npg)-1 do # Now we compute the Newton polynomial of slope number i slope:=(npg[i+1][2]-npg[i][2])/(npg[i+1][1]-npg[i][1]); res:=res,[op(npg[i]),slope,convert([seq( nthterm_laurent(coeff(f,xDF,denom(slope)*dummy+npg[i][1]), dummy*numer(slope)+npg[i][2]) *x^dummy,dummy=0..(npg[i+1][1]-npg[i][1])/denom(slope) )],`+`)] od; [res,npg[nops(npg)]] end: # Output: coprime index 1 factorizations # Input: operator f, monic in xDF # a string what: # `split over k((x))` results in a complete Newton factorization, i.e. the # broken Newton polygon and the gcd 1 reducible Newton polynomial cases # will be factored # `all right factors` gives all possible safe right factors according to # the Newton method # `all alg factors` idem # `alg factor` gives only one right factor, the one where the slope is the # least steep # `semireg` gives a coprime index 1 LCLM factorization, to be used for # computing the semi-regular parts # A Newton polynomial like (a-1)^2 will only be factored when slope=0 # (regular singular case). No algebraic extensions will be made here, they # will be made in factor_riccati factor_newton:=proc(f,what,ext) global skipped_factors; local np,v,dummy,i,j,k,d,e,res,unsafe,n_unsafe,semi; options remember; np:=Newtonpolygon(f,`include the Newton polynomials`); userinfo(7,diffop,`Newton method factorizing`,args); res:=NULL; unsafe:={}; for k from 1 to nops(np)-1 do v:=g_factors(np[k][4],x,ext); if np[k][3]<>0 then v:=[seq([g_expand(dummy[1]^dummy[2],ext,x),0],dummy=v)] else # regular singular case, we compute the unsafe factors unsafe:={}; for i in v do n_unsafe[i]:=0;semi[i]:=i[1]^i[2] od; # n_unsafe is the number of factors with an exponent that is # the exponent of this factor i minus a positive integer. # The "unsafe" means that we can not use those factors as # a right hand factor because it is not a priori known whether # such a factor can be lifted or not (hence: placing such a # factor on the right hand side would cause the risk of an # error during the lifting process). If such an "unsafe" # factor can be lifted or not depends on if there is a log in # the formal solutions. # The n_unsafe number is used in the global factorization, in # the procedure skipped_factors. for i from 1 to nops(v)-1 do for j from i+1 to nops(v) do if degree(v[i][1],x)=degree(v[j][1],x) then d:=degree(v[i][1],x)-1; e:=g_normal(coeff(v[i][1],x,d)-coeff(v[j][1],x,d)); if type(e,integer) and e<>0 and irem(e,d+1)=0 and g_expand(v[i][1]-subs(x=x+e/(d+1),v[j][1]),ext,x)=0 then if e>0 then unsafe:=unsafe union {v[i]}; n_unsafe[v[j]]:=n_unsafe[v[j]]+v[i][2]; semi[v[j]]:=semi[v[j]]*v[i][1]^v[i][2] else unsafe:=unsafe union {v[j]}; n_unsafe[v[i]]:=n_unsafe[v[i]]+v[j][2]; semi[v[i]]:=semi[v[i]]*v[j][1]^v[j][2] fi fi fi od od; v:=[op({op(v)} minus unsafe)]; if what=`semireg` then v:=[seq([g_expand(semi[i],ext,x),1],i=v)] fi fi; # Now v contains the safe right factors of slope number k for i in v do if degree(v[1][1],x)*denom(np[k][3])=degree(f,xDF) then # f allows no coprime index 1 factorization RETURN([f]) fi; j:=factor_newton2(f,i[1],np[k][3],ext); if np[k][3]=0 and what<>`semireg` then skipped_factors(j[2]):=n_unsafe[i] fi; if what=`alg factor` then RETURN([j[2]]) elif what=`split over k((x))` then RETURN([op(factor_newton(j[1],what,ext)),j[2]]) fi; res:=res,j[2] od od; [res] end: # For use with: factor_op( .. , `all right factors`, []) in the global # factorization. # It gives for each factor the number of factors we skipped, because of # the integer difference of the exponents (called the unsafe factors). # If minimum multiplicity = 1 then skipped_factors=0. If skipped_factors=0 # then not necessarily minmult=1, however, we may still treat this case in # the global factorization as if it were minmult=1, because then in this # exponential part there is only 1 unique right hand factor anyway (the others # give rise to a log in the formal solutions). # # This procedure will probably get obsolete soon. skipped_factors:=proc(f) local v; options remember; v:=[op(indets(f) intersect set_laurents)]; v:=description_laurent[v[1]]; if v[1]=lift_rightfactor then skipped_factors(v[5]) elif v[1]=lift_rsplit then skipped_factors(v[8]) elif v[1]=lift_substitute then skipped_factors(v[6]) elif member(v[1],{lift_newton,nterms_expression,truncate_x}) then # the options remember must treat the regular singular case 0 else ERROR(`?`) fi end: # Factor monic operator f, such that r is the Newton polynomial of # the right factor # Output: a list of 2 factors # r must be monic. factor_newton2:=proc(f,r,slope,ext) global rem_lift; local np,l,i,res,shift; options remember; userinfo(7,diffop,`Computing factor with Newton polynomial:`,r); np:=Newtonpolygon(f,`include the Newton polynomials`); for i from 1 to nops(np)-1 do if np[i][3]=slope then l:=i fi od; if l=evaln(l) then ERROR(`slope not found`) fi; shift:=min(0,np[l][2]+(degree(f,xDF)-np[l][1])*slope); # We will compute the Newton polynomial of the left factor np:=g_expand(x^np[l][1]*subs(x=x^denom(slope),g_quo(np[l][4],r,x,ext)),ext,x); res:=op_with_slope(degree(f,xDF)-degree(r,x)*denom(slope) ,slope,shift,[lift_newton,degree(r,x)*denom(slope),f,np, subs(x=x^denom(slope),r),shift,ext]); rem_lift[res]:=0; res end: # This procedure lifts "coprime index 1 factorizations" # Meaning of the variables: # llaur: the Laurent series we lift # acc: desired accuracy of LL. Note that the other series in left, right # are lifted too. # n_lifts: number of lifts to do # n_known: number of computed coefficients # slope # mult_args: arguments for nm_mult to determine which terms are needed # ff, left, right: exact operators f=left*right # l, r: the lowest parts of left and right # f: a middle part of ff, precisely the terms we need # lr: the product l*r # rem_lift: Those terms of the previous lift of lr, that we can use now. # l_low, r_low: the lowest line of left and right. These have gcd 1. # If they have gcd<>1 than it is not uniquely liftable. # l_extra, r_extra: the new terms we computed of l and r. # ext: the algebraic extension lift_newton:=proc(llaur,order_l,slope,order_r,ff,l_low,r_low,shift,ext,v,acc) global rem_lift,value_laurent,accuracy_laurent; local f,lr,l,r,l_extra,r_extra,s,t,n_known,i,lau ,n_lifts,mult_args,left,right,ll_low,rr_low,fe,le,re; n_lifts:=acc-accuracy_laurent[llaur]; if n_lifts<=0 then RETURN() fi; left:=v[1]; right:=v[2]; n_known:=accuracy_laurent[coeff(right,xDF,0)]+degree(right,xDF)*slope; mult_args:=ceil(n_known+shift-slope*degree(ff,xDF)), ceil(n_known+shift+n_lifts)-1,ext; f:=nm_block(ff,n_known+shift,n_known+n_lifts+shift,slope,ceil); l:=expand(subsvalueslaurents(left)); if n_known+shift<=0 then l:=l-xDF^degree(l,xDF) fi; if n_known=0 then r:=0 else r:=expand(subsvalueslaurents(right)) fi; lr:=rem_lift[v]-f+nm_mult(l,r,mult_args[2]+1-n_lifts ,mult_args[2..nops([mult_args])]); if slope=0 then # regular singular case rr_low:=subs(x=xDF,r_low); for i from n_known to n_known+n_lifts-1 do if i=0 then r_extra:=subs(x=xDF,r_low); l_extra:=expand(x^shift*subs(x=xDF,l_low)) else ll_low:=subs(x=xDF,g_expand(subs(x=x+i,l_low),ext)); s:=g_gcdex(rr_low,ll_low,1,xDF,ext); if s[1]<>1 then # ERROR(`unsafe factor`) bug() fi; t:=s[3]; s:=s[2]; r_extra:=-g_rem(t*coeff(lr,x,i+shift),rr_low,xDF,ext); l_extra:=modm(collect(-x^(i+shift)*g_quo(coeff(lr,x ,i+shift)+r_extra*ll_low,rr_low,xDF,ext),x)); r_extra:=modm(collect(x^i*r_extra,x)) fi; l:=l+l_extra; lr:=modm(lr+nm_mult(l,r_extra,mult_args)+nm_mult(l_extra,r ,mult_args)); r:=r+r_extra od else s:=g_gcdex(r_low,l_low,1,x,ext); t:=s[3]; s:=s[2]; le:=denom(slope)*shift-numer(slope)*degree(left,xDF); re:=-numer(slope)*degree(right,xDF); fe:=denom(slope)*shift-numer(slope)*degree(ff,xDF); for i from n_known*denom(slope) to (n_known+n_lifts)*denom(slope)-1 do if i=0 then r_extra:=coefs_operator(r_low,slope,re,1); l_extra:=coefs_operator(l_low,slope,le,1) else r_extra:=g_rem(t*coefs_operator(lr,slope, i+fe,0),r_low,x,ext); l_extra:=coefs_operator(modm(g_expand(-g_quo(( coefs_operator(lr,slope,i+fe,0)-r_extra*l_low),r_low,x,ext) ,ext)),slope,i+le,1); r_extra:=coefs_operator(modm(-r_extra),slope,i+re,1) fi; l:=l+l_extra; lr:=lr+nm_mult(l_extra,r,mult_args)+nm_mult(l,r_extra,mult_args); r:=r+r_extra od; fi; rem_lift[v]:=lr; for f in [[left,collect(l,xDF)],[right,collect(r,xDF)]] do for i from 0 to degree(f[1],xDF)-1 do lau:=coeff(f[1],xDF,i); value_laurent[lau]:=coeff(f[2],xDF,i); accuracy_laurent[lau]:=accuracy_laurent[lau]+n_lifts od od; NULL end: # Input: an operator or a polynomial f # Output: a polynomial if f is an operator, and an operator if f is # a polynomial. # This procedure either takes coefficients from an operator and gives them as # a polynomial, or constructs an operator from a given set of coefficients # (given as a polynomial). Something like a Newton polynomial. # slope must be > 0 coefs_operator:=proc(f,slope,i,what) local dummy,start_x,start_D; start_D:=modp(-i/numer(slope),denom(slope)); start_x:=start_D*slope+i/denom(slope); if what=0 then convert([seq(coeff(coeff(f,x,start_x+numer(slope)*dummy),xDF, start_D+dummy*denom(slope))*x^(start_D+dummy*denom(slope)), dummy=0..floor(degree(f,xDF)/denom(slope)))],`+`) else convert([seq(coeff(f,x,start_D+dummy*denom(slope))*x^(start_x +numer(slope)*dummy)*xDF^(start_D+dummy*denom(slope)) ,dummy=0..ceil(degree(f,x)/denom(slope)))],`+`) fi end: ################################################# # Factorization using a Riccati solution # ################################################# # f should be monic, have only 1 slope, and the Newton polynomial is the # power of an irreducible polynomial. The dummy argument is because of the # options remember in case of a different groundfield. # output: same as factor_op # If what=`semireg` then the output is a list of things like: # [semi-regular operator,ext,ram,exp part,`semireg`] # the exp part e is of the form: ramification_of(xDF-real exp part)=xDF-e factor_riccati:=proc(f,what,ext) local np,slope,gr,res,r,l,lr,i,v,n,dummy,ex; options remember; np:=Newtonpolygon(f,`include the Newton polynomials`)[1]; slope:=np[3]; if what=`semireg` and slope=0 then v:=op(factor_op(f,`alg factor`,ext)); if not type(v,list) then v:=[v,ext,x] fi; i:=xDF-nt(v[1],1); n:=degree(f,xDF)/degree_ext(v[2],ext); if n=1 then RETURN([[eval_laurent(v[1]+i,v[2],xDF),v[2..3],i]]) fi; # Only search for integer roots: np:=factors(numer(g_expand(subs(x=x+i,np[4]),v[2])))[2]; r:=1; for l in np do if type(x-l[1],integer) then r:=r*l[1]^l[2] fi od; if degree(r,x)<>n then bug() fi; n:=min(solve(r)); r:=expand(subs(x=x+n,r)); RETURN([[factor_newton2(substitute(xDF+i+n,f,0,0,v[2]),r,0,v[2])[2] ,v[2..3],i+n]]) fi; np:=g_factors(np[4],x,ext)[1]; if (degree(f,xDF)<=1) or (np[2]=1 and member(what ,{`split over k((x))`,`all right factors`})) then if what=`semireg` then # degree(f,xDF)=1 and slope<>0 i:=xDF-nt(f,1); RETURN([[eval_laurent(f+i,ext,xDF),ext,x,i]]) else # f is irreducible RETURN([f]) fi elif degree(np[1],x)=1 and denom(slope)=1 then # now slope<>0 np:=(x-np[1])*x^(-slope); userinfo(8,diffop,`substituting`,xDF=xDF+np); # Apply a homomorphism xDF -> xDF+np to simplify the problem. v:=factor_op(substitute(xDF+np,f,slope,0,ext),what,ext); if what=`semireg` then RETURN([seq([i[1..3],i[4]+ solve(ramification_of(xDF-np,i[3],i[2]),xDF)],i=v)]) fi; # Apply the inverse homomorphism. Note: extension could have changed RETURN(substitute(xDF-np,v,slope,0,g_ext(v))) elif what=`split over k((x))` then r:=factor_riccati(f,`alg factor`,ext); # An order 1 factor over the alg. closure of k((x)) r:=make_rightfactor(f,r[1],ext); # a righthand factor over k((x)) if r=f then RETURN([f]) fi; l:=rightdivision(f,r,ext,slope)[1]; # We factorized f here in l and r, however, these l and r lift very # slowly, so we try userinfo(7,diffop,`Computing factor using Riccati solution`,r); lr:=faster_riccati_split(f,l,r,ext,slope); RETURN([op(factor_riccati(lr[1],what,ext)),lr[2]]) elif what=`all right factors` then v:=factor_riccati(f,`all alg factors`,ext); res:=NULL; for i in v do r:=make_rightfactor(f,i,ext); if r=f then res:=f else l:=rightdivision(f,r,ext,slope)[1]; userinfo(7,diffop, `Computing factor using Riccati solution`,r); lr:=faster_riccati_split(f,l,r,ext,slope); res:=res,lr[2] fi od; RETURN([res]) elif degree(np[1],x)>1 then # Now we need an algebraic extension on k gr:=g_zero_of(np[1],x,dummy); userinfo(7,diffop,`Making an algebraic extension:`,gr); ex:=[gr,op(ext)]; r:=factor_newton2(f,g_expand((x-gr)^np[2],ex) ,slope,ex); r:=factor_riccati(r[2],what,ex); # Now we denote the algebraic extension in a list: res:=NULL; for i in r do if type(i,list) then res:=res,i else res:=res,[i,ex,x,`alg over k((x))`] fi od; RETURN([res]) else # Now we need a ramification userinfo(7,diffop,`Applying a ramification`); n:=mods(1/numer(slope),denom(slope)); r:=(x-np[1])^n*x^denom(slope); v:=factor_newton2(ramification_of(f,r,ext),g_expand((x-denom(slope) *(x-np[1])^((1-n*numer(slope))/denom(slope)))^np[2],ext) ,numer(slope),ext); v:=factor_riccati(v[2],what,ext); res:=NULL; for i in v do if type(i,list) then res:=res,[i[1],i[2],ramification_of(r,i[3],ext),i[4]] else res:=res,[i,ext,r,`alg over k((x))`] fi od; RETURN([res]) fi end: # gives a ramification, it maps xDF to xDF* 1/degree(r), and maps x to r # r is a power of x (then we have a pure ramification), or a power of x # times a constant. We allow these latter kind of ramifications, because # then we need less algebraic extensions. ramification_of:=proc(f,r,ext) local s; options remember; if r=x then RETURN(f) elif indets(f) intersect set_laurents = {} then # May 4, 1995: bug: RETURN(g_expand(subs(x=r,f),ext)) RETURN(g_expand(subs({x=r,xDF=xDF/degree(r,x)},f),ext)* degree(r,x)^degree(f,xDF)) fi; # May 4, 1995: bug: the following assumes that there is only 1 slope. # op_with_slope(degree(f,xDF),-degree(r,x)*valuation_laurent( # coeff(f,xDF,0),0)/degree(f,xDF),0,[lift_ramification,f,r,ext]) # fix: s:=Newtonpolygon(f); s:=s[nops(s)-1]; s:=s[2]/(s[1]-degree(f,xDF)); op_with_slope(degree(f,xDF),degree(r,x)*s ,0,[lift_ramification,f,r,ext]) end: lift_ramification:=proc(laur,order,slope,ff,r,ext,f,acc) global value_laurent, accuracy_laurent; local i,res,n_lifts; n_lifts:=acc-accuracy_laurent[laur]; if n_lifts<=0 then RETURN() fi; n_lifts:=ceil(n_lifts/degree(r,x))*degree(r,x); i:=(accuracy_laurent[coeff(f,xDF,0)]+order*slope)/degree(r,x); res:=g_expand(degree(r,x)^order*subs({x=r,xDF=xDF/degree(r,x)}, nm_block(ff,i,i+n_lifts/degree(r,x),slope/degree(r,x),ceil)),ext); for i from 0 to degree(f,xDF)-1 do value_laurent[coeff(f,xDF,i)]:= value_laurent[coeff(f,xDF,i)]+coeff(res,xDF,i); accuracy_laurent[coeff(f,xDF,i)]:= accuracy_laurent[coeff(f,xDF,i)]+n_lifts od; NULL end: ########################################################### # make right-factor over k((x)) using a Riccati solution # ########################################################### # This section computes LCLM(xDF-r,`and conjugates over k((x))`) where r is a # Riccati solution. In this way a local irreducible factor is obtained. This # is only necessary if the coprime index is > 1, otherwise one can factor # via the Newton method. # f is an operator # ric is a right-factor of order 1 over the algebraic closure of k((x)) # output: a right factor over k((x)) make_rightfactor:=proc(f,ric,ext) local d,v; d:=degree(ric[3],x); v:=-valuation_laurent(coeff(ric[1],xDF,0),0)/d; # Now we multiply by the degree of the algebraic extension d:=d*degree_ext(ric[2],ext); if d=degree(f,xDF) then RETURN(f) fi; op_with_slope(d,v,0,[lift_rightfactor,ric,ext,2]) end: # This procedure lifts the factor generated by make_rightfactor. It is rather # slow, therefor we also have a faster lift method. We still need this # procedure because the faster method does not work for the first few lifts. # Meaning of the variables: # ric: solution of the Riccati, i.e. 1st order factor over alg. clos. k((x)) # ram: ramification is a map x -> c * x^(ram) for a constant c. # b[0], b[1], ..: these are used as an Ansatz # This procedure generates linear equations stating that ric is a righthand # factor. This righthand factor gets computed by solving these equations. lift_rightfactor:=proc(dummy_laur,order,slope,ric,ext,acc,som,dummy_ac) global value_laurent, accuracy_laurent, description_laurent; local s,r,ram,rp,i,b,dummy,fout,l,ld,extl; ram:=degree(ric[3],x); l:=coeff(ric[1],xDF,0); r:=-nmterms_laurent(coeff(ric[1],xDF,0),acc); ld:=min(0,ldegree(r,x)); rp:=1; s:=b[0]; for i from 1 to order do l:=(i-1)*ld+acc; rp:=truncate(x*diff(rp,x)+r*rp,l,x,ext)/ram+fout[i]*x^l; s:=s+b[i]*rp od; s:=g_expand(subs(b[order]=1,s),ric[2]); extl:=NULL; if ram>1 then s:=subs(x=g_zero_of(subs(x=xDF,ric[3])-x,xDF,'extl'),s); extl:=eval(extl); #11jan: bug: g_expand does not work for negative powers of alg. expressions. s:=expand(s/extl^ldegree(s,extl)); s:=g_expand(s,[extl,op(ric[2])]) fi; s:={ext_to_coeffs(s,ext)}; assign(g_solve(s,{seq(b[dummy],dummy=0..order-1)},ext,0)); for i from 0 to order-1 do l:=coeff(som,xDF,i); b[i]:=nterms_expression(0,g_normal(b[i]),maple,ext ,accuracy_laurent[l]+3); b[i]:=collect(b[i],x,g_normal); s:=ldegree(b[i],x)-1; while indets(coeff(b[i],x,s)) intersect {seq(fout[dummy] ,dummy=1..order)} = {} and s 1 lift algorithm) # ############################################################################ # I think I'll change the global factorizer such that it does not use any # coprime index>1 factorizations anymore. When that is done, this section # and the previous section won't be important anymore (only for local # factorizations). # Input: f=left*right # Output: l and r such that l=left and r=right, and such that l and r lift # faster then left and right. faster_riccati_split:=proc(f,left,right,ext,slope) # RETURN([left,right]); op_with_slope(degree(left,xDF),slope,0,[lift_rsplit,degree(right,xDF) ,f,left,right,ext,2]) end: lift_rsplit:=proc(laur,order_l,slope,order_r,f,left,right,ext,n_ind,v,acc) global description_laurent, value_laurent, accuracy_laurent; local try,i,lau,aclau,mo,j; try:=lift_rsplit2(laur,f,slope,op(v),n_ind,ext,acc); if try=`n_ind is too small` then # I must use more indeterminates userinfo(6,diffop,`Increasing the coprime index of`,right,`to`,n_ind+1); for i in indets(v) minus {xDF} do description_laurent[i]:= [lift_rsplit,i,args[2..8],n_ind+1,v] od; elif try=`too few terms are known` then # the slow procedure lift_rightfactor will compute more terms # mo:=modulus; # if mo<>0 then # `diffop/compute_modp`(); # modulus:=0; # lau:=coeff(v[1],xDF,0); # nt(lau,accuracy_laurent[lau]+3); # modulus:=mo # else for i from 0 to degree(v[1],xDF)-1 do lau:=coeff(v[1],xDF,i); aclau:=accuracy_laurent[lau]+1; value_laurent[lau]:=nmterms_laurent(coeff(left,xDF,i) ,aclau); accuracy_laurent[lau]:=aclau od; for i from 0 to degree(v[2],xDF)-1 do lau:=coeff(v[2],xDF,i); aclau:=accuracy_laurent[lau]+1; value_laurent[lau]:=nmterms_laurent(coeff(right,xDF,i) ,aclau); accuracy_laurent[lau]:=aclau od # fi fi; NULL end: lift_rsplit2:=proc(llaur,ff,slope,left,right,n_ind,ext,acc) global value_laurent, accuracy_laurent; local f,lr,l,r,l_extra,r_extra,l_low,r_low, left_ind,right_ind, n_known,i,v,all_ind,lau,n_lifts,mult_args; n_lifts:=acc-accuracy_laurent[llaur]; # number of lifts if n_lifts<=0 then RETURN() fi; n_known:=accuracy_laurent[coeff(right,xDF,0)]+degree(right,xDF)*slope; # number of computed coefficients if n_known {} then RETURN(`n_ind is too small`) fi; lr:=lr+nm_mult(l+l_extra,r_extra,mult_args)+nm_mult(l_extra,r,mult_args); l:=l+l_extra; r:=r+r_extra; n_known:=n_known+1 od; 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; for i from 0 to degree(right,xDF)-1 do lau:=coeff(right,xDF,i); value_laurent[lau]:=coeff(r,xDF,i); accuracy_laurent[lau]:=accuracy_laurent[lau]+n_lifts od; NULL end: indeterminates_op:=proc(slope,number,f) local i,j,res; options remember; res:=0; for i from 0 to degree(f,xDF)-1 do for j from 0 to number-1 do res:=res+indeterminate()*x^(j+ceil( (i-degree(f,xDF))*slope))*xDF^i od od; res end: indeterminate:=proc() local r;r end: ############################################### # Formal solutions of differential equations #-------------------------------- ############################################### diffop_formalsol # This part of diffop computes solutions of f=xDF^n+laur*xDF^(n-1)+... # in algclos( k((x)) )[log(x),exp( .. )] # All solutions are given up to conjugation over k((x)) # The result is a list which contains a list of solutions for every exponential # part e. Such a list for 1 exp part looks like: # [ [a1,a2,..] , e] where a.i in k((x))[e,log(x)] # The corresponding solutions are: a.i*ExpInt(e) where ExpInt(e)=exp(int(e/x,x)) # If any ramifications (then x is not really x but a root of a*t^n=x) have been # used then the list is longer and looks like: # [ [a1,a2,..] , e,ext,ram,`alg over k((x))` ] where ram=a*x^n # The same form is used when algebraic extensions were needed. macro( solve_semiregular=`diffop/solve_semiregular`, integrate_logs=`diffop/integrate_logs`, lift_integral=`diffop/lift_integral`, sol_1order_eq=`diffop/sol_1order_eq`, lift_sol1=`diffop/lift_sol1`, test_result=`diffop/test_result` ): # :) # Input: an operator f of order >= 1 # Output: the formal solutions formal_sol:=proc(f,ext) local v,i,res,c,lp; options remember; if nargs=1 then RETURN(procname(f,g_ext(f))) fi; if has(f,DF) then RETURN(procname(l_p(f,0),ext)) fi; v:=factor_op(f,`semireg`,ext); res:=NULL; for i in v do c:=solve_semiregular(factor_op(i[1],`split over k((x))`,i[2]),i[2]); lp:=[c,i[4]]; if not (i[2]=ext and i[3]=x) then lp:=[op(lp),i[2..3],`alg over k((x))`] fi; test_result(formal_sol,f,ext,op(lp)[1..2],i); res:=res,lp od; [res] end: # Input: a list of factors [xDF+r1,...xDF+rn] where # the 0'th coefficients of r.i are in Z. # Output: the solutions solve_semiregular:=proc(v,ext) local s,i; if v=[] then RETURN([]) fi; s:=solve_semiregular([seq(xDF+eval_laurent(v[i]-v[nops(v)] ,ext),i=1..nops(v)-1)],ext); s:=[op(map(integrate_logs,s,ext)),1]; [seq(eval_laurent(expand(i*sol_1order_eq(coeff(v[nops(v)],xDF,0) ,ext)),ext,log(x)),i=s)] end: # used by solve_semiregular # determine g such that x*d/dx (g)=ff where ff in k((x))[log(x)] integrate_logs:=proc(ff,ext) local i,d,l,f; options remember; d:=degree(ff,log(x)); if d=0 then if has(ff,set_laurents) then f:=eval_laurent(ff,ext); # l becomes the integral of f/x: nmterms_laurent(f,1,0)*log(x)+new_laurent(1,expand(int( nmterms_laurent(f,0)/x,x)),[lift_integral,f,ext]) else g_expand(int(ff/x,x),ext) fi else # Bugfix (oct 94) l:=integrate_logs(coeff(ff,log(x),d),ext); l:=coeff(l,log(x),0)+log(x)*coeff(l,log(x),1)/(d+1); expand(l*log(x)^d+integrate_logs(convert([seq(coeff(ff,log(x),i) *log(x)^i,i=0..d-1)],`+`)-d*log(x)^(d-1)*coeff(l,log(x),0),ext)) fi end: # Used by integrate_logs # lifts LL which is int(aa/x , x) lift_integral:=proc(LL,aa,ext,acc) value_laurent[LL]+expand(int( nmterms_laurent(aa,acc,accuracy_laurent[LL])/x,x)) end: # used by solve_semiregular # Input: a Laurent series a in Z + x k[[x]] # Output: a solution L of a*L + x*d/dx L = 0 sol_1order_eq:=proc(a,ext) local n; options remember; n:=-nmterms_laurent(a,1,0); # n := coefficient of x^0 in a if not type(n,integer) then bug() fi; new_laurent(n+1,x^n,[lift_sol1,a,n,ext]) end: # used by sol_1order_eq lift_sol1:=proc(LL,aa,n,ext,acc) local L,a,t,i; L:=value_laurent[LL]; a:=nmterms_laurent(aa,acc-n); for t from accuracy_laurent[LL] to acc-1 do # divide the coefficient of x^t in a*L by n-t: L:=L+modm(g_expand(convert([seq(coeff(L,x,i)*coeff(a,x,t-i) ,i=ldegree(L,x)..degree(L,x))],`+`)/(n-t)*x^t,ext)) od; L end: ############################################### # Solutions of matrix differential equations #-------------------------------- ############################################### diffop_matsol macro( solve_mateqn=`diffop/solve_mateqn`, matrix_mult=`diffop/matrix_mult`, left_solutions=`diffop/left_solutions`, xddx_log=`diffop/xddx_log`, lift_xddx_log=`diffop/lift_xddx_log`, left_sol_rational=`diffop/left_sol_rational`, solve_matrat=`diffop/solve_matrat` ): # A matrix # n: solve modulo x^n # 3rd argument: specify `x*d/dx` SolveMat:=proc(A,n) local v,i,j,res,ram,e,k; if n=rational then RETURN(solve_matrat(A)) fi; k:=linalg[rowdim](A); if not k=linalg[coldim](A) then ERROR(`wrong arguments`) fi; res:=NULL; # This map(..,A) yields a matrix without a name. Otherwise # options remember causes errors. # v:=solve_mateqn(map(i -> convert(i,RootOf),A),args[3..nargs]); # In the latter line options remember simply did not work, try again: v:=solve_mateqn([k,k,[seq(seq(convert(A[i,j],RootOf) ,j=1..k),i=1..k)]],args[3..nargs],`not a matrix`); for i in v do if nops(i)<4 then ram:=x else ram:=i[4] fi; #bug: e:=exp(int(i[2]/degree(ram,x)/x,x)); e:=exp(int(i[2]/x,x)); for j in i[1] do res:=res,[seq(subs(x=RootOf(subs(x=_Z,ram)-x),nt(j[k,1] ,n*degree(ram,x))*e),k=1..linalg[rowdim](A))] od od; subs(g_conversion2,[res]) end: # Gives the solution of (d/dx + A)(y) = 0 in internal format. # If D=`x*d/dx` then x*d/dx is used instead of d/dx # The entries of A must be rational functions in x, Laurent series are # not yet allowed in this version. solve_mateqn:=proc(A,D) local i,f,ext,cv_mat,fs,res,v,ex,ram,m,j,k; options remember; if args[nargs]=`not a matrix` then RETURN(solve_mateqn(linalg[matrix](op(A)),args[2..nargs-1])) fi; if nargs=1 or D<>`x*d/dx` then RETURN(solve_mateqn(map(i -> x*i,A),`x*d/dx`)) fi; f:=cyclic_vector(A,D); ext:=g_ext(f); f:=subs(g_conversion1,f); cv_mat:=map((i,ex) -> eval_laurent(i,ex),linalg[transpose](evalm( f[2])),ext); f:=eval_laurent(f[1],ext,xDF); # Now try to split the problem up using an LCLM factorization: fs:=factor_newton(f,`semireg`,ext); if nops(fs)=1 then fs:={[f,cv_mat]} else # this description_laurent[][11][1] is the left hand factor fs:={seq([description_laurent[coeff(i,xDF,0)][11][1] ,matrix_mult(cv_mat,i,ext)],i=fs)} fi; res:={}; for f in fs do v:=left_solutions(f[1],ext); for i in v do if nops(i)>2 then ex:=i[3];ram:=i[4] else ex:=ext;ram:=x fi; m:=map(ramification_laur,f[2],ram,ex); m:=linalg[matrix](linalg[rowdim](m),linalg[coldim](m), [seq(seq(eval_laurent(m[j,k]*degree(ram,x)^(k-1),ex) ,k=1..linalg[coldim](m)),j=1..linalg[rowdim](m))]); res:=res union {[[seq(matrix_mult(m,j,ex) ,j=i[1])],i[2],i[3..nops(i)]]} od od; test_result(solve_mateqn,A,res,ext); res end: # r is a right factor of f in D/D*f, with it corresponds a # basis of a subspace r, xDF*r, xDF^2*r, .. and a basis # transition. This matrix transition is multiplied by m. # If a 4th argument is given then DF is used instead of xDF. matrix_mult:=proc(M,r,ext) local m,mm,i,j,d,D; if nargs=4 then D:=DF else D:=xDF fi; m:=linalg[coldim](M); mm:=m-degree(r,D); m:=linalg[multiply](M,linalg[matrix](m,mm,[seq([seq( coeff(mult(seq(D,d=2..j),r,ext),D,i-1),j=1..mm)],i=1..m)])); if nargs=4 then m else map((i,j) -> eval_laurent(i,j),m,ext) fi end: # Input: an operator f with only 1 slope. # Now consider the ring K[x*d/dx]/(K[x*d/dx] * f) # where K is an extension field of k((x)) that contains all solutions # of adjoint(f)(y)=0. # A basis of this ring is: [ (x*d/dx)^0, (x*d/dx)^1, .. ] # so elements of this ring look like polynomials in xDF (is x*d/dx) of # degree < n. (n is order of f). # The output of this procedure is a list of elements of K[x*d/dx]/(K[x*d/dx] * f) # that become 0 when multiplied by x*d/dx. left_solutions:=proc(f,ext) local i,adj_sol,sol,ff,ex,ram,g,res,tres,j; options remember; adj_sol:=formal_sol(adjoint(f,ext),ext); tres:=NULL; for sol in adj_sol do # sol[1] is a list of solutions # sol[2] is the exponential part, this appears in the ExpInt part # sol[3..] are the algebraic extensions if nops(sol)=2 then ram:=x; ex:=ext else ram:=sol[4]; ex:=sol[3] fi; ff:=ramification_of(f,ram,ex); res:=NULL; for i in sol[1] do # Now we must solve the following problem: give an operator g # such that xDF*g is zero mod ff. Use: g=i*xDF^(n-1)+.. where # i is a solution of the adjoint of ff. g:=i*xDF^(degree(f,xDF)-1); for j from degree(f,xDF)-2 by -1 to 0 do g:=g+eval_laurent(expand(i*coeff(ff,xDF,j+1)- xddx_log(coeff(g,xDF,j+1),ex) -coeff(g,xDF,j+1)*eval_laurent(sol[2],ex) ),ex,log(x))*xDF^j od; test_result(left_solutions,ExpInt(normal(sol[2]),x)*g ,ff,ram,ex); res:=res,g od; tres:=tres,[[res],sol[2..nops(sol)]] od; [tres] end: # Input: a polynomial in log(x) with laurent series coefficients. # Output: x times the derivative. xddx_log:=proc(f,ext) local m,d,i; m:=lowerbound_val(f); d:=degree(f,log(x)); m:=convert([seq(new_laurent(m,0,[lift_xddx_log,i,d,f,ext])*log(x)^i ,i=0..d)],`+`); test_result(xddx_log,f,ext,m); m end: lift_xddx_log:=proc(L,i,d,f,ext,acc) local a,res; a:=acc,accuracy_laurent[L]; res:=value_laurent[L]; if ilinalg[coldim](A) then ERROR(`wrong arguments`) fi; f:=cyclic_vector(A); ext:=g_ext(f); f:=subs(g_conversion1,f); m:=map(i -> subs(g_conversion2,i),linalg[transpose](evalm(f[2]))); v:=left_sol_rational(f[1],ext); res:=[seq(matrix_mult(m,i,ext,`use DF`),i=v)]; # Now convert this to lists: res:=[seq([seq(i[j,1],j=1..n)],i=res)]; test_result(solve_mateqn,A,res,ext); res end: # Computes g in k(x)[DF]/( k(x)[DF]f ) such that DF*g is zero mod f. left_sol_rational:=proc(f,ext) local i,sol,g,res,j; options remember; sol:=rational_solutions(adjoint(f,ext),ext); res:=NULL; for i in sol do # Now we must solve the following problem: give an operator g # such that DF*g is zero mod f. Use: g=i*DF^(n-1)+.. where # i is a solution of the adjoint of f. g:=i*DF^(degree(f,DF)-1); for j from degree(f,DF)-2 by -1 to 0 do g:=g+g_normal(i*coeff(f,DF,j+1)-diff(coeff(g,DF,j+1) ,x))*DF^j od; test_result(left_sol_rational,g,f,ext); res:=res,g od; [res] end: ############################# # debugging tools #------------------------------------------------- ############################# diffop_debug # Not included. Can be obtained by e-mail request. ############# # help text #---------------------------------------------------------------- ############# diffop_help `help/text/SolveMat`:=TEXT( ``, `FUNCTION: SolveMat - solve matrix differential equation`, ``, `CALLING SEQUENCES:`, ` SolveMat(A,n)`, ` SolveMat(A,n,``x*d/dx``)`, ` SolveMat(A,rational)`, ``, `PARAMETERS:`, ` A - a square matrix of rational functions in x`, ` n - an integer`, ``, `SYNOPSIS:`, `- The function SolveMat gives formal solutions of matrix differential`, ` equations. These are equations of the following form: A*y + diff(y,x) = 0`, ` or A*y + x*diff(y,x) = 0 if the option ``x*d/dx`` is used. Here A is a square`, ` matrix and y is a vector. The solutions are of the form: y=[y1,y2,...,ym] if`, ` A is a m by m matrix. Here y1,..ym are of the form:`, ` exp(int(element of 1/x*Q_bar[1/x],x)) times a polynomial in log(x) with`, ` coefficients Q_bar((x)) where Q_bar is the algebraic closure of Q.`, ` Or a similar expression with x replaced by x^(1/r) for some integer r.`, ``, ` SolveMat returns all solutions up to conjugation over k((x)) where k is the`, ` coefficients field. So if the equation is given over Q and`, ` exp(RootOf(_Z^2-2)/x^3)*[x,x^2] and exp(-RootOf(_Z^2-2)/x^3)*[x,x^2] are two`, ` solutions then only one of these two solutions is given. The expressions in`, ` Q_bar((x)) in the solution are given modulo x^n, i.e. the coefficient of`, ` x^(n-1) is computed but not the coefficient of x^n.`, ``, ` SolveMat(A,rational) computes all rational functions solutions of matrix`, ` differential equations. The dimension of this solution space may`, ` be smaller than the dimension of A. Currently the implementation of`, ` this option is rather slow. This may be improved in a future version`, ` of diffop.`, ``, `EXAMPLES:`, `> A:=linalg[matrix](2,2,[1,2,x,x^2]);`, ``, ` [ 1 2 ]`, ` [ ]`, ` A := [ 2 ]`, ` [ x x ]`, `> SolveMat(A,5);`, ``, ` 71 4 43 43 2 29 3 2 25 3 43 4`, `[[--- x + ---- x - ---- x + --- x - 7/36, - 1/2 + 7/72 x - --- x + --- x ]`, ` 864 36 72 216 108 288`, ``, ` ,`, ``, ` 4 2 3 2 3 4`, ` [5/12 x + 2 x - x - 1/3 x - 2, x - 2/3 x + 1/4 x ]]`, ``, `> convert( evalm(A &* "[1] + diff("[1],x)) ,list);`, ``, ` 329 4 43 5 25 4 43 6`, ` [ --- x , - --- x + --- x + --- x ]`, ` 864 288 108 288`, ``, `# Now try if the result gets more accurate if we compute more terms:`, `> v:=SolveMat(A,20):`, `> convert( evalm(A &* v[1] + diff(v[1],x)) ,list);`, ``, ` 33205214213 19`, ` [ - ------------------ x ,`, ` 162193467211776000`, ``, ` 4395898231 20 45856249 21 139185479 19`, ` ----------------- x - --------------- x - --------------- x ]`, ` 20274183401472000 519850856448000 569099884953600`, ``, `# zero modulo x^19`, ``, `> A:=linalg[matrix](2,2,[0,1/x^2,1/x,1/x]);`, `> SolveMat(A,1,``x*d/dx``):`, `> convert(",radical);`, ``, ` / 2 31 449 179807 23921791 1/2\`, `[[|- ---- - ---- - --------- - ------ + --------- x |`, ` | 3/2 32 x 1/2 786432 201326592 |`, ` \ x 4096 x /`, ``, ` 1/2 1 1 1`, ` exp(5/4 ln(x ) - ------ + --- - ------),`, ` 1/2 4 x 3/2`, ` 8 x 3 x`, ``, ` / 1 511 161281 1/2\ 1/2 1 1 1`, ` |2/x - ------- - ---- - ------ x | exp(5/4 ln(x ) - ------ + --- - ------)`, ` | 1/2 4096 786432 | 1/2 4 x 3/2`, ` \ 32 x / 8 x 3 x`, ``, ` ]]`, `# the conjugate of this expression (replace x^(1/2) by -x^(1/2)) is also`, `# a solution`, ``, ``, `> A:=linalg[matrix](2,2,[(4*x^4-8*x^3+11*x^2-12*x-3)/(x^5+6*x^3-2*x^4-8*x^2+`, ` 9*x-6),-3*x/(2*x^3+3*x^2+6*x+9),(-2*x^2-9*x-9)/(x^2+5),(2*x^5+8*x^4-17*x^3+`, ` 21*x^2-61*x-13)/(2*x^6-x^5+10*x^4-6*x^2+25*x-30)]);`, ``, `A :=`, ``, ` [ 4 3 2 ]`, ` [ 4 x - 8 x + 11 x - 12 x - 3 x ]`, ` [ --------------------------------- - 3 --------------------- ]`, ` [ 5 3 4 2 3 2 ]`, ` [ x + 6 x - 2 x - 8 x + 9 x - 6 2 x + 3 x + 6 x + 9 ]`, ` [ ]`, ` [ 2 5 4 3 2 ]`, ` [ - 2 x - 9 x - 9 2 x + 8 x - 17 x + 21 x - 61 x - 13 ]`, ` [ ---------------- --------------------------------------- ]`, ` [ 2 6 5 4 2 ]`, ` [ x + 5 2 x - x + 10 x - 6 x + 25 x - 30 ]`, ``, `# The corresponding system of differential equations is:`, `> v:={A[1,1]*y1(x)+A[1,2]*y2(x)+diff(y1(x),x)=0,`, ` A[2,1]*y1(x)+A[2,2]*y2(x)+diff(y2(x),x)=0};`, ``, `> sols:=SolveMat(A,rational);`, ``, ` x - 1 (x - 1) (2 x + 3)`, ` sols := [[- ----------, - -----------------]]`, ` 2 2`, ` x - x + 2 x - x + 2`, ``, `# Apparantly the solution space over Q(x) is one dimensional. Now check if`, `# the answer is correct:`, `> sol:=sols[1]:`, `> simplify( subs({ y1(x)=sol[1], y2(x)=sol[2]}, v));`, ``, ` {0 = 0}`, ``, `SEE ALSO: diffop, dsolve` ): `help/text/mult`:=TEXT(``, `FUNCTION: mult - Multiplication of differential operators`,``, `CALLING SEQUENCE:`,` mult(f1,f2,..,fn)`, ``,`PARAMETERS:`, ` f1 .. fn - differential operators`,``, `SYNOPSIS:`, `- This procedure computes products in the ring k(x)[DF]. The multiplication`, ` in this ring is defined by DF*a=a'+a*DF for all a in k(x).`,``, `EXAMPLE:`, `> mult(DF,x*DF);`, ` 2`, ` DF + x DF`, `` ): `help/text/diffop`:=TEXT( ``, `HELP FOR: syntax for the package diffop`, ``, ` The name diffop stands for differential operators. The main goal of this`, ` package is factorization of differential operators. A differential operator f`, ` in this package is an element of k(x)[DF] or k((x))[xDF]. Here k is the`, ` constants field, k must be Q or an algebraic extension of Q.`, ``, ` There is also an algorithm in diffop for computing formal solutions of matrix`, ` differential equations with coefficients in k(x). This algorithm uses`, ` factorization in k((x))[xDF] and a cyclic vector.`, ``, ` The symbol DF stands for differentiation w.r.t. x, and the operator x stands`, ` for multiplication by x. So DF is d/dx. Multiplication of differential`, ` operators is done by the procedure mult. We have mult(DF,x) = x * DF + 1.`, ``, ` The symbol xDF stands for x * DF. The symbol xDF is used to represent local`, ` operators (i.e. differential operators with elements of k((x)) as`, ` coefficients). The symbol DF is used to represent global operators, i.e.`, ` coefficients in k(x). We have xDF = x * d/dx and mult(xDF,x) = x * xDF + x.`, ``, ` Since in Maple DF*x gives the same result as x*DF we need to use a convention`, ` for the syntax. We will always assume that a monomial f=DF^a*x^b stands for`, ` x^b*DF^a. This representation can be obtained by sort(f,[x,DF])`, ``, `SEE ALSO: DFactor, convert[diffop], SolveMat, rational_solutions, adjoint`, ` exterior_power, symmetric_power, LCLM, GCRD, rightdivision,`, ` leftdivision, cyclic_vector, eval_laurent, nt, formal_sol,`, ` eigenring, mult, exp_solutions` ): has_liouvillian_solutions:=proc(f) if degree(f,DF)<>2 or indets(f) minus {x,DF} <> {} then ERROR(`wrong arguments`) fi; evalb(nops(DFactor(symmetric_power(f,6),`one step`))>1) end: ##################################### # Pade Hermite approximations #----------------------------------------- ##################################### diffop_pade2 # The following is a variant on Harm Derksen's (hderksen@sci.kun.nl) algorithm # for computing Pade expansions. You can find his algorithm in the share # library under the name pade2. G. Labahn and him found indepent of each other # a good method for computing these expansions. # I made the following changes: # 1) This pade2 accepts my syntax for elements of k((x)) (laur1, laur2 etc.) # 2) Harm pointed out to me that pade2 can also work with vectors of power # series. I replaced the variable z in his pade2 by z[1], z[2], etc. to # achieve this. # 3) point is always x=0 # 4) ext is the coefficients field # Note that it is useful to have the original pade2 in order to understand # how the method works. macro( lift_pade2=`diffop/lift_pade2`, pade2=`diffop/pade2`, smaller_pade2=`diffop/smaller_pade2`, smallest_pade2=`diffop/smallest_pade2`, valuation_pade2=`diffop/valuation_pade2`, wipe_pade2=`diffop/wipe_pade2` ): # functionlist is a list of lists of elements in k((x)) # pade2 looks for a k[x]-linear relation between these lists (interpret these # lists as vectors). If all these lists contain only 1 element, then this is # basically Harm Derksen's pade2 algorithm. pade2:=proc(functionlist,accuracy,ext) local i,j,n,y,z,yvars,zvars,fl,appr,degrees; n:=min(seq(seq(valuation_laurent(j,0),j=i),i=functionlist)); if n<0 then RETURN(pade2([seq([seq(eval_laurent(x^(-n)*j) ,j=i)],i=functionlist)],accuracy,ext)) fi; n:=nops(functionlist); degrees:=[seq(0,i=1..n)]; # Not used yet, for later generalizations zvars:=[seq(z[i],i=1..nops(functionlist[1]))]; yvars:=[seq(y[i],i=1..n)]; fl:=0,[seq(convert([seq(functionlist[i][j]*z[j],j=1..nops(zvars))] ,`+`),i=1..n)],n,degrees,yvars,zvars,ext,0; if subs(infinity=0,accuracy)<>accuracy then RETURN([fl]) fi; appr:=lift_pade2(fl,accuracy); smallest_pade2(op(appr)) end: # This procedure lifts the list appr # Output: new value of appr # appr, acc, n and degrees is the same as in Harm Derksen's pade2 # functionlist: same as Harm's but multiplied by z # y is a list of variables [y[1],y[2],..] are the same as in Harm's algorithm # zvars is a list of variables [z[1],z[2],..] which play the role of Harm's # variable z. # old_acc: 0 on the first call, and the previous acc otherwise # appr: 0 on the first call, and the previous result of lift_pade2 otherwise lift_pade2:=proc(appr,functionlist,n,degrees,y,zvars,ext,old_acc,extra_acc) local app,i,k,j,z,acc; userinfo(9,diffop,`lifting`,extra_acc,`steps`); acc:=old_acc+extra_acc; if old_acc=0 then # Since I can't change appr, I'll call it app here. app:=[seq(expand(x^degrees[i]*y[i]+nt( functionlist[i],acc)),i=1..n)] else app:=map(g_expand,subs({seq(y[i]=y[i]+x^(-degrees[i]) *(nt(functionlist[i],acc)-nt(functionlist[i],old_acc)) ,i=1..n)},appr),ext) fi; app:=modm(app); for i from old_acc to acc-1 do for z in zvars do k:=0; for j to n do if valuation_pade2(app[j],z,x,acc)=i and (k=0 or smaller_pade2(app[j],app[k],x,y)) then k:=j fi od; if k>0 then app:=[seq(wipe_pade2(app[j],app[k],z,x,i,ext),j=1..k-1), expand(x*app[k]),seq(wipe_pade2(app[j],app[k], z,x,i,ext),j=k+1..n)] fi od od; userinfo(9,diffop,`done lifting`); [app,functionlist,n,degrees,y,zvars,ext,acc] end: smallest_pade2:=proc(appr,dummy,n,degrees,y,zvars,ext,dummy2) local smallest,i; smallest:=1; appr; for i from 2 to n do if smaller_pade2(appr[i],appr[smallest],x,y) then smallest:=i; fi od; [seq(g_expand(coeff(appr[smallest],y[i],1)/x^degrees[i],ext),i=1..n)] end: valuation_pade2:=proc(f,z,y,acc) if coeff(f,z,1)=0 then acc else ldegree(coeff(f,z,1),y) fi end: wipe_pade2:=proc(f,g,z,y,i,ext) local ff,gg; ff:=coeff(coeff(f,z,1),y,i); gg:=g_normal(ff/coeff(coeff(g,z,1),y,i)); modm(g_expand(f-gg*g,ext)) end: smaller_pade2:=proc(f,g,x,vars) local ff,gg,i,df,dg,n; n:=nops(vars); ff:=[seq(degree(coeff(f,vars[i],1),x),i=1..n)];df:=max(op(ff)); gg:=[seq(degree(coeff(g,vars[i],1),x),i=1..n)];dg:=max(op(gg)); if dfdg then RETURN(false) fi; for i to n do if ff[i]=df and coeff(f,vars[i],1)<>0 then RETURN(false) fi; if gg[i]=dg and coeff(g,vars[i],1)<>0 then RETURN(true);fi; od; end: ####################### # mod p computation #----------------------------------------------------- ####################### diffop_modp macro( compute_modm=`diffop/compute_modm`, compute_modp=`diffop/compute_modp`, iratrecon=readlib('iratrecon'), modp_table=`diffop/modp_table`, modulus2=`diffop/modulus2`, reconstruct=`diffop/reconstruct` ): compute_modp:=proc() global accuracy_laurent, description_laurent, set_laurents, value_laurent, rem_lift, modp_table; local i,procedure,ac_old,de_old,va_old; if nargs=0 then # no modular computation accuracy_laurent:=table(modp_table[0][1]); description_laurent:=table(modp_table[0][2]); set_laurents:=modp_table[0][3]; value_laurent:=table(modp_table[0][4]); rem_lift:=table(modp_table[0][5]); RETURN() fi; userinfo(6,diffop,`computing modulo`,modulus); modp_table[0]:=[op(op(accuracy_laurent)),op(op( description_laurent)),set_laurents,op(op(value_laurent)) ,op(op(rem_lift))]; if modp_table[modulus]<>evaln(modp_table[modulus]) then accuracy_laurent:=table(modp_table[modulus][1]); description_laurent:=table(modp_table[modulus][2]); value_laurent:=table(modp_table[modulus][4]); rem_lift:=table(modp_table[modulus][5]); ac_old:=table(modp_table[0][1]); de_old:=table(modp_table[0][2]); va_old:=table(modp_table[0][4]); for i in set_laurents minus modp_table[modulus][3] do # accuracy_laurent[i]:=subs(modp_table[0][1],i); # description_laurent[i]:=subs(modp_table[0][2],i); # value_laurent[i]:=subs(modp_table[0][4],i) # these 3 lines were awfully slow, replacement: accuracy_laurent[i]:=ac_old[i]; description_laurent[i]:=de_old[i]; value_laurent[i]:=va_old[i] od; for i in modp_table[0][5] do if rem_lift[op(i)[1]] =evaln(rem_lift[op(i)[1]]) then rem_lift[op(i)[1]]:=op(i)[2] fi od else compute_modp() fi; procedure:=args[1]; i:=traperror(procedure(args[2..nargs])); if i=lasterror then userinfo(2,diffop,modulus,`gives error`,i); compute_modp(); RETURN(`new prime needed`) else modp_table[modulus]:=[op(op(accuracy_laurent)),op(op( description_laurent)),set_laurents,op(op(value_laurent)) ,op(op(rem_lift))] fi; compute_modp(); i end: # Compute modulo prime power. compute_modm:=proc(exponent) global modulus, modulus2; local v; modulus:=modulus2^exponent; v:=compute_modp(args[2..nargs]); while v=`new prime needed` do modulus2:=nextprime(modulus2*5); modulus:=modulus2^exponent; v:=compute_modp(args[2..nargs]) od; modulus:=0; v end: # Input: v is a list of polynomials, with coefficients in Z/(modl Z). # Output: a corresponding list with rational number coefficients. reconstruct:=proc(v,modl,N) local w,i,t,r; if nargs=2 then reconstruct(args,floor(evalf(sqrt(modl/5)))) # 5 instead of 2, because that way we can be almost sure that if # an answer is returned, it will be the right answer. That way we # don't compute a rightdivision in vain. elif type(v,list) then w:=NULL; for i in v do w:=w,reconstruct(i,modl,N); if w[nops([w])]=failed then RETURN(failed) fi od; [w] elif indets(v,`name`)<>{} then t:=[op(indets(v,`name`))]; t:=t[1]; w:=0; for i from ldegree(v,t) to degree(v,t) do r:=reconstruct(coeff(v,t,i),modl,N); if r=failed then RETURN(r) fi; w:=w+r*t^i od; w elif iratrecon(v,modl,N,N,'r','t') then eval(r/t) else failed fi end: ####################### # user interface #----------------------------------------------------- ####################### print_sol:=proc(v,n) local i,al,ra,ex,j; for i in v do print(`--------------------------------------------------------`); if nops(i)=2 then al:=[];ra:=x else al:=convert(subs(g_conversion2,i[3]),radical); ra:=convert(subs(g_conversion2,i[4]),radical) fi; if ra=x then ex:=i[2] else ex:=subs(x=RootOf(x-subs(x=_Z,ra)),i[2]) fi; ex:=convert(subs(g_conversion2,ex),radical)/degree(ra,x); # print(`algebraic extensions over the constants field: `,al); # print(`ramification is `,ra); print(`exponential part: `, ex); print(`Corresponding solutions are `,exp(Int(ex/`x`,x)),`times`); for j in i[1] do print(convert(subs(x=RootOf(x-subs(x=_Z,ra)) ,g_conversion2,nt(j,degree(ra,x)*n))+O(x^n),radical)) od od:NULL end: `help/text/formal_sol`:= TEXT(``,`FUNCTION: formal_sol - Compute the formal solutions of a di\ fferential operator`,``,`CALLING SEQUENCE:`,` formal_sol(f)`,``, `PARAMETERS:`,` f - a differential operator`,``,`SYNOPSIS:`, `- This procedure computes a basis of formal solutions of the differential`, ` operator f. These solutions are computed only up to conjugation over k((x))` ,` where k is the field of constants in f (i.e. k is Q extended with the`, ` RootOf's appearing in f). To obtain a full basis of solutions compute`, ` all conjugates over k((x)) of the output of formal_sol.`,``, `- The output of this procedure is rather technical. To view the solutions`, ` issue the commands:`,``,` v:=formal_sol(f):`,` print_sol(v,n);`, ``,` Then the solutions are printed in a clearer way. Here n is the accura\ cy. The`, ` power series appearing in the solutions are evaluated modulo x^n (so x^n`, ` is the first term that is not computed).`,``, `- This procedure is written for local differential operators in xDF syntax.`, ` If the operator f is in DF syntax then f will first be localized at x=0 wit\ h`, ` convert(f,diffop,x=0), and then formal_sol with return the formal solutions` ,` of this local operator. So then the formal solutions at the point x=0 are`, ` given.`,``, `- See ?convert,diffop how to localize a differential operator in DF syntax in` ,` a point p. To compute the formal solutions at the point x=p issue the`, ` command: formal_sol( convert(f,diffop,x=p) );`, ` Note that the output will be power series in x instead of x-p. So to obtain` ,` the correct answer x=x-p should be substituted in the output of print_sol.` ,``,`EXAMPLE:`,`f:=DF^2-1/x^3;`, `v:=formal_sol(f); # Note: only 1 solution is given, to get a basis of`,` \ # solutions take the conjugate (replace x^(1/2) by -x^(1/2))`, ` # as well.`,``,`print_sol(v,5);`,``, ` 1`, ` exponential part: , ---- + 3/4`, ` 1/2`, ` x`,``, ` 1`, ` ---- + 3/4`, ` / 1/2`, ` | x`, ` Corresponding solutions are , exp( | ---------- dx), times`, ` | x`, ` /`,``, ` 1/2 15 105 3/2 4725 2 72765 5/2 2837835 3` ,` 1 + 3/16 x - --- x + ---- x - ------ x + ------- x - --------- x` ,` 512 8192 524288 8388608 268435456`, ``,` 66891825 7/2 14783093325 4 468131288625 9/2 5`, ` + ---------- x - ------------ x + ------------- x + O(x )`, ` 4294967296 549755813888 8796093022208`,``,` `): ######################## # Various tools #----------------------------------------- ######################## diffop_tools # Input: f in Qbar(x)[DF] (in RootOf syntax) # Output: g such that f(y)=0 -> g(y^n)=0 symmetric_power:=proc(f,n) local i,zero,a,y,ypower,k,subs_rem,res,co; ypower:=y(x)^n; zero:=a[0]*ypower; subs_rem:=diff(y(x),seq(x,k=1..degree(f,DF)))=-convert([y(x)* coeff(f,DF,0)/lcoeff(f,DF),seq(diff(y(x),seq(x,k=1..i))* coeff(f,DF,i)/lcoeff(f,DF),i=1..degree(f,DF)-1)],`+`); co:=combinat[numbcomp](degree(f,DF)+n,degree(f,DF)); `diffop/_symext`(f,co,subs_rem,a,ypower,zero) end: symmetric_product:=proc(f,g) local i,zero,a,y,ypower,k,subs_rem,co,j; global DF,x; ypower:=y1(x)*y2(x); zero:=a[0]*ypower; subs_rem:={seq(diff(y.j(x),seq(x,k=1..degree(args[j],DF)))=-convert([y.j(x)* coeff(args[j],DF,0)/lcoeff(args[j],DF),seq(diff(y.j(x),seq(x,k=1..i))* coeff(args[j],DF,i)/lcoeff(args[j],DF),i=1..degree(args[j],DF)-1)],`+`),j=1..2)}; co:=degree(f,DF)*degree(g,DF); `diffop/_symext`([f,g],co,subs_rem,a,ypower,zero) end: `diffop/_symext`:=proc(f,co,subs_rem,a) local i,z,zero,res,k,vars; z:=args[5]; zero:=args[6]; vars:={seq(a[k],k=0..co)}; for i from 1 to co do z:=g_normal(subs(subs_rem,diff(z,x))); zero:=expand(zero+a[i]*z) od; res:=subs(solve({coeffs(zero,indets(zero) minus {x,op(vars), op(indets(f))})},vars),convert([seq(a[k]*DF^k,k=0..co)],`+`)); z:=nops(indets(res) intersect vars)-1; collect(subs(solve({coeff(res,DF,co-z)=1,seq(coeff(res,DF,i)=0 ,i=co-z+1..co)},vars),res),DF,g_normal) end: # Also called: associated equation exterior_power:=proc(f,d) local a,z,i,j,subs_rem,k,co,zero,y; if d<=0 then RETURN(`undefined`) elif d>degree(f,DF) then RETURN(1) fi; z:=linalg[det](linalg[matrix](d,d,[ seq(y[i](x),i=1..d),seq(seq(diff(y[i](x),x$j),i=1..d),j=1..d-1)])); subs_rem:=[seq( diff(y[j](x),seq(x,k=1..degree(f,DF)))=-convert([y[j](x)* coeff(f,DF,0)/lcoeff(f,DF),seq(diff(y[j](x),seq(x,k=1..i))* coeff(f,DF,i)/lcoeff(f,DF),i=1..degree(f,DF)-1)],`+`),j=1..d)]; co:=binomial(degree(f,DF),d); zero:=a[0]*z; `diffop/_symext`(f,co,subs_rem,a,z,zero) end: # Input: a differential operator d/dx + A where A is a matrix # Output: list with: matrix with columns v,(d/dx+A)v,(d/dx+A)^2v,.. # and an operator. # Specifying the cyclic vector is possible using a third argument # If the second argument is `x*d/dx` then x*d/dx is used as a derivation # instead of d/dx. cyclic_vector:=proc(A,D,vv) local v,n,i,M,m,df,MM1,MM2,AA1,AA2; if nargs=1 then RETURN(cyclic_vector(A,`d/dx`)) fi; n:=linalg[rowdim](A); # A is an n by n matrix if n<>linalg[coldim](A) then ERROR(`wrong input`) fi; if nargs<3 then v:=[1,seq(0,i=2..n)] else v:=vv fi; M:=NULL; for i from 1 to n do if i=1 then v:=evalm(v) else v:=`diffop/apply_A`(A,v,D) fi; M:=M,convert(v,list) od; m:=linalg[transpose](evalm([M])); m:=linalg[linsolve](m,`diffop/apply_A`(A,v,D)); # if indets([op(op(m))]) minus indets([x,v,op(op(A))]) <> {} then # indets doesn't work here, so try something else: MM2:=evalm(m);MM1:=op(MM2);MM2:=[op(op(MM1))];MM1:=indets(MM2,`name`); AA2:=evalm(A);AA1:=op(AA2);AA2:=[op(op(AA1))]; # Now indets should work: AA1:=indets(AA2,`name`); if MM1 minus indets([x,args[2..nargs],AA1]) <> {} then if nargs=2 then v:=[seq(randpoly(x,degree=0,coeffs=rand(-1..1)) ,i=1..n)] else v:=[seq(randpoly(x,degree=degree(vv[i],x)+1 ,coeffs=rand(-1..1)),i=1..n)] fi; RETURN(cyclic_vector(A,D,v)) fi; if D=`x*d/dx` then df:=xDF else df:=DF fi; [convert([df^n,seq(normal(-m[i+1])*df^i,i=0..n-1)],`+`),[M]] end: # used by cyclic_vector `diffop/apply_A`:=proc(A,v,D) local a; if D=`x*d/dx` then map(normal,evalm(linalg[multiply](A,v)+map(a -> x*diff(a,x),v))) else map(normal,evalm(linalg[multiply](A,v)+map(a -> diff(a,x),v))) fi end: adjoint:=proc(f) local i,ext,v,D; options remember; userinfo(7,diffop,`computing the adjoint of`,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); convert([xDF^degree(f,xDF),seq(new_laurent(v,0, [`diffop/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_expand(v,ext) else collect(v,D,g_normal) fi fi 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. `diffop/lift_adjoint`:=proc(L,n,f,ext,acc) local a,res,i; 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: # Compute the rational solutions of a differential equation by solving # linear equations. # If a third argument g is given then only the singularities of g will # be considered, not the singularities of f. rational_solutions:=proc(f,ext) local v,i,a,ansatz,bound,den,j,k,m,np,s,singularities,sol,t,zero,var,ff; if nargs=1 or not type(ext,list) then v:=g_ext(f); RETURN(subs(g_conversion2, rational_solutions(subs(g_conversion1,f),v,args[2..nargs]))) fi; if nargs>2 then den:=args[3] else den:=f fi; singularities:={seq([i[1],convert(f,diffop,x=i[1]),[op(i[3]),op(ext)]] ,i=v_ext_m(g_factors(denom(g_normal(den/lcoeff(den,DF))),x,ext),x))}; s:=0; den:=1; for k in singularities do np:=Newtonpolygon(k[2],`include the Newton polynomials`); if np[1][3]<>0 then RETURN([]) fi; v:=factors(collect(np[1][4]/lcoeff(np[1][4],x),x,g_normal))[2]; m:=-infinity; for i in v do if type(i[1]-x,integer) then m:=max(m,i[1]-x) fi od; if m=-infinity then RETURN([]) fi; t:=m/(x-k[1]); if k[3]<>ext then t:=subs(g_conversion1,evala(Trace(op(subs( g_conversion2,[t,k[3],ext]))))) fi; s:=s+t; den:=den/denom(g_normal(t))^m od; s:=g_normal(s); ff:=substitute(DF-s,f,ext); k:=convert(ff,diffop,x=infinity); np:=Newtonpolygon(k,`include the Newton polynomials`); if np[1][3]<>0 then RETURN([]) fi; v:=factors(collect(np[1][4]/lcoeff(np[1][4],x),x,g_normal))[2]; m:=-infinity; for i in v do if type(i[1]-x,integer) then m:=max(m,i[1]-x) fi od; if m=-infinity then RETURN([]) fi; bound:=m; ansatz:=convert([seq(a[i]*x^i,i=0..bound)],`+`); zero:=ansatz*coeff(ff,DF,0); t:=ansatz; for i from 1 to degree(ff,DF) do t:=diff(t,x); zero:=zero+coeff(ff,DF,i)*t od; t:=solve({coeffs(collect(numer(subs(g_conversion2,zero)),x),x)}); sol:=subs(t,ansatz); if sol=0 then RETURN([]) fi; var:={seq(a[i],i=0..bound)} intersect indets(sol); i:=[seq(g_normal(subs({seq(i=0,i=var minus {j}),j=1},sol*den)),j=var)]; test_result(rational_solutions,f,i,ext); i end: GCRD:=proc(f,g) local a; if f=0 or g=0 then f+g else a:=collect(g/lcoeff(g,DF),DF,g_normal); procname(a,rightdivision(f,a)[2]) fi end: # 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; if nargs=4 then # 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`,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; 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; if nargs=2 then a:=mult(L,1/lcoeff(L,DF)); b:=g_ext([args]); a:=rightdivision(op(subs(g_conversion1,[f,a])),b); RETURN(subs(g_conversion2,[a[1]/lcoeff(L,DF),a[2]])) fi; userinfo(5,diffop,`left division`,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: # 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() local ext,j,i,ind,sol,ansatz,a,r; ext:=g_ext([args]); if has([args],`and conjugates`) then if not type(args[nargs],list) then # ground field is Q RETURN(LCLM(args,[])) fi; if nargs>3 then RETURN(LCLM(seq(LCLM(args[j],args[nargs-1..nargs]),j=1..nargs-2))) fi; r:=args[1]; j:=subs(g_conversion1,args[3]); ansatz:=convert([seq(a[i]*DF^i,i=0..degree(r,DF) *degree_ext(ext,j))],`+`); sol:={seq(coeffs(i,DF),i=[ext_to_coeffs(g_expand(subs(g_conversion1, numer(rightdivision(ansatz,r)[2])),ext),j)])} else if type(args[nargs],list) then # ext should not be given RETURN(LCLM(args[1..nargs-1])) elif nargs=1 then RETURN(args) fi; ansatz:=convert([seq(a[i]*DF^i,i=0..convert([seq(degree(j,DF),j=args)] ,`+`))],`+`); sol:={seq(coeffs(g_expand(numer(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: `help/text/rational_solutions`:=TEXT(``, `FUNCTION: rational_solutions - compute a basis of rational solutions`,``, `CALLING SEQUENCE:`,` rational_solutions(f);`,``,`PARAMETERS:`, ` f - a differential operator`,``,`SYNOPSIS:`, `- The input f is a differential operator in DF notation, i.e. it is a`, ` polynomial in DF with rational functions in x as coefficients. The output`, ` of this procedure is a basis of the solutions in k(x) where k is the field`, ` of constants.`,``,`EXAMPLE:`,``, ` f:=DF^2+2*(185808*x^9+905961*x^8+4297020*x^7+3585946*x^6+1705844*x^5`, ` -5824272*x^4-4172916*x^3-3850014*x^2-423972*x-404928)/x/(56*x^3+49*x^2`, ` +63*x+57)^2/(55*x^5+37*x^4+35*x^3-97*x^2-50*x-79):`,``, ` rational_solutions(f);`, ` 5 4 3 2`, ` x (55 x + 37 x + 35 x - 97 x - 50 x - 79)`, ` [1/35 ---------------------------------------------]`, ` 3 2 2`, ` (56 x + 49 x + 63 x + 57)`): `help/text/exterior_power`:=TEXT(``, `FUNCTION: exterior_power - compute the n'th exterior_power of an operator`, ``,`CALLING SEQUENCE:`,` exterior_power(f,n);`,``,`PARAMETERS:`, ` f - a differential operator`,` n - an integer`,``,`SYNOPSIS:`, `- The input f is a differential operator in DF notation, i.e. it is a`, ` polynomial in DF with rational functions in x as coefficients. The output`, ` of this procedure is an operator g of minimal order with rational functions` ,` coefficients such that for all solutions y1,..,yn of f the Wronskian`, ` w = det(matrix(n,n,[y1,y1',y1'',.. , y2,y2',y2'', .. , yn,yn',yn'',..]);`, ` is a solution of g.`): `help/text/symmetric_power`:=TEXT(``, `FUNCTION: symmetric_power - compute the n'th symmetric_power of an operator`, ``,`CALLING SEQUENCE:`,` symmetric_power(f,n);`,``,`PARAMETERS:`, ` f - a differential operator`,` n - an integer`,``,`SYNOPSIS:`, `- The input f is a differential operator in DF notation, i.e. it is a`, ` polynomial in DF with rational functions in x as coefficients. The output`, ` of this procedure is an operator g of minimal order with rational functions` ,` coefficients such that for all solutions y1,..,yn of f y1*y2*..*yn is a`, ` solution of g.`, ``, `EXAMPLES:`, ``, ` Using the work of Singer and Ulmer (cf. JSC 1993, p. 9-36) one can apply`, ` DFactor to obtain information about the differential Galois group. The`, ` following procedure implements proposition 4.4, page 25, which decides if a`, ` second order differential equation has nonzero Liouvillian solutions:`,``, ` has_liouvillian_solutions:=proc(f)`, ` if degree(f,DF)<>2 or indets(f) minus {x,DF} <> {} then`, ` ERROR(``wrong arguments``)`,` fi;`, ` evalb(nops(DFactor(symmetric_power(f,6),``one step``))>1)`,` end;`,` ` ,` eqn:=diff(y(x),x,x)-(3/(16*x*(x-1))-3/(16*x^2)-3/(16*(x-1)^2))*y(x)=0;`, ``,` / 2 `, ` | d | / 3 3 3 `, ` eqn := |----- y(x)| - |------------ - ----- - -----------| y(x) = 0`, ` | 2 | |16 x (x - 1) 2 2|`, ` dx / 16 x 16 (x - 1) /`,``, ` f:=convert(eqn,diffop);`, ` # must be an equation in y(x) and x, with = 0 as righthand side`,``, ` 2`, ` 2 x - x + 1`, ` f := DF + 3/16 -----------`, ` 2 2`, ` x (x - 1)`,``, ` has_liouvillian_solutions(f);`,``, ` true`,``, ` has_liouvillian_solutions(DF^2+21/100*(x^2-x+1)/x^2/(x-1)^2);`,``, ` true`,``,`SEE ALSO: diffop`): `help/text/adjoint`:= TEXT(``,`FUNCTION: adjoint - the adjoint of a differential operator`,``, `CALLING SEQUENCE:`,` adjoint(f);`,``,`PARAMETERS:`, ` f - a differential operator`,``,`SYNOPSIS:`, `- Denote k(x)[DF] as the non-commutative ring of differential operators with`, ` rational functions coefficients. Since the commutator [x , DF] equals the`, ` commutator [-DF , x] the map A: k(x)[DF] -> k(x)[DF] defined by DF -> -DF`, ` is an anti-homomorphism, i.e. A(a*b)=A(b)*A(a),`, ` from k(x)[DF] to itself. If`,` f = a_n * (d/dx)^n + ... + a_0 * (d/dx)\ ^0 = a_n DF^n + .. + a_0 DF^0 then`, ` A(f) =(-d/dx)^n * a_n + ... + (-d/dx)^0 * a_0.`,``, `- Now the adjoint of f is defined as A(f)*(-1)^n`,``,`EXAMPLE:`, `> DF^3+1/x*DF+1;`,``,` 3 DF`, ` DF + ---- + 1`, ` x`,``,`> adjoint(");`,``, ` 2`, ` 3 DF x + 1`, ` DF + ---- - ------`, ` x 2`, ` x`,``,`> adjoint(");`,``, ` 3 DF`, ` DF + ---- + 1`, ` x`): `help/text/LCLM`:=TEXT(``, `FUNCTION: LCLM - Least Common Left Multiple of differential operators`,``, `CALLING SEQUENCE:`,` LCLM(f1,f2,..,fn)`, ` LCLM(f1,f2,..,fn,``and conjugates``,ext)`,``,`PARAMETERS:`, ` f1 .. fn - differential operators`, ` ext - (optional) a list of RootOf's describing the base field`,``, `SYNOPSIS:`, `- The least common left multiple F=LCLM(f1 .. fn) of a number of operators is` ,` defined as the operator with minimal order such that all solutions`, ` f1 .. fn are solutions of F as well.`,``, `- If the optional two arguments ``and conjugates``,ext are given then`, ` LCLM( f1 .. fn and all their conjugates over the field Q(ext))`, ` is computed. This LCLM is an element of Q(ext,x)[DF]. If ext is not`, ` specified then ext = [], i.e. Q(ext) = Q.`,``,`EXAMPLE:`, ` a:=RootOf(x^2-2);`,` b:=RootOf(x^2-3);`, ` LCLM(DF+a,DF^2+b*DF+x); # Order 3`, ` LCLM(DF+a,DF^2+b*DF+x,``and conjugates``,[b]); # Order 4`, ` LCLM(DF+a,DF^2+b*DF+x,``and conjugates``,[a]); # Order 5`,` LCLM(DF+a,DF^2+\ b*DF+x,``and conjugates``); # Order 6, nice test for DFactor`, ` 5 4 3 2 6 7`, ` 12 x - 21 x - 41 x + 56 x + 81 - 36 x + 8 x + x`, ` - 2 -----------------------------------------------------`, ` %1`,``, ` 5 4 3 2`, ` (3 x + x - 60 x + 12 x + 217 x + 207) DF`, ` + 2 --------------------------------------------`, ` %1`,``, ` 7 6 5 4 3 2 2`, ` (x + 4 x - 14 x - 21 x + 87 x - 248 x - 214 x + 65) DF`, ` + -------------------------------------------------------------`, ` %1`,``, ` 5 4 3 2 3`, ` (3 x - 9 x - 124 x - 60 x + 297 x + 233) DF`, ` - ------------------------------------------------`, ` %1`,``, ` 5 4 3 2 6 4`, ` (11 x - 16 x - 88 x + 161 x + 24 + 146 x + 2 x ) DF`, ` + --------------------------------------------------------`, ` %1`,``, ` 4 3 2 5`, ` (5 x + 32 x + 36 x - 40 x - 13) DF 6`, ` - -------------------------------------- + DF`, ` %1`,``, ` 5 4 3 2`, `%1 := x + 8 x + 12 x - 13 x - 8 - 20 x`,` `): `help/text/rightdivision`:=TEXT(``, `FUNCTION: rightdivision - division of differential operators`,``, `CALLING SEQUENCE:`,` rightdivision(a,b)`,``,`PARAMETERS:`, ` a,b - differential operators`,``,`SYNOPSIS:`,`- For operators a and \ b there exist (like for polynomials, see ?quo) operators`, ` q and r such that a = q*b + r and order(r) < order(b).`, ` The output rightdivision is the list [q,r].`,``,`EXAMPLE:`, `a:=DF^3;`,`b:=DF^2+x*DF+3;`,`qr:=rightdivision(a,b);`,``, ` 2`, ` qr := [DF - x, (- 4 + x ) DF + 3 x]`,``, `expand( mult(qr[1],b)+qr[2] );`,``, ` 3`, ` DF`): `help/text/GCRD`:=TEXT(``, `FUNCTION: GCRD - greatest common right divisor of 2 differential operators`, ``, `CALLING SEQUENCE:`,` GCRD(f,g)`,``,`PARAMETERS:`, ` a,b - differential operators`,``, `SYNOPSIS:`, `- For operators f and g there exists an operator R such that the solution`, ` space of R is the intersection of the solution spaces of f and g. The`, ` output of this procedure is this operator R. Like the gcd for polynomials,`, ` this R is a combination R=a*f+b*g for some operators a and b.`,``): `help/text/leftdivision`:= TEXT(``,`FUNCTION: leftdivision - division of differential operators` ,``,`CALLING SEQUENCE:`,` leftdivision(a,b)`,``,`PARAMETERS:`, ` a,b - differential operators`,``,`SYNOPSIS:`, `- Same as rightdivision, except now a is not q*b + r but a = b*q + r`,``, `SEE ALSO: rightdivision`): `help/text/cyclic_vector`:=TEXT(``, `FUNCTION: cyclic_vector - cyclic vector for a matrix differential equation`, ``,`CALLING SEQUENCE:`,` cyclic_vector(A);`,``,`PARAMETERS:`, ` A - an n by n matrix with entries that are rational functions in x`,``, `SYNOPSIS:`, `- For an n by n matrix A with entries A_ij in k(x) (here k is the field of`, ` constants) there exists a vector v in k(x)^n for which the vectors`,``, ` ((d/dx + A)^0)(v) , ... , ((d/dx + A)^(n-1))(v)`,``, ` are linearly independent. Then`,``, ` a_0*((d/dx + A)^0)(v) + ... + a_n((d/dx + A)^n)(v) = 0`,``,` for some\ a_0 .. a_n in k(x). By imposing the condition a_n=1 these a_0 .. a_n`, ` are uniquely determined by v. Then the adjoint of the differential operator` ,``,` f = a_n DF^n + ... a_0 DF^0`,``, ` is equivalent with the following matrix differential equation:`,``, ` (d/dx + A)(y) = 0 where y is a vector.`,``, ` This means for example that (d/dx + A)(y) = 0 has a rational solution`, ` y in k(x)^n if and only if adjoint(f) has a rational solution in k(x).`, ``,`- The output of this procedure is a list containing`, ` * The operator f.`,` * A list containing the vectors ((d/dx + A)^0)(v) .\ .. ((d/dx + A)^(n-1))(v)`, ` Note that the first element of this list is the cyclic vector v.`,``, `EXAMPLE:`, `A:=linalg[matrix](2,2,[(4*x^4-8*x^3+11*x^2-12*x-3)/(x^5+6*x^3-2*x^4-8*x^2+`,`\ 9*x-6),-3*x/(2*x^3+3*x^2+6*x+9),(-2*x^2-9*x-9)/(x^2+5),(2*x^5+8*x^4-17*x^3+` ,` 21*x^2-61*x-13)/(2*x^6-x^5+10*x^4-6*x^2+25*x-30)]);`, `c:=cyclic_vector(A);`,`f:=c[1];`,`rational_solutions(adjoint(f));`, ` 2`, ` (x + 5) (x - 1)`, ` [--------------------]`, ` 2`, ` (x + 3) (x - x + 2)`,``,``, `# So the space of rational solutions of A is 1 dimensional`, `SolveMat(A,rational);`, ` x - 1 (x - 1) (2 x + 3)`, ` [[- ----------, - -----------------]]`, ` 2 2`, ` x - x + 2 x - x + 2`): ##################################### # The mixed differential equation #---------------------------------------- ##################################### diffop_mixed_eqn # Procedures: # eigenring # endomorphism_charpoly macro( bounds_frlf0=`diffop/bounds_frlf0`, local_frlf0=`diffop/local_frlf0`, lift_frlf0=`diffop/lift_frlf0`, liftstep_frlf0=`diffop/liftstep_frlf0`, eigenring_lineq=`diffop/eigen_lineq` ): `help/text/eigenring`:=TEXT( `FUNCTION: eigenring - compute the endomorphisms of the solution space` ,``,`CALLING SEQUENCE:`,` eigenring(f,verify);`,``,`PARAMETERS:` ,` f - a differential operator`,` verify - an optional parameter`, ``,`SYNOPSIS: `, `- The input f is a differential operator. Denote V(f) as the solution space`, ` of f. This procedure computes a basis (as vector space) of the set of all`, ` operators r for which r(V(f)) is a subset of V(f).`, ` So r is an endomorphism of the solution space V(f). The characteristic`, ` polynomial of this map can be computed by the command`, ` endomorphism_charpoly(f,r);`, ` For such r we have that mult(f,r) is divisible on the right by f, so`, ` the remainder (which is rightdivision(mult(f,r),f)[2]) must be zero. If`, ` the option verify is given in the input, then this procedure will perform`, ` this check.`,``,`EXAMPLE: `, ` > f:=DF^4+6/x*DF^3+2*(x^2-1)/x^4*DF^2-2*(3*x^2-1)/x^5*DF+1/x^8:`, ` > v:=eigenring(f);`, ` 3 5 2 4 3`, ` v := [- DF x - DF x + (2 x + x) DF, 1]`,``, ` > for r in v do factor(endomorphism_charpoly(f,r)) od;`,``, ` 2 2`, ` (x - 6 x + 7)`,``, ` 4`, ` (x - 1)`): # f: monic differential operator # r: element of E_D(f,f), so r is a k linear map r: V(f) -> V(f) # RootOf syntax, parameters allowed. # Output: characteristic polynomial of the k-linear action of r on V(f) endomorphism_charpoly:=proc(f,r) local a,i,S,s,Si,j,n; global DF,x; if g_normal(subs(x=0,denom(f)*lcoeff(f,DF)))=0 then # x=0 is not a regular point RETURN(procname(seq(collect(subs(x=x+1,i),DF,g_normal),i=[f,r]))) fi; n:=degree(f,DF); # Solution of f: S:=convert([seq(a[i]*x^i,i=0..2*n-1)],`+`); # Apply f: s:=0; Si:=S; for i from 0 to n do s:=s+convert(taylor(coeff(f,DF,i),x=0,2*n),polynom)*Si; Si:=diff(Si,x) od; s:=collect(s,x); S:=subs(solve(evala({seq(coeff(s,x,i),i=0..n-1)}) ,{seq(a[i],i=n..2*n-1)}),S); # Apply r on the solution S: s:=0; Si:=S; for i from 0 to n do s:=s+convert(taylor(coeff(r,DF,i),x=0,n),polynom)*Si; Si:=diff(Si,x) od; s:=collect(s,x); S:={seq(a[i]=0,i=0..n-1)}; collect(linalg[charpoly](linalg[matrix](n,n,[seq(seq(subs(S, subs(a[j]=1,coeff(s,x,i))),i=0..n-1),j=0..n-1)]),x),x,g_normal) end: # Compute linear equations for coefficients of r, from the condition that # the r(V(f)) \subset V(f) mod x^(n+N) eigenring_lineq:=proc(f,r) local a,i,S,s,Si,n,N,var; global DF,x; if g_normal(subs(x=0,denom(f)))=0 then RETURN(`not a regular point`) fi; n:=degree(f,DF); N:=2; # We'll obtain N*n equations # Solution of f: S:=convert([seq(a[i]*x^i,i=0..2*n-1+N)],`+`); # Apply f: s:=0; Si:=S; for i from 0 to n do s:=s+convert(taylor(coeff(f,DF,i),x=0,2*n+N),polynom)*Si; Si:=diff(Si,x) od; s:=collect(s,x); S:=subs(solve(evala({seq(coeff(s,x,i),i=0..n-1+N)}) ,{seq(a[i],i=n..2*n-1+N)}),S); # Apply r on the solution S: s:=0; Si:=S; for i from 0 to n do s:=s+convert(taylor(coeff(r,DF,i),x=0,n+N),polynom)*Si; Si:=diff(Si,x) od; s:=collect(s,x); var:={seq(a[i]=0,i=0..n-1)}; s:=collect(s-convert([seq(subs(var,subs(a[i]=1,S))*coeff(s,x,i) ,i=0..n-1)],`+`),x); var:=indets(var); S:=[seq(coeffs(expand(coeff(s,x,i)),var),i=n..n-1+N)]; {seq(g_normal(numer(i)),i=S)} end: # Compute the ring of endomorphisms E_D(f), in other words the solutions r # of the mixed equation f*r+l*f=0. eigenring:=proc(F) local a,B,n,f0,L,xDFpow,i,values,zeros,r,ac,j,ne,f; global DF,x,xDF; f:=collect(F,DF,g_normal); f:=collect(f/lcoeff(f,DF),DF,g_normal); userinfo(2,diffop,`Computing endomorphism ring of`,f); if subs(x=0,denom(f))=0 then RETURN(map(collect,subs(x=x+1, procname(subs(x=x-1,f),args[2..nargs])),DF,g_normal)) fi; # Now we can assume that x=0 is regular B:=bounds_frlf0(f); userinfo(2,diffop,`bounds are:`,B); n:=degree(f,DF); # nc:=convert([seq(i[2]+1,i=B)],`+`); # Number of coeffs to be determined f0:=l_p(f,0); L:=local_frlf0(f0); # The following suffices for expressing the candidate global solutions # in n^2 variables, using the local solutions up to this accuracy. ac:=max(seq(B[i][2]+2-i,i=1..n)); # If we take ac we have the following number of equations ne:=convert([seq(ac-B[i][2]-2+i,i=1..n)],`+`); # We need at least n^2 equations: ac:=ac+max(0,ceil((n^2-ne)/n)); L:=lift_frlf0(f0,L,ac,`only r`); # Now transform to d/dx syntax: xDFpow:=[1]; for i from 1 to n-1 do xDFpow:=[op(xDFpow),mult(x*DF,xDFpow[nops(xDFpow)])] od; values:=NULL; zeros:=[seq(0,i=1..n)]; for i from 1 to n^2 do r:=collect(L[i],xDF); # convert to DF syntax as follows r:=expand(convert([seq(coeff(r,xDF,j-1) *xDFpow[j],j=1..n)],`+`)); # The accuracy of coeff(r,DF,j) is ac+j # Now divide by B[j+1][1] r:=convert([seq(convert(series( coeff(r,DF,j)/B[j+1][1],x=0,ac+j),polynom)*DF^j,j=0..n-1)],`+`); # The high part should be cancel, the low part contributes to # the value of the solution zeros:=[seq(zeros[j]+a[i]*quo(coeff(r,DF,j-1)*x^(-min(0,B[j][2]+1)),x^max(0,B[j][2]+1),x) ,j=1..n)]; values:=values,r od; values:=[values]; userinfo(3,diffop,`solving linear equations in`,n^2,`variables...`); j:={seq(a[i],i=1..n^2)}; r:=readlib(`solve/linear`)({seq(coeffs(collect(i,x),x),i=zeros)},j); userinfo(3,diffop,`done solving`); r:=collect(convert([seq(subs(r,a[i])*values[i],i=1..n^2)],`+`),DF); r:=convert([seq(g_normal(coeff(r,DF,i)*B[i+1][1])*DF^i,i=0..n-1)],`+`); i:=0; while zeros<>{0} do i:=i+1; zeros:=eigenring_lineq(op(subs(x=x-i,[f,r]))); if not member(zeros,{`not a regular point`,{0}}) then userinfo(3,diffop,`solving extra linear equations`); r:=subs(`solve/linear`(zeros,j),r) fi od; j:=indets(r) intersect j; if member(verify,[args]) then lprint(`Found `,nops(j),`linearly independent solutions. Checking correctness...`); zeros:=[coeffs(collect(rightdivision(numer(mult(f,r)),f)[2],DF),DF)]; zeros:={seq(coeffs(collect(numer(i),x),x),i=zeros)}; r:=subs(`solve/linear`(zeros,j),r); j:=indets(r) intersect j; lprint(`Verified correctness. Found `,nops(j),`independent solutions`) fi; [seq(collect(subs(i=1,solve(j),r),DF,g_normal),i=j)] end: # Bounds for the mixed equation f*r+l*f=0 bounds_frlf0:=proc(f) local v,ext,n,singularities,p,i,E,r,vp,s,m,B,D,N; global g_conversion1,g_conversion2,x,xDF; if nargs=1 then v:=g_ext(f); RETURN(subs(g_conversion2, procname(subs(g_conversion1,f),v))) fi; ext:=args[2]; n:=degree(f,DF); if n<2 then ERROR(`wrong arguments`) fi; singularities:=[seq(i[1],i=v_ext_m(g_factors(denom( g_normal(f/lcoeff(f,DF))),x,ext),x)),infinity]; userinfo(3,diffop,`Computing the bounds...`); for i from 0 to n-1 do D[i]:=1; N[i]:=0 od; for p in singularities do B[p,[]]:=0; for i from 0 to n-1 do B[p,i]:=infinity od; v:=factor_op(l_p(f,p),`semireg`); for E in v do if not type(E,list) then bug() fi; r:=degree(E[3],x); # ramification index vp:=ldegree(E[4],x)/r; # the number v' in lemma 5 s:=E[1]; # semi-regular part s:=[solve(nt(s,1),xDF)]; # list of exponents of s for i in s do if not type(i,integer) then bug() fi od; # Smallest difference in the canonical list that is an # integer divided by the ramification index: m:=(min(op(s))-max(op(s)))/r; for i from 0 to n-1 do # lemma 5, order(r) <= n-1 B[p,i]:=min(B[p,i],ceil(m+(n-1-i)*vp)) od; if p<>infinity and type(E[4],integer) and min(seq(E[4]+i,i=s))=0 then # Can improve the bound a bit for this case: i:=0; while member(i-E[4],s) do i:=i+1 od; B[p,[]]:=i fi od; for i from 0 to n-1 do if p<>infinity then if g_ext([f,p])<>ext then r:=subs(_Z=x,subs(g_conversion1,op(subs(g_conversion2,p)))) else r:=x-p fi; B[p,i]:=B[p,i]+i; # From the xDF -> DF conversion # If 0,1,..,B[p,[]]-1 are exponents and there are no negative # integer exponents then the bound can be tightened a bit, # using the following argument. # Suppose r=r_(n-1)*DF^(n-1)+..+r_0*DF^0 # Suppose that at the point x=p (where p is finite) we have lower bounds # B[i],i=0..n-1 for the valuation of r_i. Because of the xDF -> DF conversion # we always have B[0] B[0]. # However, the valuation of r(y) must be >=0 because 0 was the smallest integer # exponent and r(y) must be a (possibly 0) solution of f. If now # B[0]<0, then we can replace B[0] by B[0]+1. If then there is also an # exponent 1, apply same argument, and can increase B[1] (but then also B[0]), # if B[1]<0. # Use the property: if y in C((x))[log(x)] with valuation v>=0 and degl=0, # then the n'th derivative of y has valuation >=v-n+1 (for the +1 degl=0 is # needed). if B[p,[]]=n then B[p,i]:=0 elif i0 # then ERROR() # fi; userinfo(10,diffop,`Computing the`,x^a,`term of a local endomorphism`); if a {} then ERROR(`failed`) fi; [l+x^a*F[1],r+x^a*F[2],a+1] else gcdex(subs(xDF=xDF+a,f0),f0,-subs(g_conversion2,frlf),xDF ,'r_extra','l_extra'); subs(g_conversion1,[l+x^a*modm(l_extra),r+x^a*modm(r_extra),a+1]) fi end: #################################### # Global factorizations #------------------------------------------ #################################### diffop_global macro( compute_bound=`diffop/compute_bound`, factor_global=`diffop/factor_global`, factor_minmult1=`diffop/factor_minmult1`, factor_order1=`diffop/factor_order1`, flist=`diffop/flist`, same_charclass=`diffop/same_charclass`, try_factorization=`diffop/try_factorization`, try_factorization2=`diffop/try_factorization2`, xDFn_modr=`diffop/xDFn_modr`, try_eigenring=`diffop/try_eigenring` ): DFactor:=proc(ff) global modulus,modulus2; local i,b,f,ext,eb,opt; f:=convert(ff,RootOf); if f<>ff then RETURN(DFactor(f,args[2..nargs])) fi; if member(xDF,indets(ff)) then # Local operator as input: hence local factorization. RETURN(factor_op(ff,`split over k((x))`,g_ext([args]))) elif member(`clear all`,[args]) then `diffop/back_up`(`clear all`) fi; if degree(ff,DF)<=1 then RETURN([subs(g_conversion2,ff)]) fi; opt:=op({args[2..nargs]} minus {`clear all`}); modulus:=0; modulus2:=503; ext:=g_ext([args]); f:=collect(subs(g_conversion1,ff),DF,g_normal); eb:=NULL; for i in [args[2..nargs]] do if type(i,`=`) and op(i)[1]=bound then eb:=op(i)[2] fi od; if member(`alg factor`,[args]) then if eb=NULL then ERROR(`must specify a bound`) fi; RETURN(subs(g_conversion2,factor_global(f,ext,eb,`alg factor`))) fi; b:=factor_global(f,ext,eb); b:=subs(g_conversion2,b); if nops(b)=1 or member(`one step`,{args}) then b elif member(`right factor`,{args}) then DFactor(b[2],opt) elif member(LCLM,{args}) and degree_ext(b[2],ext)>1 then b:=subs(g_conversion1,op(DFactor(b[2],opt,`right factor`))); b:=subs(g_conversion2,[op(DFactor(rightdivision(f,LCLM(b,`and conjugates`,ext) ,ext)[1],opt)),_LCLM(b,`and conjugates`,ext)]); if b[1]=1 then b:=[b[2..nops(b)]] fi; b else [op(DFactor(b[1],opt)),op(DFactor(b[2],opt))] fi end: # Computes xDF^n mod r. xDFn_modr:=proc(n,r,ext) local a; options remember; if nfailed then userinfo(3,diffop,`possible factor found modulo prime`,modulus2 ,`of order`,i); # Computing with higher accuracy: exponent:=4; flm:=compute_modm(exponent,try_factorization2,fl,bound,eb,nstep); while flm<>failed do sr:=expand(convert([seq( flm[1][j+1]*xDF^j,j=0..nops(flm[1])-1)],`+`)); sr:=expand(reconstruct(sr/lcoeff(sr,indets([x,sr])) mod modulus2^exponent,modulus2^exponent)); if sr<>failed then if type(f,list) and f[2]=adjoint then # Note: adjoint and substitute(x*DF, ) do not commute sr:=adjoint(sr) fi; sr:=substitute(x*DF,sr,ext); if point=infinity then sr:=substitute(-x^2*DF,subs(x=1/x,sr)) else sr:=subs(x=x-point,sr) fi; if not(type(f,list)) then sr:=collect(sr/lcoeff(sr,DF),DF,g_normal); re:=rightdivision(f,sr,ext); if g_normal(re[2])=0 then userinfo(3,diffop,`found factor`,sr); RETURN([re[1],sr]) fi; userinfo(3,diffop,`unlikely failed attempt`); # to prevent always getting failed attempts: nstep:=nstep+2; # Try GCRD sr:=GCRD(f,sr); if has(sr,DF) then userinfo(3,diffop,`found: `,sr); RETURN([rightdivision(f,sr)[1],sr]) fi else # Because of f was made monic after localizing, # the adjoint of this local operator is not # the adjoint of f. To fix this sr must be # multiplied with the leading coefficient of # the localized f. if point=infinity then re:=x^(-degree(f[1],DF)) else re:=(x-point)^(-degree(f[1],DF)) fi; sr:=collect(re*lcoeff(f[1],DF)*sr,DF); sr:=mult(sr,1/lcoeff(sr,DF),ext); re:=leftdivision(f[1],sr,ext); if g_normal(re[2])=0 then userinfo(3,diffop,`found factor`,sr); RETURN([sr,collect(re[1],DF,g_normal)]) fi; userinfo(3,diffop,`unlikely failed attempt`); nstep:=nstep+2 fi fi; userinfo(4,diffop,`higher prime power needed`); exponent:=ceil(exponent*1.7); flm:=compute_modm(exponent,try_factorization2,fl,bound,eb ,nstep) od fi od; failed end: try_factorization2:=proc(fl,bound,eb,nstep) local i,v,sm,sm_old; v:=lift_pade2(op(fl),nstep); sm:=smallest_pade2(op(v)); while sm<>sm_old or sm[nops(sm)]=0 do sm_old:=sm; v:=lift_pade2(op(v),nstep); sm:=smallest_pade2(op(v)); for i in sm do if degree(i,x)>bound[nops(sm)-1]+eb then RETURN(failed) fi od od; [sm,v] end: # Input: a global operator f, the coefficients field and a bound for the # number of extra singularities. # Output: either [f] or f factored in 2 factors. factor_global:=proc(f,ext,extra_bound,what) local t,v,i,j,k,l,bound,eb,singularities,min_deg,done_s,all_one,R_left; if degree(f,DF)<=1 then RETURN([f]) fi; # point, localized f, algebraic extensions: singularities:={[infinity,convert(f,diffop,x=infinity),ext],seq([i[1], convert(f,diffop,x=i[1]),[op(i[3]),op(ext)]],i=v_ext_m(g_factors( denom(g_normal(f/lcoeff(f,DF))),x,ext),x))}; userinfo(5,diffop,`done computing localizations`); for i from 1 to degree(f,DF)-1 do bound[i]:=0 od; eb:=-1; for k in singularities do j:=Newtonpolygon(k[2]); v:=[seq(seq(j[i][2]+(l-j[i][1])*(j[i+1][2]-j[i][2]) /(j[i+1][1]-j[i][1]),l=j[i][1]..j[i+1][1]-1),i=1..nops(j)-1)]; for i from 1 to degree(f,DF)-1 do bound[i]:=bound[i]+degree_ext(k[3],ext)*v[i+1] od; eb:=eb+degree_ext(k[3],ext) od; # list of bounds for the coefficients of possible global factors bound:=[seq(-bound[degree(f,DF)-j]+j*eb,j=1..degree(f,DF)-1)]; if nargs<3 then eb:=compute_bound(singularities,ext); if eb>20 then userinfo(1,diffop, `Computed bound number of extra singularities is`, eb,`Use option bound=computed to use this bound.` ,`Now using a lower bound`); eb:=15 fi elif extra_bound=computed then eb:=compute_bound(singularities,ext) elif type(extra_bound,integer) then eb:=extra_bound else ERROR(`wrong arguments`) fi; userinfo(5,diffop,`bounds are`,bound,eb); all_one:=true; done_s:=[]; # The singularities we have considered while singularities<>{} do min_deg:=min(seq(degree_ext(i[3],ext),i=singularities)); for i in singularities do if degree_ext(i[3],ext)=min_deg then singularities:=singularities minus {i}; v:=factor_op(i[2],`all right factors`,i[3]); done_s:=[op(done_s),[op(i),v]]; for j in v do if skipped_factors(j)=0 then userinfo(3,diffop,cat(`Minimum multiplicity 1 in x=`,i[1])); RETURN(factor_minmult1(bound,i[1],f,i[2],eb,j,i[3])) fi od; all_one:=evalb(all_one and nops(v)=1 and degree(v[1],xDF)=1) fi od; od; userinfo(3,diffop,`Minimum multiplicity >1`); if all_one then # f has only 1 exponential part userinfo(3,diffop,`Looking for first order factors`); v:=factor_order1(f,done_s,ext); if v=[f] and degree(f,DF)>3 then # Play our last card: v:=try_eigenring(subs(g_conversion2,f)); if nops(v)=1 then userinfo(1,diffop,`factorization of`,f,`may be incomplete`) fi fi; RETURN(v) fi; if degree(f,DF)<=3 then bug() fi; v:=try_eigenring(subs(g_conversion2,f)); if nops(v)>1 then RETURN(v) fi; all_one:=true; for i in done_s do all_one:=evalb(all_one and nops(i[4])=1); for j in i[4] do t:=try_factorization(j,floor(degree(f,DF)/2),bound, i[1],f,eb,i[3]); if t<>failed then RETURN(t) fi od; v:=factor_op(adjoint(i[2],i[3]),`all right factors`,i[3]); for j in v do t:=try_factorization(j,degree(f,DF)-1,bound, i[1],[f,adjoint],eb,i[3]); if t<>failed then RETURN(t) fi od; # Now try the same, but searching for higher order factors: for j in i[4] do t:=try_factorization(j,degree(f,DF)-1,bound, i[1],f,eb,i[3],floor(degree(f,DF)/2)+1); if t<>failed then RETURN(t) fi od; od; if all_one then # up to conjugation only 1 exp part in every singularity userinfo(3,diffop,`Trying with algebraic extensions`); for i in done_s do if degree(i[4][1],xDF)>1 then j:=factor_op(i[4][1],`alg factor`,i[3])[1]; if degree(j[3],x)=1 then k:=j[1] elif degree_ext(j,i[3])=1 then next else k:=make_rightfactor(i[4][1],j,g_ext(j)); # the slope: l:=-valuation_laurent(coeff(i[4][1],xDF,0),0) /degree(i[4][1],xDF); k:=faster_riccati_split(i[4][1],rightdivision(i[4][1] ,k,g_ext(j),l),k,g_ext(j),l)[2] fi; t:=try_factorization(k,degree(f,DF)-1,bound,i[1],f,0,g_ext(j)); if t<>failed then RETURN(t) fi; v:=factor_op(adjoint(i[2],g_ext(j)),`all right factors`,ext); R_left:=0; for l in v while R_left=0 do if same_charclass(l ,adjoint(k,g_ext(j)),g_ext(j)) then R_left:=l fi od; if R_left=0 then bug() fi; # Now try if R gives a left hand global factor t:=try_factorization(R_left,degree(f,DF)-1,bound,i[1] ,[f,adjoint],0,g_ext(j)); if t<>failed then RETURN(t) fi fi od; fi; userinfo(1,diffop,`factorization of`,f,`may be incomplete`); [f] end: # Quick hack, this could be done somewhat more efficient. # RootOf syntax try_eigenring:=proc(f) local k,i,j,l,ext; options remember; ext:=g_ext(f); k:=eigenring(f); if nops(k)>1 then k:={seq(`if`(has(i,DF),i,NULL),i=k)}; # Remove the constant j:=min(seq(length(i),i=k)); for i in k do if length(i)=j then k:=i; break fi od; j:=g_factors(subs(g_conversion1,endomorphism_charpoly(f,k)) ,x,ext); l:=min(seq(degree(i[1],x),i=j)); for i in j do if degree(i[1],x)=l then l:=i[1]; break fi od; k:=GCRD(k-g_zero_of(l,x,'ext'),f); l:=rightdivision(f,k); if l[2]<>0 then bug() fi; [l[1],k] else [f] fi end: # R is a local factor with multiplicity 1. # The result of the procedure is a factorization. If f is returned unfactored # it is irreducible (if the specified bound is correct). factor_minmult1:=proc(bound,point,f,f_local,eb,R,ext) local t,w,fl,slope,k,l,R_left,i; if g_ext(f)<>g_ext([args]) then t:=try_eigenring(subs(g_conversion2,f)); if nops(t)>1 then RETURN(t) fi fi; fl:=floor(degree(f,DF)/2); t:=try_factorization(R,fl,bound,point,f,eb,ext); if t<>failed then RETURN(t) fi; w:=factor_op(adjoint(f_local,ext),`all right factors`,ext); R_left:=0; for i in w while R_left=0 do if same_charclass(i,adjoint(R,ext),ext) then R_left:=i fi od; if R_left=0 then bug() fi; # Now try if R gives a left hand global factor t:=try_factorization(R_left,fl,bound,point,[f,adjoint],eb,ext); if t<>failed then RETURN(t) fi; # Now try more terms: t:=try_factorization(R,degree(f,DF)-1,bound,point,f,eb,ext,fl); if t<>failed then RETURN(t) fi; t:=try_factorization(R_left,degree(f,DF)-1,bound,point,[f,adjoint],eb ,ext,fl); if t<>failed then RETURN(t) fi; userinfo(5,diffop,`trying if algebraic extensions are needed`); # if igcd(op(map(degree,w,xDF)))=1 then RETURN([f]) fi; # if igcd(seq(`if`(skipped_factors(i)=0,degree(i,xDF),0),i=w))=1 then # RETURN([f]) # fi; if igcd(degree(f,DF),degree(R,xDF))=1 then RETURN([f]) fi; slope:=-valuation_laurent(coeff(R,xDF,0),0)/degree(R,xDF); k:=op(factor_op(R,`alg factor`,ext)); if degree(k[3],x)=1 then l:=k[1] elif degree_ext(k,ext)=1 then RETURN([f]) else l:=make_rightfactor(R,k,g_ext(k)); l:=faster_riccati_split(R,rightdivision(R,l,g_ext(k),slope)[1] ,l,g_ext(k),slope)[2] fi; t:=try_factorization(l,fl,bound,point,f,0,g_ext(k)); if t<>failed then RETURN(t) fi; [f] end: # f and g irreducible. # Output: true if f and g have the same exponential parts. same_charclass:=proc(f,g,ext) local gg,ff,r,c,d; if degree(f,xDF)<>degree(g,xDF) then RETURN(false) fi; if degree(f,xDF)=1 then RETURN(type(nt(f-g,1),integer)) fi; r:=map(Newtonpolygon,[f,g],`include the Newton polynomials`); if r[1][1][3]=0 then c:=r[1][1][4]; d:=degree(c,x); c:=(coeff(c,x,d-1)-coeff(r[2][1][4],x,d-1))/d; if not type(c,integer) then RETURN(false) fi; r:=[[[r[1][1][1..3],g_expand(subs(x=x-c,r[1][1][4]),ext)] ,r[1][2..nops(r[1])]],r[2]] fi; if r[1]<>r[2] then RETURN(false) fi; gg:=factor_op(g,`alg factor`,ext)[1]; r:=nt(gg[1],1); ff:=substitute(2*xDF-r,ramification_of(f,gg[3],gg[2]),ldegree(r,x),0,gg[2]); ff:=Newtonpolygon(ff,`include the Newton polynomials`); if ff[1][3]<>0 then RETURN(false) fi; ff:=factors(collect(ff[1][4]/lcoeff(ff[1][4],x),x,g_normal)); for r in ff[2] do if degree(r[1],x)=1 and type(coeff(r[1],x,0)/coeff(r[1],x,1),integer) then RETURN(true) fi od; false end: # Search 1st left or right factors of an operator which has only 1 # exponential part in all singularities. factor_order1:=proc(ff,singularities,ext) local f,s,a,i,zero,t,var,bound,ansatz; s:=0; for i in singularities do if i[1]<>infinity then t:=subs(x=x-i[1],nt(factor_op(i[2],`split over k((x))` ,i[3])[1]-xDF,1)/x) else t:=-subs(x=1/x,x*nt(factor_op(i[2],`split over k((x))`,i[3])[1],0)) fi; if i[3]=ext then s:=s+t else s:=s+subs(g_conversion1,evala(Trace(op(subs(g_conversion2 ,[t,i[3],ext]))))) fi od; f:=substitute(DF-g_normal(s),ff,ext); bound:=nt(factor_op(convert(f,diffop,x=infinity),`split over k((x))` ,ext),1); bound:=max(seq(`if`(type(i-xDF,integer),i-xDF,-1),i=bound)); if bound<0 then RETURN([ff]) fi; ansatz:=convert([seq(a[i]*x^i,i=0..bound)],`+`); zero:=ansatz*coeff(f,DF,0); t:=ansatz; for i from 1 to degree(f,DF) do t:=diff(t,x); zero:=zero+coeff(f,DF,i)*t od; var:={seq(a[i],i=0..bound)}; t:=solve({coeffs(collect(numer(subs(g_conversion2,zero)),x),x)}); if subs(t,ansatz)=0 then if nargs=4 or degree(ff,DF)=2 then RETURN([ff]) else s:=factor_order1(adjoint(f),singularities,ext,ext); if nops(s)=1 then RETURN([ff]) fi; RETURN(map(adjoint,[s[2],s[1]])) fi fi; i:=bound; while nops(indets(subs(t,ansatz)) intersect var)>1 do t:=solve({op(t),a[i]=0},var); i:=i-1 od; ansatz:=subs(g_conversion1,subs(t,ansatz)); ansatz:=subs(op(indets(ansatz) intersect var)=1,ansatz); ansatz:=DF+g_normal(s-diff(ansatz,x)/ansatz); f:=rightdivision(ff,ansatz,ext); if g_normal(f[2])<>0 then bug() fi; [f[1],ansatz] end: # Bound for the number of extra singularities in a 1st order right factor compute_bound:=proc(singularities,ext) local v,i,res,j,Z,ma,c; global g_conversion2,x,xDF; userinfo(5,diffop,`computing number of extra singularities bound`); res:=0; for i in singularities do v:=factor_op(i[2],`semireg`,i[3]); ma:=-infinity; for j in v do c:=coeff(min(solve(nt(j[1],1)))+j[4],x,0)/degree(j[3],x); # Now compute the trace of this algebraic number divided by # the degree of this number over Q: c:=evala(Norm(Z-subs(g_conversion2,c))); c:=coeff(c,Z,degree(c,Z)-1)/lcoeff(c,Z)/degree(c,Z); if not type(c,rational) then bug() fi; ma:=max(c,ma) od; res:=res+degree_ext(i[3],ext)*ma od; res:=floor(res); userinfo(5,diffop,`bound for extra singularities is`,res); res end: `help/text/DFactor`:=TEXT( ``, `FUNCTION: DFactor - factorizes differential operators`, ``, `CALLING SEQUENCE: DFactor(f,opt)`, ``, `PARAMETERS:`, ` f - a differential operator`, ` opt - (optional) a sequence of options`, ``, `SYNOPSIS:`, ` This procedure is meant for factorizing differential operators. Global`, ` factorization is done over Q_bar(x)[DF] (Q_bar = algebraic closure of Q). The`, ` local factorization works over k((x))[xDF] where k is the coefficients field`, ` by default. A larger field k can be specified by giving a list of RootOfs, k`, ` must be Q or an algebraic extension of Q. Local operators must have the`, ` syntax that can be obtained by convert[diffop].`, ``, ` The output is a list of factors [f1..fr] such that f=mult(f1..fr)`, ``, `OPTIONS:`, ` For global factorization the following options are available:`, ` ``one step``: meant for irreducibility testing. The computation stops`, ` whenever a factorization is found.`, ` ``clear all``: clears all global variables and option remember tables`, ` used by the diffop package. Use this option to free some memory`, ` or to avoid errors when a computation was interupted. However, this`, ` option also wipes out the values of all local operators.`, ` For local factorization the following option is available:`, ` [sequence of RootOfs]: for factorization over a larger constants field.`, ``, ` Some information will be printed if infolevel[diffop] is given an integer`, ` value in the range 0 to 10. The default value of infolevel[diffop] is 1.`, ` With this value only warning messages about possibly incomplete`, ` factorizations are printed. The value 10 causes many messages, the value 0`, ` none.`, ``, `EXAMPLES:`, ` a:=RootOf(x^2-10);`, ` R:=(35/4-5/2*a)*x^2-55/2+11/2*a;`, ` f:=DF^2-R;`, ` DFactor(f); # factorization over Q_bar(x)`, ` 6 4 4 2 2`, ` [DF + (45 x - 300 x - 60 x %1 + 525 x + 150 x %1 - 110 - 34 %1)`, ``, ` / 5 3 3`, ` (- 1/6 + 1/30 %1) / (3 x - 10 x - 2 x %1 + 7 x + 2 x %1),`, ` /`, ``, ` 6 4 4 2 2`, ` DF + (45 x - 300 x - 60 x %1 + 525 x + 150 x %1 - 110 - 34 %1)`, ``, ` / 5 3 3`, ` (- 1/30 %1 + 1/6) / (3 x - 10 x - 2 x %1 + 7 x + 2 x %1) ]`, ` /`, ``, ` 2`, `%1 := RootOf(_Z - 10)`, ``, ` mult(op("));`, ``, ` 2 2 2 2`, ` DF + (15 x - 110 - 22 RootOf(_Z - 10)) (- 7/12 + 1/6 RootOf(_Z - 10))`, ``, ``, ` f:=randpoly([x,DF],degree=3);`, ` f:=mult(f,randpoly([x,DF],degree=3));`, ` f:=convert(f,diffop,x=infinity); # localize f at x=infinity`, ` DFactor(f); # Factorization over Q((x))`, ` convert(",diffop,``truncated``,5); # give the factorization modulo x^5`, ``, `See: ?LCLM and ?symmetric_power for more factorization examples`, ``, `SEE ALSO: diffop` ): # A basis for all solutions for which y'/y is a rational function. # Output: exp_solutions:=proc(f) local v,w,i,j; v:=DFactor(f,LCLM); w:=NULL; for i in v do if has (i,_LCLM) then w:=w,op(1,i) else w:=w,i fi od; v:=NULL; for i in {w} do if degree(i,DF)=1 then j:=-coeff(i,DF,0)/coeff(i,DF,1); v:=v,[j,rational_solutions( substitute(DF+j,f))] fi od; v:=[v]; if member(sorted,[args]) then v else [seq(seq(exp(Int(i[1],x))*j,j=i[2]),i=v)] fi end: `help/text/exp_solutions`:=TEXT( `FUNCTION: exp_solutions - Compute the exponential solutions`,` `, `CALLING SEQUENCE:`,` exp_solutions(f)`,` `,`PARAMETERS:`, ` f - a differential operator`,` `,`SYNOPSIS: `, `- This procedure computes a basis for all (up to conjugation) solutions`, ` which are of the form exp(int( a rational function ,x)).`,` `, `EXAMPLE: `, ` f:=DF^2-2*(8*x^5+11*x^4+8*x^3+17*x^2+14*x-1)/(x^3+3)^2/(x^2+x+1)^2;`, ` exp_solutions(f);`): ######################### # global variables #----------------------------------------------------- ######################### diffop_global_vars # `diffop/procedures`:={mult,eval_laurent,rightdivision,adjoint}: # for x in {anames()} do if type(x,procedure) and # has(remember,[op(3,op(x))]) and searchtext(diffop,x)=1 then # `diffop/procedures`:={op(`diffop/procedures`),x} # fi od: # x:=evaln(x): # Back up all the global variables `diffop/back_up`:=proc(to_do) global g_conversion1,g_conversion2,accuracy_laurent, description_laurent,set_laurents,value_laurent,rem_lift, modp_table,modulus,modulus2,`diffop/backup`; local i,P,L,j; modp_table:=evaln(modp_table); modulus:=0; modulus2:=503; if to_do=`clear all` then userinfo(3,diffop,`clearing all global variables in diffop`); g_conversion1:={}; g_conversion2:={}; accuracy_laurent:=evaln(accuracy_laurent); description_laurent:=evaln(description_laurent); set_laurents:={}; value_laurent:=evaln(value_laurent); rem_lift:=evaln(rem_lift); for i in `diffop/procedures` do assign(i=subsop(4=NULL,op(i))) od elif to_do=`back up` then if not `diffop/backup`=evaln(`diffop/backup`) then RETURN(procname(`clear all`)) fi; userinfo(3,diffop,`backing up all global variables in diffop`); `diffop/backup`:=[g_conversion1,g_conversion2, op(op(accuracy_laurent)),op(op(description_laurent)),set_laurents, op(op(value_laurent)),op(op(rem_lift)), {seq([i,[op(4,op(i))]],i=`diffop/procedures`)}]; procname(`clear all`) elif to_do=`restore` then # Necessary if the global variables were messed up by a CTRL-C userinfo(3,diffop,`restoring all global variables in diffop`); g_conversion1:=`diffop/backup`[1]; g_conversion2:=`diffop/backup`[2]; accuracy_laurent:=table(`diffop/backup`[3]); description_laurent:=table(`diffop/backup`[4]); set_laurents:=`diffop/backup`[5]; value_laurent:=table(`diffop/backup`[6]); rem_lift:=table(`diffop/backup`[7]); for i in `diffop/backup`[8] do if i[2]<>[] then P:=i[1]; L:=op(op(i[2])); for j in L do P(op(1,j)):=op(2,j) od fi od; `diffop/backup`:=evaln(`diffop/backup`) fi; NULL end: lprint(`For help type: ?diffop`);