reduce-historical/r38/packages/taylor/taylor.rlg

Tue Feb 10 12:29:11 2004 run on Linux
comment
        Test and demonstration file for the Taylor expansion package,
        by Rainer M. Schoepf.  Works with version 2.2a (01-Apr-2000);


%%% showtime;

on errcont;

 % disable interruption on errors

comment Simple Taylor expansion;


xx := taylor (e**x, x, 0, 4);


               1   2    1   3    1    4      5
xx := 1 + x + ---*x  + ---*x  + ----*x  + O(x )
               2        6        24

yy := taylor (e**y, y, 0, 4);


               1   2    1   3    1    4      5
yy := 1 + y + ---*y  + ---*y  + ----*y  + O(y )
               2        6        24


comment Basic operations, i.e. addition, subtraction, multiplication,
        and division are possible: this is not done automatically if
        the switch TAYLORAUTOCOMBINE is OFF.  In this case it is
        necessary to use taylorcombine;


taylorcombine (xx**2);


             2    4   3    2   4      5
1 + 2*x + 2*x  + ---*x  + ---*x  + O(x )
                  3        3

taylorcombine (ws - xx);


     3   2    7   3    5   4      5
x + ---*x  + ---*x  + ---*x  + O(x )
     2        6        8

taylorcombine (xx**3);


           9   2    9   3    27   4      5
1 + 3*x + ---*x  + ---*x  + ----*x  + O(x )
           2        2        8


comment The result is again a Taylor kernel;


if taylorseriesp ws then write "OK";


OK


comment It is not possible to combine Taylor kernels that were
        expanded with respect to different variables;


taylorcombine (xx**yy);


          1   2    1   3    1    4      5
(1 + x + ---*x  + ---*x  + ----*x  + O(x ))
          2        6        24

            1   2    1   3    1    4      5
**(1 + y + ---*y  + ---*y  + ----*y  + O(y ))
            2        6        24


comment But we can take the exponential or the logarithm
        of a Taylor kernel;


taylorcombine (e**xx);


             2    5*e   3    5*e   4      5
e + e*x + e*x  + -----*x  + -----*x  + O(x )
                   6          8

taylorcombine log ws;


         1   2    1   3    1    4      5
1 + x + ---*x  + ---*x  + ----*x  + O(x )
         2        6        24


comment A more complicated example;


hugo := taylor(log(1/(1-x)),x,0,5);


             1   2    1   3    1   4    1   5      6
hugo := x + ---*x  + ---*x  + ---*x  + ---*x  + O(x )
             2        3        4        5

taylorcombine(exp(hugo/(1+hugo)));


         1    4      6
1 + x + ----*x  + O(x )
         12


comment We may try to expand about another point;


taylor (xx, x, 1, 2);


 65     8             5         2            3
---- + ---*(x - 1) + ---*(x - 1)  + O((x - 1) )
 24     3             4


comment Arc tangent is one of the functions this package knows of;


xxa := taylorcombine atan ws;


             65      1536             2933040          2            3
xxa := atan(----) + ------*(x - 1) - ----------*(x - 1)  + O((x - 1) )
             24      4801             23049601


comment The trigonometric functions;


taylor (tan x / x, x, 0, 2);


     1   2      3
1 + ---*x  + O(x )
     3


taylorcombine sin ws;


          cos(1)   2      3
sin(1) + --------*x  + O(x )
            3


taylor (cot x / x, x, 0, 4);


 -2    1     1    2     2    4      5
x   - --- - ----*x  - -----*x  + O(x )
       3     45        945


comment The poles of these functions are correctly handled;


taylor(tan x,x,pi/2,0);


         pi  -1          pi
 - (x - ----)   + O(x - ----)
         2               2


taylor(tan x,x,pi/2,3);


         pi  -1    1        pi      1         pi  3           pi  4
 - (x - ----)   + ---*(x - ----) + ----*(x - ----)  + O((x - ----) )
         2         3        2       45        2               2


comment Expansion with respect to more than one kernel is possible;


xy := taylor (e**(x+y), x, 0, 2, y, 0, 2);


               1   2                            3  3
xy := 1 + y + ---*y  + x + y*x + (4 terms) + O(x ,y )
               2


taylorcombine (ws**2);


             2                                3  3
1 + 2*y + 2*y  + 2*x + 4*y*x + (4 terms) + O(x ,y )


comment We take the inverse and convert back to REDUCE's standard
        representation;


taylorcombine (1/ws);


             2                                3  3
1 - 2*x + 2*x  - 2*y + 4*y*x + (4 terms) + O(x ,y )

taylortostandard ws;


   2  2      2        2        2                    2
4*x *y  - 4*x *y + 2*x  - 4*x*y  + 4*x*y - 2*x + 2*y  - 2*y + 1


comment Some examples of Taylor kernel divsion;


xx1 := taylor (sin (x), x, 0, 4);


            1   3      5
xx1 := x - ---*x  + O(x )
            6

taylorcombine (xx/xx1);


 -1        2       1   2      3
x   + 1 + ---*x + ---*x  + O(x )
           3       3

taylorcombine (1/xx1);


 -1    1         3
x   + ---*x + O(x )
       6


tt1 := taylor (exp (x), x, 0, 3);


                1   2    1   3      4
tt1 := 1 + x + ---*x  + ---*x  + O(x )
                2        6

tt2 := taylor (sin (x), x, 0, 3);


            1   3      4
tt2 := x - ---*x  + O(x )
            6

tt3 := taylor (1 + tt2, x, 0, 3);


                1   3      4
tt3 := 1 + x - ---*x  + O(x )
                6

taylorcombine(tt1/tt2);


 -1        2         2
x   + 1 + ---*x + O(x )
           3

taylorcombine(tt1/tt3);


     1   2    1   3      4
1 + ---*x  - ---*x  + O(x )
     2        6

taylorcombine(tt2/tt1);


     2    1   3      4
x - x  + ---*x  + O(x )
          3

taylorcombine(tt3/tt1);


     1   2    1   3      4
1 - ---*x  + ---*x  + O(x )
     2        6



comment Here's what I call homogeneous expansion;


xx := taylor (e**(x*y), {x,y}, 0, 2);


                       3
xx := 1 + y*x + O({x,y} )

xx1 := taylor (sin (x+y), {x,y}, 0, 2);


                      3
xx1 := y + x + O({x,y} )

xx2 := taylor (cos (x+y), {x,y}, 0, 2);


            1   2          1   2          3
xx2 := 1 - ---*y  - y*x - ---*x  + O({x,y} )
            2              2

temp := taylorcombine (xx/xx2);


             1   2            1   2          3
temp := 1 + ---*y  + 2*y*x + ---*x  + O({x,y} )
             2                2

taylorcombine (ws*xx2);


                 3
1 + y*x + O({x,y} )


comment The following shows a principal difficulty:
        since xx1 is symmetric in x and y but has no constant term
        it is impossible to compute 1/xx1;


taylorcombine (1/xx1);


***** Not a unit in argument to invtaylor 


comment Substitution in Taylor expressions is possible;


sub (x=z, xy);


         1   2                            3  3
1 + y + ---*y  + z + y*z + (4 terms) + O(z ,y )
         2


comment Expression dependency in substitution is detected;


sub (x=y, xy);


***** Invalid substitution in Taylor kernel: dependent variables y y 


comment It is possible to replace a Taylor variable by a constant;


sub (x=4, xy);


             13   2      3
13 + 13*y + ----*y  + O(y )
             2


sub (x=4, xx1);


           3
4 + y + O(y )


sub (y=0, ws);


4


comment This package has three switches:
        TAYLORKEEPORIGINAL, TAYLORAUTOEXPAND, and TAYLORAUTOCOMBINE;


on taylorkeeporiginal;



temp := taylor (e**(x+y), x, 0, 5);


                      y         y         y
         y    y      e    2    e    3    e    4                 6
temp := e  + e *x + ----*x  + ----*x  + ----*x  + (1 term) + O(x )
                     2         6         24


taylorcombine (log (temp));


           6
y + x + O(x )


taylororiginal ws;


x + y


taylorcombine (temp * e**x);


                  y         y         y
 x   y    y      e    2    e    3    e    4                 6
e *(e  + e *x + ----*x  + ----*x  + ----*x  + (1 term) + O(x ))
                 2         6         24


on taylorautoexpand;



taylorcombine ws;


                            y           y
 y      y        y  2    4*e    3    2*e    4                 6
e  + 2*e *x + 2*e *x  + ------*x  + ------*x  + (1 term) + O(x )
                          3           3


taylororiginal ws;


 2*x + y
e


taylorcombine (xx1 / x);


   -1              2
y*x   + 1 + O({x,y} )


on taylorautocombine;



xx / xx2;


     1   2            1   2          3
1 + ---*y  + 2*y*x + ---*x  + O({x,y} )
     2                2


ws * xx2;


                 3
1 + y*x + O({x,y} )


comment Another example that shows truncation if Taylor kernels
        of different expansion order are combined;


comment First we increase the number of terms to be printed;


taylorprintterms := all;


taylorprintterms := all


p := taylor (x**2 + 2, x, 0, 10);


          2      11
p := 2 + x  + O(x  )

p - x**2;


      2      11      2
(2 + x  + O(x  )) - x

p - taylor (x**2, x, 0, 5);


       6
2 + O(x )

taylor (p - x**2, x, 0, 6);


       7
2 + O(x )

off taylorautocombine;


taylorcombine(p-x**2);


       11
2 + O(x  )

taylorcombine(p - taylor(x**2,x,0,5));


       6
2 + O(x )


comment Switch back to finite number of terms;


taylorprintterms := 6;


taylorprintterms := 6


comment Some more examples;


taylor(1/(1+y^4+x^2*y^2+x^4),{x,y},0,6);


     4    2  2    4          7
1 - y  - y *x  - x  + O({x,y} )


taylor ((1 + x)**n, x, 0, 3);


                                2
           n*(n - 1)   2    n*(n  - 3*n + 2)   3      4
1 + n*x + -----------*x  + ------------------*x  + O(x )
               2                   6


taylor (e**(-a*t) * (1 + sin(t)), t, 0, 4);


                                         3      2
                    a*(a - 2)   2     - a  + 3*a  - 1   3
1 + ( - a + 1)*t + -----------*t  + ------------------*t
                        2                   6

        3      2
    a*(a  - 4*a  + 4)   4      5
 + -------------------*t  + O(t )
           24


operator f;



taylor (1 + f(t), t, 0, 3);


                                    sub(t=0,df(f(t),t,2))   2
f(0) + 1 + sub(t=0,df(f(t),t))*t + -----------------------*t
                                              2

    sub(t=0,df(f(t),t,3))   3      4
 + -----------------------*t  + O(t )
              6


taylor(f(sqrt(x^2+y^2)),x,x0,4,y,y0,4);


         2     2                          2    2
f(sqrt(x0  + y0 )) + sub(y=y0,df(f(sqrt(x0  + y )),y))*(y - y0)

                         2    2
    sub(y=y0,df(f(sqrt(x0  + y )),y,2))          2
 + -------------------------------------*(y - y0)
                     2

                         2    2
    sub(y=y0,df(f(sqrt(x0  + y )),y,3))          3
 + -------------------------------------*(y - y0)
                     6

                         2    2
    sub(y=y0,df(f(sqrt(x0  + y )),y,4))          4
 + -------------------------------------*(y - y0)
                    24

                       2     2
 + sub(x=x0,df(f(sqrt(x  + y0 )),x))*(x - x0) + (19 terms)

             5         5
 + O((x - x0) ,(y - y0) )


clear f;



taylor (sqrt(1 + a*x + sin(x)), x, 0, 3);


                     2                     3      2
     a + 1        - a  - 2*a - 1   2    3*a  + 9*a  + 9*a - 1   3      4
1 + -------*x + -----------------*x  + -----------------------*x  + O(x )
       2                8                        48


taylorcombine (ws**2);


                 1   3      4
1 + (a + 1)*x - ---*x  + O(x )
                 6


taylor (sqrt(1 + x), x, 0, 5);


     1       1   2    1    3     5    4     7    5      6
1 + ---*x - ---*x  + ----*x  - -----*x  + -----*x  + O(x )
     2       8        16        128        256


taylor ((cos(x) - sec(x))^3, x, 0, 5);


       6
0 + O(x )


taylor ((cos(x) - sec(x))^-3, x, 0, 5);


    -6    1   -4    11    -2     347       6767    2     15377    4      6
 - x   + ---*x   + -----*x   - ------- - --------*x  - ---------*x  + O(x )
          2         120         15120     604800        7983360


taylor (sqrt(1 - k^2*sin(x)^2), x, 0, 6);


      2         2        2              2         4       2
     k    2    k *( - 3*k  + 4)   4    k *( - 45*k  + 60*k  - 16)   6      7
1 - ----*x  + ------------------*x  + ----------------------------*x  + O(x )
     2                24                          720


taylor (sin(x + y), x, 0, 3, y, 0, 3);


     1   3        1     2    1   2      1    2  3                  4  4
x - ---*x  + y - ---*y*x  - ---*y *x + ----*y *x  + (2 terms) + O(x ,y )
     6            2          2          12


taylor (e^x - 1 - x,x,0,6);


 1   2    1   3    1    4     1    5     1    6      7
---*x  + ---*x  + ----*x  + -----*x  + -----*x  + O(x )
 2        6        24        120        720


taylorcombine sqrt ws;


    1              1       2        1        3         1        4
---------*x + -----------*x  + ------------*x  + -------------*x
 sqrt(2)       6*sqrt(2)        36*sqrt(2)        270*sqrt(2)

         1         5      6
 + --------------*x  + O(x )
    2592*sqrt(2)


taylor(sin(x)/x,x,1,2);


                                       - 2*cos(1) + sin(1)         2
sin(1) + (cos(1) - sin(1))*(x - 1) + ----------------------*(x - 1)
                                               2

            3
 + O((x - 1) )


taylor((sqrt(4+h)-2)/h,h,0,5);


 1     1         1    2      5     3      7      4      21      5      6
--- - ----*h + -----*h  - -------*h  + --------*h  - ---------*h  + O(h )
 4     64       512        16384        131072        2097152


taylor((sqrt(x)-2)/(4-x),x,4,2);


    1     1               1          2            3
 - --- + ----*(x - 4) - -----*(x - 4)  + O((x - 4) )
    4     64             512


taylor((sqrt(y+4)-2)/(-y),y,0,2);


    1     1         1    2      3
 - --- + ----*y - -----*y  + O(y )
    4     64       512


taylor(x*tanh(x)/(sqrt(1-x^2)-1),x,0,3);


        7   2      4
 - 2 + ---*x  + O(x )
        6


taylor((e^(5*x)-2*x)^(1/x),x,0,2);


                   3
 3      3      73*e    2      3
e  + 8*e *x + -------*x  + O(x )
                 3


taylor(sin x/cos x,x,pi/2,3);


         pi  -1    1        pi      1         pi  3           pi  4
 - (x - ----)   + ---*(x - ----) + ----*(x - ----)  + O((x - ----) )
         2         3        2       45        2               2


taylor(log x*sin(x^2)/(x*sinh x),x,0,2);


             1   2      3
log(x)*(1 - ---*x  + O(x ))
             6


taylor(1/x-1/sin x,x,0,2);


    1         3
 - ---*x + O(x )
    6


taylor(tan x/log cos x,x,pi/2,2);


          pi  -1          pi
  - (x - ----)   + O(x - ----)
          2               2
-------------------------------
    log(pi - 2*x) - log(2)


taylor(log(x^2/(x^2-a)),x,0,3);


                2
             - x
taylor(log(--------),x,0,3)
                 2
            a - x



comment Three more complicated examples contributed by Stan Kameny;


zz2 := (z*(z-2*pi*i)*(z-pi*i/2)^2)/(sinh z-i);


                 3            2       2        3
        z*(2*i*pi  - 12*i*pi*z  - 9*pi *z + 4*z )
zz2 := -------------------------------------------
                     4*(sinh(z) - i)

dz2 := df(zz2,z);


                         3                      3               2  2
dz2 := ( - 2*cosh(z)*i*pi *z + 12*cosh(z)*i*pi*z  + 9*cosh(z)*pi *z

                      4                 3                    2
         - 4*cosh(z)*z  + 2*sinh(z)*i*pi  - 36*sinh(z)*i*pi*z

                        2                 3          2           3       3
         - 18*sinh(z)*pi *z + 16*sinh(z)*z  + 18*i*pi *z - 16*i*z  + 2*pi

                  2             2
         - 36*pi*z )/(4*(sinh(z)  - 2*sinh(z)*i - 1))

z0 := pi*i/2;


       i*pi
z0 := ------
        2

taylor(dz2,z,z0,6);


                2
           i*(pi  - 16)        i*pi      pi        i*pi  2
 - 2*pi + --------------*(z - ------) + ----*(z - ------)
                4               2        2          2

              2
    i*( - 3*pi  + 80)        i*pi  3    pi        i*pi  4
 + -------------------*(z - ------)  - ----*(z - ------)
           120                2         24         2

           2
    i*(5*pi  - 168)        i*pi  5                      i*pi  7
 + -----------------*(z - ------)  + (1 term) + O((z - ------) )
         3360               2                            2


zz3:=(z*(z-2*pi)*(z-pi/2)^2)/(sin z-1);


                  3       2            2      3
        z*( - 2*pi  + 9*pi *z - 12*pi*z  + 4*z )
zz3 := ------------------------------------------
                     4*(sin(z) - 1)

dz3 := df(zz3,z);


                   3                2  2                 3             4
dz3 := (2*cos(z)*pi *z - 9*cos(z)*pi *z  + 12*cos(z)*pi*z  - 4*cos(z)*z

                      3               2                   2              3
         - 2*sin(z)*pi  + 18*sin(z)*pi *z - 36*sin(z)*pi*z  + 16*sin(z)*z

               3        2            2       3            2
         + 2*pi  - 18*pi *z + 36*pi*z  - 16*z )/(4*(sin(z)  - 2*sin(z) + 1))

z1 := pi/2;


       pi
z1 := ----
       2

taylor(dz3,z,z1,6);


          2                                            2
        pi  - 16        pi      pi        pi  2    3*pi  - 80        pi  3
2*pi + ----------*(z - ----) + ----*(z - ----)  + ------------*(z - ----)
           4            2       2         2           120            2

                           2
    pi        pi  4    5*pi  - 168        pi  5                      pi  7
 + ----*(z - ----)  + -------------*(z - ----)  + (1 term) + O((z - ----) )
    24        2           3360            2                          2


taylor((sin tan x-tan sin x)/(asin atan x-atan asin x),x,0,6);


     5   2    1313   4    2773    6      7
1 + ---*x  + ------*x  - -------*x  + O(x )
     3        1890        11907



comment If the expansion point is not constant, it has to be taken
        care of in differentation, as the following examples show;


taylor(sin(x+a),x,a,8);


                               sin(2*a)         2    cos(2*a)         3
sin(2*a) + cos(2*a)*(x - a) - ----------*(x - a)  - ----------*(x - a)
                                  2                     6

    sin(2*a)         4    cos(2*a)         5                        9
 + ----------*(x - a)  + ----------*(x - a)  + (3 terms) + O((x - a) )
       24                   120

df(ws,a);


                               cos(2*a)         2    sin(2*a)         3
cos(2*a) - sin(2*a)*(x - a) - ----------*(x - a)  + ----------*(x - a)
                                  2                     6

    cos(2*a)         4    sin(2*a)         5                        8
 + ----------*(x - a)  - ----------*(x - a)  + (2 terms) + O((x - a) )
       24                   120

taylor(cos(x+a),x,a,7);


                               cos(2*a)         2    sin(2*a)         3
cos(2*a) - sin(2*a)*(x - a) - ----------*(x - a)  + ----------*(x - a)
                                  2                     6

    cos(2*a)         4    sin(2*a)         5                        8
 + ----------*(x - a)  - ----------*(x - a)  + (2 terms) + O((x - a) )
       24                   120



comment A problem are non-analytical terms: rational powers and
        logarithmic terms can be handled, but other types of essential
        singularities cannot;


taylor(sqrt(x),x,0,2);


 1/2      3
x    + O(x )


taylor(asinh(1/x),x,0,5);


                       1   2    3    4    5    6      7
 - log(x) + (log(2) + ---*x  - ----*x  + ----*x  + O(x ))
                       4        32        96


taylor(e**(1/x),x,0,2);


        1/x
taylor(e   ,x,0,2)


comment Another example for non-integer powers;


sub (y = sqrt (x), yy);


     1/2    1       1   3/2    1    2      5/2
1 + x    + ---*x + ---*x    + ----*x  + O(x   )
            2       6          24


comment Expansion about infinity is possible in principle...;


taylor (e**(1/x), x, infinity, 5);


     1     1   1      1   1      1    1       1    1        1
1 + --- + ---*---- + ---*---- + ----*---- + -----*---- + O(----)
     x     2    2     6    3     24    4     120    5        6
               x          x           x            x        x

xi := taylor (sin (1/x), x, infinity, 5);


       1     1   1       1    1        1
xi := --- - ---*---- + -----*---- + O(----)
       x     6    3     120    5        6
                 x            x        x


y1 := taylor(x/(x-1), x, infinity, 3);


           1     1      1        1
y1 := 1 + --- + ---- + ---- + O(----)
           x      2      3        4
                 x      x        x

z := df(y1, x);


         1        1        1        1
z :=  - ---- - 2*---- - 3*---- + O(----)
          2        3        4        5
         x        x        x        x


comment ...but far from being perfect;


taylor (1 / sin (x), x, infinity, 5);


          1
taylor(--------,x,infinity,5)
        sin(x)


clear z;



comment You may access the expansion with the PART operator;


part(yy,0);


plus

part(yy,1);


1

part(yy,4);


  3
 y
----
 6

part(yy,6);


***** Expression taylor(1 + y + 1/2*y**2 + 1/6*y**3 + 1/24*y**4,y,0,4) 
does not have part 6



comment The template of a Taylor kernel can be extracted;


taylortemplate yy;


{{y,0,4}}


taylortemplate xxa;


{{x,1,2}}


taylortemplate xi;


{{x,infinity,5}}


taylortemplate xy;


{{x,0,2},{y,0,2}}


taylortemplate xx1;


{{{x,y},0,2}}


comment Here is a slightly less trivial example;


exp := (sin (x) * sin (y) / (x * y))**2;


              2       2
        sin(x) *sin(y)
exp := -----------------
              2  2
             x *y


taylor (exp, x, 0, 1, y, 0, 1);


       2  2
1 + O(x ,y )

taylor (exp, x, 0, 2, y, 0, 2);


     1   2    1   2    1   2  2      3  3
1 - ---*x  - ---*y  + ---*y *x  + O(x ,y )
     3        3        9


tt := taylor (exp, {x,y}, 0, 2);


           1   2    1   2          3
tt := 1 - ---*y  - ---*x  + O({x,y} )
           3        3


comment An example that uses factorization;


on factor;



ff := y**5 - 1;


        4    3    2
ff := (y  + y  + y  + y + 1)*(y - 1)


zz := sub (y = taylor(e**x, x, 0, 3), ff);


                 1   2    1   3      4  4             1   2    1   3      4  3
zz := ((1 + x + ---*x  + ---*x  + O(x ))  + (1 + x + ---*x  + ---*x  + O(x ))
                 2        6                           2        6

                    1   2    1   3      4  2             1   2    1   3      4
        + (1 + x + ---*x  + ---*x  + O(x ))  + (1 + x + ---*x  + ---*x  + O(x ))
                    2        6                           2        6

                        1   2    1   3      4
        + 1)*((1 + x + ---*x  + ---*x  + O(x )) - 1)
                        2        6


on exp;



zz;


          1   2    1   3      4  5
(1 + x + ---*x  + ---*x  + O(x ))  - 1
          2        6


comment A simple example of Taylor kernel differentiation;


hugo := taylor(e^x,x,0,5);


                 1   2    1   3    1    4     1    5      6
hugo := 1 + x + ---*x  + ---*x  + ----*x  + -----*x  + O(x )
                 2        6        24        120

df(hugo^2,x);


             2    8   3    4   4      5
2 + 4*x + 4*x  + ---*x  + ---*x  + O(x )
                  3        3


comment The following shows the (limited) capabilities to integrate
        Taylor kernels. Only simple cases are supported, otherwise
        a warning is printed and the Taylor kernels are converted to
        standard representation;


zz := taylor (sin x, x, 0, 5);


           1   3     1    5      6
zz := x - ---*x  + -----*x  + O(x )
           6        120

ww := taylor (cos y, y, 0, 5);


           1   2    1    4      6
ww := 1 - ---*y  + ----*y  + O(y )
           2        24


int (zz, x);


 1   2    1    4     1    6      7
---*x  - ----*x  + -----*x  + O(x )
 2        24        720

int (ww, x);


     x   2    x    4      6
x - ---*y  + ----*y  + O(y )
     2        24

int (zz + ww, x);


  1   2    1    4     1    6      7           x   2    x    4      6
(---*x  - ----*x  + -----*x  + O(x )) + (x - ---*y  + ----*y  + O(y ))
  2        24        720                      2        24



comment And here we present Taylor series reversion.
        We start with the example given by Knuth for the algorithm;


taylor (t - t**2, t, 0, 5);


     2      6
t - t  + O(t )


taylorrevert (ws, t, x);


     2      3      4       5      6
x + x  + 2*x  + 5*x  + 14*x  + O(x )


tan!-series := taylor (tan x, x, 0, 5);


                   1   3    2    5      6
tan-series := x + ---*x  + ----*x  + O(x )
                   3        15

taylorrevert (tan!-series, x, y);


     1   3    1   5      6
y - ---*y  + ---*y  + O(y )
     3        5

atan!-series:=taylor (atan y, y, 0, 5);


                    1   3    1   5      6
atan-series := y - ---*y  + ---*y  + O(y )
                    3        5


tmp := taylor (e**x, x, 0, 5);


                1   2    1   3    1    4     1    5      6
tmp := 1 + x + ---*x  + ---*x  + ----*x  + -----*x  + O(x )
                2        6        24        120


taylorrevert (tmp, x, y);


         1         2    1         3    1         4    1         5            6
y - 1 - ---*(y - 1)  + ---*(y - 1)  - ---*(y - 1)  + ---*(y - 1)  + O((y - 1) )
         2              3              4              5


taylor (log y, y, 1, 5);


         1         2    1         3    1         4    1         5            6
y - 1 - ---*(y - 1)  + ---*(y - 1)  - ---*(y - 1)  + ---*(y - 1)  + O((y - 1) )
         2              3              4              5



comment The following example calculates the perturbation expansion
        of the root x = 20 of the following polynomial in terms of
        EPS, in ROUNDED mode;


poly := for r := 1 : 20 product (x - r);


         20        19          18            17             16               15
poly := x   - 210*x   + 20615*x   - 1256850*x   + 53327946*x   - 1672280820*x

                        14                 13                   12
         + 40171771630*x   - 756111184500*x   + 11310276995381*x

                            11                     10                      9
         - 135585182899530*x   + 1307535010540395*x   - 10142299865511450*x

                              8                       7                        6
         + 63030812099294896*x  - 311333643161390640*x  + 1206647803780373360*x

                                5                        4
         - 3599979517947607200*x  + 8037811822645051776*x

                                 3                         2
         - 12870931245150988800*x  + 13803759753640704000*x

         - 8752948036761600000*x + 2432902008176640000


on rounded;



tpoly := taylor (poly, x, 20, 4);


                                                                2
tpoly := 1.21649393692e+17*(x - 20) + 4.31564847287e+17*(x - 20)

                                      3                             4
          + 6.68609351672e+17*(x - 20)  + 6.10115975015e+17*(x - 20)

                      5
          + O((x - 20) )


taylorrevert (tpoly, x, eps);


                                                  2                        3
20 + 8.22034512178e-18*eps - 2.39726594662e-34*eps  + 1.09290580232e-50*eps

                        4        5
 - 5.97114159465e-67*eps  + O(eps )


comment Some more examples using rounded mode;


taylor(sin x/x,x,0,4);


                    2                     4      5
1 - 0.166666666667*x  + 0.00833333333333*x  + O(x )


taylor(sin x,x,pi/2,4);


                                                                   2
1 + 6.12323399574e-17*(x - 1.57079632679) - 0.5*(x - 1.57079632679)

                                        3                                      4
 - 1.02053899929e-17*(x - 1.57079632679)  + 0.0416666666667*(x - 1.57079632679)

                        5
 + O((x - 1.57079632679) )


taylor(tan x,x,pi/2,4);


                      -1
 - (x - 1.57079632679)   + 0.333333333333*(x - 1.57079632679)

                                      3                        5
 + 0.0222222222222*(x - 1.57079632679)  + O((x - 1.57079632679) )


off rounded;




comment An example that involves computing limits of type 0/0 if
        expansion is done via differentiation;



taylor(sqrt((e^x - 1)/x),x,0,15);


     1       5    2     1    3     79     4      3     5                   16
1 + ---*x + ----*x  + -----*x  + -------*x  + -------*x  + (10 terms) + O(x  )
     4       96        128        92160        40960


comment An example that involves intermediate non-analytical terms
        which cancel entirely;


taylor(x^(5/2)/(log(x+1)*tan(x^(3/2))),x,0,5);


     1       1    2    7    3    139   4     67    5      6
1 + ---*x - ----*x  - ----*x  - -----*x  + ------*x  + O(x )
     2       12        24        720        1440


comment Other examples involving non-analytical terms;


taylor(log(e^x-1),x,0,5);


           1       1    2     1     4      5
log(x) + (---*x + ----*x  - ------*x  + O(x ))
           2       24        2880


taylor(e^(1/x)*(e^x-1),x,0,5);


 1/x       1   2    1   3    1    4     1    5      6
e   *(x + ---*x  + ---*x  + ----*x  + -----*x  + O(x ))
           2        6        24        120


taylor(log(x)*x^10,x,0,5);


         10      11
log(x)*(x   + O(x  ))


taylor(log(x)*x^10,x,0,11);


         10      12
log(x)*(x   + O(x  ))


taylor(log(x-a)/((a-b)*(a-c)) + log(2(x-b))/((b-c)*(b-a))
     + log(x-c)/((c-a)*(c-b)),x,infinity,2);


           log(2)            1   1        1
 - ---------------------- - ---*---- + O(----)
                 2           2    2        3
    a*b - a*c - b  + b*c         x        x


ss := (sqrt(x^(2/5) +1) - x^(1/3)-1)/x^(1/3);


             2/5         1/3
       sqrt(x    + 1) - x    - 1
ss := ---------------------------
                  1/3
                 x

taylor(exp ss,x,0,2);


 1      1    1/15     1    2/15     1     1/5      1     4/15      1      1/3
--- + -----*x     + -----*x     + ------*x    + -------*x     + --------*x
 e     2*e           8*e           48*e          384*e           3840*e

                   31/15
 + (25 terms) + O(x     )


taylor(exp sub(x=x^15,ss),x,0,2);


 1      1         1    2      3
--- + -----*x + -----*x  + O(x )
 e     2*e       8*e


taylor(dilog(x),x,0,4);


             1   2    1   3    1   4      5
log(x)*(x + ---*x  + ---*x  + ---*x  + O(x ))
             2        3        4

       2
     pi          1   2    1   3    1    4      5
 + (----- - x - ---*x  - ---*x  - ----*x  + O(x ))
      6          4        9        16


taylor(ei(x),x,0,4);


                        1   2    1    3    1    4      5
log(x) - psi(1) + (x + ---*x  + ----*x  + ----*x  + O(x ))
                        4        18        96


comment In the following we demonstrate the possibiblity to compute the
        expansion of a function which is given by a simple first order
        differential equation: the function myexp(x) is exp(-x^2);


operator myexp,myerf;


let {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1,
     df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0};


taylor(myexp(x),x,0,5);


     2    1   4      6
1 - x  + ---*x  + O(x )
          2

taylor(myerf(x),x,0,5);


    2               2        3        1        5      6
----------*x - ------------*x  + ------------*x  + O(x )
 sqrt(pi)       3*sqrt(pi)        5*sqrt(pi)

clear {df(myexp(~x),~x) => -2*x*myexp(x), myexp(0) => 1,
       df(myerf(~x),~x) => 2/sqrt(pi)*myexp(x), myerf(0) => 0};


clear myexp,erf;



%%% showtime;

comment There are two special operators, implicit_taylor and
        inverse_taylor, to compute the Taylor expansion of implicit
        or inverse functions;



implicit_taylor(x^2 + y^2 - 1,x,y,0,1,5);


     1   2    1   4      6
1 - ---*x  - ---*x  + O(x )
     2        8


implicit_taylor(x^2 + y^2 - 1,x,y,0,1,20);


     1   2    1   4    1    6     5    8     7    10                  21
1 - ---*x  - ---*x  - ----*x  - -----*x  - -----*x   + (5 terms) + O(x  )
     2        8        16        128        256


implicit_taylor(x+y^3-y,x,y,0,0,8);


     3      5       7      9
x + x  + 3*x  + 12*x  + O(x )


implicit_taylor(x+y^3-y,x,y,0,1,5);


     1       3   2    1   3    105   4    3   5      6
1 - ---*x - ---*x  - ---*x  - -----*x  - ---*x  + O(x )
     2       8        2        128        2


implicit_taylor(x+y^3-y,x,y,0,-1,5);


        1       3   2    1   3    105   4    3   5      6
 - 1 - ---*x + ---*x  - ---*x  + -----*x  - ---*x  + O(x )
        2       8        2        128        2


implicit_taylor(y*e^y-x,x,y,0,0,5);


     2    3   3    8   4    125   5      6
x - x  + ---*x  - ---*x  + -----*x  + O(x )
          2        3        24


comment This is the function exp(-1/x^2), which has an essential
        singularity at the point 0;


implicit_taylor(x^2*log y+1,x,y,0,0,3);


***** Computation of Taylor series of implicit function failed 
      Input expression non-zero at given point 


inverse_taylor(exp(x)-1,x,y,0,8);


     1   2    1   3    1   4    1   5    1   6                  9
y - ---*y  + ---*y  - ---*y  + ---*y  - ---*y  + (2 terms) + O(y )
     2        3        4        5        6


inverse_taylor(exp(x),x,y,0,5);


         1         2    1         3    1         4    1         5            6
y - 1 - ---*(y - 1)  + ---*(y - 1)  - ---*(y - 1)  + ---*(y - 1)  + O((y - 1) )
         2              3              4              5


inverse_taylor(sqrt(x),x,y,0,5);


 2      6
y  + O(y )


inverse_taylor(log(1+x),x,y,0,5);


     1   2    1   3    1    4     1    5      6
y + ---*y  + ---*y  + ----*y  + -----*y  + O(y )
     2        6        24        120


inverse_taylor((e^x-e^(-x))/2,x,y,0,5);


     1   3    3    5      6
y - ---*y  + ----*y  + O(y )
     6        40


comment In the next two cases the inverse functions have a branch
        point, therefore the computation fails;


inverse_taylor((e^x+e^(-x))/2,x,y,0,5);


***** Computation of Taylor series of inverse function failed 


inverse_taylor(exp(x^2-1),x,y,0,5);


***** Computation of Taylor series of inverse function failed 


inverse_taylor(exp(sqrt(x))-1,x,y,0,5);


 2    3    11   4    5   5      6
y  - y  + ----*y  - ---*y  + O(y )
           12        6


inverse_taylor(x*exp(x),x,y,0,5);


     2    3   3    8   4    125   5      6
y - y  + ---*y  - ---*y  + -----*y  + O(y )
          2        3        24



%%% showtime;


comment An application is the problem posed by Prof. Stanley:
        we prove that the finite difference expression below
        corresponds to the given derivative expression;


operator diff,a,f,gg;

  % We use gg to avoid conflict with high energy
                       % physics operator.

let diff(~f,~arg) => df(f,arg);



derivative_expression :=
diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),y),x) +
diff(a(x,y)*diff(gg(x,y),x)*diff(gg(x,y),y)*diff(f(x,y),x),y) ;


derivative_expression := 2*a(x,y)*df(f(x,y),x,y)*df(gg(x,y),x)*df(gg(x,y),y)

 + a(x,y)*df(f(x,y),x)*df(gg(x,y),x,y)*df(gg(x,y),y)

 + a(x,y)*df(f(x,y),x)*df(gg(x,y),x)*df(gg(x,y),y,2)

 + a(x,y)*df(f(x,y),y)*df(gg(x,y),x,y)*df(gg(x,y),x)

 + a(x,y)*df(f(x,y),y)*df(gg(x,y),x,2)*df(gg(x,y),y)

 + df(a(x,y),x)*df(f(x,y),y)*df(gg(x,y),x)*df(gg(x,y),y)

 + df(a(x,y),y)*df(f(x,y),x)*df(gg(x,y),x)*df(gg(x,y),y)


finite_difference_expression :=
+a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
+a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
+a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
+a(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x+dx,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x,y+dy)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x+dx,y+dy)^2/(32*dx^2*dy^2)
-gg(x,y)*a(x+dx,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x+dx,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x,y+dy)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
-a(x,y)*gg(x,y)*f(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x+dx,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y+dy)/(16*dx^2*dy^2)
-gg(x+dx,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)*f(x+dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y+dy)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x+dx,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2)
-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2)
-a(x,y+dy)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x+dx,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2)
+a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2)
+a(x,y)*gg(x,y+dy)*gg(x+dx,y)*f(x+dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y+dy)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x+dx,y)*f(x+dx,y+dy)/(32*dx^2*dy^2)
-a(x,y+dy)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x,y+dy)^2*f(x+dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x,y+dy)*f(x+dx,y+dy)/(32*dx^2*dy^2)
+a(x,y)*gg(x,y)^2*f(x+dx,y+dy)/(32*dx^2*dy^2)
+f(x,y)*gg(x+dx,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y+dy)*gg(x+dx,y)*a(x+dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y+dy)^2*a(x+dx,y+dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x+dx,y+dy)/(32*dx^2*dy^2)
+a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
+a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
+a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
+a(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x+dx,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x,y-dy)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x+dx,y-dy)^2/(32*dx^2*dy^2)
-gg(x,y)*a(x+dx,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x+dx,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x,y-dy)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
-a(x,y)*gg(x,y)*f(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x+dx,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x+dx,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x,y)*gg(x+dx,y-dy)/(16*dx^2*dy^2)
-gg(x+dx,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)*f(x+dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y-dy)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x+dx,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2)
-a(x+dx,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2)
-a(x,y-dy)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x+dx,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2)
+a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2)
+a(x,y)*gg(x,y-dy)*gg(x+dx,y)*f(x+dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y-dy)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x+dx,y)*f(x+dx,y-dy)/(32*dx^2*dy^2)
-a(x,y-dy)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x,y-dy)^2*f(x+dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x,y-dy)*f(x+dx,y-dy)/(32*dx^2*dy^2)
+a(x,y)*gg(x,y)^2*f(x+dx,y-dy)/(32*dx^2*dy^2)
+f(x,y)*gg(x+dx,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y-dy)*gg(x+dx,y)*a(x+dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y-dy)^2*a(x+dx,y-dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x+dx,y-dy)/(32*dx^2*dy^2)
+f(x,y)*a(x+dx,y)*gg(x+dx,y)^2/(16*dx^2*dy^2)
+f(x,y)*a(x,y+dy)*gg(x+dx,y)^2/(32*dx^2*dy^2)
+f(x,y)*a(x,y-dy)*gg(x+dx,y)^2/(32*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x+dx,y)^2/(16*dx^2*dy^2)
-f(x,y)*gg(x,y+dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2)
-f(x,y)*gg(x,y-dy)*a(x+dx,y)*gg(x+dx,y)/(16*dx^2*dy^2)
-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x+dx,y)/(16*dx^2*dy^2)
-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x+dx,y)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y+dy)^2*a(x+dx,y)/(32*dx^2*dy^2)
+f(x,y)*gg(x,y-dy)^2*a(x+dx,y)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x+dx,y)/(16*dx^2*dy^2)
+a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
+a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
+a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
+a(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x-dx,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x,y+dy)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x-dx,y+dy)^2/(32*dx^2*dy^2)
-gg(x,y)*a(x-dx,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x-dx,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x,y+dy)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
-a(x,y)*gg(x,y)*f(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x-dx,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x,y+dy)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y+dy)/(16*dx^2*dy^2)
-gg(x-dx,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)*f(x-dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y+dy)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x-dx,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2)
-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2)
-a(x,y+dy)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x-dx,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2)
+a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2)
+a(x,y)*gg(x,y+dy)*gg(x-dx,y)*f(x-dx,y+dy)/(16*dx^2*dy^2)
-gg(x,y+dy)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x-dx,y)*f(x-dx,y+dy)/(32*dx^2*dy^2)
-a(x,y+dy)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x,y+dy)^2*f(x-dx,y+dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x,y+dy)*f(x-dx,y+dy)/(32*dx^2*dy^2)
+a(x,y)*gg(x,y)^2*f(x-dx,y+dy)/(32*dx^2*dy^2)
+f(x,y)*gg(x-dx,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y+dy)*gg(x-dx,y)*a(x-dx,y+dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y+dy)^2*a(x-dx,y+dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x-dx,y+dy)/(32*dx^2*dy^2)
+a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
+a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
+a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
+a(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x-dx,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
-f(x,y)*a(x,y-dy)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x-dx,y-dy)^2/(32*dx^2*dy^2)
-gg(x,y)*a(x-dx,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x-dx,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y)*a(x,y-dy)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
-a(x,y)*gg(x,y)*f(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x-dx,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x-dx,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y)*a(x,y-dy)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x,y)*gg(x-dx,y-dy)/(16*dx^2*dy^2)
-gg(x-dx,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)*f(x-dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y-dy)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x-dx,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2)
-a(x-dx,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2)
-a(x,y-dy)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x-dx,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2)
+a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2)
+a(x,y)*gg(x,y-dy)*gg(x-dx,y)*f(x-dx,y-dy)/(16*dx^2*dy^2)
-gg(x,y-dy)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x-dx,y)*f(x-dx,y-dy)/(32*dx^2*dy^2)
-a(x,y-dy)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2)
-a(x,y)*gg(x,y-dy)^2*f(x-dx,y-dy)/(32*dx^2*dy^2)
+gg(x,y)^2*a(x,y-dy)*f(x-dx,y-dy)/(32*dx^2*dy^2)
+a(x,y)*gg(x,y)^2*f(x-dx,y-dy)/(32*dx^2*dy^2)
+f(x,y)*gg(x-dx,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y-dy)*gg(x-dx,y)*a(x-dx,y-dy)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y-dy)^2*a(x-dx,y-dy)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x-dx,y-dy)/(32*dx^2*dy^2)
+f(x,y)*a(x-dx,y)*gg(x-dx,y)^2/(16*dx^2*dy^2)
+f(x,y)*a(x,y+dy)*gg(x-dx,y)^2/(32*dx^2*dy^2)
+f(x,y)*a(x,y-dy)*gg(x-dx,y)^2/(32*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x-dx,y)^2/(16*dx^2*dy^2)
-f(x,y)*gg(x,y+dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2)
-f(x,y)*gg(x,y-dy)*a(x-dx,y)*gg(x-dx,y)/(16*dx^2*dy^2)
-f(x,y)*a(x,y+dy)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x,y+dy)*gg(x-dx,y)/(16*dx^2*dy^2)
-f(x,y)*a(x,y-dy)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x,y-dy)*gg(x-dx,y)/(16*dx^2*dy^2)
+f(x,y)*gg(x,y+dy)^2*a(x-dx,y)/(32*dx^2*dy^2)
+f(x,y)*gg(x,y-dy)^2*a(x-dx,y)/(32*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x-dx,y)/(16*dx^2*dy^2)
+f(x,y)*a(x,y+dy)*gg(x,y+dy)^2/(16*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x,y+dy)^2/(16*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x,y+dy)/(16*dx^2*dy^2)
+f(x,y)*a(x,y-dy)*gg(x,y-dy)^2/(16*dx^2*dy^2)
+a(x,y)*f(x,y)*gg(x,y-dy)^2/(16*dx^2*dy^2)
-f(x,y)*gg(x,y)^2*a(x,y-dy)/(16*dx^2*dy^2)
-a(x,y)*f(x,y)*gg(x,y)^2/(8*dx^2*dy^2)$



comment We define abbreviations for the partial derivatives;


operator ax,ay,fx,fy,gx,gy;


operator axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy;


operator axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy;


operator axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,
         gxxxx,gxxxy,gxxyy,gxyyy,gyyyy;


operator axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy;


operator gxxxxyy,gxxxyyy,gxxyyyy;



operator_diff_rules := {
  df(a(~x,~y),~x) => ax(x,y),
  df(a(~x,~y),~y) => ay(x,y),
  df(f(~x,~y),~x) => fx(x,y),
  df(f(~x,~y),~y) => fy(x,y),
  df(gg(~x,~y),~x) => gx(x,y),
  df(gg(~x,~y),~y) => gy(x,y),

  df(ax(~x,~y),~x) => axx(x,y),
  df(ax(~x,~y),~y) => axy(x,y),
  df(ay(~x,~y),~x) => axy(x,y),
  df(ay(~x,~y),~y) => ayy(x,y),
  df(fx(~x,~y),~x) => fxx(x,y),
  df(fx(~x,~y),~y) => fxy(x,y),
  df(fy(~x,~y),~x) => fxy(x,y),
  df(fy(~x,~y),~y) => fyy(x,y),
  df(gx(~x,~y),~x) => gxx(x,y),
  df(gx(~x,~y),~y) => gxy(x,y),
  df(gy(~x,~y),~x) => gxy(x,y),
  df(gy(~x,~y),~y) => gyy(x,y),

  df(axx(~x,~y),~x) => axxx(x,y),
  df(axy(~x,~y),~x) => axxy(x,y),
  df(ayy(~x,~y),~x) => axyy(x,y),
  df(ayy(~x,~y),~y) => ayyy(x,y),
  df(fxx(~x,~y),~x) => fxxx(x,y),
  df(fxy(~x,~y),~x) => fxxy(x,y),
  df(fxy(~x,~y),~y) => fxyy(x,y),
  df(fyy(~x,~y),~x) => fxyy(x,y),
  df(fyy(~x,~y),~y) => fyyy(x,y),
  df(gxx(~x,~y),~x) => gxxx(x,y),
  df(gxx(~x,~y),~y) => gxxy(x,y),
  df(gxy(~x,~y),~x) => gxxy(x,y),
  df(gxy(~x,~y),~y) => gxyy(x,y),
  df(gyy(~x,~y),~x) => gxyy(x,y),
  df(gyy(~x,~y),~y) => gyyy(x,y),

  df(axyy(~x,~y),~x) => axxyy(x,y),
  df(axxy(~x,~y),~x) => axxxy(x,y),
  df(ayyy(~x,~y),~x) => axyyy(x,y),
  df(fxxy(~x,~y),~x) => fxxxy(x,y),
  df(fxyy(~x,~y),~x) => fxxyy(x,y),
  df(fyyy(~x,~y),~x) => fxyyy(x,y),
  df(gxxx(~x,~y),~x) => gxxxx(x,y),
  df(gxxy(~x,~y),~x) => gxxxy(x,y),
  df(gxyy(~x,~y),~x) => gxxyy(x,y),
  df(gyyy(~x,~y),~x) => gxyyy(x,y),
  df(gyyy(~x,~y),~y) => gyyyy(x,y),

  df(axxyy(~x,~y),~x) => axxxyy(x,y),
  df(axyyy(~x,~y),~x) => axxyyy(x,y),
  df(fxxyy(~x,~y),~x) => fxxxyy(x,y),
  df(fxyyy(~x,~y),~x) => fxxyyy(x,y),
  df(gxxxy(~x,~y),~x) => gxxxxy(x,y),
  df(gxxyy(~x,~y),~x) => gxxxyy(x,y),
  df(gxyyy(~x,~y),~x) => gxxyyy(x,y),
  df(gyyyy(~x,~y),~x) => gxyyyy(x,y),

  df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y),
  df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y),
  df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y)
};


operator_diff_rules := {df(a(~x,~y),~x) => ax(x,y),

   df(a(~x,~y),~y) => ay(x,y),

   df(f(~x,~y),~x) => fx(x,y),

   df(f(~x,~y),~y) => fy(x,y),

   df(gg(~x,~y),~x) => gx(x,y),

   df(gg(~x,~y),~y) => gy(x,y),

   df(ax(~x,~y),~x) => axx(x,y),

   df(ax(~x,~y),~y) => axy(x,y),

   df(ay(~x,~y),~x) => axy(x,y),

   df(ay(~x,~y),~y) => ayy(x,y),

   df(fx(~x,~y),~x) => fxx(x,y),

   df(fx(~x,~y),~y) => fxy(x,y),

   df(fy(~x,~y),~x) => fxy(x,y),

   df(fy(~x,~y),~y) => fyy(x,y),

   df(gx(~x,~y),~x) => gxx(x,y),

   df(gx(~x,~y),~y) => gxy(x,y),

   df(gy(~x,~y),~x) => gxy(x,y),

   df(gy(~x,~y),~y) => gyy(x,y),

   df(axx(~x,~y),~x) => axxx(x,y),

   df(axy(~x,~y),~x) => axxy(x,y),

   df(ayy(~x,~y),~x) => axyy(x,y),

   df(ayy(~x,~y),~y) => ayyy(x,y),

   df(fxx(~x,~y),~x) => fxxx(x,y),

   df(fxy(~x,~y),~x) => fxxy(x,y),

   df(fxy(~x,~y),~y) => fxyy(x,y),

   df(fyy(~x,~y),~x) => fxyy(x,y),

   df(fyy(~x,~y),~y) => fyyy(x,y),

   df(gxx(~x,~y),~x) => gxxx(x,y),

   df(gxx(~x,~y),~y) => gxxy(x,y),

   df(gxy(~x,~y),~x) => gxxy(x,y),

   df(gxy(~x,~y),~y) => gxyy(x,y),

   df(gyy(~x,~y),~x) => gxyy(x,y),

   df(gyy(~x,~y),~y) => gyyy(x,y),

   df(axyy(~x,~y),~x) => axxyy(x,y),

   df(axxy(~x,~y),~x) => axxxy(x,y),

   df(ayyy(~x,~y),~x) => axyyy(x,y),

   df(fxxy(~x,~y),~x) => fxxxy(x,y),

   df(fxyy(~x,~y),~x) => fxxyy(x,y),

   df(fyyy(~x,~y),~x) => fxyyy(x,y),

   df(gxxx(~x,~y),~x) => gxxxx(x,y),

   df(gxxy(~x,~y),~x) => gxxxy(x,y),

   df(gxyy(~x,~y),~x) => gxxyy(x,y),

   df(gyyy(~x,~y),~x) => gxyyy(x,y),

   df(gyyy(~x,~y),~y) => gyyyy(x,y),

   df(axxyy(~x,~y),~x) => axxxyy(x,y),

   df(axyyy(~x,~y),~x) => axxyyy(x,y),

   df(fxxyy(~x,~y),~x) => fxxxyy(x,y),

   df(fxyyy(~x,~y),~x) => fxxyyy(x,y),

   df(gxxxy(~x,~y),~x) => gxxxxy(x,y),

   df(gxxyy(~x,~y),~x) => gxxxyy(x,y),

   df(gxyyy(~x,~y),~x) => gxxyyy(x,y),

   df(gyyyy(~x,~y),~x) => gxyyyy(x,y),

   df(gxxxyy(~x,~y),~x) => gxxxxyy(x,y),

   df(gxxyyy(~x,~y),~x) => gxxxyyy(x,y),

   df(gxyyyy(~x,~y),~x) => gxxyyyy(x,y)}


let operator_diff_rules;



texp := taylor (finite_difference_expression, dx, 0, 1, dy, 0, 1);


texp := a(x,y)*fx(x,y)*gx(x,y)*gyy(x,y) + a(x,y)*fx(x,y)*gxy(x,y)*gy(x,y)

         + 2*a(x,y)*fxy(x,y)*gx(x,y)*gy(x,y) + a(x,y)*fy(x,y)*gx(x,y)*gxy(x,y)

         + a(x,y)*fy(x,y)*gxx(x,y)*gy(x,y) + ax(x,y)*fy(x,y)*gx(x,y)*gy(x,y)

                                                 2   2
         + ay(x,y)*fx(x,y)*gx(x,y)*gy(x,y) + O(dx ,dy )


comment You may also try to expand further but this needs a lot
        of CPU time.  Therefore the following line is commented out;


%texp := taylor (finite_difference_expression, dx, 0, 2, dy, 0, 2);

factor dx,dy;



result := taylortostandard texp;


result := a(x,y)*fx(x,y)*gx(x,y)*gyy(x,y) + a(x,y)*fx(x,y)*gxy(x,y)*gy(x,y)

           + 2*a(x,y)*fxy(x,y)*gx(x,y)*gy(x,y) + a(x,y)*fy(x,y)*gx(x,y)*gxy(x,y)

           + a(x,y)*fy(x,y)*gxx(x,y)*gy(x,y) + ax(x,y)*fy(x,y)*gx(x,y)*gy(x,y)

           + ay(x,y)*fx(x,y)*gx(x,y)*gy(x,y)


derivative_expression - result;


0



clear diff(~f,~arg);


clearrules operator_diff_rules;


clear diff,a,f,gg;


clear ax,ay,fx,fy,gx,gy;


clear axx,axy,ayy,fxx,fxy,fyy,gxx,gxy,gyy;


clear axxx,axxy,axyy,ayyy,fxxx,fxxy,fxyy,fyyy,gxxx,gxxy,gxyy,gyyy;


clear axxxy,axxyy,axyyy,fxxxy,fxxyy,fxyyy,gxxxx,gxxxy,gxxyy,gxyyy,gyyyy;


clear axxxyy,axxyyy,fxxyyy,fxxxyy,gxxxxy,gxxxyy,gxxyyy,gxyyyy;


clear gxxxxyy,gxxxyyy,gxxyyyy;



taylorprintterms := 5;


taylorprintterms := 5


off taylorautoexpand,taylorkeeporiginal;



%%% showtime;

comment An example provided by Alan Barnes: elliptic functions;


% Jacobi's elliptic functions 
%       sn(x,k),  cn(x,k), dn(x,k).
% The modulus and complementary modulus are denoted by K and K!'
%      usually written mathematically as k and k' respectively
%
% epsilon(x,k) is the incomplete elliptic integral of the second kind
%      usually written mathematically as E(x,k)
%
% KK(k) is the complete elliptic integral of the first kind
%       usually written mathematically as K(k)
%       K(k) = arcsn(1,k)
% KK!'(k) is the complementary complete integral
%      usually written mathematically as K'(k)
%      NB.  K'(k) = K(k')
% EE(k) is the complete elliptic integral of the second kind
%      usually written mathematically as E(k)
% EE!'(k) is the complementary complete integral
%      usually written mathematically as E'(k)
%      NB.  E'(k) = E(k')

operator sn, cn, dn, epsilon;



elliptic_rules := {
% Differentiation rules for basic functions
   df(sn(~x,~k),~x)     => cn(x,k)*dn(x,k),
   df(cn(~x,~k),~x)     => -sn(x,k)*dn(x,k),
   df(dn(~x,~k),~x)     => -k^2*sn(x,k)*cn(x,k),
   df(epsilon(~x,~k),~x)=> dn(x,k)^2,

% k-derivatives
% DF Lawden     Elliptic Functions & Applications    Springer (1989)
   df(sn(~x,~k),~k)  => (k*sn(x,k)*cn(x,k)^2
                         -epsilon(x,k)*cn(x,k)*dn(x,k)/k)/(1-k^2)
                       + x*cn(x,k)*dn(x,k)/k,
   df(cn(~x,~k),~k)  => (-k*sn(x,k)^2*cn(x,k)
                         +epsilon(x,k)*sn(x,k)*dn(x,k)/k)/(1-k^2)
                        - x*sn(x,k)*dn(x,k)/k,
   df(dn(~x,~k),~k)  => k*(-sn(x,k)^2*dn(x,k)
                            +epsilon(x,k)*sn(x,k)*cn(x,k))/(1-k^2)
                          - k*x*sn(x,k)*cn(x,k),
   df(epsilon(~x,~k),~k) => k*(sn(x,k)*cn(x,k)*dn(x,k)
                                -epsilon(x,k)*cn(x,k)^2)/(1-k^2)
                             -k*x*sn(x,k)^2,

% parity properties
   sn(-~x,~k) => -sn(x,k),
   cn(-~x,~k) => cn(x,k),
   dn(-~x,~k) => dn(x,k),
   epsilon(-~x,~k) => -epsilon(x,k),
   sn(~x,-~k) => sn(x,k),
   cn(~x,-~k) => cn(x,k),
   dn(~x,-~k) => dn(x,k),
   epsilon(~x,-~k) => epsilon(x,k),


% behaviour at zero
   sn(0,~k) => 0,
   cn(0,~k) => 1,
   dn(0,~k) => 1,
   epsilon(0,~k) => 0,


% degenerate cases of modulus
   sn(~x,0) => sin(x),
   cn(~x,0) => cos(x),
   dn(~x,0) => 1,
   epsilon(~x,0) => x,

   sn(~x,1) => tanh(x),
   cn(~x,1) => 1/cosh(x),
   dn(~x,1) => 1/cosh(x),
   epsilon(~x,1) => tanh(x)
};


elliptic_rules := {df(sn(~x,~k),~x) => cn(x,k)*dn(x,k),

   df(cn(~x,~k),~x) =>  - sn(x,k)*dn(x,k),

                           2
   df(dn(~x,~k),~x) =>  - k *sn(x,k)*cn(x,k),

                                   2
   df(epsilon(~x,~k),~x) => dn(x,k) ,

                                         2                         dn(x,k)
                        k*sn(x,k)*cn(x,k)  - epsilon(x,k)*cn(x,k)*---------
                                                                      k
   df(sn(~x,~k),~k) => -----------------------------------------------------
                                                   2
                                              1 - k

                 dn(x,k)
    + x*cn(x,k)*---------,
                    k

                                    2                                 dn(x,k)
                         - k*sn(x,k) *cn(x,k) + epsilon(x,k)*sn(x,k)*---------
                                                                         k
   df(cn(~x,~k),~k) => --------------------------------------------------------
                                                     2
                                                1 - k

                 dn(x,k)
    - x*sn(x,k)*---------,
                    k

                                    2
                           - sn(x,k) *dn(x,k) + epsilon(x,k)*sn(x,k)*cn(x,k)
   df(dn(~x,~k),~k) => k*----------------------------------------------------
                                                     2
                                                1 - k

    - k*x*sn(x,k)*cn(x,k),

   df(epsilon(~x,~k),~k)

                                                        2
          sn(x,k)*cn(x,k)*dn(x,k) - epsilon(x,k)*cn(x,k)                2
    => k*------------------------------------------------- - k*x*sn(x,k) ,
                                   2
                              1 - k

   sn( - ~x,~k) =>  - sn(x,k),

   cn( - ~x,~k) => cn(x,k),

   dn( - ~x,~k) => dn(x,k),

   epsilon( - ~x,~k) =>  - epsilon(x,k),

   sn(~x, - ~k) => sn(x,k),

   cn(~x, - ~k) => cn(x,k),

   dn(~x, - ~k) => dn(x,k),

   epsilon(~x, - ~k) => epsilon(x,k),

   sn(0,~k) => 0,

   cn(0,~k) => 1,

   dn(0,~k) => 1,

   epsilon(0,~k) => 0,

   sn(~x,0) => sin(x),

   cn(~x,0) => cos(x),

   dn(~x,0) => 1,

   epsilon(~x,0) => x,

   sn(~x,1) => tanh(x),

                   1
   cn(~x,1) => ---------,
                cosh(x)

                   1
   dn(~x,1) => ---------,
                cosh(x)

   epsilon(~x,1) => tanh(x)}


let elliptic_rules;



hugo := taylor(sn(x,k),k,0,6);


                                2                           2
                  cos(x)*(cos(x) *x + cos(x)*sin(x) + sin(x) *x - 2*x)   2
hugo := sin(x) + ------------------------------------------------------*k  + (
                                           4

                 5             4         2           4
           cos(x) *x - 2*cos(x) *sin(x)*x  + 5*cos(x) *sin(x)

                       3       2             3             2       3  2
            - 10*cos(x) *sin(x) *x + 6*cos(x) *x - 4*cos(x) *sin(x) *x

                    2       3           2         2           2
            + cos(x) *sin(x)  + 8*cos(x) *sin(x)*x  + 4*cos(x) *sin(x)

                              4                     2
            - 11*cos(x)*sin(x) *x + 30*cos(x)*sin(x) *x - 16*cos(x)*x

                      5  2           3  2             2      4
            - 2*sin(x) *x  + 8*sin(x) *x  - 8*sin(x)*x )/64*k  + (

                      7  3            7              6         2
            - 6*cos(x) *x  + 17*cos(x) *x - 99*cos(x) *sin(x)*x

                       6                   5       2  3            5       2
            + 21*cos(x) *sin(x) - 18*cos(x) *sin(x) *x  - 71*cos(x) *sin(x) *x

                       5  3            5               4       3  2
            + 36*cos(x) *x  - 18*cos(x) *x - 135*cos(x) *sin(x) *x

                        4       3             4         2             4
            - 133*cos(x) *sin(x)  + 324*cos(x) *sin(x)*x  + 172*cos(x) *sin(x)

                       3       4  3            3       4
            - 18*cos(x) *sin(x) *x  - 13*cos(x) *sin(x) *x

                       3       2  3             3       2              3  3
            + 72*cos(x) *sin(x) *x  - 156*cos(x) *sin(x) *x - 72*cos(x) *x

                        3              2       5  2             2       5
            + 160*cos(x) *x + 27*cos(x) *sin(x) *x  - 118*cos(x) *sin(x)

                        2       3             2         2            2
            + 176*cos(x) *sin(x)  - 108*cos(x) *sin(x)*x  + 32*cos(x) *sin(x)

                             6  3                   6                     4  3
            - 6*cos(x)*sin(x) *x  + 75*cos(x)*sin(x) *x + 36*cos(x)*sin(x) *x

                               4                     2  3                    2
            - 498*cos(x)*sin(x) *x - 72*cos(x)*sin(x) *x  + 888*cos(x)*sin(x) *x

                         3                           7  2             5  2
            + 48*cos(x)*x  - 384*cos(x)*x + 63*sin(x) *x  - 324*sin(x) *x

                        3  2               2        6      7
            + 540*sin(x) *x  - 288*sin(x)*x )/2304*k  + O(k )

otto := taylor(cn(x,k),k,0,6);


                                   2                           2
                  sin(x)*( - cos(x) *x - cos(x)*sin(x) - sin(x) *x + 2*x)   2
otto := cos(x) + ---------------------------------------------------------*k  + 
                                             4

                    5  2           4                    3       2  2
        ( - 2*cos(x) *x  - 5*cos(x) *sin(x)*x - 4*cos(x) *sin(x) *x

                    3       2           3  2           2       3
          - 7*cos(x) *sin(x)  + 8*cos(x) *x  + 2*cos(x) *sin(x) *x

                    2                           4  2                  4
          + 2*cos(x) *sin(x)*x - 2*cos(x)*sin(x) *x  - 3*cos(x)*sin(x)

                           2  2                  2             2           5
          + 8*cos(x)*sin(x) *x  - 4*cos(x)*sin(x)  - 8*cos(x)*x  + 7*sin(x) *x

                     3                      4               7  2
          - 22*sin(x) *x + 16*sin(x)*x)/64*k  + ( - 9*cos(x) *x

                      6         3            6                      5       2  2
            + 6*cos(x) *sin(x)*x  - 71*cos(x) *sin(x)*x + 135*cos(x) *sin(x) *x

                       5       2            5  2            4       3  3
            - 66*cos(x) *sin(x)  - 36*cos(x) *x  + 18*cos(x) *sin(x) *x

                    4       3              4         3            4
            - cos(x) *sin(x) *x - 36*cos(x) *sin(x)*x  + 18*cos(x) *sin(x)*x

                        3       4  2            3       4
            + 297*cos(x) *sin(x) *x  + 61*cos(x) *sin(x)

                        3       2  2             3       2             3  2
            - 720*cos(x) *sin(x) *x  - 208*cos(x) *sin(x)  + 252*cos(x) *x

                       2       5  3            2       5
            + 18*cos(x) *sin(x) *x  + 31*cos(x) *sin(x) *x

                       2       3  3            2       3
            - 72*cos(x) *sin(x) *x  - 24*cos(x) *sin(x) *x

                       2         3            2                             6  2
            + 72*cos(x) *sin(x)*x  + 56*cos(x) *sin(x)*x + 153*cos(x)*sin(x) *x

                              6                    4  2                    4
            + 91*cos(x)*sin(x)  - 684*cos(x)*sin(x) *x  - 212*cos(x)*sin(x)

                               2  2                   2               2
            + 900*cos(x)*sin(x) *x  - 32*cos(x)*sin(x)  - 288*cos(x)*x

                      7  3            7              5  3             5
            + 6*sin(x) *x  - 39*sin(x) *x - 36*sin(x) *x  + 318*sin(x) *x

                       3  3             3                3
            + 72*sin(x) *x  - 672*sin(x) *x - 48*sin(x)*x  + 384*sin(x)*x)/2304

          6      7
        *k  + O(k )

taylorcombine(hugo^2 + otto^2);


      2         2      7
cos(x)  + sin(x)  + O(k )


clearrules elliptic_rules;



clear sn, cn, dn, epsilon;



%%% showtime;

comment That's all, folks;


end;


Time for test: 580 ms, plus GC time: 10 ms