Finding Nemo

September 25, 2015

Nemo is a new computer algebra package for the Julia programming language. The mastermind behind Nemo is Bill Hart (lead developer of FLINT), with contributions from other people in the computer algebra group at TU Kaiserslautern (Tommy Hofmann, Claus Fieker, Oleksandr Motsak). I've also had a small hand in its development so far.

The first public release of Nemo, version 0.3, was made today. Yesterday, Bill Hart presented Nemo at the annual meeting of the German Mathematical Society: his slides (PDF) give an idea of the motivation behind Nemo. A more detailed introduction to Nemo, including a discussion of the pros and cons of Julia as a language for computer algebra, can be found in the Nemo manual (PDF).

A lot of features are still missing (notably much of the linear algebra), and there will certainly be some bugs. If you decide to give Nemo a try, please don't hesitate to leave some feedback on the GitHub issue tracker or the Nemo development mailing list.

Nemo basics

The core interface is very similar to SageMath. You first construct the mathematical ring you want to work with, and then you are free to play with elements in that ring. For example, suppose we want to generate Bernoulli polynomials $B_n$, which may be defined by the generating function identity $$\frac{t e^{xt}}{e^t-1}= \sum_{n=0}^\infty B_n(x) \frac{t^n}{n!}.$$ We first construct the ring where the Bernoulli polynomials live, $R = \mathbb{Q}[x]$, and then we construct the truncated power series ring $S = R[[t]] / \langle t^N \rangle$. We can evaluate the left-hand side of the above equation in $S$ and then read off the coefficients appearing on the right-hand side. The precision $N = 12$ gives us $B_0, \ldots, B_{10}$ (we lose one coefficient of precision in the power series division $(t + \ldots) / (t + \ldots)$).

julia> using Nemo

julia> R, x = PolynomialRing(QQ, "x");

julia> S, t = PowerSeriesRing(R, 12, "t");

# print B_10
julia> @time coeff(divexact(t * exp(x*t), exp(t)-1), 10) * fac(10)
  0.000137 seconds (369 allocations: 17.859 KB)
x^10 - 5*x^9 + 15/2*x^8 - 7*x^6 + 5*x^4 - 3/2*x^2 + 5/66

Nemo supports most of the standard mathematical domains, including integers, rationals, integers mod $n$, $p$-adic numbers, finite fields (prime and non-prime order), number fields, real and complex numbers. The base types are implemented by wrapping FLINT and other C libraries (currently Antic, Arb, Pari). Over most of these base types, Nemo supports constructing polynomials, power series, fraction fields, residue rings and matrices. Such structures can also be built recursively, allowing you to work with matrices of multivariate polynomials over residue rings over number fields, for example.

Many important "derived" domains automatically use highly optimized C implementations instead of generic code. For example, $\mathbb{Q}[x]$ and $\operatorname{Mat}_{m,n}(\mathbb{Z})$ respectively use FLINT's fmpq_poly and fmpz_mat types. However, performance in generic domains should also be quite good across the board. This has been an important consideration in the design of Nemo. A contributing factor is the Julia language's builtin JIT compilation and language design which helps writing performant code. The Nemo benchmarks page shows a few examples.

Arb transcendental functions in Julia

Nemo includes rudimentary wrapper code for Arb, which I've helped work on the last week. In Nemo's future, Arb will allow doing computer algebra over the base fields $\mathbb{R}$ and $\mathbb{C}$ in a rigorous way (for example, Arb can be used to rigorously isolate complex roots of polynomials, to compute height bounds, to provably deduce integer values from transcendental analytic formulas, etc.). At the moment real and complex numbers work, but polynomials, matrices and power series have not yet been implemented. You can currently at least do basic arithmetic and compute various transcendental functions (with rigorous error bounds).

julia> CC = AcbField(128)
Nemo.AcbField(128)   # complex field with 128 bits of precision

julia> @time besselj(CC(2,3), CC(3,4))  # J(2+3i, 3+4i)
  0.000149 seconds (7 allocations: 496 bytes)
[0.440142369123148468800270811591694988 +/- 3.48e-37] + i*[-0.117229190012375135780993149448690866 +/- 5.23e-37]

This feature could be useful for numerical Julia users who need transcendental functions with better accuracy or on a larger domain than Julia's builtins. For example, Julia's builtin Bessel functions only support real orders. We can easily write a wrapper around Arb for complex floating-point order and argument with guaranteed 53-bit accuracy for any input. Here is the function $K_{\nu}(z)$ as an example:

function besselk2(nu, z)
    nu = Complex128(nu)
    z = Complex128(z)
    prec = 64
    while true
        RR = ArbField(prec)
        CC = AcbField(prec)
        nu2 = CC(RR(real(nu)), RR(imag(nu)))
        z2 = CC(RR(real(z)), RR(imag(z)))
        t = besselk(nu2, z2)
        if accuracy_bits(t) >= 53
            return Float64(real(t)) + Float64(imag(t))*im
        end
        prec = prec * 2
    end
end

julia> @time besselk2(20+30im, 50-60im)
  0.000409 seconds (31 allocations: 1.813 KB)
8.686714193758292e-27 + 4.324369099874024e-28im

Pardon the clumsy conversion code in the above example, which is needed since some helper methods are missing. At the moment, BigFloat conversions have not been implemented. Once this is done, it will be just as easy to provide BigFloat transcendental functions along the same lines.

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