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.