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.

Sign up for GitHub

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

LinearAlgebra: Add bareiss det for BigInt Matrices (#40128) :: Take 2 #40868

Merged
merged 3 commits into from
May 21, 2021
Merged

LinearAlgebra: Add bareiss det for BigInt Matrices (#40128) :: Take 2 #40868

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
2cf1e20
LinearAlgebra: Add bareiss det for BigInt Matrices (#40128)
May 19, 2021
0d1d646
Fix bug, Add better tests
May 20, 2021
8b78248
Update NEWS.md
May 21, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
56 changes: 56 additions & 0 deletions stdlib/LinearAlgebra/src/generic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1558,6 +1558,14 @@ function det(A::AbstractMatrix{T}) where T
end
det(x::Number) = x

# Resolve Issue #40128
function det(A::AbstractMatrix{BigInt})
m, n = size(A)
if m == n || throw(DimensionMismatch("matrix is not square: dimensions are $(size(A))"))
Copy link
Member

@stevengj stevengj May 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like this test should be in det_bareiss!, and then just have:

det(A::AbstractMatrix{BigInt}) = det_bareiss(A)

Note also that you can simply call checksquare(A) in det_bareiss!

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there any reasoning behind this design choice? I'm trying to learn as I go along, so It would be helpful if you could shed some light on this.

I'll add a commit correcting this, thanks!

Copy link
Member

@stevengj stevengj May 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It shouldn't matter whether someone calls det_bareiss! directly, or through det_bareiss, or through det — in all cases, it should check whether the matrix is square, because otherwise that function doesn't make sense.

det_bareiss(A)
end
end

"""
logabsdet(M)

Expand Down Expand Up @@ -1620,6 +1628,54 @@ logdet(A) = log(det(A))

const NumberArray{T<:Number} = AbstractArray{T}

exactdiv(a, b) = a/b
exactdiv(a::Integer, b::Integer) = div(a, b)

"""
LinearAlgebra.det_bareiss!(M)

Calculates the determinant of a matrix using the
[Bareiss Algorithm](https://en.wikipedia.org/wiki/Bareiss_algorithm) using
inplace operations.

# Examples
```jldoctest
julia> M = [1 0; 2 2]
2×2 Matrix{Int64}:
1 0
2 2

julia> LinearAlgebra.det_bareiss!(M)
2
```
"""
function det_bareiss!(M)
dkarrasch marked this conversation as resolved.
Show resolved Hide resolved
n, sign, prev = size(M,1), Int8(1), one(eltype(M))
for i in 1:n-1
if iszero(M[i,i]) # swap with another col to make nonzero
swapto = findfirst(!iszero, @view M[i,i+1:end])
isnothing(swapto) && return zero(prev)
sign = -sign
Base.swapcols!(M, i, swapto)
end
for k in i+1:n, j in i+1:n
M[j,k] = exactdiv(M[j,k]*M[i,i] - M[j,i]*M[i,k], prev)
end
prev = M[i,i]
end
return sign * M[end,end]
end
"""
LinearAlgebra.det_bareiss(M)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't put the LinearAlgebra. in the docstring, I think.

Copy link
Contributor Author

@Pramodh-G Pramodh-G May 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll change it, butLinearAlgebra.jl contains docstrings like

LinearAlgebra.peakflops()

This is the one example I could find though, should I change this too?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would leave the peakflops case out of this PR.


Calculates the determinant of a matrix using the
[Bareiss Algorithm](https://en.wikipedia.org/wiki/Bareiss_algorithm).
Also refer to [`det_bareiss!`](@ref).
"""
det_bareiss(M) = det_bareiss!(copy(M))



"""
promote_leaf_eltypes(itr)

Expand Down
4 changes: 4 additions & 0 deletions stdlib/LinearAlgebra/test/generic.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -355,6 +355,10 @@ end
@test [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]] ≈ [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]]
end

@testset "Issue 40128" begin
@test LinearAlgebra.det_bareiss(BigInt[1 2; 3 4]) == -2
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fix the indenting here. (4 spaces)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also check that the result is a BigInt

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Note that det([1 2; 3 4]) == -2 returns true as well, so this is not a great test.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Specifically, you should probably test something like BigInt[1 big(2)^65+1; 3 4]) which will ensure that no Float based method would produce the right answer.

Copy link
Member

@stevengj stevengj May 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You could try det(BigInt[9 1 8 0; 0 0 8 7; 7 6 8 3; 2 9 7 7])::BigInt == -1, since for Int this is off by about 1e-12 (and even for BigInt it currently gives a non-integer result).

Copy link
Member

@stevengj stevengj May 19, 2021 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Specifically, you should probably test something like BigInt[1 big(2)^65+1; 3 4]) which will ensure that no Float based method would produce the right answer.

Actually the current LU algorithm gives the exact answer here:

julia> det(BigInt[1 big(2)^65+1; 3 4]) - (4 - 3*(big(2)^65+1))
0.0

(Remember that it uses BigFloat.)

end

# Minimal modulo number type - but not subtyping Number
struct ModInt{n}
k
Expand Down