fredrikj.net / blog /

# Bringing Calcium to Julia

*June 3, 2021*

Good news everyone! I am currently wrapping Calcium in Nemo.jl, which will provide convenient exact real and complex arithmetic in Julia. This post is a quick tour of what is already working in the current development version of Nemo (as soon as my latest pull request gets merged); more functionality will come later.

Contents

## Algebraic numbers

Nemo is getting a domain `QQBar`
representing the field of complex algebraic numbers $\overline{\mathbb{Q}}$.

julia> QQBar Field of Algebraic Numbers in minimal polynomial representation julia> x = QQBar(2) Root 2.00000 of x - 2 julia> sqrt(x) Root 1.41421 of x^2 - 2 julia> sqrt(x) ^ 2 Root 2.00000 of x - 2

Elements of `QQBar` wrap the C type `qqbar_t`.
This is not the most efficient representation for number-crunching
(see `CalciumField` below), but it
has the advantage that algebraic numbers are represented
in a canonical form and that the parent object can be a constant singleton,
just like `QQ` and `ZZ`.

### Roots and eigenvalues

The field of algebraic numbers contains all roots of polynomials over $\mathbb{Z}$ and $\mathbb{Q}$. We can solve quintic equations, for example:

julia> R, x = PolynomialRing(QQ, "x") (Univariate Polynomial Ring in x over Rational Field, x) julia> v = roots(x^5 - x - 1, QQBar) 5-element Array{qqbar,1}: Root 1.16730 of x^5 - x - 1 Root 0.181232 + 1.08395*im of x^5 - x - 1 Root 0.181232 - 1.08395*im of x^5 - x - 1 Root -0.764884 + 0.352472*im of x^5 - x - 1 Root -0.764884 - 0.352472*im of x^5 - x - 1 julia> v[1]^5 - v[1] - 1 == 0 true

Similarly, we can compute exact eigenvalues of any matrix with rational entries. Here is the 6 by 6 Hilbert matrix; we verify that adding and multiplying the eigenvalues gives the trace and determinant of the matrix:

julia> H = hilbert(MatrixSpace(QQ, 6, 6)) [ 1 1//2 1//3 1//4 1//5 1//6] [1//2 1//3 1//4 1//5 1//6 1//7] [1//3 1//4 1//5 1//6 1//7 1//8] [1//4 1//5 1//6 1//7 1//8 1//9] [1//5 1//6 1//7 1//8 1//9 1//10] [1//6 1//7 1//8 1//9 1//10 1//11] julia> v = eigenvalues(H, QQBar) 6-element Array{qqbar,1}: Root 1.61890 of 186313420339200000x^6 - 349935855575040000x^5 + 78981336366912000x^4 - 1242627237734400x^3 + 750409713900x^2 - 9316560x + 1 Root 0.242361 of 186313420339200000x^6 - 349935855575040000x^5 + 78981336366912000x^4 - 1242627237734400x^3 + 750409713900x^2 - 9316560x + 1 Root 0.0163215 of 186313420339200000x^6 - 349935855575040000x^5 + 78981336366912000x^4 - 1242627237734400x^3 + 750409713900x^2 - 9316560x + 1 Root 0.000615748 of 186313420339200000x^6 - 349935855575040000x^5 + 78981336366912000x^4 - 1242627237734400x^3 + 750409713900x^2 - 9316560x + 1 Root 1.25708e-5 of 186313420339200000x^6 - 349935855575040000x^5 + 78981336366912000x^4 - 1242627237734400x^3 + 750409713900x^2 - 9316560x + 1 Root 1.08280e-7 of 186313420339200000x^6 - 349935855575040000x^5 + 78981336366912000x^4 - 1242627237734400x^3 + 750409713900x^2 - 9316560x + 1 julia> sum(v) Root 1.87821 of 3465x - 6508 julia> prod(v) Root 5.36730e-18 of 186313420339200000x - 1 julia> tr(H) 6508//3465 julia> det(H) 1//186313420339200000

### Trigonometric functions

The field of algebraic numbers supports evaluating roots of unity and trigonometric functions (at points where the functions are algebraic):

julia> sinpi(QQBar(7) // 13) Root 0.992709 of 4096x^12 - 13312x^10 + 16640x^8 - 9984x^6 + 2912x^4 - 364x^2 + 13 julia> tanpi(atanpi(sqrt(QQBar(2)) + 1)) Root 2.41421 of x^2 - 2x - 1 julia> w = root_of_unity(QQBar, 5) Root 0.309017 + 0.951057*im of x^4 + x^3 + x^2 + x + 1 julia> is_root_of_unity(w) true julia> w ^ 5 Root 1.00000 of x - 1

### Generic structures

All the generic functionality in AbstractAlgebra/Nemo will work with
`QQBar` as a base field. For example, one can have
matrices and polynomials with algebraic coefficients:

julia> Rx, x = PolynomialRing(QQBar, "x") (Univariate Polynomial Ring in x over Field of Algebraic Numbers in minimal polynomial representation, x) julia> gcd(x^4 - 4*x^2 + 4, x^2 + sqrt(QQBar(18))*x + 4) x + (Root 1.41421 of x^2 - 2)

### Numerical evaluation and guessing

Algebraic numbers can be evaluated numerically to arbitrary precision by passing them as input to real or complex Arb fields:

julia> RR = ArbField(333) Real Field with 333 bits of precision and error bounds julia> RR(sqrt(QQBar(2))) [1.414213562373095048801688724209698078569671875376948073176679737990732478462107038850387534327641573 +/- 4.08e-100]

One of my favorite features is the reverse operation: algebraic numbers can be recovered from (sufficiently accurate) numerical values:

julia> RR = ArbField(53); guess(QQBar, RR("1.41421356 +/- 1e-6"), 2) Root 1.41421 of x^2 - 2

Let us take something nontrivial: the Dedekind eta function quotient $z = \eta(7i) / \eta(i)$. Guessing that this is an algebraic number with degree at most 30:

julia> CC = AcbField(333) Complex Field with 333 bits of precision and error bounds julia> z = modeta(CC(0,7)) / modeta(CC(0,1)) [0.208269233523258141573410018682924054967717880513667670981619809590183080146632703080166195454304592 +/- 2.06e-100] julia> guess(QQBar, z, 30) ERROR: No suitable algebraic number found

Let us try again with higher precision:

julia> CC = AcbField(1000) Complex Field with 1000 bits of precision and error bounds julia> z = modeta(CC(0,7)) / modeta(CC(0,1)) [0.208269233523258141573410018682924054967717880513667670981619809590183080146632703080166195454304591983441975797383962968667163200680700040977011368671247713335495619061581467767859763828227778622640895776232029285094042810118030389377257188940794310378412363657736045059286161683576164665929286614363 +/- 6.80e-301] julia> guess(QQBar, z, 30) Root 0.208269 of 823543x^16 + 235298x^12 + 21609x^8 + 490x^4 - 1

VoilĂ ! This is not a rigorous proof that we have the correct value, of course, but a wrong guess is very unlikely at high precision.

## Calcium exact complex field

The new `CalciumField` parent represents the field of
complex numbers.

julia> C = CalciumField() Exact Complex Field julia> x = C(pi) 3.14159 {a where a = 3.14159 [Pi]} julia> y = C(1im) 1.00000*I {a where a = I [a^2+1=0]} julia> x * y 3.14159*I {a*b where a = 3.14159 [Pi], b = I [b^2+1=0]} julia> exp(x * y) -1

A brief reminder of how it works: Calcium represents numbers as elements of extension fields of the rational numbers. For example, $\pi i$ will be represented as $a b$ in $\mathbb{Q}(a,b)$ where $a = \pi$, $b = i$. The construction of extension numbers and extension fields is fully automatic. Many simplifications are also automatic; for example, when evaluating $e^{\pi i}$, the result automatically collapses to the trivial field $\mathbb{Q}$.

The `CalciumField` parent object is essentially a wrapper
around the C type `ca_ctx_t` while the elements
wrap the C type `ca_t`.
Unlike the constant object `QQBar`, the user needs to create an instance
(`C = CalciumField()`)
of the Calcium exact complex field.
There are several reasons for this: Calcium fields cache data (potentially a lot of it),
which you don't want to stay alive forever;
Calcium fields objects also hold various evaluation options, so
you may want to create different instances for different calculations.

Here are a few more examples:

julia> log(C(10)^23) // log(C(100)) 11.5000 {23/2} julia> 4*atan(C(1)//5) - atan(C(1)//239) == C(pi)//4 true julia> Cx, x = PolynomialRing(C, "x") (Univariate Polynomial Ring in x over Exact Complex Field, x) julia> (a, b) = (sqrt(C(2)), sqrt(C(3))) (1.41421 {a where a = 1.41421 [a^2-2=0]}, 1.73205 {a where a = 1.73205 [a^2-3=0]}) julia> (x - a - b) * (x - a + b) * (x + a - b) * (x + a + b) x^4 + (-10)*x^2 + 1

As a special case, we can use `CalciumField` for efficient
number field arithmetic:

julia> phi = (sqrt(C(5)) + 1) // 2 1.61803 {(a+1)/2 where a = 2.23607 [a^2-5=0]} julia> phi^100 7.92071e+20 {(354224848179261915075*a+792070839848372253127)/2 where a = 2.23607 [a^2-5=0]} julia> (phi^100 - (1 - phi)^100) // sqrt(C(5)) 3.54225e+20 {354224848179261915075} julia> fibonacci(ZZ(100)) 354224848179261915075

### Conversions

Calcium field elements can be converted to Arb real and complex fields
for numerical evaluation.
Integer/rational/algebraic values can also be converted to `ZZ`,
`QQ` and `QQBar` elements.

julia> x = sqrt(C(2)) + sqrt(C(3)) + sqrt(C(5)) + sqrt(C(7)) 8.02808 {a+b+c+d where a = 2.64575 [a^2-7=0], b = 2.23607 [b^2-5=0], c = 1.73205 [c^2-3=0], d = 1.41421 [d^2-2=0]} julia> ArbField(333)(x) [8.0280836585063526292399244880861071066633546718813046058717184220542552442389956031063410159382160340 +/- 4.81e-102] julia> QQBar(x) Root 8.02808 of x^16 - 136x^14 + 6476x^12 - 141912x^10 + 1513334x^8 - 7453176x^6 + 13950764x^4 - 5596840x^2 + 46225

Note that `CalciumField` allows representing algebraic numbers in multivariate
form; it keeps the extension numbers $\sqrt{2}$, $\sqrt{3}$, $\sqrt{5}$, $\sqrt{7}$ separate
instead of coercing them to a single extension. A roundtrip `CalciumField` → `QQBar` → `CalciumField`
can be used to convert to a single extension:

julia> C(QQBar(x)) 8.02808 {a where a = 8.02808 [a^16-136*a^14+6476*a^12-141912*a^10+1513334*a^8-7453176*a^6+13950764*a^4-5596840*a^2+46225=0]}

Of course, the following fails:

julia> x = C(pi) 3.14159 {a where a = 3.14159 [Pi]} julia> QQBar(x) ERROR: unable to convert to an algebraic number

### Comparisons and predicates

Comparisons act on the mathematical values of Calcium field elements, not on their representations. For example, evaluating $x = \sqrt{2} \sqrt{3}$ and $y = \sqrt{6}$ results in different internal representations ($x \in \mathbb{Q}(\sqrt{3}, \sqrt{2})$ and $y \in \mathbb{Q}(\sqrt{6})$), but the numbers compare as equal:

julia> x = sqrt(C(2)) * sqrt(C(3)) 2.44949 {a*b where a = 1.73205 [a^2-3=0], b = 1.41421 [b^2-2=0]} julia> y = sqrt(C(6)) 2.44949 {a where a = 2.44949 [a^2-6=0]} julia> x == y true julia> iszero(x - y) true julia> isinteger(x - y) true

Predicate functions return `true` if the property is provably true
and `false` if the property if provably false. If Calcium is unable
to prove the truth value, an exception is thrown.
For example, with default settings, Calcium is currently
able to prove that $e^{e^{-1000}} \ne 1$,
but it fails to prove $e^{e^{-3000}} \ne 1$:

julia> x = exp(exp(C(-1000))) 1.00000 {a where a = 1.00000 [Exp(5.07596e-435 {b})], b = 5.07596e-435 [Exp(-1000)]} julia> x == 1 false julia> x = exp(exp(C(-3000))) 1.00000 {a where a = 1.00000 [Exp(1.30784e-1303 {b})], b = 1.30784e-1303 [Exp(-3000)]} julia> x == 1 ERROR: Unable to perform operation (failed deciding truth of a predicate): isequal ...

We can get an answer by creating a Calcium field that allows a higher internal working precision than the default value:

julia> C2 = CalciumField(options=Dict(:prec_limit => 10^5)); julia> exp(exp(C2(-3000))) == 1 false

### Elementary and special functions

Elementary and special functions generally create new extension numbers. In special cases, simplifications occur automatically.

julia> exp(C(1)) 2.71828 {a where a = 2.71828 [Exp(1)]} julia> exp(C(0)) 1 julia> atan(C(1)) 0.785398 {(a)/4 where a = 3.14159 [Pi]} julia> cos(C(1))^2 + sin(C(1))^2 1 julia> log(1 // exp(sqrt(C(2))+1)) == -sqrt(C(2)) - 1 true julia> gamma(C(2+3im)) -0.0823953 + 0.0917743*I {a where a = -0.0823953 + 0.0917743*I [Gamma(2.00000 + 3.00000*I {3*b+2})], b = I [b^2+1=0]} julia> gamma(C(5) // 2) 1.32934 {(3*a)/4 where a = 1.77245 [Sqrt(3.14159 {b})], b = 3.14159 [Pi]} julia> erf(C(1)) 0.842701 {a where a = 0.842701 [Erf(1)]} julia> erf(C(1)) + erfc(C(1)) 1

Calcium can prove any constant identity involving algebraic numbers, but transcendental functions are far more tricky to deal with, and the representations and simplification routines in Calcium are still a work in progress. Here is an example of two constant trigonometric identities where the first succeeds and the second fails:

julia> a = sqrt(C(2)) + 1; julia> cos(a) + cos(2*a) + cos(3*a) == sin(7*a//2)//(2*sin(a//2)) - C(1)//2 true julia> sin(3*a) == 4 * sin(a) * sin(C(pi)//3 - a) * sin(C(pi)//3 + a) ERROR: Unable to perform operation (failed deciding truth of a predicate): isequal

There is not a single universal best way to represent extension numbers. For example, depending the application, we might want to represent $\sin(x)$ either in terms of the direct function, in terms of complex exponentials, or in terms of tangents. This is partially supported:

julia> s1 = sin(C(1)) 0.841471 - 0e-24*I {(-a^2*b+b)/(2*a) where a = 0.540302 + 0.841471*I [Exp(1.00000*I {b})], b = I [b^2+1=0]} julia> s2 = sin(C(1), form=:direct) 0.841471 {a where a = 0.841471 [Sin(1)]} julia> s3 = sin(C(1), form=:exponential) 0.841471 - 0e-24*I {(-a^2*b+b)/(2*a) where a = 0.540302 + 0.841471*I [Exp(1.00000*I {b})], b = I [b^2+1=0]} julia> s4 = sin(C(1), form=:tangent) 0.841471 {(2*a)/(a^2+1) where a = 0.546302 [Tan(0.500000 {1/2})]} julia> s1 == s2 == s3 == s4 true julia> isreal(s1) && isreal(s2) && isreal(s3) && isreal(s4) true

The default representation of trigonometric functions is an example
of setting that can be changed in the parent `CalciumField`.
The exponential form is currently used by default since it tends to give
the best performance for symbolic simplification, but it is worse
for numerical evaluation when the arguments are real.
As seen in the output `0.841471 - 0e-24*I`, Calcium does not know
that the imaginary part is exactly zero at the time of printing,
but `isreal` does verify that the number is really real.

### Infinities and special values

A major design point of the Calcium field is that it
strictly represents the mathematical field of complex numbers $\mathbb{C}$ by default.
This means that infinities and NaNs are **BANNED**:

julia> 1 // C(0) ERROR: DomainError with UnsignedInfinity: Non-number result julia> 0 // C(0) ERROR: DomainError with Undefined: Non-number result julia> log(C(0)) ERROR: DomainError with -Infinity: Non-number result

We can work with such special values, but this requires creating a parent object that explicitly allows them:

julia> Cext = CalciumField(extended=true) Exact Complex Field (Extended) julia> 1 // Cext(0) UnsignedInfinity julia> log(Cext(0)) -Infinity julia> 0 // Cext(0) Undefined julia> exp(log(Cext(0))) 0 julia> 1 // (1 // Cext(0)) 0 julia> -infinity(Cext) < Cext(42) < infinity(Cext) true

See my previous blog post on infinities for more remarks on this topic.

## Future work

Things I have not yet wrapped in Nemo include Calcium types and functions for matrices, polynomials and power series. The generic implementations in AbstractAlgebra.jl will do the job, but the C implementations are faster and provide more functionality (for example, matrix functions).

There is also a lot left to do in terms of convenience functions for input, output, conversions, symbolic expressions, and random generation.

A lot of development work remains on the C side as well; some fairly basic operations are still missing or poorly implemented in Calcium. The Julia interface should make it easier to experiment with new algorithms and functions. There is already a Python interface in Calcium that I'm using for tests, but it doesn't have generic algebraic structures.

If anyone wants to volunteer, it should be possible to do a similar interface to Calcium for SageMath.