Number Theory 2 available soon!

This is pretty exciting.

I first read the first one of these about 10 years ago. The second one has finally been translated, and looks like publication is imminent. I hope the 3rd one doesn’t take as long!

What I like about these books is how poetic they are. They really expose the poetry of number theory. They don’t always give you the full details, often delegating them to 3rd party references, but they do give you that lovely beauty of the subject. Very motivating texts. I’m currently enjoying the first one.

The Japanese have an exquisite sense of aesthetics.

Complex gcd in Octave

So I wrote a couple of patches for Octave’s gcd. Jaroslav Hájek rewrote most of a bare bones implementation by David Bateman, but omitted Gaussian integers from the implementation. Octave has integral real types, but only float types for complex numbers, because for some inscrutable reason, C++ itself specifies that std::complex<T> is UB if T is not a float type (i.e., not one of float, double, or long double). Additionally, there were a couple of bugs that I also squashed with my patch.

I’m very happy because this contribution puts me in the 36th place of all of 133 contributors, but head Octave honcho John W. Eaton (jwe) also added my name to the contributor’s list. In other words, this finally puts me somewhere in the Octave map.

So let’s take commemorate the occasion by remembering how to compute the gcd of two Gaussian integers. Let’s begin with the ordinary gcd. Given two integers $$a$$ and $$b$$, their greatest common divisor, denoted $$\gcd(a,b)$$, is an integer $$d$$ such that $$d|a$$ and $$d|b$$ and if for any other $$k$$ such that $$k|a$$ and $$k|b$$, then $$k|d$$. In other words, exactly what it says on the tin, the largest number that divides both $$a$$ and $$b$$.

Now, for the Gaussian integers, i.e. complex numbers with integer real and imaginary parts, it turns out that this also makes sense if by “largest” or “greatest” you mean “largest or greatest in complex norm or absolute value” (technically, we should have also specified that we meant absolute value when we were talking about ordinary integers). The trusty Euclidean algorithm to obtain the gcd of ordinary integers also works for Gaussian integers, because given complex integers $$a$$ and $$b$$ with $$|a| < |b|$$ , you can find a quotient $$q$$ and remainder $$r$$ such that $$a = bq +r$$ with $$|r| < |b|$$. So the basic idea is that you just iterate this division, with the remainders getting smaller, $$|a| > |b| > |r| > …$$ until you get the first nonzero remainder, and that’s your gcd.

As to how to perform the actual division, it suffices to do ordinary complex division and then just round real and imaginary parts to their closest integer in order to get the quotient $$q$$ as above. The remainder $$r$$ is then found by $$r = a – qb$$. You don’t need to round to the closest integer point, since it suffices to pick $$q$$ such that $$|a – bq| < 1/2$$, but rounding is as good a choice as any, as far as I can tell. This does have the peculiarity of picking a particular gcd for the Gaussian integers, since if $$d = \gcd(a,b)$$, then $$-d, id, -id$$ are also gcd's, where $$i$$ is the imaginary unit. This actually also happens in ordinary gcd because $$d$$ and $$-d$$ are both gcds, but there you do have a positive and negative numbers, so you can make an arbitrary choice and pick the positive one. And no, there aren't positive or negative complex numbers, despite appearances. That concept just doesn't make sense for them. I am not sure the division algorithm as described above with rounding is the best way to do it, but I'm not concerned, because it's for Octave, not for a number theory package. I am now mildly curious to investigate what LiDIA does here. Oh, I should also mention, there are many many other rings where the gcd makes sense, such as for polynomials, or in $$\mathbb{Z}[\omega]$$ (Eisenstein integers, the integers plus a cube root $$\omega$$ of $$1$$), but Octave doesn't really have a convenient way to represent other rings. In fact, I'm a little surprised that Octave has a gcd at all, and I doubt many people wanted a complex gcd. Regardless, should anyone try it with complex numbers, they'll see that at least we think of everything.

A Personal Introduction to LiDIA

LiDIA is a very old C++ library for computational number theory. The copyright statements inside its source say it’s from 1994. I’m not sure if it’s any good anymore, but it probably isn’t the best one out there, although at one time, it did have the best implementation anywhere for doing arithmetic with elliptic curves.

By the way, for those not in the know, number theory is the part of mathematics that studies the integers (whole numbers) and their properties, like divisibility and distribution of primes. It’s admittedly a very broad discipline, and it does at times feel like it’s a patchwork of many different tricks, but number theory does have its own particular flavour. Pari/GP is just about the only other free software package that concentrates on number theory. In fact, the originators of both, Henri Cohen and Thomas Papanikolaou, appear to be close collaborators. Pari’s development has also slowed down but hasn’t stopped, and it helps that it’s an integral component of Sage, so it is getting attention from the mathematical free software community.

I remember reading in some Sage mailing list a very scathing review of LiDIA. It berates LiDIA for being slow, incomplete, and awkward.

LiDIA was officially abandoned in February 2010, and I had been nagging its authors to release it under the GPL for a while. It had been under some awkward makeshift noncommercial license before that. However, this was in 2006, and interest in LiDIA had already almost died by then. Although permission to release under the GPL was granted, nobody seemed to care enough about it to actually make an official free release. Finally, Christoph Ludwig, who was LiDIA’s maintainer for a long time, got around in 2010 to actually put the GPL notice in LiDIA’s sources. I nagged him about adding an “or later” clause to the version (Sage uses GPLv3, so it would need it), which he did. And then officially abandoned LiDIA.

Soon after that, I grabbed the LiDIA subversion repository, converted it to Mercurial, and uploaded it to bitbucket. So how is it doing?

Well… it compiles with a modern gcc, so that’s something. I took the liberty of squashing some bugs by compiling with -Wall and -Wextra warning flags. There were silly bits of code with signed and unsigned integer types like

u--;
if(u < 0)
  u = u + prime;

but u was an unsigned type. :-/ In fact, LiDIA overall plays fast and loose with signed and unsigned type. It has a deceptively-named lidia_size_t typedef that would make you think that it’s some unsigned integral type, like C++’s standard size_t, but no, it’s a signed integral type, and it was in fact just int. The code supposedly would let you typedef it to be long, which I did, but that broke many other things elsewhere. I also tried for about a day to make it semantically correct and make it an unsigned type, but that turned into a cascade of trouble, since a lot of code explicitly requires this typedef to be signed.

Anyways. I also got around, thanks to the suggestion of Matin Ettl, to fix more warnings in the code that Martin Ettl’s cppcheck program found. There weren’t many, only two more about dangerous calls to sprintf, and one about an uninitialised variable.

So where do I plan to go from here? This is just about getting the damn thing to compile, not yet about the actual worth of its code. Well, I plan to keep on plowing with this for a while. I want to be fairly public about my work on LiDIA, with the hope that I may catch someone’s attention (in fact, this is how I caught Martin Ettl’s attention and got a suggestion from him). My immediate plans are

  1. Move the build system to Cmake, because I just don’t get GNU autotools.
  2. Inline the documentation with Doxygen. It’s a blessing that LiDIA actually is extensively documented, even if not always accurately. The problem is that the documentation is separately maintained in LaTeX files, and I’m sure this has contributed to the problem of it going out of synch with the actual sources, so that’s the motivation for this task.
  3. Write some unit tests for it. As it is, right now all I know is that LiDIA compiles, but I have no idea if it’s any good for mathematics.
  4. Squash more bugs as I find them.

I used LiDIA myself in 2003 for an undergraduate research program related to counting points on algebraic curves. I still have that code and it still runs on the current LiDIA sources, but whether it’s accurate or fast, I wouldn’t know. It’s partly due to this youthful infatuation with LiDIA, but also because despite everything, I still like C++, and I do think it’s still the best language in which to write conceptually complex code while staying close to the machine. Furthermore, I think it’s a good opportunity to read about number theory and its techniques, that is, I’m also doing this for my own edification. Lastly, I’m doing this because it might benefit someone in ways I haven’t foreseen.

At any rate maintaining an abandoned C++ number theory library certainly won’t hurt anyone, so why not?