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.

Linked on Reddit

Oh, dear me.

A few years ago, I wrote a strongly-worded opinion piece on just why I believe so strongly on free software. Today it was linked on Reddit. While I do want free software everywhere, and I would like to completely abolish the notion of proprietary software, I see that as a distant dream. My immediate goal, my highest priority is to have free software for mathematics, because I think that hiding source and restricting distribution are the antithesis of mathematical work. Others have other priorities, but for me, free software in mathematics is highest.

In this vein, I wrote this some time ago. It’s a piece many people have read, and it always seems many have an opinion about. Many love it, but also many strongly disagree with it. They cite objections that say that the ends justifies that means, that having a good software package for doing mathematics trumps any need to have its source, that it’s ok to hide source if you’re making money, and that without the money that restricting distribution provides, you wouldn’t have quality software at all.

Needless to say, I believe all of these objections to be false if stated absolutely, as there are plenty of counterexamples to each, of varying degree. Nevertheless, I don’t have the energy to individually present them again, nor to address my detractors in the Reddit comments. I’m both sad and happy when I read those comments, and I wonder which direction the discussion overall will take. As for that beacon in the free software movement, Richard Stallman, whatever you might think of him, you have to admire him for his boundless stamina to be able to tirelessly participate in these discussions, and for being so self-consistent both in words and actions about what software freedom should be.

For the moment, I do not feel like arguing about this anymore. I feel like actions. Talk is cheap. I’ll show you the source. Perhaps we will never amount to anything, perhaps Octave will forever lag behind Matlab in one way or another, perhaps Sage will forever be a patchwork of individual programs of varying quality. Or perhaps not. Regardless, I don’t care. We must proceed. I must proceed. We must do, we must keep faith, because it’s obvious that even if the goal is unattainable, it’s the right goal. I would hardly be a hacker-errant of the rueful countenance if I didn’t think so.

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

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?

Salve, cronica electronica!

Hello, blog!

Hopefully this post is the first of many. I plan to write about several things, both personal and more public/professional.

I must say, WordPress was quite easy to set up, but I’m still not sold on PHP. I had to patch a few things here and there, in particular with this nifty $$\LaTeX$$ plugin that I plan to put to good use. Oh, and while we’re talking about meta, thanks also to Tim Saimburg for this free nifty blog design. I tweaked it a little for some spacing, colours, and font sizes, but just barely. Thanks!

It’s a Saturday afternoon, and I’m tired of tweaking and setting up this blog — which I can barely resist to call “blag”, but I will defer to established use — so I will rest for now and do other things.