Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

lcm(::Vector{<:Rational}) gives the wrong answer. #55379

Open
cyanescent opened this issue Aug 5, 2024 · 9 comments · May be fixed by #56113 or #56423
Open

lcm(::Vector{<:Rational}) gives the wrong answer. #55379

cyanescent opened this issue Aug 5, 2024 · 9 comments · May be fixed by #56113 or #56423
Labels
correctness bug ⚠ Bugs that are likely to lead to incorrect results in user code without throwing maths Mathematical functions

Comments

@cyanescent
Copy link

The extension of lcm() and gcd() to both Array and Rational was incorrect. The case of the gcd has been corrected in #34417, but not the lcm:

lcm([1//2;1//2])
1//1 # false 

lcm(1//2,1//2)
1//2 # true
@barucden
Copy link
Contributor

barucden commented Aug 5, 2024

The problem shows up even for a singleton array:

julia> lcm([1//2])
1//1

See the implementation:

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

(The problem is with the init argument; in this case, init is 1//1 and lcm(1//1, 1//2) is 1//1.)

Neither lcm nor gcd is documented to accept an array argument by the way. We only have

gcd(x, y...)

Greatest common (positive) divisor (or zero if all arguments are zero). The arguments may be integer and rational numbers.

and

lcm(x, y...)

Least common (positive) multiple (or zero if any argument is zero). The arguments may be integer and rational numbers.

What should the methods return when the input array is empty? Currently,

julia> gcd(Int[])
0

julia> lcm(Int[])
1

@cyanescent
Copy link
Author

cyanescent commented Aug 5, 2024

What should the methods return when the input array is empty? Currently,

julia> gcd(Int[])
0

julia> lcm(Int[])
1

For consistency, we should have gcd(Int[])=0 and lcm(Int[])=Inf, but Inf is not available for type Integer. Another possibility is gcd(Int[])=0//1 and lcm(Int[])=1//0, but the type is now Rational. Anyway, for such a degenerate case, I think NaN is the most sensible definition in both cases.

@barucden
Copy link
Contributor

barucden commented Aug 6, 2024

If #55318 lands, I think we could just remove the init argument of reduce for both gcd and lcm, in which case gcd([]) and lcm([]) would result in an error. I am not sure if that's (too) breaking.

@nsajko nsajko added the maths Mathematical functions label Aug 6, 2024
@mbauman mbauman added the correctness bug ⚠ Bugs that are likely to lead to incorrect results in user code without throwing label Aug 6, 2024
@cyanescent cyanescent linked a pull request Oct 11, 2024 that will close this issue
@cyanescent
Copy link
Author

cyanescent commented Oct 11, 2024

Here is a correction draft, that would make gcd(Int[]) and lcm(Int[]) error (perhaps not in the most elegant way though better now): #56113.

This is not breaking since the Array argument is still undocumented, and gcd() and lcm() also result into an error (so this is coherent). The current behaviour is also quite wrong.

@fingolfin
Copy link
Contributor

You could also define gcd(T[]) as zero(T) and lcm(T[]) could throw an error.

@cyanescent
Copy link
Author

I get your point, the inversion symmetry between gcd and lcm is broken anyway since there is no Inf in type Int.
But then shouldn't gcd() also return 0 instead of an error?
(also, I don't know any reason for choosing this value rather than NaN, my reasoning above was only deductive)

@cyanescent
Copy link
Author

cyanescent commented Oct 18, 2024

(also, I don't know any reason for choosing this value rather than NaN, my reasoning above was only deductive)

Answering this question (what are $\text{gcd}$ and $\text{lcm}$ for empty collections of numbers?):
The extension to collections of numbers is introduced recursively through associativity: $\text{gcd}(n_1,n_2,\dots,n_r)=\text{gcd}(n_1,\text{gcd}(n_2,\dots,n_r))=\text{gcd}(\text{gcd}(n_1,n_2,\dots,n_{r-1}),n_r)$ and $\text{lcm}(n_1,n_2,\dots,n_r)=\text{lcm}(n_1,\text{lcm}(n_2,\dots,n_r))=\text{lcm}(\text{lcm}(n_1,n_2,\dots,n_{r-1}),n_r)$. Thus, we should define the case $r=1$ so that this still holds.

Therefore, gcd([])=0 and perhaps gcd()=0 would result from the common (but somewhat peculiar) convention $\text{gcd}(n_1,0)=n_1$ and the associativity rule.
The consequential rule $\text{lcm}(n_1,\infty)=n_1$ is also defined in julia for type Rational:

julia> lcm(3,1//0)
3//1

leaning towards lcm([])=1//0 and perhaps lcm()=1//0, or at least lcm(Rational[])=1//0 (if the output type Rational is an issue).

You could also define gcd(T[]) as zero(T) and lcm(T[]) could throw an error.

In the end, lcm(Rational[])=1//0 should also be defined in this case, leading to the question:
Is lcm(T[])=1//zero(T) acceptable from a programmatic viewpoint?

@LilithHafner LilithHafner changed the title False array-argument lcm for non-integers lcm(::Vector{Rational}) gives the wrong answer. Oct 22, 2024
@LilithHafner LilithHafner changed the title lcm(::Vector{Rational}) gives the wrong answer. lcm(::Vector{<:Rational}) gives the wrong answer. Oct 22, 2024
@LilithHafner
Copy link
Member

lcm(Int[]) should be 1. In this case the range of the function is Int or the integers and in either case, all elements of the range are multiples of every element of the array (this is true vacuously) so we are looking for the least positive element of the range. That's 1.

In the case of lcm(Rational{Int}[]), the range is Rational{Int} or the rationals. Again, all elements of the range are multiples, but in this case there is no least positive element so we should throw.

@LilithHafner LilithHafner linked a pull request Nov 2, 2024 that will close this issue
@adienes
Copy link
Contributor

adienes commented Nov 3, 2024

I don't think it has to throw necessarily, although I'd be fine with that.

My preference would be to use the FTA definition

aka where n = prod(p_i ^ e_i) , allowing negative exponents over rationals, gcd(n1, n2) = prod(p_i ^ min(e_i1, e_i2)) and lcm(n1, n2) = prod(p_i ^ max(e_i1, e_i2)) which preserves the identity gcd(a, b) * lcm(a, b) = a * b for positive a, b

by convention gcd, lcm are defined only up to a unit, but obviously in the rationals all values are unitary so I think it would be a common enough convention to simply make the functions even. together with the common convention that gcd(a, 0) = a and lcm(a, 0) = 0 and commutativity this fully determines all outputs. I'll just put a list of what this would give over all the examples I've seen scattered over the PRs and mark if any deviate from status quo

  • lcm(1//2, 1//2) = 1//2
  • lcm([1//2, 1//2]) = 1//2
    • change from status quo as a bugfix (this issue)
  • lcm(Int[]) = 1, lcm(Rational[]) = 1//1
  • gcd(Int[]) = 1, gcd(Rational[]) = 1//1
    • change from status quo. I suggest this is more accurate as it matches prod(Int[]) = 1
  • gcd(3//7, 6//7) = 3//7

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
correctness bug ⚠ Bugs that are likely to lead to incorrect results in user code without throwing maths Mathematical functions
Projects
None yet
7 participants