Greatest Common Gaussian Divisor

Gaussian integers are complex numbers $$\x+yi\$$ such that $$\x\$$ and $$\y\$$ are both integers, and $$\i^2 = -1\$$. The norm of a Gaussian integer $$\N(x+yi)\$$ is defined as $$\x^2 + y^2 = |x+yi|^2\$$. It is possible to define a Euclidean division for Gaussian integers, which means that it is possible to define a Euclidean algorithm to calculate a greatest common divisor for any two Gaussian integers.

Unfortunately, a Euclidean algorithm requires a well-defined modulo operation on complex numbers, which most programming languages don't have (e.g. Python, Ruby), meaning that such an algorithm fails.

Gaussian division

It is possible to define the division $$\\frac a b = x+yi\$$ (where $$\a\$$ and $$\b\$$ are both Gaussian integers) as finding a quotient $$\q\$$ and a remainder $$\r\$$ such that

$$a = bq + r, \text{ and } N(r) \le \frac {N(b)} 2$$

We can further limit this to $$\q = m + ni\$$, where $$\-\frac 1 2 < x - m \le \frac 1 2\$$ and $$\-\frac 1 2 < y - n \le \frac 1 2\$$, and $$\r = b(x - m + (y - n)i)\$$

From here, a Euclidean algorithm is possible: repeatedly replace $$\(a, b)\$$ by $$\(b, r)\$$ until it reaches $$\(d, 0)\$$. $$\d\$$ can then be called the greatest common divisor of $$\a\$$ and $$\b\$$

Complex GCDs are not unique; if $$\d = \gcd(a, b)\$$, then $$\d, -d, di, -di\$$ are all GCDs of $$\a\$$ and $$\b\$$

You are to take 2 Gaussian integers $$\a, b\$$ as input and output $$\\gcd(a, b)\$$. You may take input in any convenient method, and any reasonable format, including two complex numbers, two lists of pairs [x, y], [w, z] representing $$\\gcd(x+yi, w+zi)\$$ etc. Additionally, the output format is equally lax.

You may output any of the 4 possible values for the GCD, and you don't need to be consistent between inputs.

If your language's builtin $$\\gcd\$$ function already handles Gaussian integers, and so would trivially solve this challenge by itself, please add it to the Community Wiki of trivial answers below.

This is , so the shortest code in bytes wins.

Test cases

5+3i, 2-8i    -> 1+i
5+3i, 2+8i    -> 5+3i
1-9i,  -1-7i  -> 1+i
-1+0i, 2-10i  -> 1+0i  (outputting 1 here is also fine)
4+3i,  6-9i   -> 1+0i  (outputting 1 here is also fine)
-3+2i, -3+2i  -> 2+3i
-6+6i, 3+5i   -> 1+i
4+7i, -3-4i   -> 2+i
-3+4i, -6-2i  -> 1+2i
7-7i, -21+21i -> 7+7i


Python 3, 70 bytes

f=lambda a,b:b and f(b,a-b*((t:=a/b+.5+.5j).real//1+t.imag//1*1j))or a


Please add your language to this if the builtin GCD function correctly handles Gaussian integers.

J, 2 bytes

+.


Try it online!

Wolfram Language (Mathematica), 3 bytes

GCD


Try it online!

Pari/GP, 3 bytes

gcd


Try it online!