%	diff.s
%
%	Purpose:	Symbolic differentiation and linearization
%	Copyright:	R.A. van Engelen, Leiden University, 1998

banner("Symbolic Differentiation and Linearization Package").

message("Usage: diff(Expression, Coordinate)'").
message("mode(diff, full):  full symbolic differentiation").
message("mode(diff, chain): chain-rule for symbolic differentiation").
message("mode(diff, none):  disable symbolic differentiation").

diff arity 2 is differentiation_op.
df arity 2 is differentiation_op.
int arity 2 is integration_op.
D arity 1 is differentiation_op.

:- assume df(_G2509, _G2510) has type(([T, _G2515] -> T)).
:- assume diff(_G2509, _G2510) has type(([T, _G2515] -> T)).
:- assume D(_G2509) has type(([T] -> T)).

:- insert (df(_G2509, _G2510) can           get_latex(df(U, N), none, [ "D_"                                     , N                                     , "("                                     , U                                     , ")" ]) do number(N)          ).
:- insert (df(_G2509, _G2510) can           get_latex(df(df(U, X), X), 300, [ "\textstyle{\frac{\partial^2"                                           , U                                           , "}{\partial"                                           , X                                           , "^2}}" ])          ).
:- insert (df(_G2509, _G2510) can           get_latex(df(U, X), 300, [ "\textstyle{\frac{\partial"                                    , U                                    , "}{\partial"                                    , X                                    , "}}" ])          ).
:- insert (diff(_G2509, _G2510) can           get_latex(diff(U, X), 300, [ "\textstyle{\frac{\partial"                                      , U                                      , "}{\partial"                                      , X                                      , "}}" ])          ).

:- insert (int(_G2509, _G2510) can           get_latex(int(E, X = L .. U), none, [ "\displaystyle\int_"                                               , L                                               , "^"                                               , U                                               , E                                               , "\,d"                                               , X ])          ).
:- insert (int(_G2509, _G2510) can           get_latex(int(E, X = Y), none, [ "\displaystyle\int_"                                          , Y                                          , E                                          , "\,d"                                          , X ])          ).
:- insert (int(_G2509, _G2510) can           get_latex(int(E, X), none, [ "\displaystyle\int"                                      , E                                      , "\,d"                                      , X ]) do atom(X)          ).
:- insert int(_G2509, _G2510) can get_latex(int(E, X), none, ["\displaystyle\int_", X, E]).
:- insert (integrate(_G2509, _G2510) can           get_latex(integrate(E, X = L .. U), none, [ "\displaystyle\int_"                                                     , L                                                     , "^"                                                     , U                                                     , E                                                     , "\,d"                                                     , X ])          ).
:- insert (integrate(_G2509, _G2510) can           get_latex(integrate(E, X = Y), none, [ "\displaystyle\int_"                                                , Y                                                , E                                                , "\,d"                                                , X ])          ).
:- insert (integrate(_G2509, _G2510) can           get_latex(integrate(E, X), none, [ "\displaystyle\int"                                            , E                                            , "\,d"                                            , X ]) do atom(X)          ).
:- insert integrate(_G2509, _G2510) can get_latex(integrate(E, X), none, ["\displaystyle\int_", X, E]).

diff(E, X) := proc(normalize(diff(E, X), diff)).

integrate(E, X) := proc(normalize(integrate(E, X), diff)).

(clear_D := proc((local L; L := linearized_var  get_instances; remove L))).

D(X) := proc(normalize(diff(X, nil), diff(full))).

d is inert has latex_name("\partial").
prefix d has latex_name("\partial").

d_times(_G2509, _G2510) is inert.
d_deriv(_G2509, _G2510) is inert.

d * d := d.
d ^ _G2510 := d.

d E := d_times([E], 1).
d d_times([E | Es], A) := d_times([E | Es], A).

(d_times([X | Xs], A) * d_times([Y | Ys], B) := d_times(append([Y | Ys], [X | Xs]), A * B)).

A * d_times([Y | Ys], B) := d_times([Y | Ys], A * B).
d_times([Y | Ys], A) * B := d_times([Y | Ys], A * B).

d / d_times([X | Xs], A) := d_deriv([X | Xs], 1 / A).

d_deriv([X], A) * E := A * diff(E, X).
d_deriv([X | Xs], A) * E := A * diff(d_deriv(Xs, 1) * E, X).

d_times([E], A) / d_times([X], B) := A / B * diff(E, X).
(d_times([E], A) / d_times([X | Xs], B) := A / B * diff(d_times([E], 1) / d_times(Xs, 1), X)).

(d_times([X | Xs], A) ^ N := (iff N >= 0 then d_times(flatten(fill(N, [X | Xs])), A ^ N) else error("Negative power"))).