Correction of lcm() for Array arguments

[cyanescent]

This is a draft (1st PR attempt!) to try to correct the wrong implementation of the lcm() function with Array argument, cf. bug #55379.

The init argument in the reduce() function was incorrect. It is sufficient to initialize with any value in the Array remove the init argument.
I have also added examples in the doc string, but I don't know if I did it the right way (Array argument may need a separate entry). Also, I don't know yet how to write a test.

Fixes #55379
Fixes #56166

[giordano]

@cyanescent thanks for your contribution! When you propose a change which includes a bugfix, could you please always add a relevant test, to demonstrate the bug is indeed fixed and there won't be regressions in the future? Thanks!

[cyanescent]

Why do you want to make lcm and gcd throw for empty arrays? I think the current behavior for empty input is fine. It agrees with the Fricas CAS, at least. EDIT: it agrees with Giac, too.

Interesting, I guess they have not implemented it for Rational. Or am I missing something?

In my understanding, the lcm and the gcd are functions defined on sets of rational numbers, and related by an inversion of their arguments (with neutral element 1). With this in mind, it makes no sense to define the degenerate cases for one of them as 1 and the other one as 0. If no error is thrown and that the value is not NaN (could be a nice solution too), then there are not many possibilities: either 1 for both or Inf and 0.

In practice, we have just seen that lcm(Int[])=1 as implemented here with the reduce() function yields unacceptable results for the type Array{Rational}. This would make us lean towards lcm(Int[])=Inf. Problem: there is no Inf in the type Int.
1//0 could be a strange solution, returning a Rational. But since gcd() and lcm() throw an error already, why not gcd([]) and lcm([])?

EDIT: it could be a cleaner error though. The proposed patch throws trying to evaluate an empty array at index 1.

[cyanescent]

The multiplicative identity is the lcm identity, but not the gcd identity.

That is true for integers, but not for Rationals:

julia> lcm(1//2,1)
1//1

julia> gcd(1//2,1)
1//2

I have just tried FriCAS, and I am very puzzled by it's seemingly wrong behaviour for rationals:

(1) -> lcm(1/2,1)

        1
   (1)  -
        2
                                                      Type: Fraction(Integer)
(2) -> gcd(1/2,1)

   (2)  1
                                                      Type: Fraction(Integer)

[nsajko]

Can you also add these lines to the end of the PR message, to make Github link to these issues properly:

Fixes #55379

Fixes #56166

You'll also need to add some tests, I'll try to help you with that soon.

[nsajko]

This test needs to be adjusted or removed:

julia/test/intfuncs.jl

Line 181 in 159adbf

@test lcm(T[]) ⟷ T(1)

[nsajko]

@cyanescent: I just realized you made this PR from your master branch. The correct procedure is to create a new branch for each PR (feature branches). Could you do either of these:

• Fix the branch for this PR (not sure if this is possible TBH)
• Create a new PR with the same commits but a new branch. The required git commands would be something like this, for example:

  • git checkout -b lcm_rationals_fix  # create new branch
    git checkout master  # switch back to master
    git rebase -i f3a36d74eeb1f8  # delete the PR commits from master interactively
    git checkout lcm_rationals_fix  # switch back to PR branch
    git push origin lcm_rationals_fix  # push the PR branch to your fork

[fingolfin]

I understand the change for lcm, but why also change this in gcd, it works fine, doesn't it?

Suggested change

	gcd(abc::AbstractArray{<:Real}) = reduce(gcd, abc)
	gcd(abc::AbstractArray{<:Real}) = reduce(gcd, abc; init=zero(eltype(abc)))

[cyanescent]

Right, but does it make sense?

EDIT: it makes sense if we also define:

lcm(abc::AbstractArray{<:Real}) = reduce(lcm, abc; init=1//zero(eltype(abc)))

or just lcm(T[]) = 1//zero(T), to prevent returning Vector{Rational} when the input is Vector{Int} for non-empty arrays.
cf. #55379 (comment)

[LilithHafner]

gcd(T[]) = zero(T) is correct.

From the docstring,

Greatest common (positive) divisor (or zero if all arguments are zero).

In this case, all arguments are zero (vacuously) so the answer is explicitly documented to be zero.

[LilithHafner]

Suggested change

	lcm(abc::AbstractArray{<:Real}) = reduce(lcm, abc)
	lcm(abc::AbstractArray{<:Integer}) = reduce(lcm, abc, init=one(eltype(abc)))
	lcm(abc::AbstractArray{<:Real}) = reduce(lcm, abc)

For integers, 1 is an identity of LCM (because all integers are a multiple of 1) and is also the correct answer for the empty case (see #55379 (comment)).

[LilithHafner]

Thanks for finding and fixing this bug, @cyanescent!

Right now, according to my understanding of the documented behavior, this PR makes lcm & gcd throw more than they need to which is breaking, so we'll need to address that before merging.

Specifically, this makes gcd(Rational{Int}[]) and lcm(Int[]) throw when the previously returned answers (0//1 and 1) were correct.

[nsajko]

@cyanescent are you still interested in this? Have you tried creating a new PR with a proper "feature" branch?

[LilithHafner]

It's not standard but it is perfectly fine to make a PR from a master branch. See, for example, #38823, #39002, #39147, #39312, #39331, #39654, #39710, #39737, #39850, #39897, #40253, #40391, #40456, #40633, #40635, #40668, #40806, #40830, #40942, #41045, #41276, #41622, #41758, #41847, #41924, #42140, #42261, #42385, #42500, #42501, #42543, #42661, #43111, #43346, #43427, #43523, #43674, #43699, #43966, #43986, #44376, #44412, #44566, #44709, #44760, #44805, #44999, #45228, #45392, #45793, #46033, #46183, #46410, #47070, #47390, #47564, #47872, #47875, #48073, #48289, #48358, #48536, #48730, #48814, #48973, #49158, #49192, #49343, #49808, #50123, #50126, #50688, #50774, #50843, #51055, #51547, #52252, #52441, #53378, #53592, #53767, #53807, #53861, #53873, #53954, #53981, #54149, #54206, #54283, #54358, #54751, #54955, #54983, #55173, #55315, #56000, #56008, #56122, #56202.

Three's no need to close this PR and restart with a "proper feature branch" or do any git or github fancy operations.

[nsajko]

OK, I hope this won't confuse NanoSoldier, then.

Here are the tests I wrote:

@testset "`gcd`, `lcm` binary operation properties" begin
    function test_commutativity(f, a, b)
        @test f(a, b) == f(b, a)
    end
    function test_product_property(f, g, a, b)
        @test f(a, b) * g(a, b) == abs(a * b)
    end
    function test_thomaes_function_gcd(f, a, b)
        @test f(a, b) == (abs(a) // denominator(b // a))
    end
    function gcd_array(a, b)
        gcd([a, b])
    end
    function lcm_array(a, b)
        lcm([a, b])
    end
    tested_integers = setdiff(-7:7, 0)
    for (f, g) ∈ ((gcd, lcm), (gcd_array, lcm_array))
        for a ∈ tested_integers
            for b ∈ tested_integers
                test_commutativity(f, a, b)
                test_commutativity(g, a, b)
                test_product_property(f, g, a, b)
                test_thomaes_function_gcd(f, a, b)
            end
        end
        for an ∈ tested_integers
            for ad ∈ tested_integers
                a = an // ad
                for bn ∈ tested_integers
                    for bd ∈ tested_integers
                        b = bn // bd
                        test_commutativity(f, a, b)
                        test_commutativity(g, a, b)
                        test_product_property(f, g, a, b)
                        test_thomaes_function_gcd(f, a, b)
                    end
                end
            end
        end
    end
end

Probably also something like:

@testset "`lcm` of empty" begin
    T = Rational{Int}
    for emp ∈ (Vector{T}(undef, 0), Matrix{T}(undef, 0, 0))
        @test_throws ArgumentError lcm(emp)
    end
end