fredrikj.net / blog /

Infinity in Calcium and Fungrim

January 25, 2021

When I started working on Fungrim, one of my goals was to implement an internally consistent model of infinities, limits, and values at singularities for functions of real and complex variables. Let us call this a "model of singularities" for short. A good model of singularities should 1) be powerful enough to express typical theorems in elementary real and complex analysis in a natural way, 2) allow intuitive reasoning, and 3) be suitable for automated term rewriting or formal theorem proving.

"Implement" first just meant trying to do things consistently throughout Fungrim, which is difficult to enforce, but Calcium is now in progress of implementing the same model of singularities, which can be used in computations.

IEEE 754

Let us start with something that is not a satisfactory model of singularities, for this purpose: the special values in IEEE 754 floating-point arithmetic. These special values are $-0, +0$, $-\infty, +\infty$, and $\operatorname{NaN}$ (not-a-number). On the positive side, the set $F$ of floating-point numbers (comprising both regular numbers and special values) has the useful property of being a closed system: a mathematical function can typically be implemented as a total function $F \to F$ in a way that captures intuitive consequences of taking limits (for example, $1 / \pm 0 = \pm \infty$, and $0 / 0 = \operatorname{NaN}$). Nevertheless, I have some gripes with the IEEE 754 model:

I am not saying that IEEE 754 is a bad model of computer arithmetic altogether; on the contrary, it is a fine model for maximizing the utility of 32 or 64 bits for numerical computing. For example, the distinction between $-0$ and $+0$ makes sense as a way to preserve sign information on underflow, and NaN mostly works well for representing missing data. I am, however, saying that this model is a poor fit for calculus, symbolic computation, and formal reasoning about real and complex numbers. (There are other aspects of IEEE 754 and its implementations that are more problematic generally speaking, such as stateful rounding modes, subnormals, and exception handling, but that's an entirely different topic.)

A sometimes useful feature of IEEE 754's signed zero is that it allows distinguishing between directions of approach at discontinuities, and in particular on branch cuts of multivalued complex functions. For example, it allows defining $\sqrt{-1 \pm 0 i} = \pm i$. However, this is a trick which only works when a branch cut coincides with a coordinate axis, and in any case it is questionable whether the benefits are worth the trouble of violating ring axioms for the number 0, or worse, breaking axioms of propositional logic or set theory. As will be discussed below, there are ways to deal with discontinuities without making this sacrifice.

I should mention John Gustafson's Posit arithmetic, which is an alternative to IEEE 754. In Posit arithmetic, there is only one special value (essentially a projective infinity, which also acts as a kind of NaN when needed), and 0 is unsigned. This is a far more elegant design, but perhaps a bit restrictive; especially with the lack of separate $+\infty$ and $-\infty$ values (a distinction that is useful more often than between $+0$ and $-0$).

In Arb, there is a single zero; there are signed infinities and NaN, but I gave up long ago trying to define functions consistently for special values (especially for complex infinities); it's just not worth the trouble to do complex interval arithmetic with infinities in any clever way in that model. Both internally and in applications using Arb, it tends to make sense to consider any finite interval as input worth handling, and anything non-finite (containing either NaN or infinity) as belonging to a contagious, NaN-like equivalence class. This is again fine for (rigorous) numerical computing over $\mathbb{C}$, but a better model is needed for symbolic computation involving limits.

A better model: Mathematica and Calcium

Calcium extends the real and complex numbers with the following special values, inspired by Mathematica:

The special values behave intuitively, including in the complex domain. For example:

We can now represent various commonly used extensions of the real and complex numbers, including:

Out of these sets, the most useful are arguably $\overline{\mathbb{R}}$ (for ordering things) and $\hat{\mathbb{C}}_{\infty}$ (for meromorphic functions). There are other useful specific sets containing some combination of numbers and infinities; for example, when working with modular forms, it is natural to consider $\mathbb{H} \cup \{ i \infty \}$ and $\mathbb{H} \cup \{ i \infty \} \cup \mathbb{Q}$ where $\mathbb{H}$ is the upper half-plane $\{ z : z \in \mathbb{C}, \operatorname{Im}(z) > 0 \}$.

By design, there are no rectangular signed infinities: for example, $+\infty + i \infty$ evaluates to $\mathfrak{u}$ (more generally, addition of signed infinities $c \cdot \infty + d \cdot \infty$ gives $c \cdot \infty$ if $c = d$ and $\mathfrak{u}$ otherwise). Usually, in complex analysis, either the projectively-extended or sign-extended complex numbers make more sense than $\overline{\mathbb{R}} \times \overline{\mathbb{R}} i$. Another conspicuous design choice is that signed infinities have no signed zero counterparts. There is a certain lack of symmetry in this, and one could indeed imagine having a set of signed zeros (allowing complex signs) to complement the standard unsigned zero 0. However, it would be rather inconvenient if $1 / {\tilde\infty}$ and $1 / (+\infty)$ were not both the number 0.

Combining all numbers, infinities, and undefined, we define the following singularity closures (RealSingularityClosure and ComplexSingularityClosure):

$$\overline{\mathbb{R}}_{\text{Sing}} = \mathbb{R} \cup \left\{\hat{\infty}\right\} \cup \{\pm \infty\} \cup \left\{\mathfrak{u}\right\}$$ $$\overline{\mathbb{C}}_{\text{Sing}} = \mathbb{C} \cup \left\{\hat{\infty}\right\} \cup \{[e^{i \theta}] \infty\} \cup \left\{\mathfrak{u}\right\}$$

The sets $\overline{\mathbb{R}}_{\text{Sing}}$ and $\overline{\mathbb{C}}_{\text{Sing}}$ are intended as closures of $\mathbb{R}$ and $\mathbb{C}$ for computing limits: for any function $f : \overline{\mathbb{C}}_{\text{Sing}} \to \overline{\mathbb{C}}_{\text{Sing}}$ and for any $c \in \overline{\mathbb{C}}_{\text{Sing}}$, $\lim_{z \to c} f(z)$ will have a well-defined value in $\overline{\mathbb{C}}_{\text{Sing}}$ (more on this below). Indeed, any special function on $\mathbb{C}$ can be realized as a total function $f : \overline{\mathbb{C}}_{\text{Sing}} \to \overline{\mathbb{C}}_{\text{Sing}}$ regardless of whether it has singularities of any kind (poles, branch cuts, essential singularities, boundaries of analyticity, ...). Hence the name "singularity closure".

Having such a universal domain is very convenient for describing special functions and for doing computations (as with the closure property of IEEE 754 arithmetic). For this reason, the ca_t number type in Calcium mathematically models $\overline{\mathbb{C}}_{\text{Sing}}$. Generally following Mathematica's convention, unsigned infinity represents the value of meromorphic functions at poles (for example, $1 / 0 = \hat \infty$, and $\Gamma(0) = \hat \infty$) and more generally at algebraic blow-up singularities (for example, $0 ^ {-1/2} = \hat \infty$). Signed infinity represents the value of functions at logarithmic singularities such as $\log(0) = -\infty$. Signed infinities are also used to represent directions when taking limits, and to represent limits at essential and other singularities when it is convenient to do so (for example $e^{+\infty} = +\infty$). Undefined ($\mathfrak{u}$) is generally used at essential singularities, outside boundaries of analyticity, and in other cases that are clearly, well, undefined.

It must be stressed that this choice of singularity closure is a pure matter of convention. Again, IEEE 754 makes a different choice. $\mathbb{C} \cup \{\mathfrak{u}\}$ would also serve just fine (we define any limit that is not in $\mathbb{C}$ as undefined); it just so happens that incorporating different infinities makes life easier in applications.

Adding a meta-value

The Calcium number type ca_t adds another special value: Unknown ($?$). This is a meta-value rather than a value, and represents an undetermined element of $\overline{\mathbb{C}}_{\text{Sing}}$. In IEEE 754, NaN acts both as a value and as a meta-value, leading to the problematic definition of equality ($\operatorname{NaN} \ne \operatorname{NaN}$). In Calcium, these notions are separate: undefined ($\mathfrak{u}$) represents mathematical undefinedness (evaluating a function where a number or infinity would not make sense as an answer), while unknown ($?$) represents computational indeterminacy (the operation may be defined in theory, but the implementation is unable to compute the true result). Here are some examples of calculations that produce Unknown:

Like NaN, Unknown allows nonstop computing, though some applications may want to intercept Unknown values and throw exceptions. Unknown tends to propagate contagiously, but actually less so than $\mathfrak{u}$ or NaN. Indeed, $\mathfrak{u} + ? = \mathfrak{u}$ since the result would be $\mathfrak{u}$ for any value that $?$ could represent.

What about equality? Comparisons with $?$ (including $(? = ?)$ and $(? \ne \, ?)$) are themselves undefined rather than giving false. There are two ways to realize this undefinedness. On the implementation (C) level, Calcium usually handles predicates using triple-valued logic (True, False, Unknown), allowing nonstop computation. On the application level, the intended behavior is to throw an exception when a predicate cannot be decided. To illustrate, Calcium is currently unable to determine whether $1-e^{-e^{-10000}} = 0$ at default precision:

>>> 1-exp(-exp(-10000))
0e-4342 {-a+1 where a = 1.00000 [Exp(-1.13548e-4343 {-b})], b = 1.13548e-4343 [Exp(-10000)]}

>>> 1 / (1-exp(-exp(-10000)))
Unknown

>>> 1-exp(-exp(-10000)) == 0
Traceback (most recent call last):
  ...
ValueError: unable to decide predicate: equal

Since the predicate evaluation throws an exception, the user can think in terms of ordinary mathematical logic when manipulating Calcium numbers (using if-else constructs, for example).

>>> def f(x):
...     if x == 0:
...         return "Zero"
...     else:
...         return "Nonzero"
... 
>>> 
>>> f(sin(pi))
'Zero'
>>> f(sin(pi) + 1)
'Nonzero'
>>> f(1-exp(-exp(-10000)))
Traceback (most recent call last):
  ...
ValueError: unable to decide predicate: equal

This minimizes the risk of mathematical errors. (Of course, functions implementing triple-valued logic may also be provided at the user level to allow nonstop handling of predicates for expert users.)

Again, both triple-valued logic and exceptions are "meta-level" constructs; over $\overline{\mathbb{C}}_{\text{Sing}}$, which is a closed system for arithmetic and function evaluation, ordinary logic applies. This is the important improvement over IEEE 754 semantics.

Differences between Calcium and Mathematica

As noted above, Calcium's model of singularities is closely inspired by Mathematica, but differs in some details.

There are probably other, minor differences that I'm forgetting about. A more fundamental problem with Mathematica's general semantics, of course, is that it ignores special cases when simplifying symbolic formulas, with the result that symbol substitution fails to commute with evaluation — a kind of violation of referential transparency

Different limit operators

The standard notation for limits $\lim_{x \to c} f(x)$ can be quite ambiguous, as the result may be undefined (or defined differently) depending on 1) the set through which $c$ is approached, and 2) the set of permissible limit values (and implicitly, the norms/topologies of these sets). For example:

Note that naming a variable $n$ usually makes it clear to human readers that integers are intended, but a computer algebra system or formal proof system will not understand this. In an attempt to resolve such ambiguities, I introduced several different limit operators in Fungrim:

Operator Approach set Limit set (minus $\mathfrak{u}$)
SequenceLimit $\mathbb{Z}$, $c = \pm \infty$ $\overline{\mathbb{C}}_{[e^{i \theta}] \infty} = \mathbb{C} \cup \{[e^{i \theta}] \infty\}$
RealLimit $c + (-\varepsilon, \varepsilon)$ (for finite $c$) $\overline{\mathbb{C}}_{[e^{i \theta}] \infty} = \mathbb{C} \cup \{[e^{i \theta}] \infty\}$
RightLimit $c + (0, \varepsilon)$ (for finite $c$) $\overline{\mathbb{C}}_{[e^{i \theta}] \infty} = \mathbb{C} \cup \{[e^{i \theta}] \infty\}$
LeftLimit $c + (-\varepsilon, 0)$ (for finite $c$) $\overline{\mathbb{C}}_{[e^{i \theta}] \infty} = \mathbb{C} \cup \{[e^{i \theta}] \infty\}$
ComplexLimit $D(c, \varepsilon)$ (for finite $c$) $\overline{\mathbb{C}}_{[e^{i \theta}] \infty} = \mathbb{C} \cup \{[e^{i \theta}] \infty\}$
MeromorphicLimit $D(c, \varepsilon)$ (for finite $c$) $\hat{\mathbb{C}}_{\infty} = \mathbb{C} \cup \left\{\hat{\infty}\right\}$

These operators will also be used in the Calcium symbolic expression system (though, at the moment, Calcium has no code for actually computing limits).

For all the limit operators, the result is undefined ($\mathfrak{u}$) if the approach set does not contain a sequence that converges to the limit point $c$ or if $f(z)$ does not converge to a limit point (in the intended topology). All operators give sign-extended complex numbers as limits, with the exception of MeromorphicLimit, which as the name suggests can be used to construct meromorphic functions by taking limits of meromorphic functions, with the value $\hat \infty$ at poles. That is the idea, anyway. I have omitted writing down complete definitions in the table for infinities, and actually doing so will probably be a bit tricky.

Presumably, a sequence $f(n)$ can be understood to approach $\hat \infty$ in $\hat{\mathbb{C}}_{\infty}$ if $|f(n)| \to \infty$, regardless of the sign of $f(n)$. Presumably, a sequence $f(n)$ can also be understood to approach $e^{i \theta} \cdot \infty$ in the sign-extended complex numbers if $|f(n)| \to \infty$ and $\arg f(n) \to \theta$ when $n \to \infty$. Here, $\arg f(n)$ should arguably be considered as an equivalence class modulo $2 \pi$ rather than a principal value. (This is an example of the kind of subtle distinction that makes formal definitions tricky.)

All the limit operators allow specifying an optional predicate to further restrict the set of approach. Example 1: a sequence limit can be restricted to a subset of the integers.

Equal(SequenceLimit((-1)**n, For(n, Infinity)), Undefined)
$$\lim_{n \to \infty} {\left(-1\right)}^{n} = \mathfrak{u}$$
Equal(SequenceLimit((-1)**n, For(n, Infinity), IsEven(n)), 1)
$$\lim_{n \to \infty,\,n \text{ even}} {\left(-1\right)}^{n} = 1$$

Example 2: clearly $\arg(e^x) = 0$ for real $x$, and this should be the limit as $x \to \infty$ along the real line. The limit of $\arg(e^z)$ is undefined when $z$ approaches $+\infty$ through the complex plane, but it can approach any fixed value in $(-\pi,\pi]$ by fixing the imaginary part of the approach set so that the limit is taken along a horizontal line in the complex plane:

Equal(RealLimit(Arg(Exp(x)), For(x, +Infinity)), 0)
$$\lim_{x \to +\infty} \arg(e^{x}) = 0$$
(Equal(ComplexLimit(Arg(Exp(z)), For(z, +Infinity)), Undefined)
$$\lim_{z \to +\infty} \arg(e^{z}) = \mathfrak{u}$$
(Equal(ComplexLimit(Arg(Exp(z)), For(z, +Infinity), Equal(Im(z), 1)), 1)
$$\lim_{z \to +\infty,\,\operatorname{Im}(z) = 1} \arg(e^{z}) = 1$$

Path objects

Limit operators make the distinction between $+0$ and $-0$ unnecessary, and indeed limit operators are strictly more powerful since one can encode any path of approach to any point. However, symbolic limit operators are also rather clumsy compared to numbers or number-like objects.

An alternative, or rather a complement, is to introduce path objects. In its simplest form, a path object Path(a, b) may encode a straight-line path $a \rightsquigarrow b$ between two points (complex numbers or infinities). More generally, a path object may encode a path $z_1 \rightsquigarrow z_2 \rightsquigarrow \ldots \rightsquigarrow z_n$ connecting an arbitrary number of points, or even piecewise curvilinear paths. The idea is that one should be able to evaluate a holomorphic function at a path object with a simple function call, with the interpretation that we perform analytic continuation along the path starting from a neighborhood of the initial point. We would thus have:

Equal(Log(Path(1, +NumberI, -1)), Pi * NumberI)
$$\log\left(1 \rightsquigarrow +i \rightsquigarrow -1\right) = \pi i$$
Equal(Log(Path(1, -NumberI, -1)), -Pi * NumberI)
$$\log\left(1 \rightsquigarrow -i \rightsquigarrow -1\right) = -\pi i$$

Indeed, with this interpretation, we can perform analytic continuation along the Riemann surface of the function, without being restricted to the principal branch:

Equal(Log(Path(1, NumberI, -1, -NumberI, 1)), 2 * Pi * NumberI)
$$\log\left(1 \rightsquigarrow i \rightsquigarrow -1 \rightsquigarrow -i \rightsquigarrow 1\right) = 2 \pi i$$

This is just notation for now; it remains to define concrete semantics. In terms of implementations, Marc Mezzarobba's numerical evaluation code in ore_algebra supports analytic continuation of D-finite functions along such paths, though it does not (yet) support infinite points and evaluation in irregular singularities.



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