Nemo/Hecke: Computer Algebra and Number Theory Packages for the Julia Programming Language

Julia provides two levels of types: abstract and concrete types. Concrete types are like types in Java, C or C++, etc. Abstract types can be thought of as collections of types. They are used when writing generic functions that should work for any type in a collection.

To write such a generic function, we first create an abstract type, then we create the individual concrete types that belong to that abstract type. The generic function is then specified with a type parameter, T say, belonging to the abstract type.

In Julia, the symbol ¡: is used to specify that a given type belongs to a given abstract type. For example the built-in Julia type Int64 for 64-bit machine integers belongs to the Julia abstract type Integer.

Abstract types in Julia can form a hierarchy. For example, the Nemo.Field abstract type belongs to the Nemo.Ring abstract type. An object representing a field in Nemo has type belonging to Nemo.Field. But because we define the inclusion Nemo.Field ¡: Nemo.Ring, the type of such an object also belongs to Nemo.Ring. This means that any generic function in Nemo which is designed to work with ring objects will also work with field objects.

Julia always picks the most specific function that applies to a given type. This allows one to implement a function at the most general level of a type hierarchy at which it applies. One can also write a version of a given function for specific concrete types. For example, one may wish to call a specific C implementation for a multiplication algorithm, say, when arguments with a certain very specific given type are passed.

Another way that we make use of Julia’s abstract type system in Nemo/Hecke is to distinguish between the type of elements of fields, and the fields themselves, and similarly for all other kinds of domains in Nemo.

Figure 1 shows the abstract type hierarchy in Nemo.

Naively, one may expect that specific mathematical domains in Nemo/Hecke can be modeled as types and their elements as objects of the given type. But there are various reasons why this is not a good model.

As an example, consider the ring $R=\mathbb{Z}/n\mathbb{Z}$. If we were to model the ring $R$ as a type, then the type would need to contain information about the modulus $n$. This is not possible in Julia if $n$ is an object, e.g. a multiprecision integer. Further, Julia dispatches on type, and each time we call a generic function with values of a new type, the function is recompiled by the JIT compiler for that new type. This would result in very poor performance if we were writing a multimodular algorithm, say, as recompilation would be triggered for each distinct modulus $n$. For this reason, the modulus $n$ needs to be attached to the elements of the ring, not to the type associated with those elements.

The way we deal with this is to have special (singleton) objects, known as parent objects, that act like types, but are in fact ordinary Julia objects. As ordinary objects, parents can contain arbitrary information, such as the modulus $n$ in $\mathbb{Z}/n\mathbb{Z}$. Each object representing an element of the ring then contains a pointer to this parent object.

This model of mathematical parent objects is taken from SageMath (sage, ) which in turn followed Magma (magma, ).

Julia allows ordinary objects to be made callable. We make use of the facility to write coercions and constructors for elements of mathematical domains in Nemo/Hecke. For example, the following code constructs $a=3\pmod{7}$.

R = ResidueRing(ZZ, 7)
a = R(3)

We also make use of the parent object system to encode information such as context objects needed by C libraries. As Julia objects can have precisely the same bit representation as native C objects, parent objects can be passed directly to C functions. Additional fields in these objects can be safely appended if it is desirable to retain more information at the Julia level than the C/C++ level.

Nemo wraps C types provided by libraries such as Flint for ground domains. For example, Nemo’s ZZ wraps Flint’s fmpz integer type. Nemo also uses specialised C implementations for nested structures where available. For instance, PolynomialRing(R, "x") which constructs the univariate polynomial ring $R[x]$ automatically switches to a wrapper of Flint’s fmpz_poly when $R=\mathbb{Z}$ instead of using the Julia implementation designed for generic $R$.

C libraries such as Flint allow mutating objects in-place, which improves performance especially when making incremental updates to large objects. Nemo objects are ostensibly immutable, but special mutation methods are provided for use in critical spots. Instead of writing

s += a * b

in Nemo, which creates an object for a * b and then replaces s with yet another new object, we create a reusable scratch variable t with the same type as a, b and s and use

mul!(t, a, b)
addeq!(s, t)

which creates no new objects and avoids making a copy of the data in s that is not being modified.

In Julia, in-place operators save not only memory allocation and copying overheads, but also garbage collection costs, which may be substantial in the worst case.

The Fateman benchmark tests arithmetic in $\mathbb{Z}[x_{1},\ldots,x_{k}]$ by computing $f\cdot(f+1)$ where $f=(1+x+y+z+t)^{n}$. In Nemo, this is expressed by the following Julia code:

R,x,y,z,t = PolynomialRing(ZZ, ["x","y","z","t"])
f = (1+x+y+z+t)^30
f*(f+1)

Table 1 shows timing results compared to SageMath 7.4, Magma V2.22-5 and Giac-1.2.2 on a 2.2 GHz AMD Opteron 6174. We used sparse multivariate polynomials in all cases (this benchmark could also be performed using nested dense univariate arithmetic, i.e. in $\mathbb{Z}[x][y][z][t]$).

Nemo (generic) is our Julia implementation of Johnson’s sparse multiplication algorithm which handles $R[x_{1},\ldots,x_{k}]$ for generic coefficient rings $R$. Flint (no asm) is a reimplementation of the same algorithm in C for the special case $R=\mathbb{Z}$, and Flint (asm) is an assembly optimised version of this code. Finally, Flint (array) and Giac implement an algorithm designed for dense polynomials, in which terms are accumulated into a big array covering all possible exponents.

Table 2 shows timing results on the Pearce benchmark, which consists of computing $f\cdot g$ where $f=(1+x+y+2z^{2}+3t^{3}+5u^{5})^{n}$, $g=(1+u+t+2z^{2}+3y^{3}+5x^{5})^{n}$. This problem is too sparse to use the big array method, so only the Johnson algorithm is used.

By default, Nemo uses Flint for multivariate polynomials over $R=\mathbb{Z}$ instead of the generic Julia code ²²2The new Flint type fmpz_mpoly is presently available in the git version of Flint. The mul_johnson (with and without asm enabled) and mul_array methods were timed via Nemo.. We have timed both versions here to provide a comparison between Julia and C. It is somewhat remarkable that the generic Julia code comes within a factor two of our C version in Flint (without asm), and even runs slightly faster than the C code in Magma.

The following Nemo example code computes a resultant in the ring $((\operatorname{GF}(17^{11})[y])/(y^{3}+3xy+1))[z]$, demonstrating generic recursive rings and polynomial arithmetic. Flint is used for univariate polynomials over a finite field.

R, x = FiniteField(17, 11, "x")
S, y = PolynomialRing(R, "y")
T = ResidueRing(S, y^3 + 3x*y + 1)
U, z = PolynomialRing(T, "z")

f = (3y^2+y+x)*z^2 + ((x+2)*y^2+x+1)*z + 4x*y + 3
g = (7y^2-y+2x+7)*z^2 + (3y^2+4x+1)*z + (2x+1)*y + 1
s = f^12
t = (s + g)^12
resultant(s, t)

This example takes 179907 s in SageMath 6.8, 82 s in Magma V2.21-4 and 0.2 s in Nemo 0.4.

We compute the determinant of a random $80\times 80$ matrix with entries in a cubic number field. This benchmark tests generic linear algebra over a field with coefficient blowup. The number field arithmetic is provided by Antic.

QQx, x = PolynomialRing(QQ, "x")
K, a = NumberField(x^3 + 3*x + 1, "a")
M = MatrixSpace(K, 80, 80)()

for i in 1:80
 for j in 1:80
   for k in 0:2
     M[i, j] = M[i, j] + a^k * (rand(-100:100))
   end
 end
end

det(M)

This takes 5893 s in SageMath 6.8, 21.9 s in Pari/GP 2.7.4 (pari, ), 5.3 s in Magma V2.21-4, and 2.4 s in Nemo 0.4.

For generic coefficient rings, Nemo implements dense univariate polynomials, truncated power series (supporting both relative and absolute precision models), and multivariate polynomials in sparse distributed format.

The generic polynomial resultant code in Nemo uses subresultant pseudoremainder sequences. This code can even be called for polynomials over a ring with zero divisors, so long as impossible inverses are caught. The exception can be caught using a try/catch block in Julia and a fallback method called, e.g. the resultant can be computed using Sylvester matrices. This allows an algorithm with quadratic complexity to be called generically, with an algorithm of cubic complexity only called as backup.

We make use of variants of the sparse algorithms described by Monagan and Pearce for heap-based multiplication (heapmul, ), exact division and division with remainder (heapdiv, ) and powering (heappow, ). For powering of polynomials with few terms, we make use of the multinomial formula. For multivariate polynomials over $\mathbb{Z}$, we used the generic Julia code as a template to implement a more optimised C version in Flint.

In the case of pseudodivision for sparse, multivariate polynomials, we extract one of the variables and then use a heap-based, univariate pseudoremainder algorithm inspired by an unpublished manuscript of Monagan and Pearce.

We make use of this pseudoremainder implementation to implement the subresultant GCD algorithm. To speed the latter up, we make use of a number of tricks employed by the Giac/XCAS system (giac, ), as taught to us by Bernard Parisse. The most important of these is cheap removal of the content that accumulates during the subresultant algorithm by analysing the input polynomials to see what form the leading coefficient of the GCD can take.

We also use heuristics to determine which permutation of the variables will lead to the GCD being computed in the fastest time. This heuristic favours the variable with the lowest degree as the main variable for the computation, with the other variables following in increasing order of degree thereafter. But our heuristic heavily favours variables in which the polynomial is monic.

For computing the determinant over a generic commutative ring, we implemented a generic algorithm making use of Clow sequences (clow, ) which uses only $O(n^{4})$ ring operations in the dimension $n$ of the matrix. However, two other approaches seem to always outperform it. The first approach makes use of a generic fraction free LU decomposition. This algorithm may fail if an impossible inverse is encountered. However, as a fallback, we make use of a division free determinant algorithm. This computes the characteristic polynomial and then extracts the determinant from that.

For determinants of matrices over polynomial rings, we use an interpolation approach.

We implemented an algorithm for computing the characteristic polynomial due to Danilevsky (see below) and an algorithm that is based on computing the Hessenberg form. We also make use of a generic implementation of the algorithm of Berkowitz which is division free.

As is well known, fraction free algorithms often improve the performance of matrix algorithms over an integral domain. For example, fraction free Gaussian elimination for LU decomposition, inverse, determinant and reduced row echelon form are well-known.

We have also been able to use this strategy in the computation of the characteristic polynomial. We have adapted the well-known 1937 method of Danilevsky (danilevsky, ) for computing the characteristic polynomial, into a fraction-free version.

Danilevsky’s method works by applying similarity transforms to reduce the matrix to Frobenius form. Normally such computations are done over a field, however each of the outer iterations in the algorithm introduce only a single denominator. Scaling by this denominator allows us to avoid fractions. The entries become larger as the algorithm proceeds because of the scaling, but conveniently it is possible to prove that the introduced scaling factor can be removed one step later in the algorithm. This is an exact division and does not lead to denominators.

Removing such a factor has to be done in a way that respects the similarity transforms. We achieve this by making two passes over the matrix to remove the common factor.

Flint provides arithmetic over $\mathbb{Z}$, $\mathbb{Q}$, $\mathbb{Z}/n\mathbb{Z}$, $\operatorname{GF}(q),\mathbb{Q}_{p}$ as well as matrices, polynomials and power series over most of these ground rings. Nemo implements elements of these rings as thin wrappers of the Flint C types.

Flint uses a number of specialised techniques for each domain. For example, Flint’s multiplication in $\mathbb{Z}[x]$ uses a best-of-breed algorithm which selects classical multiplication, Karatsuba multiplication, Kronecker segmentation or a Schönhage-Strassen FFT, with accurate cutoffs between algorithms, depending on the degrees and coefficient sizes.

In some cases, Flint provides separate implementations for word-size and arbitrary-size coefficients. Nemo wraps both versions and transparently selects the optimised word-size version when possible.

Additional Flint algorithms wrapped by Nemo include primality testing, polynomial factorisation, LLL, and Smith and Hermite normal forms of integer matrices.

Computing over $\mathbb{R}$ and $\mathbb{C}$ requires using some approximate model for these fields. The most common model is floating-point arithmetic. However, for many computer algebra algorithms, the error analysis necessary to guarantee correct results with floating-point arithmetic becomes impractical. Interval arithmetic solves this problem by effectively making error analysis automatic.

Nemo includes wrapper code for the C library Arb, which implements real numbers as arbitrary-precision midpoint-radius intervals (balls) $[m\pm r]$ and complex numbers as rectangular boxes $[a\pm r_{1}]$ + $[b\pm r_{2}]i$. Nemo stores the precision in the parent object. For example, R = ArbField(53) creates a field of Arb real numbers with 53-bit precision. Arb also supplies types for polynomials, power series and matrices over $\mathbb{R}$ and $\mathbb{C}$, as well as transcendental functions. Like Flint, Arb systematically uses asymptotically fast algorithms for operations such as polynomial multiplication, with tuning for different problem sizes.

Many problems can be solved using lazy evaluation: the user can try a computation with some tentative precision $p$ and restart with precision $2p$ if that fails. The precision can be set optimally when a good estimate for the minimal required $p$ is available; that is, the intervals can be used as if they were plain floating-point numbers, and the automatic error bounds simply provide a certificate.

Alternative implementations of $\mathbb{R}$ and $\mathbb{C}$ may be added to Nemo in the future. For example, it would sometimes be more convenient to use a lazy bit stream abstraction in which individual numbers store a symbolic DAG representation to allow automatically increasing their own precision.

Antic is a C library, depending on Flint, for doing efficient computations in algebraic number fields. Nemo uses Antic for number field element arithmetic. We briefly describe some of the techniques Antic uses for fast arithmetic, but refer the reader to the article (antic, ) for full details.

Antic represents number field elements as Flint polynomials thereby benefiting from the highly optimised polynomial arithmetic in Flint. However, a few more tricks are used.

Firstly, Antic makes a number of precomputations which are stored in a context object to speed up subsequent computations with number field elements. Number field parent objects in Nemo consist precisely of these context objects.

The first is a precomputed inverse of the leading coefficient of the defining polynomial of the number field. This helps speed up reduction modulo the defining polynomial.

The second is an array of precomputed powers $x^{i}$ modulo the defining polynomial $f(x)$. This allows fast reduction of polynomials whose degree exceeds that of the defining polynomial, e.g. in multiplication of number field elements $g$ and $h$ where one wants wants to compute $g(x)h(x)\pmod{f(x)}$.

The third precomputation is of Newton sums $S_{k}=\sum_{i=1}^{n}\theta_{i}^{k}$ where the $\theta_{i}$ are the roots of the defining polynomial $f(x)$. These Newton sums are precomputed using recursive formulae as described in (Cohen1993, ). They are used to speed up the computation of traces of elements of the number field.

Norms of elements of number fields are computed using resultants, for which we use the fast polynomial resultant code in Flint. Inverses of elements are computed using the fast polynomial extended GCD implementation in Flint.

Antic also offers the ability to do multiplication of number field elements without reduction. This is useful for speeding up the dot products that occur in matrix multiplication, for example. Instead of reducing after every multiplication, the unreduced products are first accumulated and their sum can be reduced at the end of the dot product.

To facilitate delayed reduction, all Antic number field elements are allocated with sufficient space to store a full polynomial product, rather than the reduction of such a product.

Singular (singular, ) is a C/C++ package for polynomial systems and algebraic geometry. Recently, we helped prepare a Julia package called Singular.jl, which is compatible with Nemo. It will be described in a future article.

A classical result is that in Dedekind domains, as for the ring of integers of a number field, any ideal can be generated using only two elements, one of which can be chosen at random in the ideal. This representation is efficient in terms of space compared to storing a $\mathbb{Z}$-basis. However, the key problem is the lack of efficient algorithms for working with two generators. We remark that one operation is always efficient using two generators: powering of ideals.

A refinement of this idea, due to Pohst and Zassenhaus (Pohst1997, , p. 400), is a normal presentation. Here a finite set $S$ of prime numbers is fixed, typically containing at least all prime divisors of the norm of the ideal. Then, based on this set, the two generators are chosen: The first as an integer having only divisors in $S$—but typically with too high multiplicities. The second generator then has the correct multiplicity at all prime ideals over primes in $S$—but is allowed to have other divisors as well.

For the remainder of this section, let $K/\mathbb{Q}$ be a number field of degree $n$ and $\mathcal{O}_{K}$ be its ring of integers.

A direct application of the Chinese remainder theorem shows the existence of such normal presentations. The algorithmic importance comes from the following result.

It should be noted that (almost all) prime ideal are naturally computed in $\{p\}$-normal presentation. The key problem in using the theorem for multiplication is that both ideals need to be given by an $S$-normal presentation for the same set $S$, which traditionally is achieved by recomputing a suitable presentation from scratch, taking random candidates for the second generator until it works. Based on the following lemma, we have developed a new algorithm that manages it at the cost of a few integer operations only.

Using this lemma on both input ideals, we can obtain compatible $S$-normal presentations at the cost of two $\gcd$ computations, a few integer multiplications and two products of integers with algebraic numbers. The final ideal multiplication is then a single integer product and a product of algebraic numbers. Thus the asymptotic cost is that of a single multiplication of two algebraic numbers. In contrast, all previous algorithms require a linear algebra step.

Finally, we improve on the computation of an $S$-normal presentation.

Together with the above lemma, this allows for the rapid computation of a normal presentation: Choose $\alpha\in A$ at random until a normal presentation is obtained. It can be shown that, unless $a$ involves many ideals of small norm, this is very efficient.

To illustrate the speed of our algorithm, we created a list of ideals of bounded norm (here: 400) and took random products of 100 ideals in the field defined by the polynomial $X^{n}+2$ for $n=16,32,64,128$. The results are presented in Table 3. Times are given in seconds.

One is often forced to study algebraic number fields in a real or complex setting, e.g. embeddings into an algebraically closed field, or computing Dedekind zeta functions. This can be due to an intrinsic property of the problem, or it may be the fastest (known) way to solve the problem. Either way, the price is working with approximations. We give a few examples of this and show how the ball arithmetic provided by Arb is employed in Hecke for this purpose.