Julia: Missing GCD methods on complex numbers

[SyxP]

In Julia 1.3.1, there is no methods gcd(::Complex{<:Integer}, ::Complex{<:Integer}).

However, it is possible to define gcd on any Euclidean domain, of which $\mathbb{Z}[i]$ is one.

[ravibitsgoa]

Hi @SyxP ,
I would like to work on this.
I think this gcd method belongs to intfuncs rather than complex. Am I correct?

[SyxP]

It seems like base/intfuncs.jl is composed of functions for Integer types. So base/complex.jl might be a better home for it. Either way, we can just write the function and some tests first and then decide where it should belong in the PR.

[ravibitsgoa]

For this, I wrote a function and a few tests, but I am wondering about these:
Should we ensure that the resultant gcd is normalized (rotated into the first quadrant)?
E.g. gcd(3+im, 1-im) = 1-im. It can be normalized to 1+im.
As discussed here: https://math.stackexchange.com/questions/1836682/calculating-the-gcd-of-complex-numbers
Also: http://mathforum.org/library/drmath/view/67068.html
(According to this mathforum post, depending on the rounding mode used, we can get 4 different GCDs (unique up to a unit rotation).)
Also, for the numerical stability, during the intermediate computation, I am using Complex{<:Rational} rather than Complex{<:AbstractFloat}. Is that correct?

[SyxP]

It would make sense that we ensure the resulting gcd will lie in the first quadrant (i.e. the argument lies within $[0, \pi/2)$). We currently ensure that gcd(::Integer, ::Integer) is non-negative.

I don't understand why one would require intermediate steps to use Complex{<:Rational}, you should be able to do everything without leaving Complex{<:Integer}

[ravibitsgoa]

Just a little clarification: Since we want the argument within [0, π/2), we should return 1 for gcd(im, 2im) = im = 1 right?

I was using Complex{<:Rational} because I was trying to do z1//z2 directly (during the remainder calculation). Now I changed it to z1*conj(z2)/abs2(z2) and used only Complex{<:Integer}.

Now I am worrying about overflow cases when abs2(z) is too large to fit within a particular type. How do we handle such cases? Currently, abs2 function does not handle overflows.
(E.g. abs2(complex(Int32(10)^5, Int32(10)^5))).

In future, perhaps we can also support gcd(::Rational, ::Rational), since implementing it using gcd(::Integer, ::Integer) is trivial.

[SyxP]

Yes, gcd(im, 2im) == 1. My feeling for this rational is that gcd(a+0im, b+0im) should equal gcd(a,b) for a, b <: Integer. The issue of overflow is indeed a problem, as of now I don't know how resolve immediately.

The support for gcd(::Rational, ::Rational) won't be very interesting. As the rational numbers form a Field, there are only two ideals $(0)$ and $\mathbb{Q}$. So it will just be whether the two inputs are both zero.

[ravibitsgoa]

I have completed the implementation and tests for gcd(::Complex{<:Integer}, ::Complex{<:Integer}).

Although I am not very familiar with Fields and ideals, according to this (https://math.stackexchange.com/questions/145054/relating-operatornamelcm-and-gcd), gcd(a/b, c/d) == gcd(a, c)/lcm(b,d). Wolfram language also gives the same result for gcd of rationals.

Also, I discovered another problem: currently, Complex{Int32} are automatically getting promoted to Complex{Int64} (to avoid overflow) while calculating the gcd. Should we just keep them as Complex{Int32}? (i.e. Should the gcd have the same type as input or should it get an automatic promotion?)

Should we separate the division and remainder calculation into methods rem(::Complex{<:Integer}, ::Complex{<:Integer}) and div(::Complex{<:Integer}, ::Complex{<:Integer})?

Also, how about adding a method checked_abs2(::Complex{<:Integer})?

[Micket]

I would be so inconvenient to say that it should probably mirror the behaviour of Integer types;

julia> typeof(gcd(Int8(4),Int16(2)))
Int16
julia> typeof(gcd(4,2000000000000000000000000000000000000000))
BigInt

But, how does promoting to Int64 really help? It can also overflow (eventually).
I'm a bit curious how your implementation is if it specifically promotes Int32->Int64. What happens with say, Int16 or BigInt?

I think separating out functions div is a good idea (because, why not?).

I experimented a bit with a naive implementation, even looked into using (which might be a terrible idea if precision becomes an issue):

function div(n::Complex{T}, d::Complex{V}) where {T<:Integer, V<:Integer}
    r = n//d
    return complex(div(r.re.num, r.re.den), div(r.im.num, r.im.den))
end

round(::Type{Complex{T}}, r::Complex) where T<:Integer = complex(round(T, r.re), round(T, r.im))

function gcd(a::Complex{T}, b::Complex{V}) where {T<:Integer,V<:Integer}
    R = promote_type(T, V)
    x, r = Complex{R}(a), Complex{R}(b)
    while r != 0
        q = round(Complex{R}, x/r)
        #q = div(x, r) # alteranative
        x, r = r, x - q*r
    end
    return x
end

Naive, not normalized, but it does preserve the return type behaviour of gcd.

[ravibitsgoa]

I agree. Actually, I had tested my code with only Int32 and Int64. Promoting to Int64 just helps with smaller data types, but I think mirroring the behaviour of gcd of Integer is the right thing to do.

For div, I had a bit of confusion because

1. It would require usage of Complex{Rational} in an intermediate step.
2. Also, mathematically division of Gaussian Integers is dependant on the rounding mode used. I think we can make it like division of Integers. div(x, y, r::RoundingMode=RoundToZero)
   In particular, for finding the remainder for GCD, it is necessary to use RoundNearest rounding mode. Otherwise, the Euclidean algorithm is not guaranteed to converge.

I have fixed the return type issue (Int64), and the current code looks like this:

function gcd(z1::Complex{T}, z2::Complex{V}) where {T<:Integer, V<:Integer}
    R = promote_type(T, V)
    a = Complex{R}(z1)
    b = Complex{R}(z2)
    if abs(a) < abs(b)
        a, b = b, a
    end
    while b != 0
        # TODO: Create rem(::Complex{<:Integer}, ::Complex{<:Integer})
        # TODO: Create div(::Complex{<:Integer}, ::Complex{<:Integer})
        b̅ = conj(b)
        # TODO: Handle overflow when calculating a*b̅
        t = a*b̅
        # TODO: Create checked_abs2(::Complex{<:Integer})
        # TODO: Handle overflow when calculating the norm of complex numbers
        abs2_b = abs2(b)
        abs2_b < 0 && __throw_gcd_overflow(z1, z2)
        r = a - b * complex(div(real(t), abs2_b, RoundNearest),
                            div(imag(t), abs2_b, RoundNearest))
        a = b
        b = r
    end
    ar, ai = reim(a)
    if ar == 0
        complex(abs(ai), zero(ar))
    elseif ai == 0
        complex(abs(ar), zero(ai))
    elseif ar>0 && ai>=0   # gcd is already in first quadrant
        a
    elseif ar<0 && ai>0     # In second quadrant
        complex(ai, -ar)
    elseif ar<0 && ai<0     # In third quadrant
        -a
    else                               # In fourth quadrant
        complex(-ai, ar)
    end
end

BTW, it seems that the gcd of integers doesn't have enough test coverage for the mixed data type tests that you proposed. e.g. gcd(Int8(4),Int16(2))

[Micket]

I don't think you actually need

    if abs(a) < abs(b)
        a, b = b, a
    end

It would just change the quadrant in the end, which you check for in the end anyway.

I'd say, definitely introduce rem, which, unless I missed something, would lend itself to a very clean gcd:

function gcd(a::Complex{T}, b::Complex{V}) where {T<:Integer,V<:Integer}
    R = promote_type(T, V)
    x, r = Complex{R}(a), Complex{R}(b)
    while r != 0
        x, r = r, rem(x, r)
    end
    return first_quadrant(x)
end

And then have to deal with problematic design choices in rem.

[StefanKarpinski]

Possibly irrelevant since you might not need to do it at all, but it would be faster to check abs2(a) < abs2(b) instead of abs(a) < abs(b) to avoid the square root operations.

[ravibitsgoa]

Actually, I was avoiding abs2 since abs2(::Complex) overflows, and hence can be slightly problematic at times. I agree that we should introduce rem function.

It seems that Julia already supports gcd(::Rational, ::Rational), but it is not well-documented and tested.

[dpsanders]

It's not clear to me that this needs to be in Base at all. It feels like it could be in a GaussianIntegers package.

[StefanKarpinski]

The issue is that but gcd and Complex{Int} are owned by Base, so that would be type piracy.

[KlausC]

The support for gcd(::Rational, ::Rational) won't be very interesting. As the rational numbers form a Field, there are only two ideals $(0)$ and $\mathbb{Q}$. So it will just be whether the two inputs are both zero.

There is a very sensible (non-tivial) implementation for rationals in place, though. See also https://math.stackexchange.com/questions/151081/gcd-of-rationals

[ravibitsgoa]

Should I add gcdx(::Complex{<:Integer}, ::Complex{<:Integer}) methods as well?

[KlausC]

If you have a concise algorithm - why not?