Oct 17, 2024
I've been using Mathematica to help with learning math (and using math to help learn Mathematica). Too often my Mathematica code has surprising behavior, which means I don't really understand the language or how it's evaluated. So I thought I'd do two things to get better.
First, write a bunch of Mathematica rules to implement a toy symbolic differentiation operator. Second, write a toy Mathematica interpreter in Typescript that's good enough to correctly run my custom differentiation code.
Writing a toy differentiator turns out to be shockingly easy. It's a near verbatim transcription of differentiation rules from any calculus textbook:
D[_?NumberQ, x_Symbol] = 0;
D[x_, x_Symbol] = 1;
D[Times[expr1_, expr2_], x_Symbol] =
D[expr1, x] expr2 + D[expr2, x] expr1;
D[Plus[expr1_, expr2_], x_Symbol] = D[expr1, x] + D[expr2, x];
D[Sin[x_], x_Symbol] = Cos[x];
D[Cos[x_], x_Symbol] = -Sin[x];
D[f_Symbol[expr_], x_Symbol] :=
(D[f[x], x] /. x -> expr) * D[expr, x];
D[Power[expr_, p_Integer], x_Symbol] := p expr^(p - 1) * D[expr, x];
My first implementation had three more rules (one for each of \(dx^k/dx\), \((c\cdot f)'(x)\), and \((1/f)'(x)\)), which I later realized you don't need. These are shortcuts for human differentiators, but they're automagically covered by the multiplication and exponentiation rules above.
Some problems I encountered while writing these rules: infinite
recursion, rules not matching, and rules matching but evaluating to a
surprising result. The hobby edition of Mathematica has some tools
for debugging these problems, but not much. Ultimately I fixed bad
rules by staring at the code and thinking really hard. I found Trace
impossible to read, and TracePrint (which is supposed to be better)
not much better. Also MatchQ is good, but somehow not as useful for
debugging as I would have liked.
I first implemented basic parsing of integers, symbols, forms, and
arithmetic operators using ts-parsec (which I wrote for this
project). In Mathematica a 2 b 3 evaluates to Times[6, a, b]
because Times has a Flat attribute. To get this behavior I
implemented attributes next-- Flat, and also HoldFirst,
HoldRest, and HoldAll which I'd eventually need. I also exposed
Attributes, SetAttributes, and ClearAttributes. These all accept
(and Attributes returns) lists, so I added those too. All this was
easy enough.
I wanted to implement assignment next so I could say a = 1. In
Mathematica even something this simple is implemented by adding
RuleDelayed[HoldPattern[a], 1] to OwnValues.1 So the next step
was to build the skeleton of a term rewriting system.2 I first
implemented a version of MatchQ that does a deep equal, extended it
to handle Blank[], extended it again to handle Blank[...], and
extended it again to handle Pattern[foo, Blank[...]]. The version
exposed to the interpreter returns True or False, but under the
hood it also returns an environment with matched variables. I built on
that next to implement term rewriting.
A really simple rewrite rule is f[1] /. f[x_] :> x + 1. In
Mathematica this parses to
ReplaceAll[f[1], RuleDelayed[f[Pattern[x, Blank[]]], Plus[x, 1]]]
With my MatchQ implementation there was now enough machinery to get
this working. I added operators /. and :> to the grammar,
implemented Replace, and built ReplaceAll on top. I tested a bunch
of rewrite rules and they all worked! From here it was also easy to
add -> and //. which I did next.
I had enough machinery to implement assignment. I added Set and
SetDelayed, and modified evaluation to apply the rules stored in
OwnValues and DownValues. This let me assign values to symbols and
define functions! I could now run code like this:
fib[1] := 1
fib[2] := 1
fib[x_] := fib[x-2] + fib[x-1]
One caveat is that Mathematica inserts rules into value lists in order
of specificity, so specific rules are tried before general rules. I
initially added rules in order of definition to get assignment
working, but then went back and added a simple specificity ordering
function3.
EDIT: I wrote specificity ordering code late at night, and just realized it was completely broken. I removed it; rules are now processed in the order they're entered. But everything else still works!
Finally, I added support for PatternTest/? and NumberQ. These
were all the features needed to do differentiation!
ts-wolframI ran D in ts-wolfram on the following examples, and cross-checked
with the results in Mathematica:
D[1, x]
D[x, x]
D[x^5, x]
D[3 x^2, x]
D[(x + 1) (x + 2), x]
D[x^2 + x^3, x]
D[Cos[x], x]
D[x^3/(x^2 + 1), x]
D[Cos[Cos[x]], x]
D[Cos[Cos[Cos[x]]], x]
D[Cos[x^2 + 1], x]
D[(x + 1)^2, x]
Somewhat miraculously, ts-wolfram got correct results on every test!
Since I didn't add any manipulation to simplify expressions, the
results weren't exactly the same. For example:
D[Cos[Cos[x]], x]
(* ts-wolfram outputs *)
Times[Minus[Sin[Cos[x]]], Minus[Sin[x]]]
(* With `FullForm`, Mathematica outputs *)
Times[Sin[x], Sin[Cos[x]]]
This would be easy to fix by adding the rule Times[Minus[x_],
Minus[y_]]=x y, but (a) Times is implemented in typescript, and the
interpreter doesn't currently support mixing kernel and userspace
rules, and (b) extending the system to simplify algebraic expressions
feels like a different undertaking.
Overall, this has been a really fun and instructive project. I built it in four days hacking on it after work, and learned a great deal about Mathematica internals. Of course this is still only scratching the surface, but now I feel a lot less lost when Mathematica doesn't behave the way I expect. I'm very happy with this outcome!
You can check this by running a = 1; OwnValues[a]. ↩
If you don't already understand how this works, reverse engineering it is tricky, even with Mathematica's great docs. For examples why, see Stack Overflow questions here, here, and here. ↩
The exact details for how Mathematica handles pattern specificity are not so easy to find. I didn't try too hard to reverse engineer it; I just did what seemed reasonable and moved on. ↩