Skip to content

Commit 8820170

Browse files
Pramodh Gopalan VAmit Shirodkar
Pramodh Gopalan V
authored and
Amit Shirodkar
committed
LinearAlgebra: Add bareiss det for BigInt Matrices (#40128) :: Take 2 (JuliaLang#40868)
1 parent fc2cad8 commit 8820170

File tree

3 files changed

+58
-0
lines changed

3 files changed

+58
-0
lines changed

NEWS.md

+1
Original file line numberDiff line numberDiff line change
@@ -108,6 +108,7 @@ Standard library changes
108108
* The shape of an `UpperHessenberg` matrix is preserved under certain arithmetic operations, e.g. when multiplying or dividing by an `UpperTriangular` matrix. ([#40039])
109109
* `cis(A)` now supports matrix arguments ([#40194]).
110110
* `dot` now supports `UniformScaling` with `AbstractMatrix` ([#40250]).
111+
* `det(M::AbstractMatrix{BigInt})` now calls `det_bareiss(M)`, which uses the [Bareiss](https://en.wikipedia.org/wiki/Bareiss_algorithm) algorithm to calculate precise values.([#40868]).
111112

112113
#### Markdown
113114

stdlib/LinearAlgebra/src/generic.jl

+52
Original file line numberDiff line numberDiff line change
@@ -1558,6 +1558,9 @@ function det(A::AbstractMatrix{T}) where T
15581558
end
15591559
det(x::Number) = x
15601560

1561+
# Resolve Issue #40128
1562+
det(A::AbstractMatrix{BigInt}) = det_bareiss(A)
1563+
15611564
"""
15621565
logabsdet(M)
15631566
@@ -1620,6 +1623,55 @@ logdet(A) = log(det(A))
16201623

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

1626+
exactdiv(a, b) = a/b
1627+
exactdiv(a::Integer, b::Integer) = div(a, b)
1628+
1629+
"""
1630+
det_bareiss!(M)
1631+
1632+
Calculates the determinant of a matrix using the
1633+
[Bareiss Algorithm](https://en.wikipedia.org/wiki/Bareiss_algorithm) using
1634+
inplace operations.
1635+
1636+
# Examples
1637+
```jldoctest
1638+
julia> M = [1 0; 2 2]
1639+
2×2 Matrix{Int64}:
1640+
1 0
1641+
2 2
1642+
1643+
julia> LinearAlgebra.det_bareiss!(M)
1644+
2
1645+
```
1646+
"""
1647+
function det_bareiss!(M)
1648+
n = checksquare(M)
1649+
sign, prev = Int8(1), one(eltype(M))
1650+
for i in 1:n-1
1651+
if iszero(M[i,i]) # swap with another col to make nonzero
1652+
swapto = findfirst(!iszero, @view M[i,i+1:end])
1653+
isnothing(swapto) && return zero(prev)
1654+
sign = -sign
1655+
Base.swapcols!(M, i, i + swapto)
1656+
end
1657+
for k in i+1:n, j in i+1:n
1658+
M[j,k] = exactdiv(M[j,k]*M[i,i] - M[j,i]*M[i,k], prev)
1659+
end
1660+
prev = M[i,i]
1661+
end
1662+
return sign * M[end,end]
1663+
end
1664+
"""
1665+
LinearAlgebra.det_bareiss(M)
1666+
1667+
Calculates the determinant of a matrix using the
1668+
[Bareiss Algorithm](https://en.wikipedia.org/wiki/Bareiss_algorithm).
1669+
Also refer to [`det_bareiss!`](@ref).
1670+
"""
1671+
det_bareiss(M) = det_bareiss!(copy(M))
1672+
1673+
1674+
16231675
"""
16241676
promote_leaf_eltypes(itr)
16251677

stdlib/LinearAlgebra/test/generic.jl

+5
Original file line numberDiff line numberDiff line change
@@ -355,6 +355,11 @@ end
355355
@test [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]] [[1,2, [3,4]], 5.0, [6im, [7.0, 8.0]]]
356356
end
357357

358+
@testset "Issue 40128" begin
359+
@test det(BigInt[9 1 8 0; 0 0 8 7; 7 6 8 3; 2 9 7 7])::BigInt == -1
360+
@test det(BigInt[1 big(2)^65+1; 3 4])::BigInt == (4 - 3*(big(2)^65+1))
361+
end
362+
358363
# Minimal modulo number type - but not subtyping Number
359364
struct ModInt{n}
360365
k

0 commit comments

Comments
 (0)