{VERSION 4 0 "IBM INTEL NT" "4.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 1 12 0 0 255 1 0 2 0 
0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 1 14 128 0 0 1 0 0 2 0 
0 0 0 0 0 1 }{CSTYLE "2D Output" 2 20 "" 1 14 0 0 128 1 0 0 0 0 0 0 0 
0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 14 0 0 0 1 
2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 3 0 3 0 2 2 0 1 }{PSTYLE "Maple Outpu
t" -1 12 1 {CSTYLE "" -1 -1 "Times" 1 16 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }
1 3 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "R3 Font 0" -1 256 1 {CSTYLE "" 
-1 -1 "Courier" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 
2 2 0 1 }{PSTYLE "R3 Font 2" -1 257 1 {CSTYLE "" -1 -1 "Times" 1 12 
128 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }}
{SECT 0 {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7973 "##    -*-Maple-*-
\n##    Title:\011gfun package (for Generating FUNctions)\n##    Creat
ed:\011Wed Mar\011 4 15:13:42 1992\n##    Authors:\011Bruno Salvy <Bru
no.Salvy@inria.fr>\n##\011\011Paul Zimmermann <Paul.Zimmermann@inria.f
r>\n##\011     and, in 1996,\n##\011\011Eithne Murray <Eithne.Murray@i
nria.fr>\n##\n## Description: converts implicit equations into differe
ntial equations,\n##\011\011differential equations into recurrences an
d vice-versa,\n##\011\011ordinary into exponential recurrences.\n##
\011\011also converts lists to linear recurrences or differential\n##
\011\011equations.\n##\n##  Some of the ideas used in this file are du
e to S. Plouffe and F. Bergeron.\n##\n## Modifications:\n##\n##   2.64
 ported to maple 6\011\011\011\011\011B.S. Apr 00\n##   2.63 fixed a p
roblem with singularities in holexprtodiffeq\011B.S. Jul 99\n##   2.62
 rectoproc generates a smaller proc in the log case\011B.S. Jun 99\n##
   2.61 removed last change due to another bug in V.5\011\011B.S. Sep \+
98\n##   2.60 added option 'normal' to convert/fullparfrac to get\011B
.S. Sep 98\n##        the expected behaviour of Maple V.4\n##   2.59 i
mproved speed of systomatrix\011\011\011\011B.S. Aug 98\n##   2.58 hol
exprtodiffeq did not work on exp(algebraic)^rat\011B.S. Aug 98\n##   2
.57 fixed a problem where algfuntoalgeq could return\011B.S. Jul 98\n#
#\011  rational function coefficients instead of polynomials.\n##   2.
56 case of non-homogeneous recurrence with constant\011B.S. Jul 98\n##
\011  coefficients in the homogeneous part now handled by\n##        r
ectoproc by the logarithic algorithm.\n##   2.55 workaround a bug in s
olve/linear in the case of\011B.S. Jul 98\n##        floating point co
efficients\n##   2.54 fixed a bug for multiple branches in algeqtoseri
es\011B.S. May 98\n##   2.53 fixed a bug in rectodiffeq in the inhomog
enous case\011B.S. May 98\n##   2.52 undeclared local in diffeqtable[h
ypergeom]\011\011B.S. Jan 98\n##   2.51 fixed a bug in diffeqtorec in \+
the inhomogeneous case\011B.S. Jan 98\n##\011  when a simplification w
as possible\n##   2.50 help pages in the new format\011\011\011\011B.S
. May 97\n##\011  added diffeqtohomdiffeq and rectohomrec\n##\011  fix
ed a problem with the order in algfuntoalgeq\n##   2.49 changed the ha
ndling of initial conditions of \011B.S. May 97\n##\011  differential \+
equations so that singular solutions\n##\011  can be considered more e
asily.\n##        added many functions to holexprtodiffeq\n##   2.48 a
dded hypergeom to holexprtodiffeq and changed the\011B.S. Feb 97\n##
\011  way it deals with functions having two arguments\n##   2.47 fixe
d a problem in algeqtoseries due to the \011\011B.S. Feb 97\n##       \+
 of solve with RootOf's\n##   2.46 faster guessing when coefficients a
re integers\011B.S. Sep 96\n##   2.45 speedup by using collect less, a
nd coeffs more\011B.S. Sep 96\n##   2.44 got rid of many D's, hence an
 important speedup\011B.S. Sep 96\n##   2.43 fixed a bug in listtorec,
 added BesselK and BesselY\011B.S. Aug 96\n##\011  to holexprtodiffeq
\n##   2.42 more careful handling of the leading coefficient in\011B.S
. Jul 96\n##\011  the recurrence given by the user\n##   2.41 added ha
ndling of pol(x)^(free(x)) to holexprtodiffeq B.S. Jul 96\n##   2.40 f
ixed a bug in guessgf for algebraic equations\011B.S. Jul 96\n##   2.3
9 fixed holexprtodiffeq on expressions not involving\011E.M. Mar 96\n#
#\011  the variable; improved error message/help page\n##   2.38 fixed
 another weakness of rectoproc\011\011\011B.S. Mar 96\n##   2.37 fixed
 weaknesses of the new options of rectoproc\011B.S. Mar 96\n##   2.36 \+
improved speed of pade2 for polynomial coefficients\011B.S. Mar 96\n##
   2.35 added options list and params to rectoproc\011\011B.S. Mar 96
\n##   2.34 improved order of equation returned by algfuntoalgeq\011E.
M. Mar 96\n##   2.33 listtorec now uses pade2 and should be much faste
r\011B.S. Mar 96\n##   2.32 improved rejection of solutions from pade2
\011\011B.S. Feb 96\n##\011  improved efficiency of pade2 by a huge fa
ctor\n##\011  improved help pages\n##\011  a few bug fixes\n##   2.31 \+
added holexprtodiffeq and algfuntoalgeq\011\011E.M. Feb 96\n##\011  E.
 Murray - Eithne.Murray@inria.fr\n##   2.3  more type checking\011\011
\011\011\011B.S. Jan 96\n##\011  extended poltodiffeq and poltorec to \+
diff and shift\n##\011  algebraicsubs now handles initial conditions\n
##\011  added algeqtoseries\n##   2.29 Fixed a bug reported by E. Murr
ay in algebraicsubs\011B.S. Jan 96\n##   2.28 same with poltorec\011
\011\011\011\011B.S. Jan 96\n##   2.27 poltodiffeq added,\011\011\011
\011\011B.S. Jan 96\n##\011  `diffeq[+*]diffeq` rewritten using poltod
iffeq\n##   2.26 Further improvements on the initial conditions\011B.S
. Jan 96\n##   2.25 Better treatment of free initial conditions\011
\011B.S. Dec 95\n##   2.24 Small improvement to gfun/isolve\011\011
\011B.S. Dec 95\n##   2.23 Bug reported by H. Prodinger. gfun/gausseli
m was not\011B.S. Nov 95\n##\011  returning the smallest possible orde
r.\n##   2.22 suppressed the use of 'easy' in the guessing part,\011B.
S. Nov 95\n##\011  this missed too many cases\n##   2.21 fixed a bug i
n rec*rec in the nonhomogeneous case\011B.S. Aug 95\n##   2.2  macro M
APLE5.4\011\011\011\011\011B.S. Apr 95\n##\011  New version for the sh
are library.\n##   2.19 algeqtodiffeq did not allow functions that wer
e\011B.S. Apr 95\n##\011   singular at the origin.\n##\011  further po
lishing of goodinitvalues/diffeq both in\n##\011   regular and singula
r case.\n##\011  same modification as version 2.13 for diffeq+diffeq\n
##   2.18 formatrec did not complain on empty recurrences\011B.S. Apr \+
95\n##\011   plus a bug fix in systomatrix for equations with\n##\011 \+
  parameters.\n##   2.17 rec+rec used diffeq+diffeq, but the order of
\011\011B.S. Mar 95\n##\011   the recurrence was not always optimal. R
ewritten\n##\011   to compute directly on recurrences.\n##   2.16 diff
eq+diffeq had a local variable leak when the\011B.S. Feb 95\n##\011   \+
initial conditions involved a linear dependency\n##\011   (like D(y)(0
)=y(0))\n##   2.15 rectoproc improved in the non-remember case\011\011
B.S. Jan 95\n##\011    following an idea of PZ.\n##\011  removed the m
acros MAPLE5.2 and MAPLE5.3\n##\011    (switching did not work anymore
)\n##\011  added a case in diffeqtorec where the order of the\n##\011 \+
   recurrence can be lowered by adding an inhomogeneous part\n##\011  \+
added Laplace as an alias to invborel\n##\011  fixed a weakness in rec
todiffeq\n##   2.14 diffeq+diffeq was returning a non-formatted equati
on\011B.S. Jan 95\n##\011    when called with two rational functions.
\n##\011  formatdiffeq reinforced to forbid y(0) in the coefficients.
\n##\011  rectodiffeq now returns a homogeneous equation when\n##\011 \+
   the inhomogeneous one contains (D@@k)(y)(0).\n##   2.13 diffeq*diff
eq did not always find all the initial\011B.S. Jan 95\n##\011  conditi
ons. This is fixed by giving an extra arg\n##\011  to goodinitvalues/d
iffeq. Same for goodinitvalues/rec.\n##   2.12 reduced the number of R
ETURN statements\011\011B.S. Dec 94\n##   2.11 fixed bug in algebraics
ubs\011\011\011\011B.S. Dec 94\n##   2.1 fixed bugs reported by Ch. Ma
llinger\011\011\011P.Z. Nov 94\n##   2.07 heavy testing, plus a refine
ment of the criterion for B.S. Aug 94\n##\011rejecting overdetermined \+
results from pade2\n##   2.06 faster reversion of power series\011\011
\011B.S. Jul 94\n##   2.05 guessing part rewritten to use pade2 for ra
tional\011B.S. Jun 94\n##\011series, differential and algebraic equati
ons\n##\011plus a bug fixed in the initial conditions of a product\n##
   2.01 a few extra checks, clarified help for listtolist\011B.S. May \+
94\n##   2.0 new version\011\011\011\011\011\011B.S. May-Nov 93\n##   \+
rectoproc, diffeqtorec, rectodiffeq, rectoproc,\n##   borel, invborel \+
rewritten.\011\011\011\011\011B.S. Dec 92\n##   formatdiffeq, formatre
c added.\011\011\011\011B.S. Dec 92\n##\n\n###########################
###############################################\n## CONVENTION: a recu
rrence is an expression of the form\n##\011sum(p[i]*u(n+i),i=0..d)\011
\011\011  (E)\n## or\n##    \{ sum(p[i]*u(n+i),i=0..d), u(0)=a0,..., u
(k)=ak \}.\n##\n## The p[i]'s are polynomials in n. The sequence(s) re
presented by such\n## a recurrence are solutions of (E) for n>=k, wher
e k is the largest\n## positive integer solution to p[d](n-d)=0, or 0 \+
if this does not cancel.\n## For n<k, the values of the sequence are g
iven by the initial conditions,\n## and it is part of the convention t
hat u(i)=0 for i<0.\n#################################################
#########################\n\ngfun := module()\n\n\011option package;\n
\n\011export \n\011\011algfuntoalgeq,\n\011\011algebraicsubs,\n\011
\011algeqtodiffeq,\n\011\011algeqtoseries,\n\011\011borel,\n" }}
{PAGEBK }{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 3259 "\011\011cauchyproduct,\n\011\011`diffeq+diffeq`,\n
\011\011`diffeq*diffeq`,\n\011\011diffeqtable,\n\011\011diffeqtohomdif
feq,\n\011\011diffeqtorec,\n\011\011guesseqn,\n\011\011guessgf,\n\011
\011hadamardproduct,\n\011\011holexprtodiffeq,\n\011\011invborel,\n
\011\011Laplace,\n\011\011listtoalgeq,\n\011\011listtodiffeq,\n\011
\011listtohypergeom,\n\011\011listtolist,\n\011\011listtoratpoly,\n
\011\011listtorec,\n\011\011listtoseries,\n\011\011listtoseriestable,
\n\011\011maxdegcoeff,\n\011\011maxdegeqn,\n\011\011maxordereqn,\n\011
\011mindegcoeff,\n\011\011mindegeqn,\n\011\011minordereqn,\n\011\011op
tionsgf,\n\011\011pade2,\n\011\011poltodiffeq,\n\011\011poltorec,\n
\011\011ratpolytocoeff,\n\011\011`rec+rec`,\n\011\011`rec*rec`,\n\011
\011rectodiffeq,\n\011\011rectohomrec,\n\011\011rectoproc,\n\011\011se
riestoalgeq,\n\011\011seriestodiffeq,\n\011\011seriestohypergeom,\n
\011\011seriestolist,\n\011\011seriestoratpoly,\n\011\011seriestorec,
\n\011\011seriestoseries,\n\011\011version;\n\n\011global\n\011\011`ty
pe/gfun/free`,\n\011\011`type/gfun/has2diffeqs`,\n\011\011`type/gfun/i
dentity`,\n\011\011`type/gfun/initeq`;\n\n\011local \n\011\011`algeqto
series/doit`,\n\011\011`algeqtoseries/prettyprint`,\n\011\011algfuntoa
lgeq2,\n\011\011`algfuntoalgeq/formpoly`,\n\011\011`algfuntoalgeq/add`
,\n\011\011`algfuntoalgeq/mul`,\n\011\011`algfuntoalgeq/RootOf`,\n\011
\011borelinvborel,\n\011\011cheapgausselim,\n\011\011`diffeqtorec/doit
`,\n\011\011firstnonzero,\n\011\011formatdiffeq,\n\011\011formatpoleq,
\n\011\011formatrec,\n\011\011funtodiffeq,\n\011\011gensolvelin,\n\011
\011getname,\n\011\011`goodinitvalues/diffeq`,\n\011\011`goodinitvalue
s/rec`,\n\011\011guessandcheck,\n\011\011indicialeq,\n\011\011inicond
\011,\n\011\011isholonomic,\n\011\011lindep,\n\011\011listprimpart,\n
\011\011`listtoseries/egf`,\n\011\011`listtoseries/Laplace`,\n\011\011
`listtoseries/ogf`,\n\011\011`listtoseries/revogf`,\n\011\011`listtose
ries/revegf`,\n\011\011`listtoseries/lgdogf`,\n\011\011`listtoseries/l
gdegf`,\n\011\011`l2r/l2r`,\n\011\011`l2h/l2h`,\n\011\011makediffeq,\n
\011\011makerec\011,\n\011\011maxindex,\n\011\011minindex,\n\011\011my
gcdex\011,\n\011\011myisolve,\n\011\011nbinicond,\n\011\011powcompose,
\n\011\011powcomposesimple,\n\011\011powcomposesimpledoit,\n\011\011po
wdivide,\n\011\011powrevert,\n\011\011powtruncate,\n\011\011pprimeknow
ingp,\n\011\011rechomrecbis,\n\011\011`rectodiffeq/doit`,\n\011\011`pa
de2/exmin`,\n\011\011`pade2/doit`,\n\011\011`ratpolytocoeff/elmt`,\n
\011\011ratsolvelin,\n\011\011rectohomrecbis,\n\011\011remove\011,\n
\011\011`s2d/s2d`,\n\011\011`s2a/s2a`,\n\011\011systomatrix,\n\011\011
typecheck\n\011\011;\n\n\n######################### Global Constants #
#############################\n\nversion:=2.64:\noptionsgf:=['ogf','eg
f']:\nmaxordereqn:=3: # default 3rd order\nminordereqn:=1: # default 1
st order\nmaxdegcoeff:=4: # default degree 4 coefficients\nmindegcoeff
:=0: # default constant coefficients\nmaxdegeqn:=3:\011# default maxim
um 3rd degree\nmindegeqn:=2:\011# default minimum 2nd degree\n\n######
################### The pade2 package ##########################\n# Th
is should not be there. It just simplifies the distribution for this\n
# version.\n## -*-Maple-*-\n##\n##    Title:\011pade2\n##    Created:
\011Sep 1993\n##    Author:\011Harm Derksen <hderksen@sci.kun.nl>\n##
\n## Description: pade-hermite approximants.\n##\n##    Modified: Oct \+
1993\n##    Author:\011Bruno Salvy <Bruno.Salvy@inria.fr>\n##    Modif
ication: rewritten for efficiency\n##\n##    June 94. Added the option
 'easy'. BS.\n##    Feb. 96. Improved efficiency.\011BS.\n\nmacro(BIG=
1000); # this is only used with option easy\n\n# This does the interfa
ce\npade2:=proc(functionlist::list(algebraic),\n\011    point::\{name,
name=algebraic\},\n\011    accuracy::\{integer,list(nonnegint)\},\n
\011    optional::identical(easy))\nlocal x, a, n, l, i, acc, m, easy,
 result;\n    if type(point,`=`) then a:=op(2,point); x:=op(1,point)\n
    else a:=0; x:=point fi;\n    n:=nops(functionlist);\n    if type(a
ccuracy,list) then\n\011m:=max(op(accuracy));\n\011l:=[seq(m-i,i=accur
acy)];\n\011acc:=convert(accuracy,`+`)+n-1\n    else\n" }}{PAGEBK }
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8524 "
\011l:=[0$n];\n\011acc:=accuracy\n    fi;\n    easy:=evalb(nargs=4);\n
    result:=`pade2/doit`(map(taylor,subs(x=x+a,functionlist),x,acc),\n
\011x,l,acc-1,easy); # it has to be taylor and not series\n    if resu
lt=FAIL then FAIL else subs(x=x-a,result) fi\nend: # pade2\n\n\n#  pad
e2/exmin\n#  Extended minimum.\n# Input: a list, a boolean function, a
nd an optional name.\n# Output: the minimum of the list according to t
he order\n# the name being assigned the index of the first occurrence \+
of\n# the minimum in the list.\n#\n# No attempt at efficiency has been
 made, since this should really be in\n# the kernel, with sort.\n`pade
2/exmin`:=proc(l,order,aux)\nlocal res;\n    res:=op(1,sort(l,order));
\n    if nargs=3 then member(res,l,aux) fi;\n    res\nend: # `pade2/ex
min`\n\n# `pade2/doit`\n# Input: a list of series in the variable x, t
he variable x, a list of degree\n#\011bounds and a bound on the number
 of iterations. This latter bound\n#\011should be at most the order of
 the series.\n# Output: a list of polynomial coefficients for the seri
es, such that\n#\011the scalar product of this list with the input has
 zero Taylor\n#\011coefficients up to a large order.\n`pade2/doit`:=pr
oc(ll,x,degs,nbiter,easy)\nlocal y, i, L, n, j, l, ind, l2, pivot, var
s, leadcoeff, k, normalize, lk,\n  rationalcoeff, tt, finished, dmin, \+
deg, lterm, jj, lterm2, inds, polycoeff,\n  ll2, dmin2, ind2, co, cco,
 cco1, cco2, cco3;\n    n:=nops(ll);\n    vars:=[seq(y[i],i=1..n),x];
\n    inds:=indets(ll) minus \{x\};\n    if inds=\{\} then\n\011normal
ize:=x->x;\n\011rationalcoeff:=type(\{op(map(op,ll))\} minus \{O(1)\},
set(rational));\n\011polycoeff:=false\011# polycoeff means non-rationa
l but polynomial\n    else\n\011normalize:=normal;\n\011rationalcoeff:
=false;\n\011polycoeff:=type(\{op(map(op,ll))\} minus \{O(1)\},\n\011 \+
   set(polynom(rational,inds)))\n    fi;\n    for i to n do\n\011L[i]:
=series(y[i]*x^degs[i],x,infinity);\n\011l[i]:=ll[i]; co[i]:=1\n    od
;\n    finished:=false;\n    for j from 0 to nbiter while not finished
 do\n\011userinfo(3,'pade2',`iteration number`,j);\n\011dmin:=infinity
;\n\011# Select the pivot\n\011for i to n do\n\011    leadcoeff[i]:=no
rmalize(coeff(l[i],x,j));\n\011    if leadcoeff[i]<>0 then\n\011\011l2
[i]:=i;\n\011\011deg:=op(nops(L[i]),L[i]);\n\011\011if deg>dmin then n
ext\n\011\011elif deg<dmin then\n\011\011    dmin:=deg; lterm:=op(nops
(L[i])-1,L[i]);\n\011\011    dmin2:=length(l[i]); ind2:=i\n\011\011els
e\n\011\011    lterm2:=op(nops(L[i])-1,L[i]);\n\011\011    if length(l
[i])<dmin2 then ind2:=i; dmin2:=length(l[i]) fi;\n\011\011    for jj t
o n while\n\011\011\011has(lterm,y[jj]) and has(lterm2,y[jj]) do od;\n
\011\011    if jj<=n and not has(lterm,y[jj]) or\n\011\011\011jj>n and
 leadcoeff[i]<>O(1) and\n\011\011\011length(leadcoeff[i])>=length(lead
coeff[ind]) then next\n\011\011    else lterm:=lterm2 fi\n\011\011fi\n
\011    else\n\011\011l2[i]:=NULL;\n\011\011if coeff(l[i],x,j)=0 then \+
next\n\011\011else # normalizing did it\n\011\011    l[i]:=map(normali
ze,l[i]);\n\011\011    if l[i]<>0 and op(1,l[i])<>O(1) then next fi\n
\011\011fi\n\011    fi;\n\011    ind:=i;\n\011    pivot:=leadcoeff[ind
];\n\011    if pivot=O(1) then break fi\n\011od;\n\011if pivot=O(1) th
en k:=ind; break fi;\n\011if dmin=infinity then next fi;\n\011ll2:=[se
q(l2[i],i=1..n)];\n\011if rationalcoeff and pivot<>1 and pivot<>-1 the
n\n\011    tt:=abs(icontent(convert([seq(leadcoeff[i]*y[i],i=ll2)],`+`
)));\n\011    if tt<>1 then\n\011\011pivot:=pivot/tt;\n\011\011for i i
n ll2 do leadcoeff[i]:=leadcoeff[i]/tt od\n\011    fi\n\011elif polyco
eff then\n\011    tt:=content(convert([seq(leadcoeff[i]*y[i],i=ll2)],`
+`),vars);\n\011    if tt<>1 then\n            divide(pivot,tt,'pivot'
);\n            for i in ll2 do divide(leadcoeff[i],tt,evaln(leadcoeff
[i])) od\n\011    fi\n\011fi;\n\011tt:=l[ind2];\n\011for i in ll2 do\n
\011    if i<>ind then\n\011\011l[i]:=# map(normalize,\n\011\011    se
ries(pivot*l[i]-leadcoeff[i]*l[ind],x,infinity);#);\n\011\011if l[i]=0
 or op(1,l[i])=O(1) then\n\011\011    L[i]:=series(pivot*co[i]*L[i]\n
\011\011\011-leadcoeff[i]*co[ind]*L[ind],x,infinity);\n\011\011    fin
ished:=true\n\011\011elif op(2,l[i])=j then l[i]:=subsop(1=0,l[i])\n
\011\011fi\n\011    fi\n\011od;\n\011if finished then # select the sma
llest one\n\011    dmin:=infinity;\n\011    for i in ll2 do\n\011\011i
f (l[i]=0 or op(1,l[i])=O(1)) and length(L[i])<dmin then\n\011\011    \+
dmin:=length(L[i]); k:=i\n\011\011fi\n\011    od;\n\011    break\n\011
fi;\n\011l[ind]:=series(x*tt,x,infinity);\n\011tt:=L[ind2]; cco:=co[in
d2];\n\011for i in ll2 do if i<>ind then\n\011    if rationalcoeff the
n\n            cco2:=pivot*co[i];\n            cco3:=co[ind]*numer(lea
dcoeff[i]);\n            cco1:=igcd(cco2,cco3);\n            cco2:=cco
2/cco1; cco3:=co[ind]*leadcoeff[i]/cco1;\n            L[i]:=series(cco
2*L[i]-cco3*L[ind],x,infinity);\n            cco2:=numer(icontent(add(
icontent(op(2*jj-1,L[i]))*y[jj],\n                jj=1..iquo(nops(L[i]
),2))));\n            if cco2<>1 and cco2<>-1 then\n                L[
i]:=series(L[i]/cco2,x,infinity) fi;\n            co[i]:=cco1*cco2\n
\011    elif polycoeff then\n            cco1:=gcd(pivot*co[i],leadcoe
ff[i]*co[ind],'cco2','cco3');\n            L[i]:=map(normal,series(cco
2*L[i]-cco3*L[ind],x,infinity));\n            cco2:=content(convert(L[
i],polynom),vars,'cco3');\n            if cco2<>1 and cco2<>-1 then\n \+
               L[i]:=series(cco3,x,infinity); co[i]:=cco1*cco2\n      \+
      else co[i]:=cco1\n            fi\n\011    else L[i]:=map(normali
ze,\n                       series(pivot*L[i]-leadcoeff[i]*L[ind],x,in
finity))\n\011    fi fi\n\011od;\n\011L[ind]:=series(x*tt,x,infinity);
 co[ind]:=cco;\n\011if easy and min(op(map(length,[seq(L[i],i=1..n)]))
)>BIG*n then\n\011    RETURN(FAIL) fi;\n\011userinfo(5,'pade2',`curren
t list`,lprint([seq(L[i],i=1..n)]));\n    od;\n    if j=nbiter+1 and n
ot finished then\n        `pade2/exmin`([seq(op(2,l[i]),i=1..n)],numer
ic,'k') fi;\n    lk:=collect(convert(L[k],polynom),vars,'distributed',
normalize);\n    if rationalcoeff then divide(lk,icontent(lk),'lk')\n \+
   elif polycoeff then divide(lk,content(lk,vars),'lk') fi;\n    map(e
xpand,[seq(coeff(lk,y[i])*x^(-degs[i]),i=1..n)])\nend: # `pade2/doit`
\n\n######################### Type Checking ##########################
####\n\n`type/gfun/identity`:=proc(x) type(x,`=`) and op(1,x)=op(2,x) \+
end:\n`type/gfun/free`:=proc(x,y) not has(x,y) end:\n\n# type(y(0),ini
teq(y))\011\011-> true\n# type(D(y)(0),initeq(y))\011-> true\n# type((
D@@k)(y)(0),initeq(y))\011-> true\n# otherwise\011\011\011-> false\n`t
ype/gfun/initeq` := proc(expr,y)\nlocal f;\n    if not type(expr,funct
ion(0)) then\n\011false\n    else\n\011f := op(0,expr);\n\011f=y or f=
'D(y)' or (type(f,function(identical(y)))) and\n\011    type(op(0,f),`
@@`(identical(D),integer))\n    fi\nend: # `type/gfun/initeq`\n\n# Thi
s procedure avoids several type checks of the same expression.\n# Besi
des, it handles the defaults.\ntypecheck:=proc (n)\nlocal i;\n    if n
=1 then\011\011# l, x, <met>\n\011if nargs>2 and type([args[2..3]],[li
st,name]) then\n\011    if nargs=3 then RETURN('stamped',args[2..3],op
(1,optionsgf))\n\011    elif nargs=4 and type(listtoseriestable[args[4
]],procedure)\n\011\011then RETURN('stamped',args[2..4])\n\011    elif
 nargs>4 then ERROR(`too many arguments`)\n\011    else ERROR(`invalid
 argument`,args[4])\n\011    fi\n\011elif nargs<3 then ERROR(`too few \+
parameters`)\n\011elif not type(args[2],list) then ERROR(`not a list`,
args[2])\n\011elif not type(args[3],name) then ERROR(`not a name`,args
[3])\n\011else ERROR(`invalid arguments`)\n\011fi\n    elif n=2 then
\011# l, y(x), <[met]>\n\011if nargs>2 and type([args[2..3]],[list,fun
ction(name)]) then\n\011    if nargs=3 then RETURN('stamped',args[2..3
],optionsgf)\n\011    elif nargs=4 and type(args[4],list) and\n\011
\011type([seq(listtoseriestable[i],i=args[4])],\n\011\011    list(proc
edure)) then RETURN('stamped',args[2..4])\n\011    elif nargs>4 then E
RROR(`too many arguments`)\n\011    else ERROR(`invalid argument`,args
[4])\n\011    fi\n\011elif nargs<3 then ERROR(`too few arguments`)\n
\011elif not type(args[2],list) then ERROR(`not a list`,args[2])\n\011
elif not type(args[3],function(name)) then\n\011    ERROR(`invalid unk
nown function`,args[3])\n\011else ERROR(`invalid arguments`)\n\011fi\n
    elif n=3 then\011# l, x, [<met>]\n\011if nargs>2 and type([args[2.
.3]],[list,name]) then\n\011    if nargs=3 then RETURN('stamped',args[
2..3],optionsgf)\n\011    elif nargs=4 and type(args[4],list) and\n
\011\011type([seq(listtoseriestable[i],i=args[4])],\n\011\011    list(
procedure)) then RETURN('stamped',args[2..4])\n\011    elif nargs>4 th
en ERROR(`too many arguments`)\n\011    elif nargs=4 then ERROR(`inval
id argument`,args[4])\n\011    fi\n\011elif nargs<3 then ERROR(`too fe
w parameters`)\n\011elif not type(args[2],list) then ERROR(`not a list
`,args[2])\n\011elif not type(args[3],name) then ERROR(`not a name`,ar
gs[3])\n\011else ERROR(`invalid argments`)\n\011fi\n    elif n=4 then
\011# s, x, y\n\011if nargs=4 and type([args[2..4]],[series,name,name]
) then\n\011    RETURN('stamped',args[2..4])\n\011elif nargs<>4 then E
RROR(`wrong number of arguments`)\n\011elif not type(args[2],series) t
hen ERROR(`not a series`,args[2])\n\011elif not type(args[3],name) the
n ERROR(`not a name`,args[3])\n\011elif not type(args[4],name) then ER
ROR(`not a name`,args[4])\n" }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" 
}}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14413 "\011else ERROR(`invalid argume
nts`)\n\011fi\n    elif n=5 then\011# l, y(x)\n\011if nargs=3 and type
([args[2..3]],[list,function(name)]) then\n\011    RETURN('stamped',ar
gs[2..3])\n\011elif nargs<>3 then ERROR(`wrong number of arguments`)\n
\011elif not type(args[2],list) then ERROR(`not a list`,args[2])\n\011
elif not type(args[3],function(name)) then\n\011    ERROR(`invalid unk
nown function`,args[3])\n\011else ERROR(`invalid arguments`)\n\011fi\n
    elif n=6 then\011# s, y(x), <[met]>\n\011if nargs>2 and type([args
[2..3]],[series,function(name)]) then\n\011    if nargs=3 then RETURN(
'stamped',args[2..3],optionsgf)\n\011    elif nargs=4 and type(args[4]
,list) and\n\011\011type([seq(listtoseriestable[i],i=args[4])],\n\011
\011    list(procedure)) then RETURN('stamped',args[2..4])\n\011    el
if nargs>4 then ERROR(`too many arguments`)\n\011    else ERROR(`inval
id argument`,args[4])\n\011    fi\n\011elif nargs<3 then ERROR(`too fe
w arguments`)\n\011elif not type(args[2],series) then ERROR(`not a ser
ies`,args[2])\n\011elif not type(args[3],function(name)) then\n\011   \+
 ERROR(`invalid unknown function`,args[3])\n\011else ERROR(`invalid ar
guments`)\n\011fi\n    elif n=7 then\011# s, <[met]>\n\011if nargs>1 a
nd type(args[2],series) then\n\011    if nargs=2 then RETURN('stamped'
,args[2],optionsgf)\n\011    elif nargs=3 and type(args[3],list) and\n
\011\011type([seq(listtoseriestable[i],i=args[3])],\n\011\011    list(
procedure)) then RETURN('stamped',args[2..3])\n\011    elif nargs>3 th
en ERROR(`too many arguments`)\n\011    else ERROR(`invalid argument`,
args[3])\n\011    fi\n\011elif nargs=1 then ERROR(`too few arguments`)
\n\011elif not type(args[2],series) then ERROR(`not a series`,args[2])
\n\011else ERROR(`invalid arguments`)\n\011fi\n    elif n=8 then\011# \+
s, <met>\n\011if nargs>1 and type(args[2],series) then\n\011    if nar
gs=2 then RETURN('stamped',args[2],'ogf')\n\011    elif nargs=3 and ty
pe(listtoseriestable[args[3]],procedure)\n\011\011then RETURN('stamped
',args[2..3])\n\011    elif nargs>3 then ERROR(`too many arguments`)\n
\011    else ERROR(`invalid argument`,args[3])\n\011    fi\n\011elif n
args=1 then ERROR(`too few arguments`)\n\011elif not type(args[2],seri
es) then ERROR(`not a series`,args[2])\n\011else ERROR(`invalid argume
nts`)\n\011fi\n    elif n=9 then\011# l, <met>\n\011if nargs>1 and typ
e(args[2],list) then\n\011    if nargs=2 then RETURN('stamped',args[2]
,'ogf')\n\011    elif nargs=3 and type(listtoseriestable[args[3]],proc
edure)\n\011\011then RETURN('stamped',args[2..3])\n\011    elif nargs>
3 then ERROR(`too many arguments`)\n\011    else ERROR(`invalid argume
nt`,args[3])\n\011    fi\n\011elif nargs=1 then ERROR(`too few argumen
ts`)\n\011elif not type(args[2],list) then ERROR(`not a list`,args[2])
\n\011else ERROR(`invalid arguments`)\n\011fi\n    elif n=10 then\011#
 recurrence, function(name), function(name)\n\011if nargs>4 then ERROR
(`too many arguments`)\n\011elif nargs<4 then ERROR(`too few arguments
`)\n\011elif not type([args[3..4]],[function(name),function(name)]) th
en\n\011    ERROR(`invalid arguments`)\n\011fi\n    else ERROR(`should
 not happen`)\n    fi;\nend: # typecheck\n\ngetname:=proc(yofz::functi
on(name), y, z)\n    y:=op(0,yofz);\n    if type(y,procedure) then ERR
OR(`not an unassigned name`,y) fi;\n    z:=op(yofz)\nend:\n\n#########
######### Modifications of existing Maple code ############\n\n## Adap
ted from solve/linear/integer\nratsolvelin := proc (equations, unknown
s)\nlocal eqn, eqns, vars, sol, sols, i, k, l, pivot, var, c, j;\n    \+
# Suppress trivial equations and\n    # normalize: convert from Q to Z
 and divide out by the content\n    eqns := \{seq(l/icontent(l)*sign(l
), l=  equations minus \{0\})\};\n    # Solve the remaining system usi
ng a sparse implementation\n    for k while eqns <> \{\} do\n\011l:=ma
p(length,[op(eqns)]);\n\011member(min(op(l)),l,'i');\n\011eqn := eqns[
i];\n\011vars := indets(eqn) intersect unknowns;\n\011var := vars[1];
\n\011pivot := coeff(eqn,var,1);\n\011for j from 2 to nops(vars) do\n
\011    # Eliminate the unknown with the smallest coefficient\n\011   \+
 c := coeff(eqn,vars[j],1);\n\011    if length(c) < length(pivot) then
 pivot := c; var := vars[j] fi\n\011od;\n\011eqns:=subs(var=-1/pivot*s
ubs(var=0,eqn),subsop(i=NULL,eqns)) minus \{0\};\n\011sol[k]:=var,eqn,
pivot\n    od;\n    sols := \{\};j:=0;\n    for i from k-1 by -1 to 1 \+
do\n\011eqn[i] := - subs( sol[i][1] = 0, sols, sol[i][2] );\n\011if eq
n[i]<>0 then j:=j+1;l:=subsop(j=i,l) fi;\n\011sols := sols union \{sol
[i][1] = eqn[i] / sol[i][3]\}\n    od;\n    # Parameterize an infinite
 solution space\n    vars := \{op(unknowns)\} minus \{seq(op(1,i),i=so
ls)\};\n    if vars=\{\} then sols\n    elif nops(vars)=1 then \{op(su
bs(op(vars)=1,sols)),op(vars)=1\}\n    elif select(type,\{seq(op(2,i),
i=sols)\},rational) minus \{0\} <> \{\} then\n\011\{op(subs(seq(i=0,i=
vars),sols)),seq(i=0,i=vars)\}\n    else FAIL # we do not want to trea
t this\n    fi\nend: # ratsolvelin\n\n# A simple interface to solve/li
near\ngensolvelin:=proc ()\nlocal s, t;\n    s:=`solve/linear`(args);
\n    t:=select(type,s,`gfun/identity`);\n    if nops(t)<>1 then FAIL
\n    else subs(op(1,op(t))=1,(1=1)=(op(1,op(t))=1),s)\n    fi\nend: #
 gensolvelin\n\n\n# this procedure does not do exactly the same as gcd
ex, but is much faster\n# it returns g' and assign s,t such that s*a+t
*b=g'\n# and g'=lambda gcd(a,b) with lambda a rational fraction indepe
ndent of y\n##########################################################
#########\n##########  ONE CAN PROBABLY IMPROVE IT MORE   ############
########\n############################################################
#######\nmygcdex := proc(a,b,y,s,t)\nlocal q,r,u,v,du,dv,alpha,beta,ol
dalpha,oldbeta,tt;\n    Order:=infinity; # We are working on polynomia
ls\n    u:=series(a,y); v:=series(b,y);\n    if u=0 then b\n    elif v
=0 then a\n    else\n\011oldalpha:=1; oldbeta:=0;\n\011alpha:=0; beta:
=1;\n\011# u = oldalpha*a + oldbeta*b\n\011# v = alpha*a + beta*b\n
\011du:=op(nops(u),u);\n\011dv:=op(nops(v),v);\n\011do\n\011    if du<
dv then\n\011\011tt:=dv; dv:=du; du:=tt;\n\011\011tt:=v; v:=u; u:=tt;
\n\011\011if nargs>3 then tt:=alpha; alpha:=oldalpha; oldalpha:=tt fi;
\n\011\011if nargs>4 then tt:=beta; beta:=oldbeta; oldbeta:=tt fi;\n
\011    fi;\n\011    userinfo(3,'gfun',`degrees of the polynomials`,du
,dv);\n\011    userinfo(5,'gfun','polynomials',u,v);\n\011    if du=0 \+
and dv=0 then break fi;\n\011    q:=op(nops(u)-1,u)/op(nops(v)-1,v)*y^
(du-dv);\n\011    r:=normal(series(u-q*v,y));\n\011    if r=0 then bre
ak fi;\n\011    if nargs>3 then tt:=oldalpha-q*alpha; oldalpha:=alpha;
 alpha:=tt fi;\n\011    if nargs>4 then tt:=oldbeta-q*beta; oldbeta:=b
eta; beta:=tt fi;\n\011    u:=v;du:=dv;\n\011    v:=r;dv:=op(nops(v),v
)\n\011od;\n    ## These rem's are much too expensive on large example
s\n\011if dv=0 then\n#\011    if nargs>3 then s:=rem(collect(alpha/coe
ff(v,y,0),y,normal),b,y) fi;\n#\011    if nargs>4 then t:=rem(collect(
beta/coeff(v,y,0),y,normal),a,y) fi;\n\011    if nargs>3 then\n\011
\011s:=convert(normal(series(alpha/coeff(v,y,0),y)),polynom) fi;\n\011
    if nargs>4 then\n\011\011t:=convert(normal(series(beta/coeff(v,y,0
),y)),polynom) fi;\n\011    1\n\011else\n#\011    if nargs>3 then s:=r
em(alpha,b,y) fi;\n#\011    if nargs>4 then t:=rem(beta,a,y) fi;\n\011
    if nargs>3 then s:=convert(normal(series(alpha,y)),polynom) fi;\n
\011    if nargs>4 then t:=convert(normal(series(beta,y)),polynom) fi;
\n\011    v\n\011fi\n    fi\nend: # mygcdex\n\n# This is really needed
 because MapleV's isolve is terrible\n# eqn is a polynomial in n\nmyis
olve:=proc (eqn, n)\nlocal sol, i, d, f, g;\n    if not has(eqn,n) the
n \{\}\n    # These two lines are due to M. Monagan:\n    elif type(eq
n,polynom(rational,n)) then\n\011select(type,\{seq(i[1],i=roots(eqn))
\},integer)\n    elif type(eqn,`*`) then\n\011`union`(op(map(myisolve,
[op(eqn)],n)))\n    else\n\011f:=collect(eqn,n);\n\011g:=primpart(f,n)
;\n\011if expand(g-f)<>0 then RETURN(myisolve(g,n)) fi;\n\011d:=degree
(f,n);\n\011if d=1 then\n\011    f:=-coeff(f,n,0)/coeff(f,n,1);\n\011 \+
   if type(f,integer) then \{f\}\n\011    elif not has(f,RootOf) and n
ot has(f,`&RootOf`) then \{\}\n\011    else f:=evala(Normal(f));\n\011
\011if type(f,integer) then \{f\} else \{\} fi\n\011    fi\n\011elif d
=2 then\n\011    f:=[seq(coeff(f,n,i),i=0..2)];\n\011    d:=sqrt(op(2,
f)^2-4*op(3,f)*op(1,f));\n\011    select(type,map(normal,\n\011\011   \+
 \{(d-op(2,f))/2/op(3,f),-(d+op(2,f))/2/op(3,f)\}),integer)\n\011elif \+
type(eqn,`*`) then `union`(op(map(myisolve,[op(eqn)],n)))\n\011else\n
\011    sol:=traperror(isolve(expand(eqn),n));\n\011    if sol=lasterr
or then \{\}\n\011    else \{seq(op(2,i),i=select(type,map(op,[sol]),
\n\011\011identical(n)=integer))\}\n\011    fi\n\011fi\n    fi\nend: #
 myisolve\n\n######################### Various Utilities #############
######\n\nfirstnonzero:=proc(pol,n) max(-1,op(myisolve(pol,n)))+1 end:
\n\n# returns the smallest i such that u(n+i) appears in a recurrence
\nminindex := proc(rec,u,n)\n    min(op(map(op,indets(rec,specfunc(lin
ear(n),u)))))-n\nend: # minindex\n\n# returns the largest i such that \+
u(n+i) appears in a recurrence\nmaxindex := proc(rec,u,n)\n    max(op(
map(op,indets(rec,specfunc(linear(n),u)))))-n\nend: # maxindex\n\n# th
e opposite of select\n# this should be in the kernel\nremove:=proc(f,a
) map(subs(_f=f,\n    proc(x) if not _f(args) then x else NULL fi end)
,\n    args[2..nargs]) end:\n\n# systomatrix\n#  Input: a system of ho
mogeneous linear equations and a list of variables V,\n#\011  and the \+
name of a vector B.\n#  Output: a matrix A such that the system is equ
ivalent to A.V=B.\n# Almost like linalg[genmatrix], but more suitable \+
for our purpose.\n# Also, if sys is a list instead of a set, then the \+
order will be preserved.\nsystomatrix:=proc (sys, V, B)\nlocal a, i, j
, eqn, zero, c, t, ind, lco;\n    if sys=\{\} or sys=[] then\n\011B:=a
rray(1..1,[0]);\n\011array(1..1,1..nops(V),[[0$nops(V)]])\n    else\n
\011zero:=[seq(i=0,i=V)];\n\011a:=array(1..nops(sys),1..nops(V),sparse
);\n\011B:=array(1..nops(sys),sparse);\n\011for i to nops(V) do ind[V[
i]]:=i od;\n\011for i to nops(sys) do\n\011    eqn:=op(i,sys);\n\011  \+
  if type(eqn,`=`) then eqn:=op(1,eqn)-op(2,eqn) fi;\n##\011    eqn:=c
ollect(eqn,V,'distributed');\n#\011    eqn:=expand(eqn);\n\011    lco:
=[coeffs(eqn,V,'t')];\n\011    t:=[t];\n\011    for j to nops(t) do\n
\011\011if t[j]<>1 then a[i,ind[t[j]]]:=lco[j] else B[i]:=-lco[j] fi\n
\011    od;\n# \011    for j to nops(V) do\n# \011\011c:=coeff(eqn,V[j
]);\n# \011\011if c<>0 then a[i,j]:=c fi\n# \011    od;\n# \011    B[i
]:=-subs(zero,eqn) od;\n\011od;\n\011op(a)\n    fi\nend: # systomatrix
\n\n#listprimpart\n# Input: a list of polynomials and a variable\n# Ou
tput: the list where common factors have been removed.\nlistprimpart:=
proc (l, z)\nlocal i, T, q;\n    if hastype(l,float) then l\n    else
\n\011content(convert([seq(op(i,l)*T^(i-1),i=1..nops(l))],`+`),T,'q');
\n\011[seq(coeff(q,T,i),i=0..nops(l)-1)]\n    fi\nend: # listprimpart
\n\n# This is used to find a linear dependency.\n# The result is the l
ast element of the last line which is not 0.\n# The pivoting must not \+
be done as gausselim does it.\n# This is a fraction free gausselim usi
ng the input matrix for the intermediate\n# computations. Also, no typ
e checking is performed.\ncheapgausselim:=proc (A, nlin, ncol)\nlocal \+
lcols, c, k, i, u, j, piv;\n    lcols:=[$1..ncol];\n    for c to nlin \+
do\n\011for k in lcols while normal(A[c,k])=0 do A[c,k]:=0 od;\n\011if
 k=ncol then RETURN(A[c,ncol]) fi;\n\011piv:=A[c,k];\n\011for i from c
+1 to nlin do\n\011    if A[i,k]=0 then next else u:=A[i,k] fi;\n#\011
    for j in lcols do A[i,j]:=piv*A[i,j]-u*A[c,j] od\n\011    for j in
 lcols do A[i,j]:=normal(piv*A[i,j]-u*A[c,j]) od\n\011od;\n\011userinf
o(5,'gfun',`line `,c,` eliminated`);\n\011# this line adds some speedu
p for very large matrices:\n\011for i in lcols do A[c,i]:=0 od;\n\011u
serinfo(6,'rsolve',`remaining matrix `,print(op(A)));\n\011lcols:=subs
(k=NULL,lcols)\n    od;\n    FAIL\nend: # cheapgausselim\n\n#lindep\n#
 Input:  a vector u = [u1 , ... , uk] and a matrix A = array(1..k,1..l
)\n#\011  such that u[i] = sum(A[i,j] e[j]) for some e[j].\n#\011  The
 coefficients are rational functions in x.\n# Output: a linear depende
ncy relation between the u[i] if there is one,\n#\011  FAIL otherwise
\nlindep := proc(A, u, x)\nlocal k, i, l, v, rel, unk, j;\n    k := no
ps(u);\n    if linalg[rowdim](A)<>k then ERROR(`incorrect number of ro
ws`,op(A)) fi;\n    l := linalg[coldim](A);\n    userinfo(2,'gfun',`lo
oking for a linear dependency in a`,k,' x ',l,'matrix');\n    unk:=[se
q(v[i],i=1..k)];\n    rel := cheapgausselim(\n\011linalg[augment](A,li
nalg[matrix](k,1,unk)),k,l+1);\n    if rel=FAIL then FAIL\n    else\n
\011rel:=primpart(numer(convert([seq(normal(subs([i=1,seq(j=0,\n\011  \+
  j=subs(i=NULL,unk))],rel))*i,i=unk)],`+`)),unk);\n\011convert([seq(s
ubs([i=1,seq(j=0,j=subs(i=NULL,unk))],rel)\n\011    *u[op(i)],i=unk)],
`+`)\n    fi\nend: # lindep\n\n#formatdiffeq\n# Input:  a list [diffeq
,y(x)] containing a linear differential equation\n#\011     with polyn
omial coefficients and its unknown\n#\011  (optional) Y and X\n#\011  \+
    names to be assigned the unknown function and its variable\n#\011 \+
 (optional) iniconds to be assigned the initial conditions\n# Output: \+
a list of polynomials in x [u(x),p_0(x),...,p_d(x)] meaning\n#\011\011
\011\011\011       (d)\n#\011     eqn = p_0(x) y(x) + ... + p_d(x) y
\011  (x) + u(x)\n#\n#   This is where the type checking is done.\n#\n
formatdiffeq:=proc (l,Y,X,iniconds)\nlocal r, y, x, i, difford, j, lva
r, Z;\n    if nops(l)<>2 then ERROR(`wrong number of arguments`,op(l))
 fi;\n    getname(op(2,l),y,x);\n    if nargs>1 then X:=x; Y:=y fi;\n \+
   if type(op(1,l),set)  then\n\011r:=select(has,op(1,l),x);\n\011if n
ops(r)>1 then ERROR(`invalid differential equation`,op(l)) fi;\n\011if
 nops(r)=0 then\n\011    ERROR(`the unknown variable does not appear i
n the equation`) fi;\n\011if nargs>3 then iniconds:=op(1,l) minus r fi
;\n\011r:=op(r)\n    else r:=op(1,l); if nargs>3 then iniconds:=\{\} f
i fi;\n    if nargs>3 and has(iniconds,D) then D(Y):=D(Y) fi;\n    if \+
type(r,`=`) then r:=op(1,r)-op(2,r) fi;\n    if has(r,D) then r:=conve
rt(r,diff) fi;\n    if indets(r,specfunc(anything,y))<>\{y(x)\} then\n
\011ERROR(`invalid differential equation`,op(l)) fi;\n    r:=expand(nu
mer(normal(r)));\n    lvar:=select(has,indets(r,specfunc(anything,diff
)),y);\n    for difford from 0 while lvar<>\{\} do\n\011lvar:=subs(dif
f(y(x),x)=y(x),lvar) minus \{y(x)\} od;\n    if not type(r,linear([y(x
),seq(diff(y(x),x$i),i=1..difford)])) then\n\011ERROR(r,`is not a line
ar differential equation in`,y(x)) fi;\n    r:=subs([y(x)=Z,seq(diff(y
(x),x$j)=Z^(j+1),j=1..difford)],r);\n    r:=[subs(Z=0,r),seq(coeff(r,Z
,j+1),j=0..difford)];\n    if not type(r,list(polynom(anything,x))) th
en\n\011ERROR(`non-polynomial coefficients`,op(l)) fi;\n    listprimpa
rt(r,x)\nend: # formatdiffeq\n\n# makediffeq\n# Input:  a differential
 equation in the format returned by formatdiffeq\n#\011  the variables
 y and x, meaning the unknown function is y(x)\n#\011  (optional) a se
t of initial conditions\n# Output: the corresponding differential equa
tion, which is in a set if\n#\011  there are initial conditions.\nmake
diffeq:=proc (deq, y, x, ini)\nlocal r, i;\n    r:=convert([seq(deq[i+
2]*diff(y(x),[x$i]),i=0..nops(deq)-2),deq[1]],`+`);\n    if nargs=3 or
 ini=\{\} then r\n    else \{r,op(ini)\} fi\nend: # makediffeq\n\n#for
matrec\n# Input:  a list [rec,u(n)] containing a linear recurrence wit
h\n#\011     polynomial coefficients and its unknown\n" }}{PAGEBK }
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17362 "
#\011  u, n names to be assigned the unknown function and variable\n#
\011  (optional) iniconds to be assigned the initial conditions\n#\011
  (optional) one: boolean saying whether one is interested in\n#\011
\011\011  only one solution.\n# Ouput:  a list of polynomials in n: [b
(n),p_0(n),...,p_d(n)] meaning\n#\n#\011\011   rec=p_0(n)u(n)+...+p_d(
n)u(n+d)+b(n)\n#\n#  This is where the type checking is done.\n#\nform
atrec:=proc (l, u, n, iniconds, one)\nlocal r, i, U, N, mi, ma, locini
, j;\n    if nops(l)<2 or nops(l)>5 then\n\011ERROR(`wrong number of a
rguments`,op(l)) fi;\n    getname(op(2,l),U,N);\n    if nargs>1 then n
:=N; u:=U fi;\n    if type(op(1,l),set) then\n\011r:=select(has,op(1,l
),N);\n\011if nops(r)>1 then ERROR(`invalid recurrence`,op(l)) fi;\n
\011if r=\{\} then ERROR(`empty recurrence`) fi;\n\011if nargs>3 then \+
locini:=op(1,l) minus r;\n\011    if not type(locini,set(U(integer)=an
ything)) then\n\011\011ERROR(`invalid initial conditions`,locini)\n
\011    fi fi;\n\011r:=op(r)\n    else r:=op(1,l); if nargs>3 then loc
ini:=\{\} fi\n    fi;\n    if nargs>4 then if nops(l)>4 and type(op(5,
l),boolean) then one:=op(5,l)\n    else one:=false fi fi;\n    if type
(r,`=`) then r:=op(1,r)-op(2,r) fi;\n    mi:=minindex(r,U,N); ma:=maxi
ndex(r,U,N);\n    r:=collect(r,[seq(U(N+i),i=mi..ma)],normal);\n    if
 not type(r,linear([seq(U(N+i),i=mi..ma)])) then\n\011ERROR(`Not a lin
ear recurrence`,op(l)) fi;\n    if mi<>0 then\n\011if nargs>3 then\n
\011    j:=-1;\n\011    for i from 0 to mi do\n\011\011if not has(loci
ni,U(i)) then\n\011\011    j:=j+1; locini:=locini union \{U(i)=_C[j]\}
\n\011\011fi\n\011    od\n\011fi;\n    # The following two lines are n
ot valid for mi>0.\n    # Example: (n+1)*u(n+1)+(n+2)*u(n+2) and recto
diffeq.\n\011if mi<0 then\n\011    r:=subs(N=N-mi,r);\n\011    ma:=ma-
mi\n\011fi\n    fi;\n    if nargs>3 then iniconds:=locini fi;\n    if \+
has(map(denom,\{op(r)\}),N) then\n\011r:=collect(numer(normal(r)),[seq
(U(N+i),i=0..ma)],normal) fi;\n    r:=[subs([seq(U(N+i)=0,i=0..ma)],r)
,seq(coeff(r,U(N+i),1),i=0..ma)];\n    if has(r,U) or not type(r,list(
polynom(anything,N))) then\n\011ERROR(`invalid recurrence or unknown`,
op(l)) fi;\n    # listprimpart should not be called here because the l
eading coefficient\n    # might vanish at some point, which has an imp
act on the initial conditions\n    #  listprimpart(r,N)\n    r\nend: #
 formatrec\n\n# makerec\n# Input:  a recurrence in the format returned
 by formatrec\n#\011  the variables u and n, meaning the unknown seque
nce is u(n)\n#\011  (optional) a set of initial conditions\n# Output: \+
the corresponding recurrence, which is in a set if there are\n#\011  i
nitial conditions.\nmakerec:=proc (rec, u, n, ini)\nlocal r, i;\n    r
:=convert([seq(rec[i+2]*u(n+i),i=0..nops(rec)-2),rec[1]],`+`);\n    if
 nargs=3 or ini=\{\} then r\n    else \{r,op(ini)\} fi\nend: # makerec
\n\n# formatpoleq\n# Input:  a list [p, y(z)] containing a polynomial \+
in two variables and\n#\011   an unknown function and possibly initial
 conditions.\n#\011  y, z names to be assigned the unknown function an
d variable\n#\011  (optional) iniconds to be assigned the initial cond
itions.\n# Output: the polynomial, type checked and without its inital
 values\nformatpoleq:=proc (l, y, z, iniconds)\nlocal Y, Z, P;\n    if
 nops(l)<2 or nops(l)>4 then\n\011ERROR(`formatpoleq: wrong number of \+
arguments`,op(l)) fi;\n    getname(op(2,l),Y,Z);\n    if nargs>1 then \+
y:=Y; z:=Z fi;\n    if type(op(1,l),`=`) then P:=op(1,op(1,l))-op(2,op
(1,l))\n    else P:=op(1,l) fi;\n    if not type(P,polynom(anything,[Y
,Z])) then ERROR(`invalid argument`,P) fi;\n    if nargs=4 then\n\011i
f nops(l)>2 then\n\011    if type(op(3,l),set) then iniconds:=op(3,l)
\n\011    else ERROR(`invalid argument`,op(3,l)) fi;\n\011else inicond
s:=\{\} fi fi;\n    collect(primpart(P,Y),Y)\nend: # formatpoleq\n\n# \+
`goodinitvalues/rec`\n# Input:  a recurrence in the format returned by
 formatrec\n#\011  the unknown sequence and its variable\n#\011  some \+
initial conditions\n#\011  a boolean flag cleanup\n#\011  (optional) a
n integer p\n# Output: a set of equalities u(k)=v_k, from which all th
e other values\n#\011  can be deduced by solving the recurrence for it
s maximal index.\n#\011  This set is continued up to the pth term if p
 is given.\n#\011  New variables _C[i] are introduced when a value is \+
arbitrary.\n#\011  When the flag is true, the unnecessary entries u(k)
=_C[i] are removed\n#\011  The result is an ERROR when no initial cond
ition can be found.\n`goodinitvalues/rec`:=proc (rec, u, n, ini, flag,
 p)\nlocal n0, order, i, inds, minind, maxind, sys, r, sol, b, a,minin
d2,lmin, j, k, rej, inds2, dorej;\n    order:=nops(rec)-2;\n    maxind
:=order-1;\n    if type(ini,set) then inds:=map(op,indets(ini,u(intege
r)))\n    else inds:=\{\}; maxind:=max(maxind,ini) fi;\n    if nargs=6
 then maxind:=max(maxind,p) fi;\n    # The initial conditions are u(mi
nind)..u(maxind)\n    n0:=firstnonzero(subs(n=n-order,rec[nops(rec)]),
n);\n    # u(n0-1) cannot be deduced from the previous ones\n    maxin
d:=max(maxind,op(inds),n0-1);\n    minind:=min(op(inds),0);\n    # min
ind2 is used when the recurrence is not valid for small n\n    minind2
:=max(order-1,n0-1,op(inds))-order+1;\n    if minind2>minind then lmin
:=[minind,minind2] else lmin:=[minind] fi;\n    r:=makerec(rec,u,n);\n
    for j to nops(lmin) do\n\011sys:=\{op(ini),seq(subs(n=i,r),i=op(j,
lmin)..maxind-order)\};\n\011if sys=\{\} then\n\011    if not flag the
n RETURN(\{seq(u(i)=_C[i],i=minind..maxind)\})\n\011    else RETURN(\{
seq(u(i)=_C[i],\n\011\011i=\{$minind..maxind\} minus \{$0..order-1\})
\})\n\011    fi\n\011fi;\n\011a:=systomatrix(sys,[seq(u(i),i=minind..m
axind)],'b');\n\011# This traperror is due to a bug in solve/linear in
 V.5 (and earlier ?)\n\011# in the case of floating point coefficients
\n\011sol:=traperror(linalg[linsolve](a,b));\n\011if (sol=NULL or sol=
lasterror) and j=nops(lmin) then\n\011    ERROR(`no valid initial cond
itions`)\n\011elif sol<>NULL and sol<>lasterror then break fi\n    od;
\n    sol:=convert(sol,list);\n    inds:=indets(sol,_t[anything]);\n  \+
  # replace the _t[anything] by _C[anything] depending on flag\n    do
rej:=flag and (j=1);\n    j:=max(op(map(op,indets([rec,ini],_C[anythin
g]))));\n    if j=-infinity then j:=-1 fi;\n    for i in inds do\n\011
if member(i,sol,'k') and (not dorej or nops(select(has,sol,i))>1\n\011
\011or k<minind-1 or k+minind>order) then\n\011    j:=j+1;\n\011    so
l:=subs(i=_C[minind+j],sol);\n\011    rej[i]:=NULL\n\011else rej[i]:=k
\n\011fi\n    od;\n    sys:=\{seq(u(i+minind-1)=sol[i],i=\{$1..nops(so
l)\} minus\n\011\{seq(rej[i],i=inds)\})\};\n    # clean the _C[] for i
nitial conditions of the type _C[1]+_C[2],_C[1]-_C[2]\n    if hastype(
remove(type,sys,u(anything)=name),_C[anything])\n\011# limitations due
 to the cost of Groebner basis computation:\n\011and max(op(map(degree
,[seq(op(2,i),i=sys)],\n\011    indets(sys,_C[anything]))))<3 and not \+
has(sys,RootOf)\n\011# also it has to be a system of polynomials\n\011
and type([seq(op(2,i),i=sys)],\n\011    list(polynom(rational,indets(s
ys,_C[anything]))))\n    then\n\011inds:=indets(sys,_C[anything]);\n
\011inds2:=sort([op(map(op,indets(sys,u(anything))))]);\n\011sys:=subs
([seq(u(i)=u[i],i=inds2)],\{seq(op(1,i)-op(2,i),i=sys)\}):\n\011# this
 should be replaced by an elimination order\n\011sol:=subs(solve(\{op(
remove(hastype,grobner[gbasis](sys,[op(inds),\n\011    seq(u[i],i=inds
2)],'plex'),_C[anything]))\},\n\011\011\{seq(u[i],i=inds2)\}),[seq(u[i
],i=inds2)]);\n\011j:=-1;\n\011for i to nops(sol) do\n\011    if sol[i
]=u[inds2[i]] then\n\011\011if dorej and nops(select(has,sol,sol[i]))=
1 then\n\011\011    rej[i]:=i\n\011\011else\n\011\011    rej[i]:=NULL;
\n\011\011    j:=j+1;\n\011\011    sol:=subs(u[inds2[i]]=_C[j],sol);\n
\011\011fi\n\011    fi\n\011od;\n\011\{seq(u(inds2[i])=sol[i],i=\{$1..
nops(sol)\} minus\n\011\011\011\011    \{seq(rej[i],i=1..nops(sol))\})
\}\n    else sys\n    fi\nend: # `goodinitvalues/rec`\n\n# indicialeq
\n# Input: a differential equation in the format returned by formatdif
feq\n#\011 the unknown function and its variable\n#\011 a variable for
 the resulting polynomial\n#\011 (optional) the shift of exponent\n# O
utput: the indicial equation (a polynomial in the variable vanishing\n
#\011 for possible valuations of series solutions)\nindicialeq:=proc (
deq, z, alpha, val)\nlocal ldeg, i, v, res, j, dd;\n    # this is beca
use of the new convention for degree(0) in V.5\n    dd:=map(proc(x) if
 x=0 then 1 else x fi end,deq);\n    ldeg:=[seq(ldegree(dd[i],z)-i+2,i
=2..nops(deq))];\n    v:=min(op(ldeg));\n    if nargs=4 then val:=v fi
;\n    for i to nops(ldeg) do\n\011if ldeg[i]=v then\n\011    res[i]:=
tcoeff(deq[i+1],z)*convert([seq(alpha-j,j=0..i-2)],`*`)\n\011else res[
i]:=0 fi\n    od;\n    collect(convert([seq(res[i],i=1..nops(ldeg))],`
+`),alpha)\nend: # indicialeq\n\n# nbinicond\n# Input:  a differential
 equation in the format returned by formatdiffeq\n#\011  the unknown f
unction and its variable\n# Output: the number of initial conditions n
ecessary to specify a solution\n#\011  an error if no series solution \+
exist\n#\011  this is mostly useful in the singular case\nnbinicond:=p
roc(deq,y,z)\nlocal pol, val;\n    pol:=indicialeq(deq,z,z,val);\n    \+
if deq[1]=0 then firstnonzero(pol,z)-1\011\011# homogeneous case\n    \+
elif subs(z=ldegree(deq[1],z)-val,pol)=0 then\011\n\011if subs(z=0,deq
[nops(deq)])=0 then 0\n\011else ERROR(`no valid initial conditions`) f
i\n    elif nops(deq)=2 then -1\n    else ldegree(deq[1],z)-val-1 fi\n
end: # nbinicond\n\n# `goodinitvalues/diffeq`\n# Input:  a differentia
l equation in the format returned by formatdiffeq\n#\011  the unknown \+
function and its variable\n#\011  (optional) some initial conditions\n
#\011  (optional) an integer p\n# Output: a set of equalities (D@@k)(y
)(0)=v_k, from which all the others\n#\011   can be computed. This set
 is continued up to order p if p is given.\n#\011  The result is an ER
ROR when no initial condition can be found,\n#\011  except in the case
 when the origin is a singular point. Then it\n#\011  returns \{\}.\n`
goodinitvalues/diffeq`:=proc (deq, y, z, ini, p)\nlocal u, init, i, so
l, maxorder, ord, inds, j, rej, k, inds2, sys, v, gb, dorej;\n    maxo
rder:=nbinicond(deq,y,z);\n    if nargs=5 then maxorder:=max(maxorder,
p) fi;\n    ord:=maxorder;\n    if nargs>3 and type(ini,set) then\n
\011init:=ini;\n\011maxorder:=max(maxorder,seq(op(2,op(0,op(0,i))),\n
\011    i=indets(init,`gfun/initeq`(y)) minus \{y(0),D(y)(0)\}));\n
\011if maxorder<=0 and has(init,D(y)(0)) then maxorder:=1 fi\n    else
 init:=\{\} fi;\n    u:=series(eval(subs(y(z)=subs(init,convert([seq((
D@@i)(y)(0)*z^i/i!,\n\011i=0..maxorder),O(1)*z^(maxorder+1)],`+`)),\n
\011makediffeq(deq,y,z))),z,infinity);\n    u:=\{seq(coeff(u,z,i),i=0.
.maxorder)\} minus \{O(1)\} union init;\n    sol:=solve(u,\{seq((D@@i)
(y)(0),i=0..maxorder)\}\n\011intersect indets(u,function(0)));\n    if
 sol=NULL then  # it might be more singular than it looks\n\011if ord<
>-1 and subs(z=0,op(nops(deq),deq))<>0 then\n\011    ERROR(`no valid i
nitial conditions`)\n\011elif ord=-1 then \{\}\n\011else init\n\011fi
\n    elif ord=-1 then \{\}\n    else\n\011dorej:=evalb(subs(z=0,deq[n
ops(deq)])<>0);\n\011if dorej then sol:=subs(sol,[seq((D@@i)(y)(0),i=0
..ord)])\n\011else sol:=subs(sol,[seq((D@@i)(y)(0),i=0..maxorder)]) fi
;\n\011inds:=select(has,indets(sol,function),\n\011    \{seq((D@@i)(y)
(0),i=0..maxorder)\});\n\011j:=max(op(map(op,indets([deq,sol],_C[anyth
ing]))));\n\011if j=-infinity then j:=-1 fi;\n\011for i in inds do\n
\011    if member(i,sol,'k') and (# not dorej or\n\011\011nops(select(
has,sol,i))>1) then\n\011\011j:=j+1;\n\011\011sol:=subs(i=_C[j],sol);
\n\011\011rej[i]:=NULL\n\011    else rej[i]:=k\n\011    fi\n\011od;\n
\011inds2:=\{$1..nops(sol)\} minus \{seq(rej[i],i=inds)\};\n\011# clea
n the _C[] for initial conditions of the type\n\011# _C[1]+_C[2],_C[1]
-_C[2]\n\011if not hastype(deq,_C[anything]) and\n\011    hastype(remo
ve(type,[seq(sol[i],i=inds2)],name),_C[anything])\n\011    # limitatio
ns due to the cost of Groebner basis computation:\n\011    and max(op(
map(degree,sol,indets(sol,_C[anything]))))<3\n\011    and not has(sol,
RootOf)\n\011    # also it has to be a system of polynomials\n\011    \+
and type([seq(sol[i],i=inds2)],\n\011\011list(polynom(rational,indets(
sol,_C[anything]))))\n\011then\n\011    sys:=\{seq(v[i]-sol[i],i=inds2
)\};\n\011    inds:=indets(sys,_C[anything]);\n\011    inds2:=sort([op
(inds2)]);\n\011    # this should be replaced by an elimination order
\n\011    gb:=\{op(remove(hastype,grobner[gbasis](sys,[op(inds),\n\011
\011seq(v[i],i=inds2)],'plex'),_C[anything]))\};\n\011    if type(gb,s
et(linear([seq(v[i],i=inds2)]))) then\n\011\011sol:=subs(solve(gb,\{se
q(v[i],i=inds2)\}),[seq(v[i],i=inds2)]);\n\011\011j:=-1;\n\011\011for \+
i to nops(sol) do\n\011\011    if sol[i]=v[inds2[i]] then\n\011\011
\011if nops(select(has,sol,sol[i]))=1 and i<=maxorder then\n\011\011
\011    rej[i]:=i\n\011\011\011else\n\011\011\011    rej[i]:=NULL;\n
\011\011\011    j:=j+1;\n\011\011\011    sol:=subs(v[inds2[i]]=_C[j],s
ol);\n\011\011\011fi\n\011\011    fi\n\011\011od;\n\011    fi;\n\011  \+
  \{seq((D@@(inds2[i]-1))(y)(0)=sol[i],i=\{$1..nops(sol)\} minus\n\011
\011\011\011\011\{seq(rej[i],i=1..nops(sol))\})\}\n\011else \{seq((D@@
(i-1))(y)(0)=sol[i],i=inds2)\}\n\011fi\n    fi;\nend: # `goodinitvalue
s/diffeq`\n\n######################### Conversion Routines ###########
########\n\nseriestolist:=proc()\nlocal s, meth, l, x, i;\n    if args
[1]<>'stamped' then RETURN(seriestolist(typecheck(8,args))) fi;\n    s
:=args[2]; meth:=args[3]; x:=op(0,s);\n    if op(nops(s)-1,s)=O(1) the
n l:=[seq(coeff(s,x,i),i=0..op(nops(s),s)-1)]\n    else l:=[seq(coeff(
s,x,i),i=0..op(nops(s),s))] fi;\n    if meth='ogf' then l\n    else se
riestolist('stamped',listtoseries('stamped',l,x,meth),'ogf')\n    fi\n
end: # seriestolist\n\nlisttolist:=proc()\nlocal x;\n    if args[1]<>'
stamped' then listtolist(typecheck(9,args))\n    elif args[3]='ogf' th
en args[2]\n    else seriestolist('stamped',\n\011listtoseries('stampe
d',args[2],x,args[3]),'ogf')\n    fi\nend: #listtolist\n\nseriestoseri
es:=proc ()\n    if args[1]<>'stamped' then seriestoseries(typecheck(8
,args))\n    else listtoseries('stamped',seriestolist('stamped',args[2
],'ogf'),\n\011op(0,args[2]),args[3])\n    fi\nend: # seriestoseries\n
\nlisttoseries:=proc ()\n    if args[1]<>'stamped' then listtoseries(t
ypecheck(1,args))\n    else map(normal,listtoseriestable[args[4]](args
[2],args[3])) fi\nend: # listtoseries\n\nlisttoseriestable[egf]:=proc(
l,x)\nlocal i;\n    series(convert([seq(op(i,l)*x^(i-1)/(i-1)!,i=1..no
ps(l)),\n\011O(x^(nops(l)))],`+`),x,nops(l))\nend:\n\nlisttoseriestabl
e[Laplace]:=proc(l,x)\nlocal i;\n    series(convert([seq(op(i,l)*x^(i-
1)*(i-1)!,i=1..nops(l)),\n\011O(x^(nops(l)))],`+`),x,nops(l))\nend:\n
\nlisttoseriestable[ogf]:=proc(l,x)\nlocal i;\n    series(convert([seq
(op(i,l)*x^(i-1),i=1..nops(l)),\n\011O(x^(nops(l)))],`+`),x,nops(l))\n
end: # `listtoseries/ogf`\n\nlisttoseriestable[revogf]:=proc(L,x)\nloc
al l, i, nl;\n    l:=L;\n    while l<>[] and l[1]=0 do l:=subsop(1=NUL
L,l) od;\n    if l=[] then ERROR(`cannot revert 0 series`) fi;\n    nl
:=nops(l);\n    powrevert(series(convert([seq(l[i]*x^i,i=1..nl),x^(nl+
1)],\n\011`+`),x,nl+1),x,nl)\nend: # `listtoseries/revogf`\n\nlisttose
riestable[revegf]:=proc(L,x)\nlocal l, i, nl;\n    l:=L;\n    while l<
>[] and l[1]=0 do l:=subsop(1=NULL,l) od;\n    if l=[] then ERROR(`can
not revert 0 series`) fi;\n    nl:=nops(l);\n    powrevert(series(conv
ert([seq(l[i]*x^i/i!,i=1..nl),x^(nl+1)]\n\011,`+`),x,nl+1),x,nl)\nend:
 # `listtoseries/revegf`\n\nlisttoseriestable[lgdogf]:=proc(L,x)\nloca
l l, i, nl;\n    l:=L;\n    while l<>[] and l[1]=0 do l:=subsop(1=NULL
,l) od;\n    if l=[] then ERROR(`cannot revert 0 series`) fi;\n    nl:
=nops(l);\n    series(convert([seq(i*l[i+1]*x^(i-1),i=1..nl-1)],`+`)/
\n\011    convert([seq(l[i]*x^(i-1),i=1..nl)],`+`),x,nl-1)\nend: # `li
sttoseries/lgdogf`\n\nlisttoseriestable[lgdegf]:=proc(L,x)\nlocal l, i
, nl;\n    l:=L;\n    while l<>[] and l[1]=0 do l:=subsop(1=NULL,l) od
;\n    if l=[] then ERROR(`cannot revert 0 series`) fi;\n    nl:=nops(
l);\n    series(convert([seq(l[i]*x^(i-2)/(i-2)!,i=2..nl)],`+`)/\n\011
convert([seq(l[i]*x^(i-1)/(i-1)!,i=1..nl)],`+`),x,nl-1)\nend: # `listt
oseries/lgdegf`\n\nlisttodiffeq:=proc()\nlocal result, ex, methods, me
thod, y, x, s, unkn, expr;\n    if args[1]<>'stamped' then RETURN(list
todiffeq(typecheck(2,args))) fi;\n    expr:=args[2];unkn:=args[3];meth
ods:=args[4];\n    y:=op(0,unkn);x:=op(unkn);ex:=expr;\n    for method
 in methods do\n\011s:=traperror(listtoseries('stamped',ex,x,method));
\n\011if s=lasterror then next fi;\n\011userinfo(3,'gfun',`Trying the \+
`,method,s);\n\011result:=`s2d/s2d`(s,x,y);\n\011if result<>FAIL then
\n\011    userinfo(2,'gfun','The',method,'`seems to satisfy`',result);
\n\011    RETURN([inicond(s,result,y,x),method])\n\011fi\n    od;\n   \+
 FAIL\nend: # listtodiffeq\n\nseriestodiffeq:=proc ()\n    if args[1]<
>'stamped' then seriestodiffeq(typecheck(6,args))\n    else listtodiff
eq('stamped',\n\011seriestolist('stamped',args[2],'ogf'),args[3],args[
4]) fi\nend: # seriestodiffeq\n\nlisttoalgeq:=proc()\nlocal result, ex
, methods, method, y, x, s, unkn, expr;\n    if args[1]<>'stamped' the
n RETURN(listtoalgeq(typecheck(2,args))) fi;\n    expr:=args[2];unkn:=
args[3];methods:=args[4];\n    y:=op(0,unkn);x:=op(unkn);ex:=expr;\n  \+
  for method in methods do\n\011s:=traperror(listtoseries('stamped',ex
,x,method));\n\011if s=lasterror then next fi;\n\011userinfo(3,'gfun',
`Trying the `,method,s);\n\011result:=`s2a/s2a`(s,x,y);\n\011if result
<>FAIL then\n\011    userinfo(2,'gfun','The',method,'`seems to satisfy
`',result);\n\011    RETURN([result,method])\n\011fi\n    od;\n    FAI
L\nend: # listtoalgeq\n\nseriestoalgeq:=proc ()\n    if args[1]<>'stam
ped' then seriestoalgeq(typecheck(6,args))\n    else listtoalgeq('stam
ped',\n\011seriestolist('stamped',args[2],'ogf'),args[3],args[4]) fi\n
end: # seriestoalgeq\n\n# guessandcheck\n# Input:\n#   listseries\011a
 list of series in x\n#   x\011\011name of the variable\n#   ord\011
\011order of the Pade-Hermite approximant\n#   nbnonzero\011max number
 of non zero coefficients expected in the result\n#   tryfactors\011bo
olean (try whether factoring produces less nonzero coefficients)\n#   \+
ordcheck\011number of zero terms in the final series\nguessandcheck:=p
roc(listseries,x,ord,nbnonzero,tryfactors,ordcheck)\nlocal lpol, k, j,
 check;\n    lpol:=map(collect,pade2(listseries,x,ord),x);\n    if lpo
l=FAIL then RETURN(FAIL) fi;\n    if nops(subs(0=NULL,map(coeffs,lpol,
x)))>nbnonzero then\n\011if not tryfactors then RETURN(FAIL) fi;\n\011
lpol:=map(factors,lpol);\n\011if nops(subs(0=NULL,[seq(coeffs(op(1,k),
x),k=seq(op(op(2,k)),\n" }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 25188 "\011\011k=lpol))]))+nops(remove(
has,lpol,x))>nbnonzero then\n\011    RETURN(FAIL) fi;\n\011lpol:=map(c
onvert,[seq([k[1],seq(j[1]^j[2],j=k[2])],k=lpol)],`*`)\n    fi;\n    c
heck:=map(normal,series(convert([seq(lpol[j]*listseries[j],\n\011j=1..
nops(listseries))],`+`),x,ordcheck));\n    if check=0 or op(1,check)=O
(1) then RETURN(lpol) fi;\n    FAIL\nend: # guessandcheck\n\nlisttorat
poly:=proc()\nlocal result, ex, methods, method, s, x, bigO, ord, nbz,
 tryfactors, nbzeros;\n    if args[1]='stamped' then ex:=args[2];x:=ar
gs[3];methods:=args[4]\n    else RETURN(listtoratpoly(typecheck(3,args
))) fi;\n    tryfactors:=type(ex,list(rational));\n    bigO:=nops(ex);
\n    nbzeros:=nops(ex)-nops(subs(0=NULL,ex));\n    for method in meth
ods do\n\011s:=traperror(listtoseries('stamped',ex,x,method));\n\011if
 s=lasterror then next fi;\n\011userinfo(3,'gfun',`Trying the `,method
,s);\n\011if type(s,series) then ord:=order(s); nbz:=max(op(2,s),nbzer
os)\n\011else ord:=bigO-1; nbz:=nbzeros fi;\n\011result:=guessandcheck
([1,s],x,ord,ord-nbz-1,tryfactors,ord+3);\n\011if result<>FAIL then\n
\011    userinfo(2,'gfun','The',method,'`seems to be`',result);\n\011 \+
   RETURN([-result[1]/result[2],method])\n\011fi\n    od;\n    FAIL\ne
nd: # listoratpoly\n\nseriestoratpoly:=proc () # yes, it's stupid to c
onvert it to a list now.\n    if args[1]<>'stamped' then seriestoratpo
ly(typecheck(7,args))\n    else listtoratpoly('stamped',\n\011    seri
estolist('stamped',args[2],'ogf'),op(0,args[2]),args[3]) fi\nend: # se
riestoratpoly\n\n# s2d/s2d\n# Input: a series s, its variable x, a nam
e y for the function y(x) whose\n#\011 series it is.\n# Output: a line
ar differential equation satisfied by y(x), if possible.\n#\011  The o
rder of the differential equation is bounded by the global\n#\011  var
iable maxordereqn. The degree of the coefficients is bounded by\n#\011
  the global variable maxdegcoeff.\n`s2d/s2d`:=proc (s, x, y)\nlocal s
erie, i, bigO, lpol, j, maxord, nbz, k, tryfactors;\n    if op(nops(s)
-1,s)<>O(1) then bigO:=op(nops(s),s)\n    else bigO:=op(nops(s),s)-1 f
i; nbz:=op(2,s);\n    tryfactors:=type(subs(O(1)=0,[op(s)]),list(ratio
nal));\n    maxord:=min(bigO-1,maxordereqn);\n    for i from minordere
qn to maxord do\n\011# diff has option remember\n\011serie:=[1,s,seq(d
iff(s,x$k),k=1..i)];\n\011userinfo(4,'gfun',`Looking for differential \+
equation of order`,i);\n\011lpol:=guessandcheck(serie,x,min(bigO-i+1,(
maxdegcoeff+1)*(i+2)),\n\011    bigO-i-nbz,tryfactors,bigO+2);\n\011if
 lpol<>FAIL then\n\011    RETURN(convert([lpol[1],lpol[2]*y(x),\n\011
\011    seq(lpol[j]*diff(y(x),x$(j-2)),j=3..i+2)],`+`)) fi;\n\011if i=
maxord then RETURN(FAIL) fi\n    od;\n    FAIL\nend: # `s2d/s2d`\n\n# \+
s2a/s2a\n# Input: a series s, its variable x, a name y for the functio
n y(x) whose\n#\011 series it is.\n# Output: a polynomial equation sat
isfied by y(x), if possible.\n#\011  The degree in y of the equation i
s bounded by the global\n#\011  variable maxdegeqn. The degree of the \+
coefficients is bounded by\n#\011  the global variable maxdegcoeff.\n`
s2a/s2a`:=proc (s, x, y)\nlocal serie, i, bigO, lpol, j, nbz, tryfacto
rs;\n    if op(nops(s)-1,s)<>O(1) then bigO:=op(nops(s),s)\n    else b
igO:=op(nops(s),s)-1 fi; nbz:=op(2,s);\n    tryfactors:=type(subs(O(1)
=0,[op(s)]),list(rational));\n    serie:=map(series,[seq(s^i,i=0..mind
egeqn-1)],x,bigO+2);\n    for i from mindegeqn to maxdegeqn do\n\011us
erinfo(4,'gfun',`Looking for polynomial equation of degree`,i);\n\011s
erie:=[op(serie),series(serie[nops(serie)]*s,x,bigO+2)];\n\011lpol:=gu
essandcheck(serie,x,min(bigO+1,(maxdegcoeff+1)*(i+1)),\n\011    bigO-n
bz,tryfactors,bigO+2);\n\011if lpol<>FAIL then\n\011    RETURN(convert
([seq(lpol[j]*y(x)^(j-1),j=1..i+1)],`+`))\n\011fi\n    od;\n    FAIL\n
end: # `s2a/s2a`\n\nlisttorec:=proc()\nlocal result, methods, method, \+
u, n, s, unkn, expr;\n    if args[1]<>'stamped' then RETURN(listtorec(
typecheck(2,args))) fi;\n    expr:=args[2];unkn:=args[3];methods:=args
[4];\n    u:=op(0,unkn);n:=op(unkn);\n    for method in methods do\n
\011s:=traperror(listtolist('stamped',expr,method));\n\011if s=lasterr
or then next fi;\n\011userinfo(3,'gfun',`Trying the `,method,s);\n\011
result:=`l2r/l2r`(s,n,u);\n\011if result<>FAIL then\n\011    userinfo(
2,'gfun','The',method,'`seems to satisfy`',result);\n\011    RETURN([r
esult,method])\n\011fi\n    od;\n    FAIL\nend: # listtorec\n\nseriest
orec:=proc ()\n    if args[1]<>'stamped' then seriestorec(typecheck(6,
args))\n    else listtorec('stamped',\n\011    seriestolist('stamped',
args[2],'ogf'),args[3],args[4]) fi\nend: # seriestorec\n\n`l2r/l2r`:=p
roc (l,n,u)\nlocal i, s, z, j, nb, lpol, ordereqn, result, k, lvar, nb
inicond, inicoeff;\n    nb:=nops(l);\n    s:=[listtoseries(l,z),\n\011
seq(listtoseries([seq(l[j]*(j-1)^i,j=1..nb)],z),i=1..mindegcoeff-1)];
\n    for i from max(1,mindegcoeff) to min(maxdegcoeff,nops(l)-1) do\n
\011s:=[op(s),listtoseries([seq(l[j]*(j-1)^i,j=1..nb)],z)];\n\011useri
nfo(4,'gfun',\n\011    `Looking for recurrence with coefficients of de
gree`,i);\n\011lpol:=guessandcheck(s,z,min(nb,(i+1)*(maxordereqn+1)),n
b,false,nb+2);\n\011if lpol<>FAIL then\n\011    ordereqn:=max(op(map(d
egree,lpol,z)));\n\011    lvar:=[seq(u(n+j),j=0..ordereqn)];\n\011    \+
result:=convert([seq(\n\011\011convert([seq(coeff(lpol[k],z,j)*(n+orde
reqn-j)^(k-1),\n\011\011   k=1..i+1)],`+`)*u(n+ordereqn-j),j=0..ordere
qn)],`+`);\n\011    inicoeff:=subs(n=n-ordereqn,coeff(result,u(n+order
eqn),1));\n\011    nbinicond:=max(ordereqn,firstnonzero(inicoeff,n))-1
;\n\011    result:=collect(primpart(result,lvar),lvar);\n\011    RETUR
N(\{result,seq(u(j)=l[j+1],j=0..min(nbinicond,nb-1))\})\n\011fi\n    o
d;\n    FAIL\nend: # `l2r/l2r`\n\ninicond:=proc (s, eqn, y ,x)\nlocal \+
order, deq, i;\n    deq:=select(has,indets(eqn,'diff(anything,identica
l(x))'),y(x));\n    if deq=\{\} then eqn\n    else\n\011for order whil
e deq<>\{y(x)\} do deq:=subs(diff(y(x),x)=y(x),deq) od;\n\011\{eqn,seq
((D@@i)(y)(0)=coeff(s,x,i)*i!,i=0..order-2)\}\n    fi\nend: # inicond
\n\nlisttohypergeom:=proc()\nlocal result, methods, method, s, unkn, e
xpr;\n    if args[1]<>'stamped' then RETURN(listtohypergeom(typecheck(
3,args))) fi;\n    expr:=args[2];unkn:=args[3];methods:=args[4];\n    \+
for method in methods do\n\011s:=traperror(listtolist('stamped',expr,m
ethod));\n\011if s=lasterror then next fi;\n\011userinfo(3,'gfun',`Try
ing the `,method,s);\n\011result:=`l2h/l2h`(s,unkn);\n\011if result<>F
AIL then\n\011    userinfo(2,'gfun','The',method,'`seems to satisfy`',
result);\n\011    RETURN([result,method])\n\011fi\n    od;\n    FAIL\n
end: # listtohypergeom\n\nseriestohypergeom:=proc ()\n    if args[1]<>
'stamped' then seriestohypergeom(typecheck(7,args))\n    else listtohy
pergeom('stamped',\n\011    seriestolist('stamped',args[2],'ogf'),op(0
,args[2]),args[3]) fi\nend: # seriestohypergeom\n\n`l2h/l2h`:=proc (l,
 x)\nlocal a, a0, k, eqn, u, v, w, den, i, z, c;\n    a:=l;\n    for k
 while op(1,a)=0 do a:=subsop(1=NULL,a) od;\n    a0:=op(1,a);k:=k-1;\n
    if nops(a)<5 then RETURN(FAIL) fi;\n    a:=[seq(op(i,a)/a0,i=2..no
ps(a))];\n    eqn:=normal((6*a[4]*a[1]**2*a[2]+9*a[2]*a[3]**2+6*a[4]*a
[1]*a[3]-6*\n\011a[3]**2*a[1]**2+a[2]**2*a[3]*a[1]-16*a[4]*a[2]**2)*x*
*2+(-32*a[4]\n\011*a[2]**2+5*a[2]**2*a[3]*a[1]+6*a[4]*a[1]*a[3]+18*a[4
]*a[1]**2*\n\011a[2]+27*a[2]*a[3]**2-24*a[3]**2*a[1]**2)*x+6*a[2]**2*a
[3]*a[1]-\n\01118*a[3]**2*a[1]**2+12*a[4]*a[1]**2*a[2]);\n    if eqn=0
 then v:=1\n    elif degree(eqn,x)=0 then RETURN(FAIL)\n    else v:=op
(1,[solve(eqn,x)])\n    fi;\n    den:=normal(4*a[2]**2*v**2-3*a[3]*v**
2*a[1]-a[2]*v**2*a[1]**2-3*a[2]*v*\n\011a[1]**2+8*a[2]**2*v-3*a[3]*v*a
[1]-2*a[2]*a[1]**2);\n    if den=0 then RETURN(FAIL) fi;\n    w:=norma
l(-2*(2*a[2]**2*v+4*a[2]**2-3*a[3]*v*a[1]-3*a[3]*a[1])*v)/den;\n    if
 type(w,negint) or w=0 then RETURN(FAIL) fi;\n    z:=normal(-3*a[3]*v*
a[1]**2+a[1]*a[2]**2*v+3*a[3]*v*a[2]-3*a[3]*\n\011a[1]**2+2*a[2]**2*a[
1]);\n    if z=0 then RETURN(FAIL) fi;\n    u:=-normal(2*a[2]**2*v+4*a
[2]**2-3*a[3]*v*a[1]-3*a[3]*a[1])*a[1]/z;\n    z:=2*z/den;\n    userin
fo(3,'gfun',`candidate: hypergeom(`,[u,v],[w],z*x,`)`);\n    c:=u*(u+1
)*(u+2)*(u+3)*v*(v+1)*(v+2)*(v+3)/w/(w+1)/(w+2)/(w+3)*z^4/24;\n    for
 i from 5 to nops(a) do\n\011c:=c*(u+i-1)*(v+i-1)*z/(w+i-1)/i;\n\011if
 c<>op(i,a) then RETURN(FAIL) fi;\n    od;\n    userinfo(2,'gfun',`hyp
ergeom found, parameters:`,[u,v],[w],z*x);\n    RETURN(simplify(a0*x^k
*hypergeom([u,v],[w],z*x),hypergeom))\nend: # `l2h/l2h`\n\n#ratpolytoc
oeff\n# Input: a rational function of x\n#\011 x its variable\n#\011 n
 a name\n# Output: the nth coefficient of the Taylor expansion at the \+
origin of f.\nratpolytocoeff:=proc(f,x,n)\nlocal g;\n    g:=convert(f,
'fullparfrac',x,'sqrfree');\n    if type(g,`+`) then g:=[op(g)] else g
:=[g] fi;\n    convert(map(`ratpolytocoeff/elmt`,g,x,n),`+`)\nend: # r
atpolytocoeff\n\n`ratpolytocoeff/elmt`:=proc(g,x,n)\nlocal k, a, c, i;
\n    if type(g,function) then\n\011op(0,g)(`ratpolytocoeff/elmt`(op(1
,g),x,n),op(2..nops(g),g))\n    elif type(g,polynom(anything,x)) then \+
0\n    else\n\011# g must be c(a)*(x-a)^(-k)\n\011k:=select(has,indets
(g,`^`),x);\n\011if nops(k)<>1 then ERROR(`report this as a bug`,g,x,n
) fi;\n\011k:=op(k);\n\011a:=x-op(1,k);\n\011c:=g/k;\n\011k:=-op(2,k);
\n\011c/(-a)^k*a^(-n)*convert([seq(n+i,i=1..k-1)],`*`)/(k-1)!\n    fi
\nend: # `ratpolytocoeff/elmt`\n\nguesseqn:=proc ()\nlocal y, result, \+
l, x, methods, ll, i;\n    if args[1]<>'stamped' then RETURN(guesseqn(
typecheck(2,args))) fi;\n    l:=args[2];y:=op(0,args[3]);x:=op(args[3]
);methods:=args[4];\n    # First try to find a rational function\n    \+
userinfo(1,'gfun',`Trying to find a rational generating function`);\n \+
   for i in methods do\n\011ll[i]:=traperror(listtolist(l,i));\n\011if
 ll[i]=lasterror then next fi;\n\011result:=listtoratpoly('stamped',ll
[i],x,['ogf']);\n\011if result=FAIL then next fi;\n\011RETURN([denom(r
esult[1])*y(x)-numer(result[1]),i])\n    od;\n    # Then an algebraic \+
equation\n    userinfo(1,'gfun',`Trying to find an algebraic generatin
g function`);\n    for i in methods do\n\011result:=listtoalgeq('stamp
ed',ll[i],y(x),['ogf']);\n\011if result<>FAIL then RETURN([result[1],i
]) fi\n    od;\n    # Then a linear differential equation\n    userinf
o(1,'gfun',`Trying to find a linear differential equation`);\n    for \+
i in methods do\n\011result:=listtodiffeq('stamped',ll[i],y(x),['ogf']
);\n\011if result<>FAIL then RETURN([result[1],i]) fi\n    od;\n    FA
IL\nend: # guesseqn\n\nguessgf:=proc ()\nlocal interres, y, result, l,
 x, methods, inds, s, i, ll, j, ord, sol, mini, tmp;\n    if args[1]<>
'stamped' then RETURN(guessgf(typecheck(3,args))) fi;\n    l:=args[2];
x:=args[3];methods:=args[4];\n    # First try to find a rational funct
ion\n    userinfo(1,'gfun',`Trying to find a rational generating funct
ion`);\n    for i in methods do\n\011ll[i]:=traperror(listtolist(l,i))
;\n\011if ll[i]=lasterror then next fi;\n\011result:=listtoratpoly('st
amped',ll[i],x,['ogf']);\n\011if result<>FAIL then RETURN([result[1],i
]) fi\n    od;\n    # Then trap easy hypergeometrics\n    userinfo(1,'
gfun',`Trying to find an hypergeometric generating function`);\n    fo
r i in methods do\n\011result:=listtohypergeom('stamped',ll[i],x,['ogf
']);\n\011if result<>FAIL then RETURN([result[1],i]) fi\n    od;\n    \+
# Then algebraic functions\n    userinfo(1,'gfun',`Trying to find an a
lgebraic generating function`);\n    for i in methods do\n\011result:=
listtoalgeq('stamped',ll[i],y(x),['ogf']);\n\011if result=FAIL then ne
xt fi;\n\011userinfo(1,'gfun',`Trying to solve the equation`);\n\011so
l:=[solve(result[1],y(x))];\n\011if nops(sol)=1 then RETURN([sol,i]) f
i;\n\011ord:=0; mini:=0; s:=listtoseries(ll[i],x);\n\011for j to nops(
sol) do\n\011    tmp:=series(sol[j]-s,x,nops(ll[i]));\n\011    if tmp=
0 then RETURN([sol[j],i])\n\011    elif type(tmp,series) and op(2,tmp)
>ord then\n\011\011ord:=op(2,tmp); mini:=j\n\011    fi\n\011od;\n\011i
f mini=0 then next fi;\n\011RETURN([sol[j],i])\n    od;\n    # Then a \+
linear differential equation\n    userinfo(1,'gfun',`Trying to find a \+
linear differential equation`);\n    for i in methods do\n\011interres
:=listtodiffeq('stamped',ll[i],y(x),['ogf']);\n\011if interres=FAIL th
en next fi;\n\011userinfo(1,'gfun',`Trying to solve it`);\n\011result:
=dsolve(op(1,interres),y(x));\n\011if result<>FAIL and result<>NULL th
en\n\011    inds:=(indets(op(2,result),name) minus indets(l,name))\n
\011\011minus \{x,constants\};\n\011    if inds=\{\} then RETURN([op(2
,result),i]) fi;\n\011    s:=series(op(2,result),x,nops(l)+1);\n\011  \+
  s:=solve(\{seq(coeff(s,x,j-1)-op(j,l),j=1..nops(l))\},inds);\n\011  \+
  if s<>NULL and type(s,set) then\n\011\011RETURN(subs(s,[op(2,result)
,i])) fi\n\011fi\n    od;\n    FAIL\nend: # guessgf\n\n###############
######## Power Series Reversion #######################\n\n# p should \+
be a series with no constant term and a non-zero linear term.\npowreve
rt := proc (s::series,x,o)\nlocal p, v, k, pv, ppv;\n    v:=x/op(1,s);
\n    p:=convert(s,polynom);\n    k:=1;\n    while 2*k+1<=o do\n\011pv
:=powcompose(p,v,x,2*k+1);\n\011ppv:=pprimeknowingp(pv,v,x,2*k+1);\n
\011v:=v-powdivide(pv-x,ppv,x,2*k+1);\n\011k:=2*k+1\n    od;\n    if k
<o then\n\011pv:=powcompose(p,v,x,o);\n\011ppv:=pprimeknowingp(pv,v,x,
o);\n\011v:=v-powdivide(pv-x,ppv,x,o)\n    fi;\n    series(v+O(x^(o+1)
),x,o+1)\nend: # powrevert\n\n# pol, pol -> pol\npowcompose:=proc (Q,P
,x,n) # this assumes P(0)=0.\nlocal m, pm, pr, pr1, l, i, s, p, q;\n  \+
  p:=powtruncate(P,x,n);\n    q:=powtruncate(Q,x,n);\n    if n<9 then
\n\011s[0]:=coeff(Q,x,0);\n\011for i to degree(Q,x) do\n\011    s[i]:=
collect(coeff(Q,x,i)*powtruncate(P,x,n-i+1)^i,x)\n\011od;\n\011powtrun
cate(convert([seq(s[i],i=0..degree(Q,x))],`+`),x,n)\n    else\n\011m:=
isqrt(trunc(3.32192809*n/length(n)));\n\011pm:=powtruncate(p,x,m);\n
\011pr:=1;pr1:=p-pm;\n\011l:=powcomposesimple(q,pm,x,n);\n\011s[0]:=l;
\n\011for i to iquo(n,m)+1 do\n\011    l:=pprimeknowingp(l,pm,x,n-i);
\n\011    pr:=powtruncate(collect(pr*pr1,x),x,n);\n\011    s[i]:=colle
ct(l*pr/i!,x)\n\011od;\n\011powtruncate(convert([seq(s[i],i=0..iquo(n,
m)+1)],`+`),x,n)\n    fi\nend: # powcompose\n\n# pol, pol -> pol\npowc
omposesimple:=proc (q,p,x,n)\nlocal s, j, pk, i;\n    j:=degree(q,x);
\n    s:=1; while s<j do s:=s*2 od;\n    pk[1]:=powtruncate(p,x,n);\n \+
   i:=1; while i<s do\n\011pk[2*i]:=powtruncate(collect(pk[i]^2,x),x,n
); i:=2*i od;\n    powcomposesimpledoit(q,s,pk,x,n)\nend: # powcompose
simple\n\n# pol, pol -> pol\npowcomposesimpledoit:=proc (q, s, pk, x, \+
n)\nlocal q1, q2;\n    if s=8 then\n\011powtruncate(collect(coeff(q,x,
0)+pk[1]*(coeff(q,x,1)+\n\011coeff(q,x,3)*pk[2]+pk[4]*(coeff(q,x,5)+co
eff(q,x,7)*pk[2]))+\n\011pk[2]*(coeff(q,x,2)+coeff(q,x,6)*pk[4])+coeff
(q,x,4)*pk[4]+\n\011coeff(q,x,8)*pk[8],x),x,n)\n    elif s>8 then\n
\011q1:=powtruncate(q,x,s/2-1);\n\011q2:=collect((q-q1)/x^(s/2),x);\n
\011powtruncate(powcomposesimpledoit(q1, s/2, pk, x, n)+\n\011    coll
ect(pk[s/2]*\n\011\011powcomposesimpledoit(q2, s/2, pk, x, n-s/2),x),x
,n)\n    elif s=4 then\n\011powtruncate(collect(coeff(q,x,0)+pk[1]*(co
eff(q,x,1)+\n\011    coeff(q,x,3)*pk[2])+coeff(q,x,2)*pk[2]+coeff(q,x,
4)*pk[4],x),x,n)\n    elif s=2 then\n\011collect(coeff(q,x,0)+coeff(q,
x,1)*pk[1]+coeff(q,x,2)*pk[2],x)\n    else # ASSERTION s=1\n\011subs(x
=pk[1],q)\n    fi\nend: # powcomposesimpledoit\n\n# pol, pol -> pol\np
primeknowingp:=proc (poff, f, x, n)\n    powdivide(diff(poff,x),diff(f
,x),x,n)\nend: # pprimeknowingp\n\n# pol, pol -> pol\npowdivide:=proc \+
(p, q, x, n)\nlocal u, b, i, j;\n    # assuming q[0]<>0, otherwise a s
hift of this version is all we need\n    # This is slower that convert
(series(p/q,x,n+1),polynom), but\n    # uses less memory.\n    if subs
(x=0,q)=0 then ERROR(`not implemented`) fi;\n    for i from 0 to n do
\n\011b[i]:=coeff(q,x,i);\n\011u[i]:=(coeff(p,x,i)-convert([seq(u[j]*b
[i-j],j=0..i-1)],`+`))/b[0]\n    od;\n    convert([seq(u[i]*x^i,i=0..n
)],`+`)\nend: # powdivide\n\npowtruncate:=proc (pol, x, n)\nlocal i;\n
    if degree(pol,x)<=n then pol\n    elif ldegree(pol,x)>n then 0\n  \+
  else convert([seq(coeff(pol,x,i)*x^i,i=0..n)],`+`)\n    fi\nend: # p
owtruncate\n\n######################## Puiseux Expansions ############
##############\n\nalgeqtoseries:=proc(Pol::polynom,x::name,y::name,ord
::nonnegint,optional_positive_slopes)\n    map(`algeqtoseries/prettypr
int`,`algeqtoseries/doit`(args),x)\nend: # algeqtoseries\n\n`algeqtose
ries/doit`:=proc (Pol,x,y,ord,opt)\nlocal pol, a, u, i, j, pts, alpha,
 mini, deg, theta, jmin, sl, lastpt, p, pu, nb2, res, r, nb, q, a0, no
rmalizer, eq;\noption `Copyright Bruno Salvy, INRIA Rocquencourt, Fran
ce`;\n    _EnvExplicit:=false; # otherwise some RootOf's will be pertu
rbed\n    pol:=collect(Pol,[y,x],evala);\n    # Newton polygon\n    de
g:=degree(pol,y);\n    pts:=select(proc(x) x[2]<>0 end,[seq([i,coeff(p
ol,y,i)],i=0..deg)]);\n    pts:=[seq([i[1],ldegree(i[2],x)],i=pts)];\n
    nb:=0;\n    lastpt:=pts[1];\n    for i from 2 to nops(pts) do\n
\011mini:=infinity;\n\011for j from i to nops(pts) do\n\011    theta:=
(pts[j][2]-lastpt[2])/(pts[j][1]-lastpt[1]);\n\011    if theta<mini th
en mini:=theta; jmin:=j fi\n\011od;\n\011nb:=nb+1; alpha[nb]:=-mini; i
:=jmin; lastpt:=pts[i]\n    od;\n    # Treat each slope\n    nb2:=0;\n
    alpha:=\{seq(alpha[i],i=1..nb)\};\n    if nargs=5 then alpha:=sele
ct(type,alpha,nonnegative) fi;\n    for sl in alpha do\n\011r:=denom(s
l);\n\011p:=collect(subs(x=x^r,y=x^numer(sl)*y,pol),x);\n\011if ldegre
e(p,x)<>0 then p:=collect(p/x^ldegree(p,x),x) fi;\n\011q:=collect(coef
f(p,x,0),y);\n\011if ldegree(q,y)<>0 then q:=collect(q/y^ldegree(q,y),
y) fi;\n\011for u in sqrfree(q,y)[2] do\n\011    a0:=traperror(evala(R
ootOf(u[1],y)));#These 5 lines are not\n\011    if a0<>lasterror then \+
a0:=[a0]\011# strictly necessary but they\n\011    elif op(1,[a0])<>`r
educible RootOf detected.  Substitutions are`\n\011\011then ERROR(a0)
\011\011\011\011# save much trouble in\n\011    else a0:=map(subs,a0[2
],RootOf(u[1],y)) fi; # later computations.\n\011    for a0 in a0 do\n
\011\011if u[2]=1 then\011# regular case\n\011\011    if type(a0,RootO
f) then normalizer:=evala\n\011\011    else normalizer:=normal fi;\n
\011\011    a[0]:=a0; pu:=p;\n\011\011    for i to ord-1 do # ord rela
tively small\n\011\011\011pu:=collect(subs(y=a[i-1]+x*y,pu),x);\n\011
\011\011# do not use solve: it does strange things with RootOf's\n\011
\011\011# a[i]:=normalizer(solve(coeff(pu,x,i),y))\n\011\011\011eq:=co
eff(pu,x,i); # this is assumed to be linear in y\n\011\011\011if degre
e(eq,y)=1 then\n\011\011\011    a[i]:=normalizer(-coeff(eq,y,0)/coeff(
eq,y,1))\n\011\011\011else ERROR(`unforecast case`) fi\n\011\011    od
; # in this loop, normal saves memory\n\011\011    nb2:=nb2+1;\n\011
\011    res[nb2]:=[seq([a[i],sl+i/r],i=0..ord-1),[O(1),sl+ord/r]]\n
\011\011elif ord=0 then # several branches, O() term\n\011\011    to u
[2] do\n\011\011\011nb2:=nb2+1;\n\011\011\011res[nb2]:=[[O(1),sl]]\n
\011\011    od\n\011\011else\011\011# several branches\n\011\011    fo
r i in `algeqtoseries/doit`(\n\011\011\011    subs([x=x^u[2],y=a0+x*y]
,p),x,y,ord-1,1) do\n\011\011\011nb2:=nb2+1;\n\011\011\011res[nb2]:=[[
a0,sl],seq([j[1],sl+(j[2]+1)/u[2]/r],j=i)]\n\011\011    od\n\011\011fi
\n\011    od\n\011od\n    od;\n    [seq(res[i],i=1..nb2)]\nend: # gfun
/algeqtoseries/doit\n\n`algeqtoseries/prettyprint`:=proc (l, x)\nlocal
 i;\n    series(convert([seq(i[1]*x^i[2],i=l)],`+`),x,max(4,ceil(l[nop
s(l)][2])))\nend: # `algeqtoseries/prettyprint`\n\n\n#################
####### Holonomic Functions #########################\n\n# diffeqtohom
diffeq\n# Input:  deq: differential equation\n#\011  y(z): its unknown
 function\n# Output: a homogeneous linear differential equation having
 the solutions\n#         of deq as solutions\n#\ndiffeqtohomdiffeq:=p
roc (Deq, yofx)\nlocal deq, y, x, ini, c, dc, n, i;\n    deq:=formatdi
ffeq([args],y,x,ini);\n    if deq[1]=0 then Deq\n    else\n\011n:=nops
(deq);\n\011if ini<>\{\} then ini:=`goodinitvalues/diffeq`(deq,y,x,ini
,n-2) fi;\n\011c:=deq[1];\n\011dc:=-diff(c,x);\n\011makediffeq(map(col
lect,[0,dc*deq[2]+c*diff(deq[2],x),\n\011    seq(dc*deq[i]+c*diff(deq[
i],x)+c*deq[i-1],i=3..n),c*deq[n]],x),\n\011    y,x,ini)\n    fi\nend:
 # diffeqtohomdiffeq\n\n#algeqtodiffeq\n# Input:  a polynomial in two \+
variables y and z\n#\011  the unknown function y(z)\n#\011  (optional)
 a set of initial conditions\n# Output: a linear differential equation
 verified by RootOf(P,y)\n#\011  or FAIL if the initial conditions can
not be satisfied\nalgeqtodiffeq := proc()\nlocal y, z, P, inits, g, u,
 d, i, Y, deq, j, r, inity, P0;\n    P:=formatpoleq([args],'y','z','in
its');\n    g:=mygcdex(diff(P,y),P,y,'u');\n    if has(g,y) then RETUR
N(algeqtodiffeq(normal(P/g),y(z),inits)) fi;\n    d := degree(P,y); us
erinfo(3,'gfun',`degree is `,d);\n    if d<=1 then deq:=subs(y=y(z),P)
\n    elif not has(P,z) then deq:=y(z)-RootOf(P,y)\n    else\n\011Y[1]
:=rem(-u/g*diff(P,z),P,y);\n\011for i from 2 to d-1 do # compute Y[i] \+
= diff(y,z$i) mod P\n\011    Y[i]:=rem(collect(diff(Y[i-1],z)+diff(Y[i
-1],y)*Y[1],y),P,y) od;\n\011deq:=lindep(array(1..d+1,1..d,normal([[1,
0$(d-1)],[0,1,0$(d-2)],\n\011    seq([seq(coeff(Y[i],y,j),j=0..d-1)],i
=1..d-1)])),\n\011    [1,seq(diff(y(z),[z$i]),i=0..d-1)],z)\n    fi;\n
    userinfo(1,'gfun',`differential equation is`,deq);\n    P0:=normal
(subs(z=0,y=y(0),inits,P));\n    if not has(P0,y) and P0<>0 then # thi
s means that the origin\n\011\011\011\011    # is a singular point\n
\011if inits=\{\} then RETURN(deq)\n\011else ERROR(`invalid initial co
nditions`) fi\n    fi;\n    if P0<>0 then\n#     This will be the corr
ect way when Maple does not insist on the argument\n#     of RootOf be
ing irreducible.\n#\011inits:=inits union \{y(0)=RootOf(subs(y(0)=_Z,P
0))\} fi;\n\011if traperror(evala(RootOf(subs(y(0)=_Z,P0))))<>lasterro
r then\n\011    inits:=inits union \{y(0)=RootOf(subs(y(0)=_Z,P0))\} f
i;\n    fi;\n    inity:=y=subs(inits,y(0));\n    inits:=\{y(0)=op(2,in
ity)\} union subs(y(0)=op(2,inity),inits);\n    for i to d-1 do\n\011r
:=traperror((D@@i)(y)(0)=subs([z=0,inity],Y[i]));\n\011if r<>lasterror
 then inits:=inits union \{r\} else break fi\n    od;\n    inits:=`goo
dinitvalues/diffeq`(formatdiffeq([deq,y(z)]),y,z,\n\011select(type,ini
ts,anything=`gfun/free`(y(0))));\n    if inits=\{\} then RETURN(deq) e
lse RETURN(\{deq,op(inits)\}) fi\nend: # algeqtodiffeq\n\n#diffeqtorec
\n# Input:  eqn: differential equation (for example output of algeqtod
iffeq)\n#\011  y(z): its unknown function\n#\011  u(n): the name of th
e sequence of Taylor coefficients at the origin\n# Output: the linear \+
recurrence satisfied by u(n)\n#\ndiffeqtorec:=proc (eqn,yofz,uofk)\nlo
cal iniconds, f, y, z, u, k, Y, Z;\n    if nargs<>3 then ERROR(`wrong \+
number of arguments`) fi;\n    getname(uofk,u,k);\n    f:=formatdiffeq
([eqn,yofz],'y','z','iniconds');\n    # Avoid problems when u=y or u=z
,...\n    f:=subs([y=Y,z=Z],f); iniconds:=subs(y=Y,iniconds);\n    `di
ffeqtorec/doit`(f,Y,Z,u,k,iniconds)\nend: # diffeqtorec\n\n`diffeqtore
c/doit` := proc(R,y,z,u,k,iniconds)\nlocal l, ini, i, rec, j, minordre
c, maxordrec, m, r, dr1, inhdeg, inhpart, p, rr;\n    if has(R,k) then
 ERROR(k,`cannot appear in the differential equation`)\n    elif has(R
,u) then ERROR(u,`cannot appear in the differential equation`)fi;\n   \+
 # initial conditions\n    l:=\{seq(op(2,op(0,op(0,i))),i=indets([inic
onds,R],`gfun/initeq`(y))\n\011minus \{y(0),D(y)(0)\})\};\n    ini:=[y
(0)=u(0),D(y)(0)=u(1),seq((D@@i)(y)(0)=u(i)*i!,i=l)];\n    r:=subs(ini
,R); ini:=subs(ini,iniconds);\n    # In very special cases, this loop \+
makes it possible to return\n    # an inhomogeneous equation of lower \+
order.\n    # Ex: z*(-1+z)^3*(D@@2)(y)(z)+(-1+z)^3*D(y)(z)-(-1+z)^3*y(
z)-z*(z-3)\n    if r[1]<>0 then for inhdeg from 0 do\n\011p:=1-z;\n
\011for i from 2 to nops(r) while degree(p,z)=1 do p:=gcd(p,r[i]) od;
\n\011if not has(p,z) then break fi;\n\011r:=[r[1],seq(quo(r[i],1-z,z)
,i=2..nops(r))]\n    od fi;\n    # main loop\n    ## this fixes a prob
lem with the change in the definition of degree(0)\n    ## in V.5\n   \+
 rr:=map(proc(x) if x=0 then 1 else x fi end,r);\n    minordrec:=min(s
eq(i-degree(op(i+2,rr),z),i=0..nops(r)-2));\n    maxordrec:=max(seq(i-
ldegree(op(i+2,rr),z),i=0..nops(r)-2));\n    rec:=array(sparse,minordr
ec..maxordrec);\n    for i from 2 to nops(r) do\n\011for j from ldegre
e(op(i,r),z) to degree(op(i,r),z) do\n\011    rec[i-2-j]:=rec[i-2-j]+c
oeff(op(i,r),z,j)*\n\011\011expand(convert([seq(k+m,m=1-j..i-2-j)],`*`
)) od od;\n    # inhomogeneous part of the differential equation\n    \+
if r[1]=0 then dr1:=-1; inhpart:=0\n    else\n\011dr1:=degree(r[1],z);
\n\011if inhdeg<>0 then # r[1] stands for r[1]/(1-z)^inhdeg\n\011    i
nhpart:=expand(convert([seq(coeff(r[1],z,i)*convert(\n\011\011[seq(k+j
,j=1-i..inhdeg-i-1)],`*`),i=0..dr1)],`+`))\n\011\011/(inhdeg-1)!;\n
\011    r:=subsop(1=series(r[1]/(1-z)^inhdeg,z,\n\011\011max(dr1,maxor
drec-minordrec)+1),r)\n\011else\n\011    inhpart:=0\n\011fi\n    fi;\n
    ini:=\{op(ini),op(map(convert,[seq(subs(k=i,[coeff(r[1],z,i),seq(
\n\011rec[j]*u(i+j),j=max(minordrec,-i)..maxordrec)]),i=0..dr1)],`+`))
\};\n    for i from dr1+1 while i<-minordrec or\n\011subs(k=i,\{seq(re
c[maxordrec-j],j=0..(i-dr1-1))\})=\{0\} do\n\011ini:=\{op(ini),convert
(subs(k=i,[seq(rec[j]*u(i+j),\n\011    j=max(minordrec,-i)..maxordrec)
,inhpart]),`+`)\} od;\n    rec:=listprimpart(\n\011subs(k=k-minordrec,
[inhpart,seq(rec[i],i=minordrec..maxordrec)]),k);\n    while nops(rec)
>2 and rec[nops(rec)]=0 do\n\011rec:=subsop(nops(rec)=NULL,rec) od;\n \+
   ini:=`goodinitvalues/rec`(rec,u,k,ini,true);\n    makerec(map(colle
ct,rec,k),u,k,ini)\nend: # `diffeqtorec/doit`\n\n#rectohomrec\n# Input
: a recurrence and its unknown function u(n)\n# Output: a homogeneous \+
recurrence cancelling the solutions of the original one\n#\nrectohomre
c:=proc (Rec,uofk)\nlocal rec, u, k, ini, n, i, c, dc;\n    rec:=forma
trec([args],u,k,ini);\n    if rec[1]=0 then Rec\n    else\n\011n:=nops
(rec);\n\011if ini<>\{\} then ini:=`goodinitvalues/rec`(rec,u,k,ini,tr
ue,n-2) fi;\n\011c:=rec[1];\n\011dc:=-subs(k=k+1,c);\n\011makerec(map(
collect,[0,dc*rec[2],\n\011    seq(dc*rec[i]+c*subs(k=k+1,rec[i-1]),i=
3..n),\n\011    c*subs(k=k+1,rec[n])],k),u,k,ini)\n    fi\nend: # rect
ohomrec\n\n# This one is for the case of constant coefficients\n# exce
pt possibly in the inhomogeneous part. It returns\n# a homogeneous rec
urrence with constant coefficients and order\n# the order of the origi
nal one plus the degree of the inhomogeneous term.\nrectohomrecbis:=pr
oc(Rec,uofk)\nlocal rec, u, k, ini, n, i, c, dc, co, j;\n    rec:=form
atrec([args],u,k,ini);\n    if rec[1]=0 then Rec\n    elif has(subsop(
1=NULL,rec),k) then ERROR(`invalid recurrence`,Rec)\n    else\n\011n:=
nops(rec);\n\011c:=collect(rec[1],k);\n\011dc:=degree(c,k);\n\011if in
i<>\{\} then ini:=`goodinitvalues/rec`(rec,u,k,ini,true,n-2+dc) fi;\n
\011co[-1]:=0;\n\011for i from 0 to n-2 do co[i]:=rec[i+2] od;\n\011fo
r i from 0 to dc do\n\011    co[n+i-1]:=co[n+i-2];\n\011    for j from
 n+i-2 by -1 to 0 do co[j]:=co[j-1]-co[j] od\n\011od;\n\011makerec([0,
seq(co[j],j=0..n+dc-1)],u,k,ini)\n    fi\nend: # rectohomrecbis\n\n" }
}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }{MPLTEXT 1 0 1002 "\n#rectopr
oc\n# Input:  a recurrence and its unknown function u(n).\n#\011  (opt
ional) 'remember'\n#\011  (optional) 'list'\n#\011  (optional) 'params
'=[n1,n2,...]\n# Output: a procedure of n,n1,n2,... computing values o
f the sequence up to\n#\011  the index n or a list of these vaues if t
he option 'list' is given.\n#\n#   If 'remember' is given, then the pr
ocedure uses option remember to\n# compute values of the sequence in l
inear time and linear space.\n# Otherwise, if the recurrence has const
ant coefficients, then binary\n# powering is used to make the procedur
e logarithmic time/logarithmic space.\n# If the coefficients are not c
onstants, then the procedure is a simple loop\n# that will compute the
 values in linear time/constant space. The same loop\n# is used in the
 'list' case.\n#\n# The variables that are reserved for the produced c
ode are lowercase\n# all the other variables are uppercase.\nrectoproc
 := proc(expr,yofn)\nlocal T,Y,N,R,INITS,ARGLIST,RECSTAT,A,U,IIg,REMBR
,ORDER,N0,n,a,l,res,i, INIVECT,INIMAT, RES, " }{TEXT -1 1 "\n" }
{MPLTEXT 1 0 4200 "RECLOOP, THRESHOLD, LIST, PARAMS, u, TOSUBS, J, SAR
GS,j,k,m;\n    if nargs<2 then ERROR(`wrong number of parameters`) fi;
\n    R:=formatrec([args[1..2]],'Y','N','INITS');\n    if nargs>2 then
\n\011SARGS:=\{args[3..nargs]\};\n\011if member('remember',SARGS) then
\n\011    REMBR:=['remember']; SARGS:=SARGS minus \{'remember'\}\n\011
else REMBR:=[] fi;\n\011LIST:=member('list',SARGS);\n\011if LIST then \+
SARGS:=SARGS minus \{'list'\} fi;\n\011PARAMS:=op(select(type,SARGS,id
entical('params')=list(name)));\n\011if PARAMS<>NULL then\n\011    SAR
GS:=SARGS minus \{PARAMS\};\n\011    TOSUBS:=seq(op(IIg,op(2,PARAMS))=
b||IIg,IIg=1..nops(op(2,PARAMS)));\n\011    PARAMS:=seq(b||IIg,IIg=1..
nops(op(2,PARAMS)))\n\011else TOSUBS:=NULL\n\011fi;\n\011if SARGS<>\{
\} then ERROR(`invalid arguments`,op(SARGS)) fi\n    else REMBR:=[]; L
IST:=false; PARAMS:=NULL; TOSUBS:=NULL fi;\n    ORDER:=nops(R)-2;\n   \+
 INITS:=subs([TOSUBS],`goodinitvalues/rec`(R,Y,N,INITS,false));\n    N
0:=nops(INITS);\n    R:=subs([N=N-ORDER,TOSUBS],listprimpart(R,N));\n \+
   ARGLIST:=[n,PARAMS];\n    RECSTAT:=-convert(convert([R[1],seq(R[IIg
+2]*U[IIg],IIg=0..ORDER-1)],`+`),\n\011horner,N)/convert(R[nops(R)],ho
rner,N);\n    if REMBR<>[] and not LIST then ##### Easy case: linear t
ime/linear space###\n\011if nops(ARGLIST)=1 then\n\011    T:=subsop(4=
NULL,procmake(`&proc`(ARGLIST,[],['remember'],\n\011\011subs([N=n,seq(
U[IIg]=`&args`[-2](n-ORDER+IIg),IIg=0..ORDER-1)],\n\011\011RECSTAT))))
\n\011else\n\011    T:=subsop(4=NULL,procmake(`&proc`(ARGLIST,[],['rem
ember'],\n\011\011subs(INITS,`&if`(seq(op([n=IIg,Y(IIg)]),IIg=0..N0-1)
,\n\011\011subs([N=n,seq(U[IIg]=`&args`[-2](n-ORDER+IIg,PARAMS),\n\011
\011    IIg=0..ORDER-1)],\n\011\011RECSTAT))))))\n\011fi\n    elif LIS
T then ######## Again linear time/linear space###\n\011RECLOOP:=`&stat
seq`(op(subs(INITS,\n\011    [seq(`&:=`(u[IIg],Y(IIg)),IIg=0..N0-1)]))
,\n\011    `&for`(i,N0,1,n,true,\n\011\011`&:=`(u[i],\n\011\011    sub
s([N=i,seq(U[IIg]=u[i-ORDER+IIg],IIg=0..ORDER-1)],RECSTAT))),\n\011   \+
 ['seq'(u[i],i=0..n)]);\n\011if nops(ARGLIST)=1 then\n\011    T:=procm
ake(`&proc`(ARGLIST,[i,u],REMBR,RECLOOP))\n\011else\n\011    T:=procma
ke(`&proc`(ARGLIST,[i,u],REMBR,\n\011\011subs(INITS,`&if`(seq(op([n=II
g,[seq(Y(J),J=0..IIg)]]),IIg=0..N0-1),RECLOOP))))\n\011fi\n    else\n
\011RECLOOP:=`&statseq`(op(subs(INITS,\n\011    [seq(`&:=`(evaln(u||II
g),Y(N0-ORDER+IIg)),IIg=0..ORDER-1)])),\n\011    `&for`(i,N0,1,n-1,tru
e,`&statseq`(\n\011\011`&:=`(evaln(u||ORDER),\n\011\011    subs([N=i,s
eq(U[IIg]=evaln(u||IIg),IIg=0..ORDER-1)],RECSTAT)),\n\011\011seq(`&:=`
(evaln(u||IIg),evaln(u||(IIg+1))),IIg=0..ORDER-1))),\n\011    subs([N=
n,seq(U[IIg]=evaln(u||IIg),IIg=0..ORDER-1)],RECSTAT));\n\011if nops(AR
GLIST)>1 then\n\011    RECLOOP:=subs(INITS,\n\011\011`&if`(seq(op([n=I
Ig,Y(IIg)]),IIg=0..N0-1),RECLOOP)) fi;\n\011if has(subsop(1=NULL,R),N)
 then##linear time/constant space###########\n\011    T:=procmake(`&pr
oc`(ARGLIST,\n\011\011[i,seq(evaln(u||IIg),IIg=0..ORDER)],[],RECLOOP))
\n\011else\011\011 ######## logarithmic time/constant space#######\n
\011    if R[1]<>0 then\n\011\011RETURN(procname(rectohomrecbis(expr,y
ofn),yofn)) fi;\n\011    THRESHOLD:=round(evalf(-2*ORDER^2/ln(2)*\n
\011\011LambertW(-1,-1/2*ln(2)/ORDER^2)));\n\011    INIMAT:=[[seq(-R[n
ops(R)-IIg]/R[nops(R)],IIg=1..ORDER)],\n\011\011    seq([0$(IIg-1),1,0
$(ORDER-IIg)],IIg=1..ORDER-1)];\n\011    INIVECT:=subs(INITS,[seq(Y(N0
-IIg),IIg=1..ORDER)]);\n\011    A:=array(1..ORDER,1..ORDER);RES:=array
(1..ORDER);\n\011    T:=procmake(`&proc`(ARGLIST,\n\011\011[a,res,i,j,
k,m,l,seq(evaln(u||IIg),IIg=0..ORDER)],[],\n\011\011`&if`(n<=THRESHOLD
,\n\011\011    RECLOOP,\n\011\011`&statseq`(\n\011\011    `&:=`(a,'arr
ay'(1..ORDER,1..ORDER,INIMAT)),\n\011\011    `&:=`(res,'array'(1..ORDE
R,INIVECT)),\n\011\011    `&:=`(l,'convert'(n-N0+1,base,2)),\n\011\011
    `&if`(l[1]=1,\n\011\011\011`&:=`(res,'array'(1..ORDER,op(3,linalg[
'multiply'](\n\011\011\011    array(INIMAT),array(1..ORDER,INIVECT))))
)),\n\011\011    `&if`(l=[1],\n\011\011\011res[1],\n\011\011    `&stat
seq`(\n\011\011\011`&for`(i,'subsop'(1=NULL,'nops'(l)=NULL,l),true,\n
\011\011\011`&statseq`(\n\011\011\011    `&:=`(a,'array'(1..ORDER,1..O
RDER,[\n\011\011\011\011'seq'(['seq'('add'(a[j,m]*a[m,k],m=1..ORDER),
\n\011\011\011\011    k=1..ORDER)],j=1..ORDER)])),\n\011\011\011    `&
if`(i=1,`&:=`(res,'array'(1..ORDER,\n\011\011\011\011['seq'('add'(a[j,
m]*res[m],m=1..ORDER),\n\011\011\011\011    j=1..ORDER)]))))),\n\011
\011\011'add'('add'(a[1,j]*a[j,i],j=1..ORDER)*res[i],\n\011\011\011   \+
 i=1..ORDER)))))))\n\011fi\n    fi;\n    # put the initial conditions \+
in the remember table\n    if not LIST and nops(ARGLIST)=1 then\n\011f
or IIg in INITS do T(op(op(1,IIg))):=op(2,IIg) od;\n    fi;\n    op(T)
\nend: # rectoproc\n\n" }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 1 "\n" }
{MPLTEXT 1 0 11353 "#rectodiffeq\n# Input:  expr: a linear recurrence \+
(with or without initial conditions)\n#\011  a(n): its unknown functio
n\n#\011  f(t): the function sum(a(n)*t^n,n=0..infinity)\n# Output: th
e linear differential equation satisfied by f(t).\n#\nrectodiffeq := p
roc(expr,aofn,foft)\nlocal r, a, n, f, t, iniconds, A, N;\n    if narg
s<>3 then ERROR(`wrong number of arguments`) fi;\n    getname(foft,f,t
);\n    r:=formatrec([expr,aofn],'a','n','iniconds');\n    iniconds:=`
goodinitvalues/rec`(r,a,n,iniconds,false);\n    # Avoid problems when \+
a=f or a=t or n=f or ...\n    r:=subs([a=A,n=N],r); iniconds:=subs([a=
A,n=N],iniconds);\n    `rectodiffeq/doit`(listprimpart(r,N),A,N,f,t,in
iconds)\nend: # rectodiffeq\n\n`rectodiffeq/doit`:=proc (r,u,n,f,z,ini
conds)\nlocal order, diffeq, P, k, p, a, rr, ini, i, k0, inds, res, l,
 c, aa, j;\n    if has(r,f) then ERROR(f,`cannot appear in the recurre
nce`)\n    elif has(r,z) then ERROR(z,`cannot appear in the recurrence
`) fi;\n    order:=max(op(map(degree,r,n)));\n    diffeq:=array(sparse
,-1..order);\n    k0:=nops(r)-2;\n    D(f):=D(f);\n    # To keep polyn
omial coefficients, multiply by z^(nops(r)-2)\n    for k from 0 to k0 \+
do\n\011P:=op(k+2,r);\n\011for p from 0 to degree(P,n) do\n\011    a:=
subs(n=p-k,P);\n\011    P:=quo(P-a,n+k-p,n);\n\011    diffeq[p]:=diffe
q[p]+a*z^(k0+p-k);\n\011    if p>0 then\n\011\011diffeq[-1]:=diffeq[-1
]-collect(a*z^(k0+p-k)*diff(\n\011\011    convert([seq((D@@i)(f)(0)*z^
i/i!,i=p..k-1)],`+`),z$p),z)\n\011    else\n\011\011diffeq[-1]:=diffeq
[-1]-collect(a*z^(k0+p-k)*convert(\n\011\011[seq((D@@i)(f)(0)*z^i/i!,i
=p..k-1)],`+`),z) fi od od;\n    P:=op(1,r);\n    for p from 0 to degr
ee(P,n) do\n\011rr[p]:=subs(n=-p-1,P);\n\011P:=quo(P-rr[p],(n+p+1)/(p+
1),n) od;\n    diffeq:=[diffeq[-1]*(1-z)^p+z^k0*convert([seq(rr[k]*(1-
z)^(p-k-1),\n\011k=0..p-1)],`+`),seq((1-z)^p*diffeq[k],k=0..order)];\n
    # remove apparent singularity at the origin\n    k:=min(op(map(lde
gree,diffeq,z)));\n    if k>0 then diffeq:=map(quo,diffeq,z^k,z) fi;\n
    # initial conditions\n    inds:=map(op,indets(iniconds,u(anything)
));\n    ini:=solve(subs([seq(u(i)=(D@@i)(f)(0)/i!,i=inds)],iniconds),
\n\011\{seq((D@@i)(f)(0),i=inds)\});\n    diffeq:=subs(ini,diffeq);\n \+
   # some initial conditions may correspond to polynomial inhomogeneit
ies\n    # at least when the equation is not singular at the origin.\n
    inds:=max(op(inds));\n    if # subs(z=0,diffeq[nops(diffeq)])<>0 a
nd\n    inds<>-infinity and inds>=order\n\011then\n\011diffeq:=subsop(
1=collect(diffeq[1],z)-convert(series(eval(subs(f(z)=\n\011    convert
([seq((D@@i)(f)(0)*z^i/i!,i=0..inds),O(1)*z^(inds+1)],`+`)\n\011    ,i
ni,makediffeq(diffeq,f,z))),z,infinity),polynom),diffeq)\n    fi;\n   \+
 # Case when the inhomogeneous part contains (D@@k)(f)(0)\n    # Then \+
make the equation homogeneous.\n    if has(diffeq[1],f) then\n\011l:=s
elect(has,indets(diffeq[1],function),f);\n\011if not type(diffeq[1],li
near(l)) then\n\011    ERROR(`invalid inhomogeneous part`) fi;\n\011di
ffeq:=subs(\{seq(l[i]=aa[i],i=1..nops(l))\},diffeq);\n\011l:=[seq(aa[i
],i=1..nops(l))];\n\011diffeq:=subsop(1=collect(diffeq[1],l),diffeq);
\n\011for i in l do\n\011    if has(diffeq[1],i) then\n\011\011c:=coef
f(diffeq[1],i,1);\n\011\011res:=-diff(c,z);\n\011\011diffeq:=[res*diff
eq[1]+c*diff(diffeq[1],z),\n\011\011\011 res*diffeq[2]+c*diff(diffeq[2
],z),\n\011\011\011 seq(res*diffeq[j]+c*diff(diffeq[j],z)+\n\011\011
\011    c*diffeq[j-1],j=3..nops(diffeq)),\n\011\011\011 c*diffeq[nops(
diffeq)]];\n\011\011diffeq:=subsop(1=collect(diffeq[1],l,normal),diffe
q)\n\011    fi\n\011od;\n\011if has(diffeq,aa) then ERROR(`some assert
ion was wrong`) fi;\n\011diffeq:=map(collect,diffeq,z)\n    fi;\n    d
iffeq:=listprimpart(diffeq,z);\n    ini:=`goodinitvalues/diffeq`(diffe
q,f,z,ini);\n    res:=makediffeq(subs(ini,diffeq),f,z);\n    if ini<>
\{\} then \{res,op(ini)\} else res fi\nend: # `rectodiffeq/doit`\n\n# \+
poltodiffeq\n#Input: a polynomial P in n+1 variable z, y1, y2, ..., yn
\n#\011a list of n differential equations\n#\011a list of n variables \+
y1, y2, ..., yn\n#\011a variable y(z)\n#Output: a differential equatio
n in y(z) satisfied by P(z,y1(z),...,yn(z))\n#        where yi is a so
lution to the ith differential equation.\n#\n# It is important to avoi
d using D, otherwise too much time is spent\n# converting from diff to
 D, and back.\n#\npoltodiffeq:=proc (p, ldeq::list, ly::list, yofz)\nl
ocal i, deqs, ord, ini, y, z, lco, vars, n, h, u, tosubs, origdiff, j,
 k, l2, l, yy, dorder, jmax, inds, pp, linds, v, mon, lind, lcoef, lma
x, basis, back, uu, sing;\n    if nops(ldeq)<>nops(ly) then\n\011ERROR
(`not the same number of elements`,ldeq,ly) fi;\n    getname(yofz,y,z)
;\n    n:=nops(ldeq);\n    tosubs:=subsop(4=NULL,proc() args end):\n  \+
  for i to n do\n\011deqs[i]:=formatdiffeq([ldeq[i],ly[i]],'yy','z',in
i[i]);\n\011ord[i]:=nops(deqs[i])-2;\n\011lco:=deqs[i][nops(deqs[i])];
\n\011mon:=diff(y[i](z),[z$ord[i]]);\n\011v:=-deqs[i][1]/lco-\n\011   \+
 add(deqs[i][j]/lco*diff(y[i](z),[z$(j-2)]),j=2..ord[i]+1);\n\011tosub
s(mon):=v;\n\011lcoef[mon]:=[coeffs(v,[seq(diff(y[i](z),[z$j]),j=0..or
d[i]-1)],lind[mon])]\n    od;\n    if has(p,D) then pp:=convert(p,diff
) else pp:=p fi;\n    u:=subs([seq(op(0,ly[i])=y[i],i=1..n)],pp);\n   \+
 linds:=\{seq(diff(y[i](z),[z$ord[i]]),i=1..n)\};\n    origdiff:=[seq(
i=tosubs(i),i=linds)];\011# initial rewriting rules\n    # treat the d
iff's in u\n    pp:=select(has,indets(u,specfunc(anything,diff)),\{seq
(y[i],i=1..n)\});\n    for lmax from 0 while pp<>\{\} do pp:=select(ha
s,map(op,pp),diff) od;\n    while has(u,linds) do u:=eval(subs(origdif
f,u)) od;\n    pp:=u;\n    vars:=\{seq(seq(diff(y[i](z),[z$j]),j=0..or
d[i]-1),i=1..n)\};\n    if not type(u,polynom(ratpoly(anything,z),vars
)) then\n\011if nops(ly)=1 then\n\011    ERROR(p,`not a polynomial in \+
`,op(ly),`and its derivatives`)\n\011else ERROR(p,`not a polynomial in
 `,op(ly),`and their derivatives`) fi\n    fi;\n    basis:=vars;\n    \+
linds:=\{1\} union linds;\n    for dorder do # this loop includes an i
ncremental Gaussian elimination\n\011# first the derivatives of high o
rder are rewritten\n\011# and the polynomial is collected in the small
er order derivatives\n\011u:=collect(subs(origdiff,u),vars,'distribute
d',normal);\n\011u:=[coeffs(u,vars,'l')];l:=[l];\n\011# then those mon
omials corresponding to previous lines of the Gaussian\n\011# eliminat
ion are eliminated\n\011u:=collect(convert([seq(u[j]*tosubs(l[j]),j=1.
.nops(u))],`+`),vars,\n\011    'distributed',normal);\n\011# the list \+
of remaining monomials will be stored in l:\n\011lco:=[coeffs(u,vars,'
l')];\n\011# linds is the set of monomials corresponding to previous l
ines\n\011l2:=\{l\} minus linds;\n\011# if no other monomial remains, \+
then elimination is finished\n\011if l2=\{\} then break fi;\n\011l:=[l
];\n\011# the monomial which will be selected as a pivot might be invo
lved in\n\011# previous lines. In this case, we also compute the back \+
substitution.\n\011back:=evalb(l2 minus basis=\{\});\n\011if not back \+
then mon:=op(1,l2 minus basis) else mon:=l2[1] fi;\n\011member(mon,l,'
k');\n\011# new line\n\011tosubs(mon):=collect((diff(h(z),[z$(dorder-1
)])-convert([seq(lco[j]*\n\011    tosubs(l[j]),j=1..k-1),seq(lco[j]*to
subs(l[j]),j=k+1..nops(l))],\n\011    `+`))/lco[k],vars,'distributed')
;\n\011lcoef[mon]:=[coeffs(tosubs(mon),vars,lind[mon])];\n\011basis:=b
asis union \{lind[mon]\};\n\011# back substitution\n\011if back then\n
\011    for i in linds do\n\011\011if member(mon,[lind[i]],'k') then\n
\011\011    lind[i]:=[lind[i]];\n\011\011    tosubs(i):=collect(conver
t([\n\011\011\011seq(lcoef[i][j]*lind[i][j],j=1..k-1),\n\011\011\011  \+
  lcoef[i][k]*tosubs(mon),\n\011\011\011seq(lcoef[i][j]*lind[i][j],j=k
+1..nops(lcoef[i]))\n\011\011\011],`+`),vars,'distributed');\n\011\011
    lcoef[i]:=[coeffs(tosubs(i),vars,evaln(lind[i]))]\n\011\011fi\n
\011    od;\n\011fi;\n\011linds:=linds union \{mon\};\n\011u:=diff(u,z
)\n    od;\n    # final equation\n    u:=subs(h=y,collect(primpart(num
er(diff(h(z),[z$(dorder-1)])-u)),\n\011    [seq(diff(h(z),[z$j]),j=0..
dorder-1)]));\n    ## initial conditions\n    # if one of the differen
tial equations was not given with initial\n    # conditions and it is \+
singular, it may mean that a singular solution\n    # is considered. T
hen don't try to return initial conditions.\n    sing:=false;\n    for
 i to n while not sing do\n\011if ini[i]=\{\} and subs(z=0,deqs[i][nop
s(deqs[i])])=0 then\n\011    sing:=true fi od;\n    if sing then ini:=
\{\}\n    else\n\011jmax:=-1;\n\011for i to n do\n\011    D(y[i]):=D(y
[i]);\n\011    ini[i]:=`goodinitvalues/diffeq`(deqs[i],op(0,ly[i]),op(
ly[i]),\n\011\011ini[i],dorder-2+lmax);\n\011    inds:=indets(ini[i],_
C[anything]);\n\011    if inds<>\{\} then\n\011\011ini[i]:=subs([seq(i
nds[j]=_C[jmax+j],j=1..nops(inds))],ini[i]);\n\011\011jmax:=jmax+nops(
inds)\n\011    fi\n\011od;\n\011ini:=`union`(seq(subs(op(0,ly[i])=y[i]
,ini[i]),i=1..n));\n\011jmax:=max(seq(op(2,op(0,op(0,op(1,i)))),i=sele
ct(has,ini,`@@`)),\n\011    seq(1,i=select(has,ini,\{seq(D(y[i]),i=1..
n)\})),0);\n\011uu:=formatdiffeq([u,y(z)]);\n\011jmax:=max(jmax,nbinic
ond(uu,y,z));\n\011v[0]:=traperror(limit(pp,z=0,right));\n\011if v[0]=
lasterror or has(v[0],infinity) then jmax:=-1 fi;\n\011for j to jmax d
o\n\011    pp:=convert(diff(pp,z),D);\n\011    v[j]:=subs(z=0,pp);\n
\011od;\n\011ini:=`goodinitvalues/diffeq`(uu,y,z,\n\011    remove(has,
subs(ini,\{seq((D@@j)(y)(0)=v[j],j=0..jmax)\}),\n\011\011    \{seq(y[i
],i=1..n)\}))\n    fi;\n    if ini=\{\} then u else \{u,op(ini)\} fi\n
end: # poltodiffeq\n\n# diffeq+diffeq\n#Input: two differential equati
ons Eq1 and Eq2 in the variable y(z)\n#Output: a differential equation
 satisfied by the sum of a solution of Eq1 and\n#\011 a solution of Eq
2.\n`diffeq+diffeq` := proc(eq1,eq2,yofz)\nlocal y1,y2,y,z;\n    getna
me(yofz,y,z);\n    poltodiffeq(y1(z)+y2(z),[subs(y=y1,eq1),subs(y=y2,e
q2)],[y1(z),y2(z)],yofz)\nend:\n\n# diffeq*diffeq\n#Input: two differe
ntial equations Eq1 and Eq2 in the variable y(z)\n#Output: a different
ial equation satisfied by the product of a solution of Eq1\n#\011 and \+
a solution of Eq2.\n`diffeq*diffeq` := proc(eq1,eq2,yofz)\nlocal y1,y2
,y,z;\n    getname(yofz,y,z);\n    poltodiffeq(y1(z)*y2(z),[subs(y=y1,
eq1),subs(y=y2,eq2)],[y1(z),y2(z)],yofz)\nend:\n\n# cauchyproduct\n#In
put: two linear recurrences rec1 and rec2 in the variable uofn (u(n))
\n#Output: a linear recurrence satisfied by \\sum_\{k=0\}^n\{u_kv_\{n-
k\}\}, where\n#\011\011 u is a solution of rec1 and v is a solution of
 rec2.\ncauchyproduct := proc(rec1,rec2,uofn)\nlocal y1, y2, y, z, d1,
 d2, inds, j, i;\n    d1:=rectodiffeq(rec1,uofn,y1(z));\n    d2:=recto
diffeq(rec2,uofn,y2(z));\n    inds:=indets(d1,_C[anything]) intersect \+
indets(d2,_C[anything]);\n    if inds<>\{\} then\n\011j:=max(op(map(op
,\n\011    indets(d1,_C[anything]) union indets(d2,_C[anything]))));\n
\011d2:=subs([seq(inds[i]=_C[j+i],i=1..nops(inds))],d2)\n    fi;\n    \+
diffeqtorec(poltodiffeq(y1(z)*y2(z),[d1,d2],[y1(z),y2(z)],y(z)),y(z),u
ofn)\nend: # cauchyproduct\n\n# poltorec\n#Input: a polynomial P in k+
1 variable n, u1, u2, ..., uk\n#\011a list of k recurrence equations\n
#\011a list of k variables u1, u2, ..., uk\n#\011a variable u(n)\n#Out
put: a recurrence equation in u(n) satisfied by P(n,u1(n),...,uk(n)) w
here\n# ui is a solution to the ith recurrence equation.\npoltorec:=pr
oc (p, lrec::list, lu::list, uofn)\nlocal u, n, k, tosubs, i, rec, uu,
 v, lcoef, lind, j, vars, linds, w, rorder, l, l2, m, jmax, inds, ini,
 lco, mon, ord, n0, pp, n1, lmax, basis, back;\n    if nops(lrec)<>nop
s(lu) then\n\011ERROR(`not the same number of elements`,lrec,lu) fi;\n
    getname(uofn,u,n);\n    k:=nops(lrec);\n    tosubs:=subsop(4=NULL,
proc() args end):\n    for i to k do\n\011rec[i]:=formatrec([lrec[i],l
u[i]],'uu','n',ini[i]);\n\011ord[i]:=nops(rec[i])-2;\n\011lco:=rec[i][
nops(rec[i])];\n\011mon:=u[i](n+ord[i]);\n\011if ord[i]>0 then\n\011  \+
  v:=convert(map(normal,[-rec[i][1]/lco,\n\011\011seq(-rec[i][j]/lco*u
[i](n+j-2),j=2..ord[i]+1)]),`+`)\n\011else v:=-normal(rec[i][1]/lco)\n
\011fi;\n\011tosubs(mon):=v;\n\011lcoef[mon]:=[coeffs(v,[seq(u[i](n+j)
,j=0..ord[i]-1)],lind[mon])]\n    od;\n    pp:=subs([seq(op(0,lu[i])=u
[i],i=1..k)],p);\n    w:=pp; lmax:=0;\n    # a few checks on the polyn
omial\n    for i to k do\n" }}{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9237 "\011l2:=maxindex(w,u[i],n);\n
\011if l2<>-infinity-n then\n\011    lmax:=max(lmax,l2);\n\011    for \+
j from l2 by -1 to ord[i] do\n\011\011w:=subs(u[i](n+j)=subs(n=n+j-ord
[i],\n\011\011    tosubs(u[i](n+ord[i]))),w)\n\011    od\n\011fi\n    \+
od;\n    vars:=\{seq(seq(u[i](n+j),j=0..ord[i]-1),i=1..k)\};\n    if n
ot type(w,polynom(ratpoly(anything,n),vars)) then\n\011ERROR(`Not a po
lynomial`, p) fi;\n    basis:=\{seq(seq(u[i](n+j),j=0..ord[i]-1),i=1..
k)\};\n    linds:=\{1,seq(u[i](n+ord[i]),i=1..k)\};\n    for rorder do
 # this loop includes an incremental Gaussian elimination\n\011w:=coll
ect(subs([seq(u[i](n+ord[i])=tosubs(u[i](n+ord[i])),\n\011    i=1..k)]
,w),vars,'distributed');\n\011w:=[coeffs(w,vars,'l')];l:=[l];\n\011w:=
collect(convert([seq(w[j]*tosubs(l[j]),j=1..nops(w))],`+`),vars,\n\011
    'distributed',normal);\n\011lco:=[coeffs(w,vars,'l')];\n\011l2:=\{
l\} minus linds;\n\011if l2=\{\} then break fi;\n\011l:=[l];\n\011back
:=evalb(l2 minus basis=\{\});\n\011if not back then mon:=op(1,l2 minus
 basis) else mon:=l2[1] fi;\n\011member(mon,l,'m');\n\011tosubs(mon):=
collect((h(n+rorder-1)-convert([seq(lco[j]*\n\011    tosubs(l[j]),j=1.
.m-1),seq(lco[j]*tosubs(l[j]),j=m+1..nops(l))],\n\011    `+`))/lco[m],
vars,'distributed');\n\011lcoef[mon]:=[coeffs(tosubs(mon),vars,lind[mo
n])];\n\011basis:=basis union \{lind[mon]\};\n\011# back substitution
\n\011if back then\n\011    for i in linds do\n\011\011if member(mon,[
lind[i]],'m') then\n\011\011    lind[i]:=[lind[i]];\n\011\011    tosub
s(i):=collect(convert([\n\011\011\011seq(lcoef[i][j]*lind[i][j],j=1..m
-1),\n\011\011\011    lcoef[i][m]*tosubs(mon),\n\011\011\011seq(lcoef[
i][j]*lind[i][j],j=m+1..nops(lcoef[i]))\n\011\011\011],`+`),vars,'dist
ributed');\n\011\011    lcoef[i]:=[coeffs(tosubs(i),vars,evaln(lind[i]
))]\n\011\011fi\n\011    od\n\011fi;\n\011linds:=linds union \{mon\};
\n\011w:=subs(n=n+1,w)\n    od;\n    # final equation\n    w:=subs(h=u
,collect(primpart(numer(h(n+rorder-1)-w)),\n\011    [seq(h(n+j),j=0..r
order-1)]));\n    ## initial conditions\n    n0:=max(rorder-2,firstnon
zero(subs(n=n-rorder+2,\n\011coeff(w,u(n+rorder-1),1)),n));\n    jmax:
=-1;\n    for i to k do\n\011ini[i]:=`goodinitvalues/rec`(rec[i],op(0,
lu[i]),\n\011    op(lu[i]),ini[i],true,n0);\n\011inds:=indets(ini[i],_
C[anything]);\n\011if inds<>\{\} then\n\011    ini[i]:=subs([seq(inds[
j]=_C[jmax+j],j=1..nops(inds))],ini[i]);\n\011    jmax:=jmax+nops(inds
)\n\011fi\n    od;\n    n1:=max(seq(op(op(1,i)),i=seq(op(ini[i]),i=1..
k)),lmax+n0);\n    if n1>n0 then # this happens when one of the recurr
ence equations\n\011\011  # is valid only for large values of the inde
x\n\011jmax:=-1; # this has to be redone in case one of the ini was \{
\}\n\011for i to k do\n\011    ini[i]:=`goodinitvalues/rec`(rec[i],op(
0,lu[i]),\n\011\011op(lu[i]),ini[i],true,n1);\n\011    inds:=indets(in
i[i],_C[anything]);\n\011    ini[i]:=subs([seq(inds[j]=_C[jmax+j],j=1.
.nops(inds))],ini[i]);\n\011    jmax:=jmax+nops(inds)\n\011od\n    fi;
\n    ini:=`union`(seq(subs(op(0,lu[i])=u[i],ini[i]),i=1..k));\n    in
i:=`goodinitvalues/rec`(formatrec([w,u(n)],'u','n'),u,n,\n\011remove(h
as,subs(ini,\{seq(u(i)=subs(n=i,pp),\n\011    i=\{seq(op(op(1,i)),i=in
i)\})\}),\{seq(u[i],i=1..k)\}),true);\n    if ini=\{\} then w else \{w
,op(ini)\} fi\nend: # poltorec\n\n# rec+rec\n#Input: two linear recurr
ences rec1 and rec2 in the variable uofn (u(n))\n#Output: a linear rec
urrence satisfied by the sum of a solution of rec1 and\n#\011 a soluti
on of rec2.\n`rec+rec`:=subs(poltodiffeq=poltorec,eval(`diffeq+diffeq`
)):\n\n# rec*rec\n#Input: two linear recurrences rec1 and rec2 in the \+
variable uofn (u(n))\n#Output: a linear recurrence satisfied by the pr
oduct of a solution of\n#\011\011rec1 by a solution of rec2.\n`rec*rec
` := subs(poltodiffeq=poltorec,eval(`diffeq*diffeq`)):\n\n# borel, inv
borel\n# Input:  a linear recurrence or differential equation\n#\011  \+
u(n) or y(x) the variable\n#\011  (optional) a flag 'diffeq' saying th
at it's a differential equation\n#\011      by default it is a recurre
nce\n# Output: the linear recurrence or differential equation in u(n) \+
(or y(x))\n#\011  satisfied by the sequence u(n)/n! in the borel case,
 u(n)*n! in the\n#\011  invborel case. For differential equations, the
 equation is the\n#\011  equation satisfied by the generating function
 of the borel/invborel\n#\011  transform of the sequence of Taylor coe
fficients.\ninvborel:=proc() borelinvborel(false,args) end:\nLaplace:=
invborel:\nborel:=\011proc() borelinvborel(true, args) end:\n\nborelin
vborel := proc(borel,expr,aofn)\nlocal a, b, n, rec2;\n    if nargs=3 \+
then\n\011getname(aofn,a,n);\n\011if borel then rec2:=\{n*b(n)=b(n-1),
b(0)=1\}\n\011else rec2:=\{b(n)=n*b(n-1),b(0)=1\} fi;\n\011poltorec(a(
n)*b(n),[expr,rec2],[a(n),b(n)],a(n))\n    elif args[4]<>'diffeq' then
 ERROR(`invalid argument`,args[4])\n    else\n\011rectodiffeq(procname
(borel,\n\011    diffeqtorec(expr,aofn,a(n)),a(n)),a(n),aofn)\n    fi
\nend: # borelinvborel\n\n# hadamardproduct\n#Input: two linear differ
ential equations eq1 and eq2 in the variable yofz\n#Output: a linear d
ifferential equation satisfied by the Hadamard product\n# of any solut
ion of eq1 with any solution of eq2.\nhadamardproduct := proc(eq1,eq2,
yofz)\nlocal u, u1, u2, n, r1, r2, inds, i, j;\n    r1:=diffeqtorec(eq
1,yofz,u1(n));\n    r2:=diffeqtorec(eq2,yofz,u2(n));\n    inds:=indets
(r1,_C[anything]) intersect indets(r2,_C[anything]);\n    if inds<>\{
\} then\n\011j:=max(op(map(op,\n\011    indets(r1,_C[anything]) union \+
indets(r2,_C[anything]))));\n\011r2:=subs([seq(inds[i]=_C[j+i],i=1..no
ps(inds))],r2)\n    fi;\n    rectodiffeq(poltorec(u1(n)*u2(n),[r1,r2],
[u1(n),u2(n)],u(n)),u(n),yofz)\nend: # hadamardproduct\n\n# algebraics
ubs\n#Input: a linear differential equation Deq in the variable yofz (
which is y(z))\n#\011a polynomial eq in y and z\n#\011the variable y(z
)\n#\011(optional) initial conditions for y(z) in the polynomial\n#Out
put: a linear differential equation satisfied by f(y(z)) for any solut
ion\n# f of Deq and y of eq.\nalgebraicsubs := proc(Deq,eq,yofz,inipol
)\nlocal y, z, deq, P, u, d, i, d1, k, A, C, Dg, g, j, ord_eqn, c, f, \+
inhomog, eqn, inip, inid, deq0, ini, P0, v, deq1, P2, tosubs, lvars, r
educe,Y;\n    P:=formatpoleq([args[2..nargs]],'y','z',inip); d:=degree
(P,y);\n    deq0:=subs(z=y,formatdiffeq([Deq,yofz],'y','z',inid)); d1:
=nops(deq0)-2;\n    if d=0 then ERROR(`not a polynomial in the right v
ariable`,eq,y)\n    elif d=1 then # special case to speedup things\n
\011P2:=-coeff(P,y,0)/coeff(P,y,1)-z;\n\011if P2=0 then RETURN(Deq)\n
\011elif not has(P2,z) then RETURN(makediffeq(map(collect,subs(y=z+P2,
deq0),z),y,z))\n\011elif d1=0 then RETURN(makediffeq(subs(y=-coeff(P,y
,0)/coeff(P,y,1),deq0),y,z))\n\011fi\n    fi;\n    g:=mygcdex(diff(P,y
),P,y,'u');\n    if has(g,y) then RETURN(algebraicsubs(Deq,normal(P/g)
,yofz)) fi;\n    Dg:=rem(-u/g*diff(P,z),P,y);\n    g:=mygcdex(deq0[d1+
2],P,y,'u');\n    if has(g,y) then RETURN(algebraicsubs(Deq,normal(P/g
),yofz)) fi;\n    deq:=map(rem,[seq(-deq0[i]*u/g,i=1..d1+1)],P,y);\n  \+
  inhomog:=evalb(op(1,deq)<>0);\n    # if inhomog then ord_eqn:=d*d1 e
lse ord_eqn:=d*(d1+1) fi;\n    # Previous order was wrong. Fixed BS De
c 94.\n    if inhomog then ord_eqn:=d*(d1+1) else ord_eqn:=d*d1 fi;\n \+
   tosubs:=diff(f(y),[y$d1])=makediffeq(deq,f,y);\n    lvars:=[seq(dif
f(f(y),[y$k]),k=0..d1-1)];\n    eqn[0]:=f(y);\n    reduce:=subs([_P=P,
_Y=y],proc(Q) rem(Q,_P,_Y) end);\n    for k to ord_eqn do\n\011eqn[k]:
=collect(diff(eqn[k-1],z)+subs(tosubs,diff(eqn[k-1],y))*Dg,\n\011    l
vars,'distributed',reduce)\n    od;\n    tosubs:=[seq(diff(f(y),[y$k])
=Y^(k+1),k=0..d1-1)];\n    for k from 0 to ord_eqn do eqn[k]:=subs(tos
ubs,eqn[k]) od;\n    A:=array(1..ord_eqn+1,1..ord_eqn,sparse);\n    C:
=array(1..ord_eqn+1,sparse);\n    # f^\{(i)\}(g).g^j (i.e. F^i y^j) ->
 column j*d1+i+1, 0<=i<d1, 0<=j<d\n    A[1,1]:=1;\n    for k to ord_eq
n+1 do for i from 0 to d1-1 do\n\011c:=collect(coeff(eqn[k-1],Y,i+1),y
);\n\011for j from 0 to d-1 do A[k,j*d1+i+1]:=coeff(c,y,j) od od\n    \+
od;\n    if inhomog then # g^j -> column d*d1+j, 0<j<d; g^0 -> C\n\011
for k to ord_eqn+1 do\n\011    c:=coeff(eqn[k-1],Y,0);\n\011    for j \+
to d-1 do A[k,d*d1+j]:=normal(coeff(c,y,j)) od;\n\011    C[k]:=coeff(c
,y,0)\n\011od\n    fi;\n    deq:=subs([seq(v[i]=diff(y(z),[z$i]),i=0..
ord_eqn)],collect(primpart(\n\011numer(lindep(A,[seq(v[i]-C[i+1],i=0..
ord_eqn)],z)),\n\011[seq(v[i],i=0..ord_eqn)]),[seq(v[i],i=0..ord_eqn)]
));\n    # initial conditions\n    deq1:=formatdiffeq([deq,y(z)]);\n  \+
  ord_eqn:=firstnonzero(indicialeq(deq1,z,z),z);\n    if ord_eqn=0 the
n ini:=\{\} # singular case, don't look for initial conditions\n    el
se\n\011ini:=\{\};\n\011for i from 0 while has(inip,(D@@i)(y)(0)) do o
d;\n\011P0:=convert(subs(inip,[seq((D@@j)(y)(0)/j!*z^j,j=0..i-1)]),`+`
);\n\011if i<ord_eqn then\n\011    P:=subs(y=P0+z^i*y,P);\n\011    ini
p:=algeqtoseries(P,z,y,ord_eqn-i,true);\n\011    if nops(inip)=0 then
\n\011\011if i>0 then inip:=P0\n\011\011else inip:=infinity fi\n\011  \+
  elif nops(inip)>1 then inip:=P0 # do not choose\n\011    else\n\011
\011inip:=P0+z^i*inip[1];\n\011\011# this is due to a bug in V.3\n\011
\011inip:=subs(seq(i=O(z^ceil(op(2,op(i)))),i=indets(inip,\n\011\011  \+
  specfunc(identical(z)^rational,O))),inip)\n\011    fi\n\011else inip
:=P0 fi;\n\011P0:=eval(subs([O=0,z=0],inip));\n\011if P0=0 then\n\011 \+
   inid:=`goodinitvalues/diffeq`(subs(y=z,deq0),y,z,\n\011\011inid,ord
_eqn-1)\n\011elif P0=infinity then\n\011    inid:=\{\}\n\011else # for
get the initial conditions which were given\n\011    inid:=`goodinitva
lues/diffeq`(map(collect,subs(y=P0+z,deq0),z),\n\011\011y,z,\{\},ord_e
qn-1)\n\011fi;\n\011if inid<>\{\} and inip<>infinity then\n\011    ini
:=series(subs(v=inip-P0,convert(subs(inid,\n\011\011[seq((D@@j)(y)(0)*
v^j/j!,j=0..ord_eqn-1),O(v^ord_eqn)])\n\011\011    ,`+`)),z,infinity);
\n\011    if not has(ini,O(1)) then ini:=convert(ini,polynom) fi\n\011
fi;\n\011if type(ini,series) then\n\011    ini:=`goodinitvalues/diffeq
`(deq1,y,z,\n\011\011\{seq((D@@j)(y)(0)=coeff(ini,z,j)*j!,\n" }}
{PAGEBK }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 
0 18504 "\011\011j=0..min(ord_eqn-1,op(nops(ini),ini)-1))\})\n\011elif
 type(ini,polynom) then\n\011    ini:=`goodinitvalues/diffeq`(deq1,y,z
,\n\011\011\{seq((D@@j)(y)(0)=coeff(ini,z,j)*j!,j=0..ord_eqn-1)\})\n
\011fi\n    fi;\n    if ini=\{\} then deq else \{deq,op(ini)\} fi\nend
: # algebraicsubs\n\n# the following code is used by the holexprtodiff
eq function\n# added by E. Murray, Feb, 1996\n\n# a table referenced b
y the various holonomic functions giving the\n# differential equations
 they satisfy and singularity information.\n# Entries are of the form
\n# 'toto' = ( (y,z)->\n# [diffeq in y(z), \{singular points\}, \{[sin
g. point of diffeq,\n# order of first derivative to take for initial c
ondions at the point],...\}] )\n#\n# the diffeq does NOT include initi
al conditions\n# the singular points are the singular points of the fu
nction\n# the entry for BesselJ is an example of a function whose diff
eq has\n# singularity\n#\n# To access this information\n# you must cal
l diffeqtable['toto'](y,z). This eliminates problems with\n# y and z b
eing local/global variables.\n\n# for the functions which take 2 or mo
re arguments (eg. the Bessel's),\n# the assumption is that these extra
 argument comes first in the\n# function call (i.e. toto(nu,x)), and i
n the entries of this table\n# these extra values come last.\ndiffeqta
ble := table([\n'exp' = ((y,z)->[diff(y(z),z) - y(z),\{infinity\},\{\}
]),\n'ln' = ((y,z)->[diff(y(z),z)*z-1,\{0,infinity\},\{\}]),\n'sin' = \+
((y,z)->[diff(y(z),z,z) + y(z),\{infinity\},\{\}]),\n'cos' = ((y,z)->[
diff(y(z),z,z) + y(z),\{infinity\},\{\}]),\n'sinh' = ((y,z)->[y(z) - d
iff(y(z),z,z),\{infinity\},\{\}]),\n'cosh' = ((y,z)->[y(z) - diff(y(z)
,z,z),\{infinity\},\{\}]),\n'arctan' = ((y,z)->[(1+z^2)*diff(y(z),z) -
 1,\{I,-I\},\{\}]),\n'arctanh' = ((y,z)->[1+(z^2-1)*diff(y(z),z),\{1,-
1\},\{\}]),\n'arccot' = ((y,z)->[(1+z^2)*diff(y(z),z) + 1, \{I,-I\},\{
\}]),\n'arccoth' = ((y,z)->[1+(z^2-1)*diff(y(z),z),\{1,-1\},\{\}]),\n'
arcsin' = ((y,z)->[ diff(y(z),z,z)*(1-z^2) - diff(y(z),z)*z,\n\{1,-1,i
nfinity\},\{\}]),\n'arcsinh' = ((y,z)->[z*diff(y(z),z)+(z^2+1)*diff(di
ff(y(z),z),z),\n    \{I,-I,infinity\},\{\}]),\n'arccos' = ((y,z)->[z*d
iff(y(z),z)+(z^2-1)*diff(diff(y(z),z),z),\n    \{1,-1,infinity\},\{\}]
),\n'arccosh' = ((y,z)->[z*diff(y(z),z)+(z^2-1)*diff(diff(y(z),z),z),
\n    \{1,-1,infinity\},\{\}]),\n'arccsch' = ((y,z)-> [(2*z^2+1)*diff(
y(z),z)+(z^3+z)*diff(diff(y(z),z),z),\n    \{0,I,-I\},\{\}]),\n'arcsec
h' = ((y,z)-> [(2*z^2-1)*diff(y(z),z)+(z^3-z)*diff(diff(y(z),z),z),\n \+
   \{0,1,-1\},\{\}]),\n'arccsc' = ((y,z)-> [ (2*z^2-1)*diff(y(z),z)+(z
^3-z)*diff(diff(y(z),z),z),\n    \{0,1,-1\},\{\}]),\n'arcsec' = ((y,z)
-> [ (2*z^2-1)*diff(y(z),z)+(z^3-z)*diff(diff(y(z),z),z),\n    \{0,1,-
1\},\{\}]),\n'dilog' = ((y,z)->[diff(y(z),z)*z+(-z+z^2)*diff(diff(y(z)
,z),z)+1,\n    \{0,infinity\},\{1,0\}]),\n'erf' = ((y,z)->[2*z*diff(y(
z),z)+diff(y(z),z,z), \{infinity\},\{\}]),\n'erfc' = ((y,z)->[2*z*diff
(y(z),z)+diff(y(z),z,z),\{infinity\},\{\}]),\n'erfc2' = ((y,z,nu)->[di
ff(diff(y(z),z),z)+2*z*diff(y(z),z)-2*nu*y(z),\n    \{infinity\},\{\}]
),\n'erfi' = ((y,t)->[-2*t*diff(y(t),t)+diff(diff(y(t),t),t),\{infinit
y\},\{\}]),\n'hypergeom' = proc(y,z,lnum,lden) local j, n, deq, sing, \+
u;\n    deq:=rectodiffeq(u(n+1)*(n+1)*mul(n+j,j=lden)-u(n)*mul(n+j,j=l
num),\n\011u(n),y(z));\n    if type(deq,set) then deq:=op(select(has,d
eq,z)) fi;\n    if nops(lnum)>nops(lden)+1 then sing:=\{0,infinity\}\n
    elif nops(lnum)=nops(lden)+1 then sing:=\{1,infinity\}\n    else s
ing:=\{infinity\} fi;\n    [deq,sing,\{\}] end,\n'ln' = ((y,z) -> [z*d
iff(y(z),z)-1,\{0,infinity\},\{\}]),\n'AiryAi' = ((y,z)->[diff(y(z),z,
z)-z*y(z),\{infinity\},\{\}]),\n'AiryAi2' = proc(y,z,k::nonnegint) loc
al deq;\n    deq:=poltodiffeq(diff(y(z),[z$k]),[diff(y(z),z,z)-z*y(z)]
,[y(z)],y(z));\n    if type(deq,set) then deq:=op(select(has,deq,y(z))
) fi;\n    [deq,\{infinity\},\{\}] end,\n'AiryBi' = ((y,z)->[diff(y(z)
,z,z)-z*y(z),\{infinity\},\{\}]),\n'AiryBi2' = proc(y,z,k::nonnegint) \+
local deq;\n    deq:=poltodiffeq(diff(y(z),[z$k]),[diff(y(z),z,z)-z*y(
z)],[y(z)],y(z));\n    if type(deq,set) then deq:=op(select(has,deq,y(
z))) fi;\n    [deq,\{infinity\},\{\}] end,\n'AngerJ' = ((y,z,nu)->\n  \+
  [z^2*diff(y(z),z,z)+z*diff(y(z),z)+(z^2-nu^2)*y(z)=(z-nu)*sin(nu*Pi)
/Pi,\n    \{infinity\},\{\}]),\n# these should work for nu an integer \+
(both positive and negative) because\n# maple will automatically evalu
ate BesselI(-nu,x)->BesselI(nu,x) and\n# BesselJ(-nu,x) -> -BesselJ(nu
,x) or BesselJ(nu,x) as required\n# if there can be cases where this e
valuation does -not- take place,\n# then in those cases this code will
 barf.\n'BesselJ' = ((y,z,nu)->[z^2*(D@@2)(y)(z) + z*D(y)(z)+(z^2-nu^2
)*y(z),\n    \{infinity\},\{[0,nu]\}]),\n'BesselI' = ((y,z,nu)->[z^2*(
D@@2)(y)(z) + z*D(y)(z)-(z^2+nu^2)*y(z),\n    \{infinity\},\{[0,nu]\}]
),\n'BesselY' = ((y,z,nu)->[z^2*(D@@2)(y)(z) + z*D(y)(z)+(z^2-nu^2)*y(
z),\n    \{0,infinity\},\{\}]),\n'BesselK' = ((y,z,nu)->[z^2*(D@@2)(y)
(z) + z*D(y)(z)-(z^2+nu^2)*y(z),\n    \{0,infinity\},\{\}]),\n'Chi' = \+
((y,z) -> [z*diff(diff(diff(y(z),z),z),z)+2*diff(diff(y(z),z),z)\n    \+
-diff(y(z),z)*z,\{0,infinity\},\{\}]),\n'Ci' = ((y,z) -> [diff(y(z),z)
*z+2*diff(diff(y(z),z),z)+\n    z*diff(diff(diff(y(z),z),z),z),\{0,inf
inity\},\{\}]),\n'FresnelC' = ((y,z)->[Pi^2*z^3*diff(y(z),z)-diff(diff
(y(z),z),z)+z*diff(diff(diff(y(z),z),z),z),\{infinity\},\{[0 ,0]\}]),
\n'FresnelS' = ((y,z)->[Pi^2*z^3*diff(y(z),z)-diff(diff(y(z),z),z)+z*d
iff(diff(diff(y(z),z),z),z),\{infinity\},\{[0 ,0]\}]),\n'HankelH1' = (
(y,z,nu)->[z^2*(D@@2)(y)(z) + z*D(y)(z)+(z^2-nu^2)*y(z),\n    \{0,infi
nity\},\{\}]),\n'HankelH2' = ((y,z,nu)->[z^2*(D@@2)(y)(z) + z*D(y)(z)+
(z^2-nu^2)*y(z),\n    \{0,infinity\},\{\}]),\n'KelvinBer' = ((w,x,nu)-
>[x^4*diff(w(x),x$4)+2*x^3*diff(w(x),x$3)-(1+2*nu^2)*(x^2*diff(w(x),x$
2)-x*diff(w(x),x) )+(nu^4-4*nu^2+x^4)*w(x),\{infinity\},\{[0,nu]\}]),
\n'KelvinBei' = ((w,x,nu)->[x^4*diff(w(x),x$4)+2*x^3*diff(w(x),x$3)-(1
+2*nu^2)*(x^2*diff(w(x),x$2)-x*diff(w(x),x) )+(nu^4-4*nu^2+x^4)*w(x),
\{infinity\},\{[0,nu]\}]),\n'KelvinKer' = ((w,x,nu)->[x^4*diff(w(x),x$
4)+2*x^3*diff(w(x),x$3)-(1+2*nu^2)*(x^2*diff(w(x),x$2)-x*diff(w(x),x) \+
)+(nu^4-4*nu^2+x^4)*w(x),\{0,infinity\},\{\}]),\n'KelvinKei' = ((w,x,n
u)->[x^4*diff(w(x),x$4)+2*x^3*diff(w(x),x$3)-(1+2*nu^2)*(x^2*diff(w(x)
,x$2)-x*diff(w(x),x) )+(nu^4-4*nu^2+x^4)*w(x),\{0,infinity\},\{\}]),\n
'KelvinHer' = ((y,z,nu)->[(-4*nu^2+nu^4+z^4)*y(z)+(-2*z^2*nu^2-z^2)*di
ff(diff(y(z),z),z)+(2*z*nu^2+z)*diff(y(z) ,z)+diff(diff(diff(diff(y(z)
,z),z),z),z)*z^4+2*diff(diff(diff(y(z),z),z),z)*z^3,\{0,infinity\},\{
\}]) ,\n'KelvinHei' = ((y,z,nu)->[(-4*nu^2+nu^4+z^4)*y(z)+(-2*z^2*nu^2
-z^2)*diff(diff(y(z),z),z)+(2*z*nu^2+z)*diff(y(z) ,z)+diff(diff(diff(d
iff(y(z),z),z),z),z)*z^4+2*diff(diff(diff(y(z),z),z),z)*z^3,\{0,infini
ty\},\{\}]) ,\n'KummerM' = ((y,z,mu,nu)->[z*diff(y(z),z,z)+(nu-z)*diff
(y(z),z)-mu*y(z),\{infinity\},\{[0,0]\}]),\n'KummerU' = ((y,z,mu,nu)->
[z*diff(y(z),z,z)+(nu-z)*diff(y(z),z)-mu*y(z),\{0,infinity\},\{\}]),\n
'LommelS1'=proc(y,z,mu::integer,nu)\n    [z^2*diff(y(z),z,z)+z*diff(y(
z),z)+(z^2-nu^2)*y(z)-z^(mu+1),\{infinity\},\{[0,0]\}] end,\n'LommelS2
'=proc(y,z,mu::integer,nu)\n    [z^2*diff(y(z),z,z)+z*diff(y(z),z)+(z^
2-nu^2)*y(z)-z^(mu+1),\{0,infinity\},\{\}] end,\n'Shi' = ((y,z) -> [z*
diff(diff(diff(y(z),z),z),z)+2*diff(diff(y(z),z),z)\n    -diff(y(z),z)
*z,\{infinity\},\{[0,0]\}]),\n'Si' = ((y,z) -> [diff(y(z),z)*z+2*diff(
diff(y(z),z),z)+\n    z*diff(diff(diff(y(z),z),z),z),\{infinity\},\{[0
,0]\}]),\n'Ssi' = ((y,z) -> [diff(y(z),z)*z+2*diff(diff(y(z),z),z)+\n \+
   z*diff(diff(diff(y(z),z),z),z),\{infinity\},\{[0,0]\}]),\n'StruveH'
=((y,z,nu::integer)->[z^2*diff(y(z),z,z)+z*diff(y(z),z)+(z^2-nu^2)*y(z
)=4*(z/2)^(nu+1)/Pi^(1/2)/GAMMA(nu+1/2),\{0,infinity\},\{\}]),\n'Struv
eL'=((y,z,nu::integer)->[z^2*diff(y(z),z,z)+z*diff(y(z),z)-(z^2+nu^2)*
y(z)=4*(z/2)^(nu+1)/Pi^(1/2)/GAMMA(nu+1/2),\{0,infinity\},\{\}]),\n'We
berE' = ((y,z,nu)->\n    [z^2*diff(y(z),z,z)+z*diff(y(z),z)+(z^2-nu^2)
*y(z)=\n\011((nu-z)*cos(nu*Pi)-(nu+z))/Pi,\{infinity\},\{\}]),\n'Whitt
akerM'=((y,z,mu,nu)->[diff(y(z),z,z)+(-1/4+mu/z+(1/4-nu^2)/z^2)*y(z),
\n    \{0,infinity\},\{\}]),\n'WhittakerW'=((y,z,mu,nu)->[diff(y(z),z,
z)+(-1/4+mu/z+(1/4-nu^2)/z^2)*y(z),\n    \{0,infinity\},\{\}])\n]):\n
\n`type/gfun/has2diffeqs` := proc(x)\n    member(x,\{'AiryAi','AiryBi'
,'erfc'\})\nend:\n\n# given a holonomic expression, return the differe
ntial equation it\n# satifies, giving initial conditions when we can\n
holexprtodiffeq := proc(expr,yofx)\n# global diffeqtable;\nlocal y,x,d
eqs, newexpr, funs,nf,i,j,sortbylen,eq,inits, pows, np, ini, deq,exps;
\n    getname(yofx, y,x);\n    if has(expr, y) then\n\011ERROR(`first \+
argument should not contain`, y, expr)\n    elif not isholonomic(expr,
x) then\n\011ERROR(`expression is not holonomic (or involves a functio
n which is not implemented yet)`, expr)\n    fi;\n    if not has(expr,
 x) then\n\011RETURN(yofx - expr)\n    fi;\n\n    newexpr:=expr;\n    \+
# special case for exp(anything)^rational\n    exps:=indets(newexpr,sp
ecfunc(anything,exp)^rational);\n    if exps<>\{\} then\n\011newexpr:=
subs([seq(i=exp(op(2,i)*op(op(1,i))),i=exps)],newexpr)\n    fi;\n\n   \+
 # deal with holonomic functions f(expr)\n    funs := [op(remove(type,
indets(newexpr,function), RootOf))];\n    nf := nops(funs);\n    deqs \+
:= map(funtodiffeq, funs, y(x));\n    deqs := table([seq(subs(y=y[i],d
eqs[i]),i=1..nf)]);\n    newexpr := subs([seq(funs[i]=y[i](x),i=1..nf)
], newexpr);\n\n    # special case for polynomial(x)^(freeof(x))\n    \+
pows := remove(type,indets(newexpr,polynom(anything,x)^`gfun/free`(x))
,\n\011anything^rational);\n    np:=nops(pows);\n    for i to np do\n
\011ini:=subs(x=0,op(1,pows[i]));\n\011deq:=op(1,pows[i])*diff(y[i+nf]
(x),x)\n\011    -op(2,pows[i])*diff(op(1,pows[i]),x)*y[i+nf](x);\n\011
if ini=0 then deqs[i+nf]:=deq\n\011else deqs[i+nf]:=\{deq,y[i+nf](0)=i
ni^op(2,pows[i])\}\n\011fi\n    od;\n    newexpr:=subs([seq(pows[i]=y[
i+nf](x),i=1..np)],newexpr);\n\n    # pick out the non-polynomial alge
braic expressions,\n    # the bits that are polynomial in x will be ha
ndled by poltodiffeq\n    funs  := remove(type,\n      indets(newexpr,
 \{radfun(anything,x),algfun(anything,x)\}),\n      polynom(anything,x
));\n    sortbylen := (a,b) -> evalb(length(a)>=length(b));\n    for i
 while funs <> \{\} do\n\011# find the longest one - guaranteed to be \+
at the outermost level\n\011funs := sort([op(funs)],sortbylen);\n\011e
q := algfuntoalgeq(funs[1],y(x),'inits','algebraic');\n\011deqs[i+nf+n
p] := subs(y=y[i+nf+np],algeqtodiffeq(eq,y(x),inits));\n\011newexpr :=
 subs(funs[1]=y[i+nf+np](x),newexpr);\n\011funs  := remove(type,\n\011
    indets(newexpr, \{radfun(anything,x),algfun(anything,x)\}),\n\011 \+
   polynom(anything,x));\n    od;\n\n    poltodiffeq(newexpr, [seq(deq
s[j],j=1..nf+np+i-1)],\n\011[seq(y[j](x),j=1..nf+np+i-1)],y(x));\nend:
 # holexprtodiffeq\n\n# determine if an expression is holonomic\nishol
onomic := proc(expr,x)\n# global diffeqtable;\nlocal i;\n    if not ha
s(expr, x) then true\n    elif type(expr, function) then\n\011if not t
ype(expr, RootOf) then\n\011    # holonomic function with argument alg
ebraic in x\n\011    evalb(assigned(diffeqtable[op(0,expr)]) and\n\011
    type(op(nops(expr),expr),\{radfun(anything,x),algfun(anything,x)\}
))\n\011else\n\011    # this will reject things like RootOf(arcos(_Z))
\n\011    type(op(1,expr),\{algfun(anything,_Z),radfun(anything,_Z)\})
;\n\011fi;\n    elif type(expr, \{radfun(anything,x),algfun(anything,x
)\}) then\n\011true\n    elif type(expr, \{`+`,`*`\}) then\n\011for i \+
to nops(expr) while procname(op(i,expr),x) do od;\n\011evalb(i=nops(ex
pr)+1);\n    elif type(expr, `^`) then\n\011# polynomial^nice\n\011typ
e(op(1,expr),polynom(anything,x)) and\n\011    not has(op(2,expr),x) o
r\n\011# holonomic^posint (i.e. multiplication)\n\011type(op(2,expr), \+
posint) and procname(op(1,expr),x) or\n\011# exp(algebraic)^rational\n
\011type(op(2,expr),rational) and\n\011  type(op(1,expr),specfunc(\{ra
dfun(anything,x),algfun(anything,x)\},exp))\n    else\n\011false\n    \+
fi\nend: # isholonomic\n\n# Given a holonomic expression of the form f
un(arg(s)), returns the\n# differential equation it satisfies.\n# Remo
ved the following restriction (BS Aug.96).\n## Returns an error if it \+
encounters a singularity.\n\nfuntodiffeq := proc(expr, yofx)\nlocal x,
y,fun,deq,lim,deq2,ord,initconds,i,initpts,cond,tabref,nu,info,eq,\n  \+
pinit, limpt, newx, sing, issing;\n# global diffeqtable;\n    getname(
yofx, y,x);\n    fun := op(0,expr);\n    if not has(expr, x) then RETU
RN(yofx - expr) fi;\n    tabref := fun;\n    # things like erfc have 2
 diff eqs, depending on the\n    # number of arguments.  The second on
e is called erfc2.\n    if type(tabref, `gfun/has2diffeqs`) and nops(e
xpr)=2 then\n\011tabref := cat(tabref,2)\n    fi;\n    # worry about h
aving 2 or 3 arguments\n    nu:=op(1..nops(expr)-1,expr);\n    info :=
 diffeqtable[tabref](y,x,nu);\n    deq := info[1];\n\n    lim := limit
(op(nops(expr),expr),x=0, right);\n    if has(lim,infinity) then lim :
= infinity fi;\n    issing:=member(lim,info[2]);\n    if lim = 0 then \+
deq2:=deq\n    elif lim <> infinity then\n\011# find diffeq for fun(li
m+x)\n\011deq2 := algebraicsubs(deq, y-(lim+x),y(x));\n    else\n\011#
 find diffeq for fun(1/x)\n\011deq2 := algebraicsubs(deq, y*x-1, y(x))
;\n    fi;\n\n    # how many initial conditions do we need?\n    if is
sing or nops(expr)=2 and not type(nu, integer) then\n\011ord := 0  # w
e can't give any\n    else\n\011ord := nops(formatdiffeq([deq2,yofx]))
-2;\n    fi;\n\n    # where do we take them? - is the diff eq singular
 at this point?\n    sing := select((a,l)->evalb(op(1,a)=l), info[3],l
im);\n    if sing = \{\} then\n\011initpts := \{seq(i,i=0..ord-1)\}\n \+
   else\n\011initpts := \{seq(i+op(2,op(sing)),i=0..ord-1)\};\n    fi;
\n\n    if lim = infinity or nops(expr)>1 or sing<>\{\} then\n\011# ca
n't use the D notation - take limits of diff\n\011if lim = infinity th
en\n\011    limpt := 0;\n\011    newx := 1/x\n\011else\n\011    limpt \+
:= lim;\n\011    newx := x;\n\011fi;\n\n\011if member(0,initpts) then
\n\011    cond := \{y(0) = limit(fun(nu,newx),x=limpt)\};\n\011    ini
tpts := initpts minus \{0\}\n\011else\n\011    cond := \{\};\n\011fi;
\n\011initconds := cond union\n\011   \{seq((D@@i)(y)(0)=limit(diff(fu
n(nu,newx),x$i),x=limpt),i=initpts)\};\n    else\n\011# use the D nota
tion\n\011initconds := \{seq((D@@i)(y)(0)=(D@@i)(fun)(lim),i=initpts)
\};\n    fi;\n\n    if not type(deq2,set) then\n\011deq2 := \{deq2\};
\n    fi;\n    deq2 := deq2 union simplify(initconds,infinity);\n\n   \+
 if lim <> infinity then\n\011# find eq for op(expr), diff eq for fun(
lim + (op(expr)-lim))\n\011eq := algfuntoalgeq(op(nops(expr),expr),yof
x,'pinit','algebraic');\n\011# and modify the initial conditions of eq
\n\011# y(0)=lim  in op(expr) becomes y(0)=0 in op(expr)-lim\n\011alge
braicsubs(deq2,\n\011    subs(y=y+lim,eq),y(x),\{y(0)=0\} union remove
(has,pinit,y(0)));\n    else\n\011# find eq for 1/ op(expr), diff eq f
or fun( 1/ 1/(op(expr)) )\n\011eq := algfuntoalgeq(1/op(nops(expr),exp
r),yofx,'algebraic');\n\011# y(0)=infinity in op(expr) becomes y(0)=0 \+
in 1/op(expr)\n\011algebraicsubs(deq2,eq,y(x),\{y(0)=0\});\n    fi;\ne
nd: # funtodiffeq\n\n# given an algebraic expression expr in the varia
ble x, return a polynomial\n# (possibly with algebraic number coeffici
ents) in\n# x and y which has that expression as a root\n# iniconds (o
ptional) - will be assigned the initial condition