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
Calcium exact complex field
Future work

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.

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