This covers applications where manipulation of mathematical expressions, in general containing variables, is emphasized at the expense of actual numerical calculation.
There are several successful computer algebra systems, such as Axiom, Maple, and Mathematica, which can do certain symbolic operations that are useful in mathematics.
Examples including factorizing polynomials and differentiating and integrating expressions.
We will show how ML can be used for such applications. Our example will be symbolic differentiation.
This will illustrate all the typical components of symbolic computation systems: data representation, internal algorithms, parsing, and prettyprinting.
We will allow mathematical expressions to be built up from variables and constants by the application of n-ary operators.
Therefore we define a recursive type as follows:
- datatype term =
Var of string
| Const of string
| Fn of string * (term list);
For example the expression sin(x+y)/cos(x-exp(y))-ln(1+x) is represented by:
Fn("-",
[Fn("/",
[Fn("sin", [Fn("+", [Var "x", Var "y"])]),
Fn("cos", [Fn("-", [Var "x", Fn("exp", [Var "y"])])])]),
Fn("ln", [Fn("+", [Const "1", Var "x"])])]);
Reading and writing expressions in their raw form is rather unpleasant. This is a general problem in all symbolic computation systems. Typically one wants:
We will use parsing as another example of functional programming, so we will defer that for now.
However, we will now write a simple printer for our expressions.
How do we want to print expressions?
We can have a list of binary operators together with their precedences, e.g.
- val infixes = [("+",10), ("-",10),
("*",20), ("/",20),
("^",30)];
We can use this as an association list, which is just a simple form of dictionary as we've seen before. To get the precedence associated with a key we use:
- fun assoc ((x,y)::rest) a =
if a = x then y else assoc rest a;
In our case, we define:
- val get_precedence = assoc infixes;
Linear searches are inefficient for large dictionaries, but they are simple and have low overhead which makes them ideal for small lists like we use here.
Note that get_precedence captures a specific list of operators; there is no way to change the list of operators dynamically without going back to the assoc function and the list directly.
We will treat an operator as infix exactly when it appears in the list of infixes with a precedence. We can do:
- fun is_infix s =
(get_precedence s; true)
handle _ => false;
because get_precedence fails on non-infixes. An alternative coding is to use an auxiliary function can which finds out whether the application of its first argument to its second succeeds:
- fun can f x =
(f x; true)
handle _ => false;
> val can = fn ('a -> 'b) -> 'a -> bool
- val is_infix = can get_precedence;
> val is_infix = fn : string -> bool
The printer consists of two mutually recursive functions.
The function string_of_term takes two arguments. The first indicates the `currently active precedence', and the second is the term.
For example, in printing the right-hand argument of x * (y + z), the currently active precedence is the precedence of *. If the function prints an application of an infix operator (here +), it puts brackets around it unless its own precedence is higher.
We have a second mutually recursive function to print a list of terms separated by commas. This is for the argument list of non-unary and non-infix functions of the form f(x1,…,xn).
fun string_of_term prec =
fn (Var s) => s
| (Const c) => c
| (Fn(f,args)) =>
if length args = 2 andalso is_infix f
then let val prec' = get_precedence f
val s1 = string_of_term prec' (hd args)
val s2 = string_of_term prec' (hd (tl args))
val ss = s1^" "^f^" "^s2
in if prec' <= prec then "("^ss^")"
else ss
end
else f^"("^(string_of_terms args)^")"
and string_of_terms [] = ""
| string_of_terms [t] = string_of_term 0 t
| string_of_terms (x::xs) = (string_of_term 0 x)^","^
(string_of_terms xs);
Moscow ML has special facilities from its Caml Light heritage for installing user-defined printers in the toplevel read-eval-print loop.
Once our printer is installed, anything of type term will be printed using it.
- load "PP";
> val it = () : unit
- fun print_term pps s =
let open PP
in begin_block pps INCONSISTENT 0;
add_string pps
("'"^(string_of_term 0 s)^"'");
end_block pps
end;
> val print_term = fn : ppstream -> term -> unit
- installPP print_term;
> val it = () : unit
- t;
> val it =
Fn("-",
[Fn("/",
[Fn("sin", [Fn("+", [Var "x", Var "y"])]),
Fn("cos", [Fn("-", [Var "x", Fn("exp", [Var "y"])])])]),
Fn("ln", [Fn("+", [Const "1", Var "x"])])]) : term
- installPP print_term;
> val it = () : unit
- t;
> val it = 'sin(x + y) / cos(x - exp(y)) - ln(1 + x)' : term
- val x = t;
> val x = 'sin(x + y) / cos(x - exp(y)) - ln(1 + x)' : term
- [x,t,x];
> val it =
['sin(x + y) / cos(x - exp(y)) - ln(1 + x)',
'sin(x + y) / cos(x - exp(y)) - ln(1 + x)',
'sin(x + y) / cos(x - exp(y)) - ln(1 + x)'] : term list
There is a well-known method (taught in schools) to differentiate complicated expressions.
This translates very easily into a recursive computer algorithm.
fun differentiate x =
fn Var y => if y = x then Const "1" else Const "0"
| Const c => Const "0"
| Fn("-",[t]) => Fn("-",[differentiate x t])
| Fn("+",[t1,t2]) => Fn("+",[differentiate x t1,
differentiate x t2])
| Fn("*",[t1,t2]) =>
Fn("+",[Fn("*",[differentiate x t1, t2]),
Fn("*",[t1, differentiate x t2])])
| Fn("inv",[t]) => chain x t
(Fn("-",[Fn("inv",[Fn("^",[t,Const "2"])])]))
| Fn("^",[t,n]) => chain x t
(Fn("*",[n, Fn("^",[t, Fn("-",[n, Const "1"])])]))
| (tm as Fn("exp",[t])) => chain x t tm
| Fn("ln",[t]) => chain x t (Fn("inv",[t]))
| Fn("sin",[t]) => chain x t (Fn("cos",[t]))
| Fn("cos",[t]) => chain x t (Fn("-",[Fn("sin",[t])]))
| Fn("/",[t1,t2]) => differentiate x
(Fn("*",[t1, Fn("inv",[t2])]))
| Fn("tan",[t]) => differentiate x
(Fn("/",[Fn("sin",[t]), Fn("cos",[t])]))
and chain x t u = Fn("*",[differentiate x t, u]);
Let's try a few examples:
- val t1 = Fn("sin",[Fn("*",[Const "2",Var "x"])]);
> val t1 = 'sin(2 * x)' : term
- differentiate "x" t1;
> val it = '(0 * x + 2 * 1) * cos(2 * x)' : term
- val t2 = Fn("tan",[Var "x"]);
> val t2 = 'tan(x)' : term
- differentiate "x" t2;
> val it =
'(1 * cos(x)) * inv(cos(x)) +
sin(x) * ((1 * -(sin(x))) *
-(inv(cos(x) ^ 2)))' : term
- differentiate "y" t2;
> val it =
'(0 * cos(x)) * inv(cos(x)) +
sin(x) * ((0 * -(sin(x))) *
-(inv(cos(x) ^ 2)))' : term
It seems to work okay, but it isn't making certain obvious simplifications, like 0 * x = 0. These arise partly because we apply the recursive rules like the chain rule even in trivial cases. We'll write a simplifier:
- val simp =
fn (Fn("+",[Const "0", t])) => t
| (Fn("+",[t, Const "0"])) => t
| (Fn("-",[t, Const "0"])) => t
| (Fn("-",[Const "0", t])) => Fn("-",[t])
| (Fn("+",[t1, Fn("-",[t2])])) => Fn("-",[t1, t2])
| (Fn("*",[Const "0", t])) => Const "0"
| (Fn("*",[t, Const "0"])) => Const "0"
| (Fn("*",[Const "1", t])) => t
| (Fn("*",[t, Const "1"])) => t
| (Fn("*",[Fn("-",[t1]), Fn("-",[t2])])) =>
Fn("*",[t1, t2])
| (Fn("*",[Fn("-",[t1]), t2])) =>
Fn("-",[Fn("*",[t1, t2])])
| (Fn("*",[t1, Fn("-",[t2])])) =>
Fn("-",[Fn("*",[t1, t2])])
| (Fn("-",[Fn("-",[t])])) => t
| t => t;
We need to apply simp recursively, bottom-up:
- fun dsimp (Fn(f,args)) =
simp(Fn(f,map dsimp args))
| dsimp t = simp t;
Now we get better results:
- dsimp (differentiate "x" t1);
> val it = '2 * cos(2 * x)' : term
- dsimp (differentiate "x" t2);
> val it = 'cos(x) * inv(cos(x)) +
sin(x) *
(sin(x) * inv(cos(x) ^ 2))' : term
- dsimp (differentiate "y" t2);
> val it = '0' : term
There are still other simplifications to be made, e.g., cos(x)*inv(cos(x))=1. Good algebraic simplification is a difficult problem.