Grim symbolic computation

January 31, 2020

First blog post of the new decade! I am currently working on symbolic computation in Pygrim, the Python library I'm developing as the backend for Fungrim. It is now possible to use Pygrim for basic expression simplification (or theorem proving). Pygrim is not anywhere near a fully-fledged computer algebra system, but it starting to be useful for one of my initial goals: randomized testing of formulas in Fungrim itself.

The *grim projects

Just to clear up any confusion, here are the projects I will be talking about in this post:

Fungrim (the Mathematical Functions Grimoire) is a website and database of formulas for mathematical functions, currently comprising some 2600 entries.
Pygrim is a Python library for symbolic mathematics, used to build Fungrim.
Fungrim and Pygrim use the Grim formula language to represent mathematical formulas and theorems semantically. It looks like this: Implies(Element(n, ZZ), Equal(Sin(Mul(Pi, n)), 0)) represents the logical formula $\left(n \in \mathbb{Z}\right) \implies \left(\sin\!\left(\pi n\right) = 0\right)$.

Simple symbolics with Pygrim

I will illustrate some symbolic operations using the git version of Pygrim. As a first step, we may import Pygrim (import * imports the classes as well as all predefined symbol names) and construct a simple Grim expression:

>>> from pygrim import *
>>> 
>>> 1 + Pi + 1
Add(Add(1, Pi), 1)
>>> _.latex()
'1 + \\pi + 1'

Grim expressions (instances of the pygrim.Expr class) are preserved verbatim by default. In order to simplify an expression, we will have to create a pygrim.Brain instance (this is essentially a symbolic computation engine) and call the brain's simple method:

>>> brain = Brain()
>>> brain.simple(1 + Pi + 1)
Add(2, Pi)
>>> brain.simple(Pi**2 / Pi**3 - 1 / Pi)
0

The simple method returns an expression that is mathematically equivalent to the input expression, simplified if possible. In particular, logical formulas may be simplified to True or False. The brain can prove simple theorems about real and complex numbers, using a combination of symbolic simplification, domain inference and numerical evaluation (it proves inequalities between explicit complex numbers using Arb numerical evaluation):

>>> brain.simple(Element(Pow(10, Pow(10, Pow(10, 10000000000000000000000000000000))), ZZ))   # a BIG integer
True_
>>> brain.simple(Element(Pi, SetMinus(RR, QQ)))   # pi is real but not rational
True_
>>> brain.simple(Equal(Pi, Log(640320**3 + 744) / Sqrt(163)))  # a good approximation but not equal to pi!
False_
>>> brain.simple(Less(Abs(Pi - Log(640320**3 + 744) / Sqrt(163)), Pow(10,-30)))
True_
>>> brain.simple(Element(Exp(3*ConstI*Pi/7), AlgebraicNumbers))
True_
>>> brain.simple(Equal(Abs(Exp(3*ConstI*Pi/7)), 1))
True_

Now let us try a simplification with a variable instead of a constant:

>>> brain.simple(1 + x + 1)
Add(Add(1, x), 1)

Nothing happens! This brings us to a fundamental principle: Pygrim will only perform simplification that are correct for all possible values of the variables. A priori we know nothing about x in this example; we do not even know that addition is associative. To enable familiar simplifications, we have to create a brain with knowledge about the variables we are going to use; for example, we can specify that x is a complex number:

>>> brain = Brain(variables=[x], assumptions=Element(x, CC))
>>> brain.simple(1 + x + 1)
Add(2, x)
>>> brain.simple(2*x - 3*x + x)
0

The expression $x / x$ does not simplify to 1, because the identity $x / x = 1$ would not be valid when $x = 0$.

>>> brain.simple(x / x)
Div(x, x)

To fix this, we can exclude $x = 0$ from the domain:

>>> brain = Brain(variables=[x], assumptions=And(Element(x, CC), NotEqual(x, 0)))
>>> brain.simple(x / x)
1

Another fun example: showing that a whole class of numbers described by a parameter are transcendental:

>>> brain = Brain([x], Element(x, SetMinus(QQ, Set(0))))
>>> brain.simple(Element(Exp(x), AlgebraicNumbers))
False_

The examples I'm showing are somewhat cherry-picked. The symbolic engine is still rather limited; for example, it does not yet have any builtin facilities for polynomial arithmetic and will thus fail to prove the following identity:

>>> brain.simple(Equal(x*(x+1), x**2+x))
Equal(Mul(x, Add(1, x)), Add(x, Pow(x, 2)))

It also has no ability to perform complex logical simplifications (e.g. SAT solving), reduce systems of inequalities, etc. But there is at least a foundation to build on now, and we will see how far it develops in the future.

Testing formulas

One of the uses for the symbolic simplification code in Pygrim is to test formulas. As an example, let us construct the symbolic equation $\sqrt{x^2} = x$:

>>> formula = Equal(Sqrt(x**2), x)
>>> formula
Equal(Sqrt(Pow(x, 2)), x)

Is it a correct formula? We can run the test method, which works by automatically generating "random" specific values for the free variables, plugging those values into the formula, and simplifying. If a logical formula simplifies to False, we have found a counterexample. If it simplifies to True, we have proved that the formula is valid for that particular choice of specific values (not a proof that it is true for all values, of course). If it cannot be simplified to True or False, the truth value is unknown; the symbolic algorithms at least failed to prove that the equation is False, so we should gain some confidence that the formula is valid.

When calling the method, we need to list the free variables in the formula (here x) and specify assumptions (domain of the variables), for example $x \in \mathbb{R}$:

>>> formula.test([x], Element(x, RR))
{x: 0}    ...  True
{x: Div(1, 2)}    ...  True
{x: Sqrt(2)}    ...  True
{x: Pi}    ...  True
{x: 1}    ...  True
{x: Neg(Div(1, 2))}    ...  False
Traceback (most recent call last):
  ...
ValueError

The formula is valid for the values $x = 0, 1/2, \sqrt{2}, \pi, 1$ but the test code found that the negative value $x = -1/2$ is a counterexample. We can add the condition $x \ge 0$ to the assumptions to get a correct formula:

formula.test([x], And(Element(x, RR), GreaterEqual(x, 0)))
...
Passed 69 instances (67 True, 2 Unknown, 0 False)

More generally, this formula is valid for complex numbers, provided correct conditions for the real and imaginary parts:

formula.test([x], And(Element(x, CC), Or(Greater(Re(x), 0), And(Equal(Re(x), 0), Greater(Im(x), 0)))))
...
Passed 91 instances (66 True, 25 Unknown, 0 False)

The test code can handle several variables and more complex assumptions. For example, $((x+y)^n)^{1/n} = x+y$ will be valid for complex $x, y$ and integers $n \ge 1$ if $|\arg(x+y)| < \pi / |n|$:

>>> # incorrect assumptions
>>> Equal(((x+y)**n)**(1/n), x+y).test([x,y,n], And(Element(x,CC), Element(y,CC), Element(n,ZZGreaterEqual(1))))
{x: 0, y: Div(1, 2), n: 1}    ...  True
{x: 0, y: Sqrt(2), n: 1}    ...  True
{x: Div(1, 2), y: Div(1, 2), n: 1}    ...  True
{x: 0, y: ConstI, n: 2}    ...  Unknown
{x: 0, y: ConstI, n: 1}    ...  True
{x: Div(1, 2), y: Sqrt(2), n: 1}    ...  True
{x: Sqrt(2), y: Div(1, 2), n: 1}    ...  True
{x: 0, y: Pi, n: 2}    ...  Unknown
{x: 0, y: Pi, n: 1}    ...  True
{x: Div(1, 2), y: ConstI, n: 2}    ...  Unknown
{x: Div(1, 2), y: ConstI, n: 1}    ...  True
{x: Sqrt(2), y: Sqrt(2), n: 1}    ...  True
{x: ConstI, y: Div(1, 2), n: 1}    ...  True
{x: 0, y: Mul(Pi, ConstI), n: 3}    ...  False
Traceback (most recent call last):
  ...
ValueError

>>> # correct assumptions
>>> Equal(((x+y)**n)**(1/n), x+y).test([x,y,n], And(Element(x,CC), Element(y,CC), Element(n,ZZGreaterEqual(1)), Less(Abs(Arg(x+y)),Pi/Abs(n))))
...
Passed 100 instances (58 True, 42 Unknown, 0 False)

By default, the test code tries up to 100 iterations (assignments of specific values satisfying the assumptions of the formula). It is reassuring that it finds a counterexample quickly. The proportion of True/Unknown results gives some indication of how well the internal symbolic simplification code performs; ideally, we would get True all the time.

Testing Fungrim

It is now very simple to test all the entries in Fungrim: pygrim.test_fungrim() runs the test described above for each entry in Fungrim. This currently takes nearly two hours (the code is far from optimized). The verdict? When I did so for the first time yesterday, the test code found 24 entries that had errors out of the 2618 entries in Fungrim (I have now fixed all of them). The errors can be categorized as follows:

Semantically nonsensical expression (Equal() called with a single argument because of misplaced parentheses, etc.): 9
Wrong list of variables in metadata: 5
Mathematically wrong assumptions: 6
Mathematically wrong formula: 4

The first two kinds of errors are not that interesting; detecting them is essentially a matter of type-checking. What about the 9 entries with clearly wrong mathematical content? I show the corrected versions below and explain what was wrong.

Entry 30bd5b:

$$W_{0}\!\left(x \log(x)\right) = \log(x)$$

Assumptions: $x \in \left[\frac{1}{e}, \infty\right)$

The error? The assumptions were incorrectly given as $x \in \left[-\frac{1}{e}, \infty\right)$. I guess when I entered this formula in Fungrim, I just thought about the domain of the Lambert W function and did not stop to think through the fact that the argument is being modified!

Entry 52ea5f:

$$\operatorname{Li}_{s}\!\left(z\right) = \frac{\Gamma\!\left(1 - s\right)}{{\left(2 \pi\right)}^{1 - s}} \left({i}^{1 - s} \zeta\!\left(1 - s, \frac{1}{2} + \frac{\log\!\left(-z\right)}{2 \pi i}\right) + {i}^{s - 1} \zeta\!\left(1 - s, \frac{1}{2} - \frac{\log\!\left(-z\right)}{2 \pi i}\right)\right)$$

Assumptions: $s \in \mathbb{C} \,\mathbin{\operatorname{and}}\, z \in \mathbb{C} \,\mathbin{\operatorname{and}}\, z \notin \left\{0, 1\right\} \,\mathbin{\operatorname{and}}\, s \notin \mathbb{Z}_{\ge 0}$

The error? The argument to the second Hurwitz zeta function was given as $s$ instead of $1-s$ (a transcription error).

Entry d10873:

$${\left(-1\right)}^{n} B_{2 n + 2} > 0$$

Assumptions: $n \in \mathbb{Z}_{\ge 0}$

The error? The prefactor was $(-1)^{n+1}$ instead of $(-1)^n$, so the sign was exactly wrong. Doh!

Entries b25089, 504717 and 90ac58 (these are very similar; only showing b25089 here): $$\,{}_2{\textbf F}_1\!\left(a, b, c, z\right) = {\left(1 - z\right)}^{-a} \,{}_2{\textbf F}_1\!\left(a, c - b, c, \frac{z}{z - 1}\right)$$

Assumptions: $a \in \mathbb{C} \,\mathbin{\operatorname{and}}\, b \in \mathbb{C} \,\mathbin{\operatorname{and}}\, c \in \mathbb{C} \,\mathbin{\operatorname{and}}\, z \in \mathbb{C} \,\mathbin{\operatorname{and}}\, z \notin \left[1, \infty\right)$

The error? The last assumption (excluding the branch cut) was missing. I thought this was correct because I had used it unconditionally in Arb (it happens to be correct when used in Arb because of additional convergence requirements).

Entry a68e0e:

$$\operatorname{sgn}(k) \operatorname{Im}\!\left(W_{k}\!\left(z\right)\right) \in \left(\left(2 \left|k\right| - 2\right) \pi, \left(2 \left|k\right| + 1\right) \pi\right)$$

Assumptions: $z \in \mathbb{C} \setminus \left\{0\right\} \,\mathbin{\operatorname{and}}\, k \in \mathbb{Z} \setminus \left\{-1, 0, 1\right\}$

The error? The factor $\operatorname{sgn}(k)$ was originally in the endpoints of the interval on the right, so the endpoints were in the wrong order for $k < 0$ and described an empty interval.

Entry d1a0ec:

$$\log G\!\left(1 - x\right) = \log G\!\left(1 + x\right) + x \log\!\left(\frac{\left|\sin\!\left(\pi x\right)\right|}{\pi}\right) + \frac{1}{2 \pi} \operatorname{Im}\!\left(\operatorname{Li}_{2}\!\left({e}^{2 \pi i x}\right)\right) + \operatorname{sgn}(x) n \left(n + 1\right) \frac{\pi i}{2}\; \text{ where } n = \left\lfloor x \right\rfloor$$

Assumptions: $x \in \mathbb{R} \,\mathbin{\operatorname{and}}\, x \notin \mathbb{Z}$

The error? This formula had two transcription errors: a missing factor $x$ and a multiplication instead of an addition. Sloppy!

Entry e722ca:

$$\operatorname{Im}\!\left(\sqrt{z}\right) = \operatorname{sgn}\!\left(\operatorname{Im}(z)\right) \sqrt{\frac{\left|z\right| - \operatorname{Re}(z)}{2}}$$

Assumptions: $z \in \mathbb{C} \setminus \left(-\infty, 0\right)$

The error? The assumptions did not exclude negative real numbers. A dumb oversight.

Entry 185efc:

$$\sqrt{\frac{z}{c - z}} = \sqrt{z} \sqrt{\frac{1}{c - z}}$$

Assumptions: $z \in \mathbb{R} \,\mathbin{\operatorname{and}}\, c \in \left[0, \infty\right) \,\mathbin{\operatorname{and}}\, c - z \ne 0$

The error? The original assumptions allowed $z \in \mathbb{C}$, but this is not actually valid. The current assumptions could be relaxed to allow complex numbers with the right conditions; I have not had time to look at this.

It was fun to see that the test code really managed to find errors. An error rate around 1% (24 / 2618) for the formulas I have entered manually into Fungrim sounds plausible. Granted, the symbolic engine has limited simplification powers and 100 inputs isn't enough to probe the parameter space effectively for some formulas, so it's safe to assume that at least as many errors still remain undetected! I hope I will manage to improve the symbolic simplification code so that it finds more problems!

Automatic simplification with the Fungrim formula database

What is the purpose of the Brain class, as opposed to just having a function for symbolic simplification? From the implementation point of view, a class is a convenient way to organize all the related subroutines for simplification. Some form of context object is also needed to store information about variables. Moreover, you could subclass Brain and overloads its internal methods to define additional simplification heuristics.

The code organization I have in mind for Pygrim is something like this:

Symbolic expressions are very simple objects (just used as a data structure). The Brain implements a set of heuristics for simplifying or rewriting expressions; this may include delegating specific subproblems to dedicated implementations (for example, it currently uses Python-FLINT for rational arithmetic and Arb numerical evaluation). The third potential component is the library of formulas in Fungrim; the mathematical knowledge encoded in this database could be used for simplifications.

The experimental FungrimBrain class does precisely this. Fundamentally, for every subexpression $a$ in the given expression, it checks if there is an entry of the form $a = b$ in Fungrim where $b$ is simpler than $a$ (by some metric), and in that case it replaces $a$ with $b$. Here are some examples:

>>> brain = FungrimBrain()
>>> brain.simple(RiemannZeta(2))
Div(Pow(Pi, 2), 6)
>>> brain.simple(RiemannZeta(2) / Pi**2)  # combining the above with arithmetic
Div(1, 6)
>>> brain.simple(Element(RiemannZeta(3), QQ))
False_
>>> brain.simple(Zeros(Gamma(z), ForElement(z, CC)))
Set()
>>> brain.simple(Integral(1/Sinc(x), For(x, 0, Pi/2)))  # an integral representation for Catalan's constant
Mul(2, ConstCatalan)

Fungrim entries involving variables are applied using pattern matching. For example, Fungrim has an entry for the identity $\sin(z)^2 + \cos(z)^2 = 1$ (for any complex number $z$), as well as an entry for $\operatorname{erf}(z) + \operatorname{erfc}(z) = 1$:

>>> brain.simple(Sin(1)**2 + Cos(1)**2)
1
>>> brain.simple(Erf(Sin(1)**2 + Cos(1)**2) + Erfc(1))
1

Here is an even more interesting example: entry 463077 in Fungrim gives the integration rule $$\int_{z}^{\infty} \frac{1}{{\left(a x + b\right)}^{c}} \, dx = \frac{1}{a \left(c - 1\right) {\left(a z + b\right)}^{c - 1}}$$ with the stated assumptions $a \in \mathbb{R} \,\mathbin{\operatorname{and}}\, b \in \mathbb{R} \,\mathbin{\operatorname{and}}\, c \in \mathbb{R} \,\mathbin{\operatorname{and}}\, z \in \mathbb{R} \,\mathbin{\operatorname{and}}\, a > 0 \,\mathbin{\operatorname{and}}\, a z + b > 0 \,\mathbin{\operatorname{and}}\, c > 1$. Sure enough, it works:

>>> brain.simple(2 * Integral(1/(2*x+3)**Div(3,2), For(x, 1, Infinity)))
Mul(2, Pow(5, Div(-1, 2)))

Note that assumptions are handled correctly: the rule is not applied when the integral does not converge:

>>> brain.simple(2 * Integral(1/(2*x+3)**Div(1,2), For(x, 1, Infinity)))
Mul(2, Integral(Div(1, Pow(Add(Mul(2, x), 3), Div(1, 2))), For(x, 1, Infinity)))

It is also possible to apply Fungrim formulas one by one. For example, entry ad6c1c in Fungrim is the trigonometric identity $$\sin(a) \sin(b) = \tfrac{1}{2}\left(\cos\left(a - b\right) - \cos\left(a + b\right)\right)$$ which transforms a product of sine functions to a sum of cosines. This is arguably not a simplification (so it does not get applied automatically), but we can invoke this rewrite rule manually:

>>> brain.rewrite_fungrim(Sin(2)*Sin(Sqrt(2)), "ad6c1c")
Div(Sub(Cos(Sub(2, Sqrt(2))), Cos(Add(2, Sqrt(2)))), 2)

If any of this seems impressive, I'm afraid I have to temper your expectations!

The current pattern matching code is naive and requires an exact match; for example, it will not match $1 + x$ with $x + 1$. The Fungrim brain will not be useful for any real applications until the pattern matching algorithm can recognize associativity and commutativity and solve for (at least) linear transformations of the matching variables.
The current code is also rather buggy; it does not handle bound variables correctly and sometimes gives bogus matches when the same symbol is used by the user and in a Fungrim entry. (Not too hard to fix.)
A useful global simplification algorithm needs some heuristics for how to combine sequences of rewrite rules. It's relatively straightforward to apply directed rules for a specific purpose (for example, converting formulas involving integrals to closed-form solutions without integrals), but achieving general formula simplification based purely (or mostly) on declarative knowledge is a very, very hard problem which I don't expect to solve!

Considering these problems, this is just a proof of concept, and a lot needs to be done before Pygrim can compete with other computer algebra systems.

Performance

Right now, I'm only concerned about making the rewriting and simplification algorithms powerful and correct. I haven't looked at performance optimization at all. Indeed, there are lots of obvious things that can be optimized. The symbolic computation code is still quite stupid; it does a lot of redundant work (recomputes the same thing over and over) and sometimes uses algorithms with poor complexity. On top of that, a lot could be gained by implementing core symbolic operations in something other than Python. Such optimizations are a goal for the more distant future.

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor