fredrikj.net / blog /

# Things I would like to see in a computer algebra system

*April 20, 2022*

If I were to design a computer algebra system (CAS) from scratch today, I would try to achieve the following goals. Several ideas have already been implemented in one system or another, but certainly not all at once. Some ideas could be adopted by existing CASes; others would require fundamentally different designs or new programming languages. Some ideas are possibly unrealistic or naive.

Contents

## Free and open source

This immediately rules out using or building on certain non-free systems. No matter how many other boxes they check.

## Inert expressions

I would like a CAS to distinguish clearly between syntactical symbolic expressions and values. Most CASes perform some kind of automatic rewriting or "canonicalization": for example, if I input the expression $2(b-a)$, the CAS might return $-2a + 2b$. This is fine if I want an expanded normal form of a value in $\mathbb{Q}[a,b]$, but sometimes I really want to preserve the exact expression or work with some canonical form other than the one the CAS author thought of. Working around builtin canonicalization is a nuisance! A far better solution is to have syntactical symbolic expressions that are completely inert by default and require explicit commands for rewriting them. (This is what I did for the Fungrim backend.) This is simpler to implement too, so it's a double win.

Here, I really consider symbolic expressions separate from data structures for explicit computation in particular algebraic domains like $\mathbb{Q}[a,b]$, where canonicalization indeed is crucial for performance.

## Based on math (even in analysis!)

I would like to have a CAS that is based on math. Real numbers should actually be (isomorphic to) Cauchy sequences over $\mathbb{Q}$, rings should actually satisfy the ring axioms, sets should actually be able to hold infinitely many elements (but not themselves), holomorphic functions should not occasionally have branch cuts on their domain, "solving an equation" should give me exactly the set of solutions (neither more nor less) over some well-defined domain, $x = 0$ should be true if and only if $x$ is zero...

Pretty much the only systems that more or less get math right are the mathematical libraries for proof assistants (Lean, Coq, etc.). The mathematical structures and interfaces in Lean's MathLib are actually pretty close to what I'd like to see in a CAS.

The algebraically oriented CASes like Sage, Magma and GAP are pretty good as long as you stick to pure algebra involving finite objects. They generally start to cheat when it comes to mathematical analysis. A CAS worth its name ought to be able to represent structures like transcendental number fields, rings of holomorphic functions, differential fields, transseries, and so on. Unfortunately, to the extent current-generation CASes do mathematical analysis at all, it is either numerically or using ad-hoc "symbolics" without solid algebraic foundations. (More on this in a future blog post, if I have time.)

## Based on typed math

The programming language of a CAS should have a type system suitable for doing mathematics. This presumably means having fully dependent types in some form. You can implement mathematical structures on top of a poor type system (Sage does this, for example), but proper language support would be infinitely better. I have no idea what the ideal CAS type system might look like, and I believe this is a hard research problem; the following are just some scattered thoughts on the matter.

I will note that "type system" does not necessarily mean "static type system", nor "decidable", although statically decidable properties are useful (especially for reusable library code). A core ability of a CAS is to handle untyped expressions and map them to different algebraic structures, often requiring type information to be computed at runtime, and sometimes depending on properties that are not guaranteed to be computable. (Make everything decidable and you get a formal theorem prover, but a practical CAS also requires the ability to work with unproven math!)

There are tradeoffs involved in how much information you want to encode in types. On one extreme, you have a system where the type signature of a function completely specifies the function. On the other extreme (most mainstream programming languages), the type system does not even guarantee that functions are total functions. I do think that the type system in a CAS should be expressive enough to encode arbitrary mathematical relations as invariants/preconditions/postconditions (perhaps as "type restrictions"), even if you don't end up using this all the time in practice.

When it comes to actually modelling mathematical structures with types, it is important to understand that mathematics requires the ability to consider the same abstract value as an element of different algebraic structures. The number 6 can be an element of the ring of integers $\mathbb{Z}$, but it can also be an element of the semiring $\mathbb{N}$, or the set $\{n \in \mathbb{Z} : n \text{ composite}\}$, the rational field $\mathbb{Q}$, the complex field $\mathbb{C}$, the polynomial ring $\mathbb{Z}[[x]]$, and so on. You might want to distinguish between "$\mathbb{Z}$" as a type representing the integers, as a set containing all the values of this type, or as a subset or subtype of something else (like the set of rational numbers or the type of rational numbers).

The most useful type or domain in a given situation typically depends on what kind
of closure properties you need rather than on the values you start with.
This means, e.g. rings tend to be more useful than semirings.
For example, 6 has an additive inverse in $\mathbb{Z}$ but not in $\mathbb{N}$,
so if you ever want to do a subtraction,
chances are that you want to represent 6 as an element of
$\mathbb{Z}$ rather than as an element of $\mathbb{N}$
(this is the reason why `int` often makes more sense than `unsigned int` in C, even when
you have an unsigned value).
In this situation, I'd still like to be able to attach restrictions like "$n \ge 1$" in some form.

Classically "analysis-oriented" CASes like Mathematica just implicitly work over a single mathematical universe containing the field of the complex numbers and also different objects like formal indeterminates and infinities. This maximizes closure properties (every function can be applied to everything!) and avoids a lot of hairy type conversion issues, but gives very weak invariants (you cannot assume properties like associativity... though the system often does so anyway, as a practical necessity, resulting in bugs).

Modern "algebra-oriented" CASes like Sage give elements "parents" corresponding explicit algebraic structures. This has pros and cons, but on balance I think it is the better approach. The main danger is that the user gets bogged down in types, but this is perhaps a language/interface issue that can be solved. You certainly need a robust system for automatic coercions; Sage has the right idea here with its coercion system based on category theory, even if the implementation isn't perfect.

I will end this section with a more concrete example.
Consider a function that divides two integers: `a / b -> c`. What is its type
signature? Here are some possibilities:

`(Integer, Integer) -> Integer`. Inexact division floors the result, and division by zero returns a dummy value like 0. (This is the preferred behavior in proof assistants.)`(Integer, Integer) -> Union(Integer, DivisionByZeroError)`. As above, but division by zero signals an error.`(Integer, Integer) -> ProjectiveExtendedInteger`or`(Integer, Integer) -> Union(Integer, Infinity)`. As above, but division by zero returns an infinity in a structure corresponding to $\mathbb{Z} \cup \{ \infty \}$. (What about $0 / 0$ though?)`(Integer, Integer) -> Union(Integer, InexactError, DivisionByZeroError)`. Division by zero as well as inexact division signal an error.`(Integer, Integer) -> Union(Integer, Rational, DivisionByZeroError)`. Returns an`Integer`when the division is exact,`Rational`when it is a fraction, and signals an error on division by zero.`(Integer, Integer) -> Union(Rational, DivisionByZeroError)`. Always returns a fraction except when dividing by zero.`(Integer, Integer) -> Rational`. Always returns a fraction, with a dummy value like 0 on division by zero.`(Integer, Integer) -> ProjectiveExtendedRational`. Returns a fraction, or an infinity on division by zero.`(Integer, NonZeroInteger) -> Integer`. Floor division; division by zero excluded by type-level contract.`(Integer, NonZeroInteger) -> Rational`. Division producing fraction; division by zero excluded by type-level contract.`(Integer x, NonZeroInteger y [y divides x]) -> Integer`. Inexact division also excluded by type-level contract.`(Integer, Integer) -> Union(Integer, DivisionByZeroError)`. Division presumed to be exact; inexact division returns a garbage value. (This is GMP's`divexact`.)`(Integer, Integer) -> IEEEFloat64`. This one makes no sense whatsoever but it is the way some mainstream programming languages define division.`(Integer, Integer) -> Union(IEEEFloat64, DivisionByZeroError)`. Ditto.

I can think of situations where any of these would be useful.
One possible takeway is that it is
*not sufficient for a CAS to infer the target type of an abstract mathematical function from the
types of the arguments, and certainly not from their values*.
A CAS should, accordingly:

- Make it easy to choose both the domain
*and*the codomain of a mathematical function that can have multiple "concretizations". - Base any automatic coercions on mathematically sensible structure (e.g. preserving ring homomorphisms).
- Provide robust error handling (mismatched "concretizations" of abstract operations are a leading source of bugs).

## Type-integrated symbolics and enclosures

A CAS should be able to deal with objects that are known precisely ("the integer 3") as well as objects that are known partially ("some integer $n \ge 3$") or approximately ("some real number near 3.14159"). My earlier post computing with metavalues discusses this among other issues.

There are several poor ways to deal with this. The Mathematica approach is to make everything a symbolic expression, where symbolic expressions can return unevaluated or partially unevaluated. For example, if you attempt to evaluate $x = 1 + 1$ where nothing is known about $x$, you simply get the expression $x = 2$ back. This is an elegant solution, but you don't get the benefits of a proper type system as discussed above. Another issue is that symbolic expressions easily blow up in size when fed through algorithms.

The Sage approach is to have a nice type system designed for precise objects but then throw inexact objects into "inexact rings" (which don't really behave like rings) and all symbolics into a Mathematica-like "symbolic ring" (which is certainly not a mathematical ring). (You then add a buggy "assumptions system" as a band-aid to make the symbolic ring vaguely useful.)

I believe a better approach is to have a type system with a clean separation between "mathematical types" and "implementation types/data" with uniform support for alternative representations (canonical, exact, symbolic, partial, inexact) on the implementation level. For example:

- An
`Integer`can hold a binary blob representing the integer 3 in canonical form, but it can alternatively hold a non-canonical, symbolic or partial description of an integer like "$10^{10^{1000}}$" (unexpanded!), "$n \ge 3$" or "$n \equiv 15 \bmod 60$" or "$n$ is-prime", or "$n$ not equal to this other symbolic integer $m$". Or some combination of such attributes. -
A
`Matrix`instance can hold a canonical array of coefficients, but it might also just hold information about the number of rows and columns (and these numbers may themselves be symbolic integers), the values on the main diagonal, the symbolic description $A = B C$ where $B$ and $C$ are other (potentially symbolic) matrices, etc. A "matrix over $\mathbb{Z}$" can be $[[1, 1], [0, 1]]$, but it can also be $[[1, n], [0, 1]]$ where $n$ is a symbolic integer, an unknown $2 \times 2$ integer matrix known to have determinant 1, or a completely unspecified integer matrix. - An "element of a a ring" can be an element of the specific ring $\mathbb{Z}$, but it can also be an element of a symbolic ring $R$. This symbolic ring $R$ itself can have additional information associated with it, e.g. the property "commutative", or "characteristic $n$" (where $n$ can be symbolic). I should be able to represent the specific element $1 \in R$ of a symbolic ring $R$.

Enclosures in the sense of interval arithmetic
($x \in [3.14, 3.15]$) also fit into this framework.
A key point here is that a "symbolic integer", "symbolic function"
"symbolic matrix", "symbolic element of a ring"
or "real number represented inexactly by an interval"
should not be a second-class citizen;
I should be able to create a symbolic instance or value enclosure of *any* mathematical type,
and these instances should work with the usual
mathematical type system (unlike Sage's symbolics).
I should also be able (indeed, encouraged)
to write code that is agnostic
about the underlying representation or level of specification of mathematical objects,
at least until
I *specifically* need access to the internals.

This design does lead to two difficulties, though I think they can be solved. The first difficulty is that even basic programming elements (logic, collections and control flow) need to fit into the same framework:

- Predicates must be able to return booleans that can be unknown or symbolic.
- A loop or recursion with respect to a symbolic integer $n$ must be evaluated abstractly or symbolically.
- A sorted list or associative array needs to deal correctly with elements whose representations are not comparable or hashable.

`x == 1 + 1`returns False in SymPy; Sage returns the symbolic expression

`x == 2`, but then

`bool(x == 1 + 1)`still returns False. In Calcium, I use triple-valued booleans at the C level (with great care) and throw exceptions in the Python wrapper; these are workable solutions but not perfect.

The second issue is the potential implementation
complexity and unpredictable performance characteristics.
It would be difficult to hardcode data structures
and algorithms for all combinations
of internal representations.
You'd rather want to define a mathematical type in declarative form
as a set of possible properties (e.g. for a `Matrix`: domain,
shape, actual elements, elements on the diagonal, determinant,
trace, eigenvalues,
has-$x$-as-eigenvalue,
is-invertible, is-symmetric, is-real, ...)
with a schema for types and invariants (e.g. shape is a pair of natural numbers; is-symmetric implies is-square, ...),
have a compiler that automatically generates compact binary data representations,
and then provide user-configurable rules for evaluation
(e.g. "reduce all booleans to True/False/Unknown",
or "only track the shapes of matrices" or "enclose real numbers within 53 bits of precision").

Some formal theorem provers have the right ideas here already; it is mostly a question of making things efficient and versatile for calculation.

## Good for inequalities

This is a rather specific feature request. Any CAS will let me compute some terms of Stirling's series for $\Gamma(z)$. How come I can't get a truncation of Stirling's series with an explicit error term? Analysis-oriented CASes are generally good at manipulating equalities and limits, but strangely poor at manipulating inequalities. This is bizarre since inequalities are at least as important as equalities in analysis!

## Large but lean

The most advanced CASes are multi-gigabyte behemoths,
and even more specialized systems
start to weigh a bit heavy when you count all dependencies.
They take time to download, they are complex to compile, and they are difficult to develop.
Even if we get around the build/installation hurdles
(let's just put everything into clouds! or containers! or containers in the clouds!),
we have startup time and interactive latency to deal with.
There are times where I want to run a multi-hour job in a CAS,
but most of the time I just want to do a very quick calculation.
I will not launch `sage` in the terminal (≈ 3 seconds startup)
if I can get away with
`python` or `gp` (≈ 0.01 seconds startup).
I will certainly not reach for a cloud notebook (perhaps more than
a minute to
navigate a website, log in,
wait for some JS framework to load, and spin up a remote VM) if I have
a better option.
Add to this the latency for loading any extra libraries/modules
and compiling code at runtime (Julia currently leaves a lot
to be desired here).

Now, size isn't something you can escape from, at least entirely.
A useful CAS needs to know a great deal of mathematics, how to translate
that mathematics to machine instructions, and how to interact with the
user and other software;
this knowledge will not fit on a floppy disk.
However, I think the user-unfriendliness of large CASes is
mostly of question of design.
Consider Pari/GP for a counterpoint: it is self-contained (except for an optional dependency on GMP),
the entire binary is something like 10 MB, starts instantly (it even loads pretty
quickly in a web browser; try it here),
and it does a *ton* of advanced number theory very efficiently.

Pie-in-the-sky: I would like to see a CAS deliberately designed
around having
a small language/interface kernel
together with a large database of mathematics (types,
functions, constants, theorems, properties, algorithms),
with
*automatic, efficient lazy loading*.

I say "database" rather than "library"
because I think the usual
library/module way to organize software leaves something to be desired.
Computational mathematics is inherently *highly integrated*; a numerical
calculation might invoke some number-theoretic subroutine
that transforms the problem to a calculation over finite fields
aided by a lookup in a database of Conway polynomials.
Even if a computer algebra system is neatly
organized into libraries and modules,
a single calculation may end up using a sparse subset of functionality
of several components. This shouldn't require downloading,
possibly compiling and then loading into memory
all the contents of several huge libraries.
(And as a user, I should certainly not have to do
any of these things manually!)

I propose that the mathematical library of a CAS
should be some kind of semantic database with code-as-data,
fine-grained automatic dependency tracking,
efficient loading/installation/caching,
and tiered interpretation/JIT compilation for dynamic code.
With this design, you could have a kernel that fits in
a small WASM or native executable,
making it possible to launch a local instance
in a terminal or a web browser where you can do
`factor(12345)` or
`plot(sin(x))`
instantly while also having convenient access to
efficient implementations of
modular forms on quaternion algebras (or whatever)
that will be loaded quickly on demand.

I could see computer algebra moving server-side entirely, but it will be a tall order to to solve the usability and performance problems. Even the free-as-in-beer web services offered by the multimillion dollar company Wolfram (Alpha and Cloud) have notable issues. Surprisingly, SymPy Live and Gamma are much more usable; that is, at least until you try to calculate something nontrivial and the backend chokes.

## Integer performance

The core datatype in numerical computing is the machine-precision floating-point number. High-performance computing is virtually synonymous with turning tasks into cache-friendly and vectorized operations on homogeneous arrays of floats. Hardware, compilers, programming languages and numerical libraries are strongly designed around this idea.

The core datatype in computer algebra is the
arbitrary-size integer.
Hardware, compilers, programming languages
are *not at all* designed for performant computation
with this datatype.
Indeed, applications that support arbitrary-size integers
typically just wrap GMP's `mpz_t` (sometimes behind
extra layers of indirection) and implement mathematical operations
naively on top.
This results in massive amounts of overhead for small integers
and often in algorithms with poor asymptotic complexity,
giving users the misleading impression that exact integer arithmetic is
inherently inefficient (often 10 to 1000 times slower than necessary),
which in turn leads to a vicious cycle
of neglect.

Performance-oriented computer algebra libraries like FLINT take pains packing integers and arrays of integers into arrays of bits or words, using modular arithmetic where appropriate to avoid intermediate coefficient swell and facilitate vectorization. This requires a lot of manual coding: you have one standard representation for the coefficients of polynomials, several other internal representations for the coefficients when multiplying or computing GCDs, others for exponents of multivariate polynomials... with lots of fiddling for conversions, resizing, overflow detection, modular arithmetic, precomputed inverses, delayed reduction...

I imagine that a CAS created by an advanced extraterrestrial civilization would be written in a programming language with a native notion of arbitrary-size integer with the ability to generate context-optimized low-level packed representations, bit-level magic and parallelization/vectorization from high-level descriptions of data structures, morphisms and algorithms. Barring alien technology, a good start would be some kind of compiler/language infrastructure where doing all these things manually is easier than what it is right now!

## Good math display

Superficial but important: a CAS should be able to produce publication-quality formulas, graphics and tables without special effort by the user. In practice, the results are often less than impressive. These should not even be particularly hard problems to solve nowadays; they are just not prioritized, or don't get worked on by persons with a sufficiently developed eye for aesthetics or detail (sorry).

## Text-friendly (and human-friendly)

Tools that work well with plain human-readable ASCII input and output are wonderful. I would like every object in a CAS to print to a plain text expression in the syntax of the language which evaluates back to precisely the same object, independent of context. Serialization will, as a consequence, be trivial. Other forms (pretty-printed, condensed, binary) can by all means be supported as well, but a verbose and faithful text representation should be trivial to access (if not the default form of output). This also calls for having a simple, regular syntax so that it is easy to extract and compose subexpressions, and so that the syntax is comprehensible for new users.

(I'm on the fence about Unicode math. It seems to occupy an uncanny valley between easily-typed ASCII math and proper typesetting.)

Concerns about readability and serializability may also influence semantic design decisions. It seems desirable to have strongly functional semantics that avoid implicit context and state. I also think a CAS should have good support for decimal floating-point numbers, without forced conversions to binary. (Lossy binary-decimal conversions are a frequent interface pain point with Arb.)

## Well-named

Since mathematics is highly integrated,
a CAS should provide a unified namespace
with a consistent and convenient naming convention.
You can by all means have sub-namespaces for specialized functionality,
and the library implementation can be hierarchical,
but for the
love of things don't make users guess which modules they need to `import`.
Granted, naming things is the hardest problem in computer science,
so this is perhaps wishful thinking!

Actually, I lied: the hardest problem in computer science is backward compatibility. A simple, feature-frozen kernel with a versioned standard library could be a way forward (but easier said than done).

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