Try it on GitHub: https://github.com/coffeemug/ts-wolfram

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.

Toy differentiation

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.

Toy interpreter

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 (You can check this by running a = 1; OwnValues[a]). So the next step was to build the skeleton of a term rewriting system (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.). 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 function (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.).

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!

Differentiating in ts-wolfram

I 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!