Taylor expansion wrt. several variables; taylor, taytorat bugs?

The linearization of the function f(x,y) at a point [a1, a2]
should be something like this, I think:

mytaylorn(f(x,y), [x,y], [a1,a2], 1);
''%;
(%o11)
(at('diff(f(x,y),y,1),[x=a1,y=a2]))*(y-a2)+(x-a1)*(at('diff(f(x,y),x,1),[x=a1,y=a2]))+at(f(x,y),[x=a1,y=a2])

taylor() doesn't seem to specify the point [x,y] sufficiently
(only either x or y are set):

ft1: taylor(f(x,y),[x,y], [a1,a2], [1,1]);
(%o12)
f(a1,a2)+((at('diff(f(x,y),x,1),x=a1))*(x-a1)+(at('diff(f(a1,y),y,1),y=a2))*(y-a2))+...

If I understand correctly, taytorat() can be used to get rid
of the ellipsis, however it produces an unknown symbol
(probably gensym()ed):

taytorat(ft1);
(%o13)
((at('diff(f(x,y),x,1),x=a1))*(x-a1)+(at('diff(f(a1,y),y,1),y=a2))*(y-a2))*g34365+f(a1,a2)

The documentation says taytorat() does the same as
rat(ratdisrep()), but it yields a different result:

rat(ratdisrep(ft1));
(%o14)
(at('diff(f(a1,y),y,1),y=a2))*y+(x-a1)*(at('diff(f(x,y),x,1),x=a1))-(at('diff(f(a1,y),y,1),y=a2))*a2+f(a1,a2)

Best regards,
Martin



some weird code to get the above result:


wxMaxima 0.8.3a http://wxmaxima.sourceforge.net
Maxima 5.19.1 http://maxima.sourceforge.net
Using Lisp GNU Common Lisp (GCL) GCL 2.6.8 (aka GCL)
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.



(%i1) kill(all);
(%o0) done



  aux. functions

(%i1) nestdiff(f,xv,nv) :=
if length(xv)=1
then diff(f, first(xv), first(nv))
else nestdiff(
     diff(f, first(xv), first(nv)),
     rest(xv),
     rest(nv)
)$



(%i2) nestdiff(f(x,y), [x,y], [n1,n2]);
(%o2) 'diff(f(x,y),x,n1,y,n2)



(%i3) nestsum(e, vv, n) := block([vvv],
if length(vv)=1
then buildq(
     [e, vvv:[first(vv)], n],
     sum(e,splice(vvv), 0, n)
)
else nestsum(
     buildq(
         [e, vvv:[first(vv)], n],
         sum(e,splice(vvv), 0, n)),
     rest(vv),
     n
))$



(%i4) nestsum(a, [gg,hh], n2);
(%o4) sum(sum(a,gg,0,n2),hh,0,n2)



  symbol version

(%i5) mytaylor(f, x, a, n) := block(
     [
         i,
         len: length(x),
         expr,
         varlist
     ],
     varlist: makelist(?gensym(),i,1,len),
     expr: product(
         ((part(x,i)-part(a,i))^part(varlist,i))
         /
         (part(varlist,i)!),
         i,1,len
     )
     *
     (at(nestdiff(f, x, varlist), map("=", x, a))),
     expr: nestsum(expr, varlist, n)
)$


(How could I get the buildq()ed expression evaluated _inside_
of mytaylor, instead of calling ev(%) or ''% after the
mytaylor() call?)

(%i6) mytaylor(f(x,y), [x,y], [a1,a2], n2);
''%;
(%o6)
sum(sum(((x-a1)^g34175*(at('diff(f(x,y),x,g34175,y,g34176),[x=a1,y=a2]))*(y-a2)^g34176)/(g34175!*g34176!),g34175,0,n2),g34176,0,n2)
(%o7)
sum(((y-a2)^g34176*sum(((x-a1)^g34175*(at('diff(f(x,y),x,g34175,y,g34176),[x=a1,y=a2])))/g34175!,g34175,0,n2))/g34176!,g34176,0,n2)



  number version

(%i8) nestdiffn(f, xv, nv, n) :=
(if apply("+", nv)>n then 0 else 1)*(if length(xv)=1
then diff(f, first(xv), first(nv))
else nestdiff(
     diff(f, first(xv), first(nv)),
     rest(xv),
     rest(nv)
))$



(%i9) mytaylorn(f, x, a, n) := block(
     [
         i,
         len: length(x),
         expr,
         varlist
     ],
     varlist: makelist(?gensym(),i,1,len),
     expr: product(
         ((part(x,i)-part(a,i))^part(varlist,i))
         /
         (part(varlist,i)!),
         i,1,len
     )
     *
     (at(nestdiffn(f, x, varlist, n), map("=", x, a))),
     expr: nestsum(expr, varlist, n)
)$



(%i10) mytaylorn(f(x,y), [x,y], [a1,a2], 1);
''%;
(%o10)
sum(sum(((x-a1)^g34358*(at('diff(f(x,y),x,g34358,y,g34359),[x=a1,y=a2]))*(y-a2)^g34359*(if 


g34359+g34358>1 then 0 else
1))/(g34358!*g34359!),g34358,0,1),g34359,0,1)
(%o11)
(at('diff(f(x,y),y,1),[x=a1,y=a2]))*(y-a2)+(x-a1)*(at('diff(f(x,y),x,1),[x=a1,y=a2]))+at(f(x,y),[x=a1,y=a2])


  What we get from taylor() is


(%i12) ft1: taylor(f(x,y),[x,y], [a1,a2], [1,1]);
taytorat(ft1);
rat(ratdisrep(ft1));
(%o12)
f(a1,a2)+((at('diff(f(x,y),x,1),x=a1))*(x-a1)+(at('diff(f(a1,y),y,1),y=a2))*(y-a2))+...
(%o13)
((at('diff(f(x,y),x,1),x=a1))*(x-a1)+(at('diff(f(a1,y),y,1),y=a2))*(y-a2))*g34365+f(a1,a2)
(%o14)
(at('diff(f(a1,y),y,1),y=a2))*y+(x-a1)*(at('diff(f(x,y),x,1),x=a1))-(at('diff(f(a1,y),y,1),y=a2))*a2+f(a1,a2)