From fa493b356c3032ac97e965ef7178c9b147aa5aad Mon Sep 17 00:00:00 2001 From: Tommy Hofmann Date: Mon, 12 Feb 2024 10:17:56 +0100 Subject: [PATCH] feat: Change solving interface - Use the `Solve.*` functionality - Hide all previous `*solve*` and `*kernel*` functions using underscores [breaking] --- docs/src/linear_solving.md | 16 +- docs/src/matrix.md | 170 +--------- src/AbstractAlgebra.jl | 25 +- src/Generic.jl | 6 +- src/Matrix-Strassen.jl | 28 +- src/Matrix.jl | 516 +++++++++++++++--------------- src/Module.jl | 4 +- src/ModuleHomomorphism.jl | 4 +- src/Solve.jl | 2 +- src/generic/MatRing.jl | 8 +- src/generic/ModuleHomomorphism.jl | 2 +- src/generic/Submodule.jl | 2 +- test/generic/MatRing-test.jl | 4 +- test/generic/Matrix-test.jl | 124 +++---- test/runtests.jl | 3 +- 15 files changed, 373 insertions(+), 541 deletions(-) diff --git a/docs/src/linear_solving.md b/docs/src/linear_solving.md index aefa68ff25..ed04c6852a 100644 --- a/docs/src/linear_solving.md +++ b/docs/src/linear_solving.md @@ -2,11 +2,7 @@ CurrentModule = AbstractAlgebra.Solve ``` -# Linear solving - -!!! note - - This functionality is experimental and subject to change. +# [Linear solving](@id solving_chapter) ## Overview of the functionality @@ -17,23 +13,21 @@ The module `AbstractAlgebra.Solve` provides the following four functions for sol * `can_solve_with_solution_and_kernel` All of these take the same set of arguments, namely: -* a matrix $A$ of type `MatElem{T}`; -* a vector or matrix $B$ of type `Vector{T}` or `MatElem{T}`; +* a matrix $A$ of type `MatElem`; +* a vector or matrix $B$ of type `Vector` or `MatElem`; * a keyword argument `side` which can be either `:left` (default) or `:right`. If `side` is `:left`, the system $xA = B$ is solved, otherwise the system $Ax = B$ is solved. -For matrices defined over a field, the functions internally rely on `rref`. -If the matrices are defined over a ring, the function `hnf_with_transform` is required internally. The functionality of the functions can be summarized as follows. * `solve`: return a solution, if it exists, otherwise throw an error. * `can_solve`: return `true`, if a solution exists, `false` otherwise. * `can_solve_with_solution`: return `true` and a solution, if this exists, and `false` and an empty vector or matrix otherwise. -* `can_solve_with_solution_and_kernel`: like `can_solve_with_solution` and additionally return a matrix whose columns (respectively rows) give a basis of the kernel of $A$. +* `can_solve_with_solution_and_kernel`: like `can_solve_with_solution` and additionally return a matrix whose rows (respectively columns) give a basis of the kernel of $A$. ## Solving with several right hand sides -Systems $Ax = b_1,\dots, Ax = b_k$ with the same matrix $A$, but several right hand sides $b_i$ can be solved more efficiently, by first initializing a "context object" `C` of type `SolveCtx`. +Systems $xA = b_1,\dots, xA = b_k$ with the same matrix $A$, but several right hand sides $b_i$ can be solved more efficiently, by first initializing a "context object" `C`. ```@docs solve_init ``` diff --git a/docs/src/matrix.md b/docs/src/matrix.md index ecb4de6db4..dbfa768870 100644 --- a/docs/src/matrix.md +++ b/docs/src/matrix.md @@ -959,143 +959,7 @@ julia> pfaffians(M, 2) ### Linear solving -```@docs -solve{T <: FieldElem}(::MatElem{T}, ::MatElem{T}) -``` - -```@docs -solve_rational{T <: RingElem}(::MatElem{T}, ::MatElem{T}) -``` - -```@docs -can_solve_with_solution{T <: RingElement}(::MatElem{T}, ::MatElem{T}) -``` - -```@docs -can_solve{T <: RingElement}(::MatElem{T}, ::MatElem{T}) -``` - -```@docs -solve_left{T <: RingElem}(::MatElem{T}, ::MatElem{T}) -``` - -```@docs -solve_triu{T <: FieldElem}(::MatElem{T}, ::MatElem{T}, ::Bool) -``` - -```@docs -can_solve_left_reduced_triu{T <: RingElement}(::MatElem{T}, ::MatElem{T}) -``` - -**Examples** - -```jldoctest -julia> R, x = polynomial_ring(QQ, "x") -(Univariate polynomial ring in x over rationals, x) - -julia> K, = residue_field(R, x^3 + 3x + 1); a = K(x); - -julia> S = matrix_space(K, 3, 3) -Matrix space of 3 rows and 3 columns - over residue field of univariate polynomial ring modulo x^3 + 3*x + 1 - -julia> U = matrix_space(K, 3, 1) -Matrix space of 3 rows and 1 column - over residue field of univariate polynomial ring modulo x^3 + 3*x + 1 - -julia> A = S([K(0) 2a + 3 a^2 + 1; a^2 - 2 a - 1 2a; a^2 + 3a + 1 2a K(1)]) -[ 0 2*x + 3 x^2 + 1] -[ x^2 - 2 x - 1 2*x] -[x^2 + 3*x + 1 2*x 1] - -julia> b = U([2a, a + 1, (-a - 1)]) -[ 2*x] -[ x + 1] -[-x - 1] - -julia> x = solve(A, b) -[ 1984//7817*x^2 + 1573//7817*x - 937//7817] -[ -2085//7817*x^2 + 1692//7817*x + 965//7817] -[-3198//7817*x^2 + 3540//7817*x - 3525//7817] - -julia> A = matrix(ZZ, 2, 2, [1, 2, 0, 2]) -[1 2] -[0 2] - -julia> b = matrix(ZZ, 2, 1, [2, 1]) -[2] -[1] - -julia> can_solve(A, b, side = :right) -false - -julia> A = matrix(QQ, 2, 2, [3, 4, 5, 6]) -[3//1 4//1] -[5//1 6//1] - -julia> b = matrix(QQ, 1, 2, [2, 1]) -[2//1 1//1] - -julia> can_solve_with_solution(A, b; side = :left) -(true, [-7//2 5//2]) - -julia> A = S([a + 1 2a + 3 a^2 + 1; K(0) a^2 - 1 2a; K(0) K(0) a]) -[x + 1 2*x + 3 x^2 + 1] -[ 0 x^2 - 1 2*x] -[ 0 0 x] - -julia> bb = U([2a, a + 1, (-a - 1)]) -[ 2*x] -[ x + 1] -[-x - 1] - -julia> x = solve_triu(A, bb, false) -[ 1//3*x^2 + 8//3*x + 13//3] -[-3//5*x^2 - 3//5*x - 12//5] -[ x^2 + 2] - -julia> R, x = polynomial_ring(ZZ, "x") -(Univariate polynomial ring in x over integers, x) - -julia> S = matrix_space(R, 3, 3) -Matrix space of 3 rows and 3 columns - over univariate polynomial ring in x over integers - -julia> U = matrix_space(R, 3, 2) -Matrix space of 3 rows and 2 columns - over univariate polynomial ring in x over integers - -julia> A = S([R(0) 2x + 3 x^2 + 1; x^2 - 2 x - 1 2x; x^2 + 3x + 1 2x R(1)]) -[ 0 2*x + 3 x^2 + 1] -[ x^2 - 2 x - 1 2*x] -[x^2 + 3*x + 1 2*x 1] - -julia> bbb = U(transpose([2x x + 1 (-x - 1); x + 1 (-x) x^2])) -[ 2*x x + 1] -[ x + 1 -x] -[-x - 1 x^2] - -julia> x, d = solve_rational(A, bbb) -([3*x^4-10*x^3-8*x^2-11*x-4 -x^5+3*x^4+x^3-2*x^2+3*x-1; -2*x^5-x^4+6*x^3+2*x+1 x^6+x^5+4*x^4+9*x^3+8*x^2+5*x+2; 6*x^4+12*x^3+15*x^2+6*x-3 -2*x^5-4*x^4-6*x^3-9*x^2-4*x+1], x^5 + 2*x^4 + 15*x^3 + 18*x^2 + 8*x + 7) - -julia> S = matrix_space(ZZ, 3, 3) -Matrix space of 3 rows and 3 columns - over integers - -julia> T = matrix_space(ZZ, 3, 1) -Matrix space of 3 rows and 1 column - over integers - -julia> A = S([BigInt(2) 3 5; 1 4 7; 9 2 2]) -[2 3 5] -[1 4 7] -[9 2 2] - -julia> B = T([BigInt(4), 5, 7]) -[4] -[5] -[7] -``` +See [Linear Solving & Kernel](@ref solving_chapter) ### Inverse @@ -1162,38 +1026,6 @@ julia> X, d = pseudo_inv(A) nullspace{T <: FieldElem}(::MatElem{T}) ``` -### Kernel - -```@docs -kernel{T <: RingElem}(::MatElem{T}) -left_kernel{T <: RingElem}(::MatElem{T}) -right_kernel{T <: RingElem}(::MatElem{T}) -``` - -**Examples** - -```jldoctest -julia> S = matrix_space(ZZ, 4, 4) -Matrix space of 4 rows and 4 columns - over integers - -julia> M = S([1 2 0 4; - 2 0 1 1; - 0 1 1 -1; - 2 -1 0 2]) -[1 2 0 4] -[2 0 1 1] -[0 1 1 -1] -[2 -1 0 2] - -julia> nr, Nr = kernel(M) -(1, [-8; -6; 11; 5]) - -julia> nl, Nl = left_kernel(M) -(1, [0 -1 1 1]) - -``` - ### Hessenberg form ```@docs diff --git a/src/AbstractAlgebra.jl b/src/AbstractAlgebra.jl index 199c6de967..c1b0f61680 100644 --- a/src/AbstractAlgebra.jl +++ b/src/AbstractAlgebra.jl @@ -696,6 +696,11 @@ import .Generic: sturm_sequence include("Solve.jl") +import ..Solve: solve +import ..Solve: solve_init +import ..Solve: can_solve +import ..Solve: can_solve_with_solution + # Do not export inv, div, divrem, exp, log, sqrt, numerator and denominator as we define our own export _check_dim export _checkbounds @@ -727,10 +732,10 @@ export basis export block_diagonal_matrix export cached export can_solve -export can_solve_left_reduced_triu -export can_solve_with_kernel +#export can_solve_left_reduced_triu +#export can_solve_with_kernel export can_solve_with_solution -export can_solve_with_solution_interpolation +#export can_solve_with_solution_interpolation export canonical_unit export change_base_ring export change_coefficient_ring @@ -920,7 +925,7 @@ export leading_exponent_vector export leading_exponent_word export leading_monomial export leading_term -export left_kernel +#export left_kernel export leglength export length export lift @@ -1057,7 +1062,7 @@ export reverse_cols export reverse_cols! export reverse_rows export reverse_rows! -export right_kernel +#export right_kernel export rising_factorial export rising_factorial2 export rowlength @@ -1088,11 +1093,11 @@ export snf_kb_with_transform export snf_kb! export snf_with_transform export solve -export solve_ff -export solve_left -export solve_rational -export solve_triu -export solve_with_det +#export _solve_ff +#export _solve_left +#export _solve_rational +#export _solve_triu +#export _solve_with_det export sort_terms! export SparsePolynomialRing export Strassen diff --git a/src/Generic.jl b/src/Generic.jl index 6184e543e6..b6c3eedad1 100644 --- a/src/Generic.jl +++ b/src/Generic.jl @@ -116,7 +116,9 @@ import Base: // import Base: / import Base: != -# The type and helper function for the dictionaries for hashing +import ..AbstractAlgebra: _can_solve_with_solution_fflu +import ..AbstractAlgebra: _can_solve_with_solution_lu +import ..AbstractAlgebra: _left_kernel import ..AbstractAlgebra: @attributes import ..AbstractAlgebra: @enable_all_show_via_expressify import ..AbstractAlgebra: add! @@ -125,8 +127,6 @@ import ..AbstractAlgebra: addmul! import ..AbstractAlgebra: base_ring import ..AbstractAlgebra: base_ring_type import ..AbstractAlgebra: CacheDictType -import ..AbstractAlgebra: can_solve_with_solution_fflu -import ..AbstractAlgebra: can_solve_with_solution_lu import ..AbstractAlgebra: canonical_unit import ..AbstractAlgebra: change_base_ring import ..AbstractAlgebra: characteristic diff --git a/src/Matrix-Strassen.jl b/src/Matrix-Strassen.jl index bcf62330ab..f679c1109c 100644 --- a/src/Matrix-Strassen.jl +++ b/src/Matrix-Strassen.jl @@ -1,9 +1,9 @@ """ Provides generic asymptotically fast matrix methods: - mul and mul! using the Strassen scheme - - solve_tril! + - _solve_tril! - lu! - - solve_triu + - _solve_triu Just prefix the function by "Strassen." all 4 functions support a keyword argument "cutoff" to indicate when the base case should be used. @@ -171,9 +171,9 @@ end # B C Y V # => X = solve(A, U) # Y = solve(C, V - B*X) -function solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0; cutoff::Int = 2) where T +function _solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0; cutoff::Int = 2) where T if nrows(A) < cutoff || ncols(A) < cutoff - return AbstractAlgebra.solve_tril!(A, B, C, f) + return AbstractAlgebra._solve_tril!(A, B, C, f) end n = nrows(B) n2 = div(n, 2) @@ -184,10 +184,10 @@ function solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0; cu X2 = view(A, n2+1:n, 1:ncols(A)) C1 = view(C, 1:n2, 1:ncols(A)) C2 = view(C, n2+1:n, 1:ncols(A)) - solve_tril!(X1, B11, C1, f; cutoff) + _solve_tril!(X1, B11, C1, f; cutoff) x = B21 * X1 # strassen... sub!(X2, C2, x) - solve_tril!(X2, B22, X2, f; cutoff) + _solve_tril!(X2, B22, X2, f; cutoff) end function apply!(A::MatElem, P::Perm{Int}; offset::Int = 0) @@ -249,8 +249,8 @@ function lu!(P::Perm{Int}, A; cutoff::Int = 300) if r1 > 0 #Note: A00 is a view of A0 thus a view of A # A0 is lu!, thus implicitly two triangular matrices giving the - # lu decomosition. solve_tril! looks ONLY at the lower part of A00 - solve_tril!(A01, A00, A01, 1) + # lu decomosition. _solve_tril! looks ONLY at the lower part of A00 + _solve_tril!(A01, A00, A01, 1) X = A10 * A01 sub!(A11, A11, X) end @@ -271,11 +271,11 @@ function lu!(P::Perm{Int}, A; cutoff::Int = 300) return r1 + r2 end -function solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff) +function _solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff) #b*inv(T), thus solves Tx = b for T upper triangular n = ncols(T) if n <= cutoff - R = AbstractAlgebra.solve_triu(T, b) + R = AbstractAlgebra._solve_triu(T, b) return R end @@ -292,16 +292,16 @@ function solve_triu(T::MatElem, b::MatElem; cutoff::Int = cutoff) B = view(T, 1:n2, 1+n2:n) C = view(T, 1+n2:n, 1+n2:n) - S = solve_triu(A, U; cutoff) - R = solve_triu(A, X; cutoff) + S = _solve_triu(A, U; cutoff) + R = _solve_triu(A, X; cutoff) SS = mul(S, B; cutoff) sub!(SS, V, SS) - SS = solve_triu(C, SS; cutoff) + SS = _solve_triu(C, SS; cutoff) RR = mul(R, B; cutoff) sub!(RR, Y, RR) - RR = solve_triu(C, RR; cutoff) + RR = _solve_triu(C, RR; cutoff) return [S SS; R RR] end diff --git a/src/Matrix.jl b/src/Matrix.jl index fca9e99a18..7d0a135920 100644 --- a/src/Matrix.jl +++ b/src/Matrix.jl @@ -1948,7 +1948,7 @@ function rref!(A::MatrixElem{T}) where {T <: FieldElement} V[j, i] = A[j, pivots[np + i]] end end - V = solve_triu(U, V, false) + V = _solve_triu(U, V, false) for i = 1:rnk for j = 1:i A[j, pivots[i]] = i == j ? one(R) : R() @@ -2631,13 +2631,13 @@ end # permutations implicitly and add units along the diagonal of the lower # triangular matrix L instead of removing columns (and corresponding rows of # the upper triangular matrix U). We also set free variables to zero. -function can_solve_with_solution_fflu(A::MatElem{T}, b::MatElem{T}) where {T <: RingElement} - base_ring(A) != base_ring(b) && error("Base rings don't match in can_solve_with_solution_fflu") - nrows(A) != nrows(b) && error("Dimensions don't match in can_solve_with_solution_fflu") +function _can_solve_with_solution_fflu(A::MatElem{T}, b::MatElem{T}) where {T <: RingElement} + base_ring(A) != base_ring(b) && error("Base rings don't match in _can_solve_with_solution_fflu") + nrows(A) != nrows(b) && error("Dimensions don't match in _can_solve_with_solution_fflu") FFLU = deepcopy(A) p = one(SymmetricGroup(nrows(A))) rank, d = fflu!(p, FFLU) - flag, y = solve_fflu_precomp(p, FFLU, b) + flag, y = _solve_fflu_precomp(p, FFLU, b) n = nrows(A) if flag && rank < n b2 = p*b @@ -2654,7 +2654,7 @@ end # exists. If not, `flag` may be set to `false`, however it is not required to # be. If `r` is the rank of `A` then the first `r` rows of `p(A)y = p(b)d` will # hold iff `flag` is `true`. The remaining rows must be checked by the caller. -function solve_fflu_precomp(p::Perm, FFLU::MatElem{T}, b::MatElem{T}) where {T <: RingElement} +function _solve_fflu_precomp(p::Perm, FFLU::MatElem{T}, b::MatElem{T}) where {T <: RingElement} x = p * b n = nrows(x) m = ncols(x) @@ -2748,7 +2748,7 @@ end # permutations implicitly and add units along the diagonal of the lower # triangular matrix L instead of removing columns (and corresponding rows of # the upper triangular matrix U). We also set free variables to zero. -function can_solve_with_solution_lu(A::MatElem{T}, b::MatElem{T}) where {T <: FieldElement} +function _can_solve_with_solution_lu(A::MatElem{T}, b::MatElem{T}) where {T <: FieldElement} base_ring(A) != base_ring(b) && error("Base rings don't match in can_solve_with_solution_lu") nrows(A) != nrows(b) && error("Dimensions don't match in can_solve_with_solution_lu") @@ -2764,7 +2764,7 @@ function can_solve_with_solution_lu(A::MatElem{T}, b::MatElem{T}) where {T <: Fi p = one(SymmetricGroup(nrows(A))) rank = lu!(p, LU) - y = solve_lu_precomp(p, LU, b) + y = _solve_lu_precomp(p, LU, b) n = nrows(A) flag = true @@ -2782,7 +2782,7 @@ end # this function will return `y` such that the first `r` rows of `p(A)y = p(b)` # hold, where `r` is the rank of `A`. The remaining rows must be checked by # the caller. -function solve_lu_precomp(p::Perm, LU::MatElem{T}, b::MatElem{T}) where {T <: FieldElement} +function _solve_lu_precomp(p::Perm, LU::MatElem{T}, b::MatElem{T}) where {T <: FieldElement} x = p * b n = nrows(x) m = ncols(x) @@ -2855,12 +2855,12 @@ function solve_lu_precomp(p::Perm, LU::MatElem{T}, b::MatElem{T}) where {T <: Fi return y end -function solve_ff(M::MatrixElem{T}, b::MatrixElem{T}) where {T <: FieldElement} +function _solve_ff(M::MatrixElem{T}, b::MatrixElem{T}) where {T <: FieldElement} base_ring(M) != base_ring(b) && error("Base rings don't match in solve") nrows(M) != nrows(b) && error("Dimensions don't match in solve") m = nrows(M) - flag, x, d = can_solve_with_solution_fflu(M, b) - !flag && error("System not solvable in solve_ff") + flag, x, d = _can_solve_with_solution_fflu(M, b) + !flag && error("System not solvable in _solve_ff") for i in 1:nrows(x) for j in 1:ncols(x) x[i, j] = divexact(x[i, j], d) @@ -2869,8 +2869,8 @@ function solve_ff(M::MatrixElem{T}, b::MatrixElem{T}) where {T <: FieldElement} return x end -function can_solve_with_solution_with_det(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} - # We cannot use solve_fflu directly, since it forgot about the (parity of +function _can_solve_with_solution_with_det(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} + # We cannot use _solve_fflu directly, since it forgot about the (parity of # the) permutation. R = base_ring(M) FFLU = deepcopy(M) @@ -2888,7 +2888,7 @@ function can_solve_with_solution_with_det(M::MatElem{T}, b::MatElem{T}) where {T c += 1 end end - flag, x = solve_fflu_precomp(p, FFLU, b) + flag, x = _solve_fflu_precomp(p, FFLU, b) n = nrows(M) if flag && rank < n b2 = p*b @@ -2913,31 +2913,31 @@ function can_solve_with_solution_with_det(M::MatElem{T}, b::MatElem{T}) where {T return true, rank, p, pivots, x, d end -function can_solve_with_solution_with_det(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} - flag, r, p, pivots, x, d = can_solve_with_solution_interpolation_inner(M, b) +function _can_solve_with_solution_with_det(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} + flag, r, p, pivots, x, d = _can_solve_with_solution_interpolation_inner(M, b) return flag, r, p, pivots, x, d end -# This can be removed once Nemo implements can_solve_with_solution_with_det +# This can be removed once Nemo implements _can_solve_with_solution_with_det # It's here now only because Nemo overloads it -function solve_with_det(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} - flag, r, p, piv, x, d = can_solve_with_solution_with_det(M, b) - !flag && error("System not solvable in solve_with_det") +function _solve_with_det(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} + flag, r, p, piv, x, d = _can_solve_with_solution_with_det(M, b) + !flag && error("System not solvable in _solve_with_det") return x, d end -function solve_ff(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} +function _solve_ff(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} m = nrows(M) n = ncols(M) if m == 0 return zero_matrix(base_ring(M), ncols(M), ncols(b)), base_ring(M)() end if n == 0 - b != 0 && error("System not soluble in solve_ff") + b != 0 && error("System not soluble in _solve_ff") return zero_matrix(base_ring(M), ncols(M), ncols(b)), base_ring(M)() end - flag, S, d = can_solve_with_solution_fflu(M, b) - !flag && error("System not soluble in solve_ff") + flag, S, d = _can_solve_with_solution_fflu(M, b) + !flag && error("System not soluble in _solve_ff") return S, d end @@ -3015,7 +3015,7 @@ function can_solve_with_solution_interpolation_inner(M::MatElem{T}, b::MatElem{T if bad_evaluation error("Bad evaluation point") end - flag, r, p, pv, Vl, dl = can_solve_with_solution_with_det(X, Y) + flag, r, p, pv, Vl, dl = _can_solve_with_solution_with_det(X, Y) if !flag return flag, r, p, pv, zero(x), zero(R) end @@ -3093,13 +3093,13 @@ function can_solve_with_solution_interpolation_inner(M::MatElem{T}, b::MatElem{T return true, rnk, prm, pivots, x, interpolate(R, y, d) end -function can_solve_with_solution_interpolation(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} +function _can_solve_with_solution_interpolation(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} flag, r, p, pv, x, d = can_solve_with_solution_interpolation_inner(M, b) return flag, x, d end @doc raw""" - solve_rational(M::MatElem{T}, b::MatElem{T}) where T <: RingElement + _solve_rational(M::MatElem{T}, b::MatElem{T}) where T <: RingElement Given a non-singular $n\times n$ matrix over a ring and an $n\times m$ matrix over the same ring, return a tuple $x, d$ consisting of an @@ -3107,57 +3107,57 @@ $n\times m$ matrix $x$ and a denominator $d$ such that $Ax = db$. The denominator will be the determinant of $A$ up to sign. If $A$ is singular an exception is raised. """ -function solve_rational(M::MatElem{T}, b::MatElem{T}) where T <: RingElement - return solve_ringelem(M, b) +function _solve_rational(M::MatElem{T}, b::MatElem{T}) where T <: RingElement + return _solve_ringelem(M, b) end -function solve_ringelem(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} +function _solve_ringelem(M::MatElem{T}, b::MatElem{T}) where {T <: RingElement} base_ring(M) != base_ring(b) && error("Base rings don't match in solve") nrows(M) != nrows(b) && error("Dimensions don't match in solve") - return solve_ff(M, b) + return _solve_ff(M, b) end -function solve_rational(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} +function _solve_rational(M::MatElem{T}, b::MatElem{T}) where {T <: PolyRingElem} base_ring(M) != base_ring(b) && error("Base rings don't match in solve") nrows(M) != nrows(b) && error("Dimensions don't match in solve") flag = true try - flag, x, d = can_solve_with_solution_interpolation(M, b) - !flag && error("No solution in solve_rational") + flag, x, d = _can_solve_with_solution_interpolation(M, b) + !flag && error("No solution in _solve_rational") return x, d catch e if !isa(e, ErrorException) rethrow(e) end - !flag && error("No solution in solve_rational") - return solve_ff(M, b) + !flag && error("No solution in _solve_rational") + return _solve_ff(M, b) end end -@doc raw""" - solve(a::MatElem{S}, b::MatElem{S}) where {S <: RingElement} - -Given an $m\times r$ matrix $a$ over a ring and an $m\times n$ matrix $b$ -over the same ring, return an $r\times n$ matrix $x$ such that $ax = b$. If -no such matrix exists, an exception is raised. -See also [`solve_left`](@ref). -""" -function solve(a::MatElem{S}, b::MatElem{S} - ) where S <: RingElement - can, X = can_solve_with_solution(a, b, side = :right) - can || throw(ArgumentError("Unable to solve linear system")) - return X -end +# @doc raw""" +# solve(a::MatElem{S}, b::MatElem{S}) where {S <: RingElement} +# +# Given an $m\times r$ matrix $a$ over a ring and an $m\times n$ matrix $b$ +# over the same ring, return an $r\times n$ matrix $x$ such that $ax = b$. If +# no such matrix exists, an exception is raised. +# See also [`_solve_left`](@ref). +# """ +# function solve(a::MatElem{S}, b::MatElem{S} +# ) where S <: RingElement +# can, X = can_solve_with_solution(a, b, side = :right) +# can || throw(ArgumentError("Unable to solve linear system")) +# return X +# end @doc raw""" - solve_left(a::MatElem{S}, b::MatElem{S}) where S <: RingElement + _solve_left(a::MatElem{S}, b::MatElem{S}) where S <: RingElement Given an $r\times n$ matrix $a$ over a ring and an $m\times n$ matrix $b$ over the same ring, return an $m\times r$ matrix $x$ such that $xa = b$. If no such matrix exists, an exception is raised. See also [`solve`](@ref). """ -function solve_left(a::MatElem{S}, b::MatElem{S} +function _solve_left(a::MatElem{S}, b::MatElem{S} ) where S <: RingElement (flag, x) = can_solve_with_solution(a, b; side = :left) flag || error("Unable to solve linear system") @@ -3190,13 +3190,13 @@ end # ############################################################################### -function can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: FieldElement +function _can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: FieldElement @assert base_ring(A) == base_ring(B) if side === :right @assert nrows(A) == nrows(B) - return _can_solve_with_kernel(A, B) + return __can_solve_with_kernel(A, B) elseif side === :left - b, C, K = _can_solve_with_kernel(transpose(A), transpose(B)) + b, C, K = __can_solve_with_kernel(transpose(A), transpose(B)) @assert ncols(A) == ncols(B) if b return b, transpose(C), transpose(K) @@ -3208,7 +3208,7 @@ function can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) wher end end -function _can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}) where T <: FieldElement +function __can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}) where T <: FieldElement R = base_ring(A) mu = [A B] rk, mu = rref(mu) @@ -3258,14 +3258,14 @@ and $K$ is the kernel. Otherwise returns `false, S, K` where $S$ and $K$ are undefined, though of the right type for type stability. Tries to solve $Ax = B$ for $x$ if `side = :right` or $xA = B$ if `side = :left`. """ -function can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: RingElement +function _can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) where T <: RingElement @assert base_ring(A) == base_ring(B) if side === :right @assert nrows(A) == nrows(B) - b, c, K =_can_solve_with_kernel(transpose(A), B) + b, c, K = __can_solve_with_kernel(transpose(A), B) return b, transpose(c), K elseif side === :left - b, C, K = _can_solve_with_kernel(A, transpose(B)) + b, C, K = __can_solve_with_kernel(A, transpose(B)) @assert ncols(A) == ncols(B) if b return b, C, transpose(K) @@ -3278,7 +3278,7 @@ function can_solve_with_kernel(A::MatElem{T}, B::MatElem{T}; side = :right) wher end # Note that _a_ must be supplied transposed and the solution is transposed -function _can_solve_with_kernel(a::MatElem{S}, b::MatElem{S}) where S <: RingElement +function __can_solve_with_kernel(a::MatElem{S}, b::MatElem{S}) where S <: RingElement H, T = hnf_with_transform(a) z = zero_matrix(base_ring(a), ncols(b), nrows(a)) l = min(nrows(a), ncols(a)) @@ -3368,7 +3368,7 @@ function is_upper_triangular(M::MatrixElem) end @doc raw""" - solve_triu(U::MatElem{T}, b::MatElem{T}, unit::Bool = false) where {T <: FieldElement} + _solve_triu(U::MatElem{T}, b::MatElem{T}, unit::Bool = false) where {T <: FieldElement} Given a non-singular $n\times n$ matrix $U$ over a field which is upper triangular, and an $n\times m$ matrix $b$ over the same field, return an @@ -3376,7 +3376,7 @@ $n\times m$ matrix $x$ such that $Ux = b$. If $U$ is singular an exception is raised. If unit is true then $U$ is assumed to have ones on its diagonal, and the diagonal will not be read. """ -function solve_triu(U::MatElem{T}, b::MatElem{T}, unit::Bool = false) where {T <: FieldElement} +function _solve_triu(U::MatElem{T}, b::MatElem{T}, unit::Bool = false) where {T <: FieldElement} n = nrows(U) m = ncols(b) R = base_ring(U) @@ -3413,16 +3413,16 @@ function solve_triu(U::MatElem{T}, b::MatElem{T}, unit::Bool = false) where {T < end @doc raw""" - solve_triu(U::MatElem{T}, b::MatElem{T}) where {T <: RingElement} + _solve_triu(U::MatElem{T}, b::MatElem{T}) where {T <: RingElement} Given a non-singular $n\times n$ matrix $U$ over a field which is upper triangular, and an $n\times m$ matrix $b$ over the same ring, return an $n\times m$ matrix $x$ such that $Ux = b$. If this is not possible, an error will be raised. -See also [`AbstractAlgebra.solve_triu_left`](@ref) +See also [`AbstractAlgebra.__solve_triu_left`](@ref) """ -function solve_triu(U::MatElem{T}, b::MatElem{T}) where {T <: RingElement} +function _solve_triu(U::MatElem{T}, b::MatElem{T}) where {T <: RingElement} n = nrows(U) m = ncols(b) R = base_ring(U) @@ -3450,17 +3450,17 @@ function solve_triu(U::MatElem{T}, b::MatElem{T}) where {T <: RingElement} end @doc raw""" - solve_triu_left(b::MatElem{T}, U::MatElem{T}) where {T <: RingElement} + __solve_triu_left(b::MatElem{T}, U::MatElem{T}) where {T <: RingElement} Given a non-singular $n\times n$ matrix $U$ over a field which is upper triangular, and an $m\times n$ matrix $b$ over the same ring, return an $m\times n$ matrix $x$ such that $xU = b$. If this is not possible, an error will be raised. -See also [`solve_triu`](@ref) or [`can_solve_left_reduced_triu`](@ref) when +See also [`_solve_triu`](@ref) or [`can__solve_left_reduced_triu`](@ref) when $U$ is not square or not of full rank. """ -function solve_triu_left(b::MatElem{T}, U::MatElem{T}) where {T <: RingElement} +function __solve_triu_left(b::MatElem{T}, U::MatElem{T}) where {T <: RingElement} n = ncols(U) m = nrows(b) R = base_ring(U) @@ -3489,8 +3489,8 @@ end #solves A x = B for A intended to be lower triangular #only the lower part is used. if f is true, then the diagonal is assumed to be 1 #used to use lu! -#can be combined with Strassen.solve_tril! -function solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0) where T +#can be combined with Strassen._solve_tril! +function _solve_tril!(A::MatElem{T}, B::MatElem{T}, C::MatElem{T}, f::Int = 0) where T # a x u ax = u # b c * y = v bx + cy = v @@ -3603,7 +3603,7 @@ normal form or reduced row echelon form and $r$ and $x$ are row vectors with $m$ columns (i.e. $1 \times m$ matrices). If there is no solution, flag is set to `false` and $x$ is set to zero. """ -function can_solve_left_reduced_triu(r::MatElem{T}, +function _can_solve_left_reduced_triu(r::MatElem{T}, M::MatElem{T}) where T <: RingElement ncols(r) != ncols(M) && error("Incompatible matrices") r = deepcopy(r) # do not destroy input @@ -3646,164 +3646,164 @@ function can_solve_left_reduced_triu(r::MatElem{T}, return true, x end -function can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: FracElem{T} where T <: PolyRingElem - if side == :left - (f, x) = can_solve_with_solution(transpose(a), transpose(b); side=:right) - return (f, transpose(x)) - elseif side == :right - d = numerator(one(base_ring(a))) - for i = 1:nrows(a) - for j = 1:ncols(a) - d = lcm(d, denominator(a[i, j])) - end - end - for i = 1:nrows(b) - for j = 1:ncols(b) - d = lcm(d, denominator(b[i, j])) - end - end - A = matrix(parent(d), nrows(a), ncols(a), [numerator(a[i, j]*d) for i in 1:nrows(a) for j in 1:ncols(a)]) - B = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) - flag = false - x = similar(A, ncols(A), nrows(B)) - den = one(base_ring(a)) - try - flag, x, den = can_solve_with_solution_interpolation(A, B) - catch - flag, x, den = can_solve_with_solution_fflu(A, B) - end - X = change_base_ring(base_ring(a), x) - X = divexact(X, base_ring(a)(den)) - return flag, X - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -# The fflu approach is the fastest over a fraction field (see benchmarks on PR 661) -function can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: Union{FracElem, Rational{BigInt}} - if side == :left - (f, x) = can_solve_with_solution(transpose(a), transpose(b); side=:right) - return (f, transpose(x)) - elseif side == :right - d = numerator(one(base_ring(a))) - for i = 1:nrows(a) - for j = 1:ncols(a) - d = lcm(d, denominator(a[i, j])) - end - end - for i = 1:nrows(b) - for j = 1:ncols(b) - d = lcm(d, denominator(b[i, j])) - end - end - A = matrix(parent(d), nrows(a), ncols(a), [numerator(a[i, j]*d) for i in 1:nrows(a) for j in 1:ncols(a)]) - B = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) - flag, x, den = can_solve_with_solution_fflu(A, B) - X = change_base_ring(base_ring(a), x) - X = divexact(X, base_ring(a)(den)) - return flag, X - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -@doc raw""" - can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - -Given two matrices $a$ and $b$ over the same ring, try to solve $ax = b$ -if `side` is `:right` or $xa = b$ if `side` is `:left`. In either case, -return a tuple `(flag, x)`. If a solution exists, `flag` is set to true and -`x` is a solution. If no solution exists, `flag` is set to false and `x` -is arbitrary. If the dimensions of $a$ and $b$ are incompatible, an exception -is raised. -""" -function can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - if side == :right - (f, x) = can_solve_with_solution(transpose(a), transpose(b); side=:left) - return (f, transpose(x)) - elseif side == :left - @assert ncols(a) == ncols(b) - H, T = hnf_with_transform(a) - b = deepcopy(b) - z = zero(a, nrows(b), nrows(a)) - l = min(ncols(a), nrows(a)) - t = base_ring(a)() - for i = 1:nrows(b) - for j = 1:l - k = 1 - while k <= ncols(H) && is_zero_entry(H, j, k) - k += 1 - end - if k > ncols(H) - continue - end - q, r = divrem(b[i, k], H[j, k]) - if r != 0 - return (false, zero(a, 0, 0)) - end - z[i, j] = q - q = -q - for h = k:ncols(H) - t = mul!(t, q, H[j, h]) - b[i, h] = addeq!(b[i, h], t) - end - end - end - if b != 0 - return (false, zero(a, 0, 0)) - end - return (true, z*T) - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -function can_solve_with_solution(A::MatElem{T}, B::MatElem{T}; side::Symbol = :right) where T <: FieldElement - if side == :right - (f, x) = can_solve_with_solution(transpose(A), transpose(B), side = :left) - return (f, transpose(x)) - elseif side == :left - R = base_ring(A) - ncols(A) != ncols(B) && error("Incompatible matrices") - mu = zero_matrix(R, ncols(A), nrows(A) + nrows(B)) - for i = 1:ncols(A) - for j = 1:nrows(A) - mu[i, j] = A[j, i] - end - for j = 1:nrows(B) - mu[i, nrows(A) + j] = B[j, i] - end - end - rk, mu = rref(mu) - p = find_pivot(mu) - if any(i -> i > nrows(A), p) - return (false, zero(A, 0, 0)) - end - sol = zero_matrix(R, nrows(B), nrows(A)) - for i = 1:length(p) - for j = 1:nrows(B) - sol[j, p[i]] = mu[i, nrows(A) + j] - end - end - return (true, sol) - else - error("Unsupported argument :$side for side: Must be :left or :right.") - end -end - -@doc raw""" - can_solve(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - -Given two matrices $a$ and $b$ over the same ring, check the solubility -of $ax = b$ if `side` is `:right` or $xa = b$ if `side` is `:left`. -Return true if a solution exists, false otherwise. If the dimensions -of $a$ and $b$ are incompatible, an exception is raised. If a solution -should be computed as well, use `can_solve_with_solution` instead. -""" -function can_solve(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement - return can_solve_with_solution(a, b; side=side)[1] -end + function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: FracElem{T} where T <: PolyRingElem + if side == :left + (f, x) = _can_solve_with_solution(transpose(a), transpose(b); side=:right) + return (f, transpose(x)) + elseif side == :right + d = numerator(one(base_ring(a))) + for i = 1:nrows(a) + for j = 1:ncols(a) + d = lcm(d, denominator(a[i, j])) + end + end + for i = 1:nrows(b) + for j = 1:ncols(b) + d = lcm(d, denominator(b[i, j])) + end + end + A = matrix(parent(d), nrows(a), ncols(a), [numerator(a[i, j]*d) for i in 1:nrows(a) for j in 1:ncols(a)]) + B = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) + flag = false + x = similar(A, ncols(A), nrows(B)) + den = one(base_ring(a)) + try + flag, x, den = _can_solve_with_solution_interpolation(A, B) + catch + flag, x, den = _can_solve_with_solution_fflu(A, B) + end + X = change_base_ring(base_ring(a), x) + X = divexact(X, base_ring(a)(den)) + return flag, X + else + error("Unsupported argument :$side for side: Must be :left or :right.") + end + end + + # The fflu approach is the fastest over a fraction field (see benchmarks on PR 661) +function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: Union{FracElem, Rational{BigInt}} + if side == :left + (f, x) = _can_solve_with_solution(transpose(a), transpose(b); side=:right) + return (f, transpose(x)) + elseif side == :right + d = numerator(one(base_ring(a))) + for i = 1:nrows(a) + for j = 1:ncols(a) + d = lcm(d, denominator(a[i, j])) + end + end + for i = 1:nrows(b) + for j = 1:ncols(b) + d = lcm(d, denominator(b[i, j])) + end + end + A = matrix(parent(d), nrows(a), ncols(a), [numerator(a[i, j]*d) for i in 1:nrows(a) for j in 1:ncols(a)]) + B = matrix(parent(d), nrows(b), ncols(b), [numerator(b[i, j]*d) for i in 1:nrows(b) for j in 1:ncols(b)]) + flag, x, den = _can_solve_with_solution_fflu(A, B) + X = change_base_ring(base_ring(a), x) + X = divexact(X, base_ring(a)(den)) + return flag, X + else + error("Unsupported argument :$side for side: Must be :left or :right.") + end + end + + @doc raw""" + _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement + + Given two matrices $a$ and $b$ over the same ring, try to solve $ax = b$ + if `side` is `:right` or $xa = b$ if `side` is `:left`. In either case, + return a tuple `(flag, x)`. If a solution exists, `flag` is set to true and + `x` is a solution. If no solution exists, `flag` is set to false and `x` + is arbitrary. If the dimensions of $a$ and $b$ are incompatible, an exception + is raised. + """ +function _can_solve_with_solution(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement + if side == :right + (f, x) = _can_solve_with_solution(transpose(a), transpose(b); side=:left) + return (f, transpose(x)) + elseif side == :left + @assert ncols(a) == ncols(b) + H, T = hnf_with_transform(a) + b = deepcopy(b) + z = zero(a, nrows(b), nrows(a)) + l = min(ncols(a), nrows(a)) + t = base_ring(a)() + for i = 1:nrows(b) + for j = 1:l + k = 1 + while k <= ncols(H) && is_zero_entry(H, j, k) + k += 1 + end + if k > ncols(H) + continue + end + q, r = divrem(b[i, k], H[j, k]) + if r != 0 + return (false, zero(a, 0, 0)) + end + z[i, j] = q + q = -q + for h = k:ncols(H) + t = mul!(t, q, H[j, h]) + b[i, h] = addeq!(b[i, h], t) + end + end + end + if b != 0 + return (false, zero(a, 0, 0)) + end + return (true, z*T) + else + error("Unsupported argument :$side for side: Must be :left or :right.") + end + end + +function _can_solve_with_solution(A::MatElem{T}, B::MatElem{T}; side::Symbol = :right) where T <: FieldElement + if side == :right + (f, x) = _can_solve_with_solution(transpose(A), transpose(B), side = :left) + return (f, transpose(x)) + elseif side == :left + R = base_ring(A) + ncols(A) != ncols(B) && error("Incompatible matrices") + mu = zero_matrix(R, ncols(A), nrows(A) + nrows(B)) + for i = 1:ncols(A) + for j = 1:nrows(A) + mu[i, j] = A[j, i] + end + for j = 1:nrows(B) + mu[i, nrows(A) + j] = B[j, i] + end + end + rk, mu = rref(mu) + p = find_pivot(mu) + if any(i -> i > nrows(A), p) + return (false, zero(A, 0, 0)) + end + sol = zero_matrix(R, nrows(B), nrows(A)) + for i = 1:length(p) + for j = 1:nrows(B) + sol[j, p[i]] = mu[i, nrows(A) + j] + end + end + return (true, sol) + else + error("Unsupported argument :$side for side: Must be :left or :right.") + end + end + + @doc raw""" + _can_solve(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement + + Given two matrices $a$ and $b$ over the same ring, check the solubility + of $ax = b$ if `side` is `:right` or $xa = b$ if `side` is `:left`. + Return true if a solution exists, false otherwise. If the dimensions + of $a$ and $b$ are incompatible, an exception is raised. If a solution + should be computed as well, use `can_solve_with_solution` instead. + """ + function _can_solve(a::MatElem{S}, b::MatElem{S}; side::Symbol = :right) where S <: RingElement + return can_solve_with_solution(a, b; side=side)[1] + end ############################################################################### # @@ -3822,14 +3822,14 @@ is raised. """ function pseudo_inv(M::MatrixElem{T}) where {T <: RingElement} is_square(M) || throw(DomainError(M, "Can not invert non-square Matrix")) - flag, X, d = can_solve_with_solution_fflu(M, identity_matrix(M)) + flag, X, d = _can_solve_with_solution_fflu(M, identity_matrix(M)) !flag && error("Singular matrix in pseudo_inv") return X, d end function Base.inv(M::MatrixElem{T}) where {T <: FieldElement} is_square(M) || throw(DomainError(M, "Can not invert non-square Matrix")) - flag, A = can_solve_with_solution_lu(M, identity_matrix(M)) + flag, A = _can_solve_with_solution_lu(M, identity_matrix(M)) !flag && error("Singular matrix in inv") return A end @@ -4018,14 +4018,14 @@ end ############################################################################### @doc raw""" - left_kernel(A::MatElem{T}) where T <: RingElement + _left_kernel(A::MatElem{T}) where T <: RingElement Return a tuple $(n, M)$ where $M$ is a matrix whose rows generate the kernel of $A$ and $n$ is the rank of the kernel. The transpose of the output of this function is guaranteed to be in flipped upper triangular format (i.e. upper triangular format if columns and rows are reversed). """ -function left_kernel(A::MatElem{T}) where T <: RingElement +function _left_kernel(A::MatElem{T}) where T <: RingElement !is_domain_type(elem_type(base_ring(A))) && error("Not implemented") R = base_ring(A) H, U = hnf_with_transform(A) @@ -4042,45 +4042,45 @@ function left_kernel(A::MatElem{T}) where T <: RingElement end end -function left_kernel(A::MatElem{T}) where T <: FieldElement +function _left_kernel(A::MatElem{T}) where T <: FieldElement n, M = nullspace(transpose(A)) return n, transpose(M) end @doc raw""" - right_kernel(A::MatElem{T}) where T <: RingElement + _right_kernel(A::MatElem{T}) where T <: RingElement Return a tuple $(n, M)$ where $M$ is a matrix whose columns generate the kernel of $A$ and $n$ is the rank of the kernel. """ -function right_kernel(A::MatElem{T}) where T <: RingElement - n, M = left_kernel(transpose(A)) +function _right_kernel(A::MatElem{T}) where T <: RingElement + n, M = _left_kernel(transpose(A)) return n, transpose(M) end -function right_kernel(A::MatElem{T}) where T <: FieldElement +function _right_kernel(A::MatElem{T}) where T <: FieldElement return nullspace(A) end -@doc raw""" - kernel(A::MatElem{T}; side::Symbol = :right) where T <: RingElement - -Return a tuple $(n, M)$, where $n$ is the rank of the kernel of $A$ and $M$ is a -basis for it. If side is `:right` or not specified, the right kernel is -computed, i.e. the matrix of columns whose span gives the right kernel -space. If side is `:left`, the left kernel is computed, i.e. the matrix -of rows whose span is the left kernel space. -""" -function kernel(A::MatElem{T}; side::Symbol = :right) where T <: RingElement - if side == :right - return right_kernel(A) - elseif side == :left - return left_kernel(A) - else - error("Unsupported argument: :$side for side: must be :left or :right") - end -end - +# @doc raw""" +# _kernel(A::MatElem{T}; side::Symbol = :right) where T <: RingElement +# +# Return a tuple $(n, M)$, where $n$ is the rank of the kernel of $A$ and $M$ is a +# basis for it. If side is `:right` or not specified, the right kernel is +# computed, i.e. the matrix of columns whose span gives the right kernel +# space. If side is `:left`, the left kernel is computed, i.e. the matrix +# of rows whose span is the left kernel space. +# """ +# function _kernel(A::MatElem{T}; side::Symbol = :right) where T <: RingElement +# if side == :right +# return _right_kernel(A) +# elseif side == :left +# return _left_kernel(A) +# else +# error("Unsupported argument: :$side for side: must be :left or :right") +# end +# end + ################################################################################ # # Kernel over different rings diff --git a/src/Module.jl b/src/Module.jl index 4f530c49cc..a1e0dee1c1 100644 --- a/src/Module.jl +++ b/src/Module.jl @@ -228,14 +228,14 @@ function ==(M::FPModule{T}, N::FPModule{T}) where T <: RingElement mat2 = reduced_form(mat2) # Check containment of rewritten gens of M in row space of mat2 for v in G1 - flag, r = can_solve_left_reduced_triu(Generic._matrix(v), mat2) + flag, r = _can_solve_left_reduced_triu(Generic._matrix(v), mat2) if !flag return false end end # Check containment of rewritten gens of N in row space of mat1 for v in G2 - flag, r = can_solve_left_reduced_triu(Generic._matrix(v), mat1) + flag, r = _can_solve_left_reduced_triu(Generic._matrix(v), mat1) if !flag return false end diff --git a/src/ModuleHomomorphism.jl b/src/ModuleHomomorphism.jl index 26c56b209c..eb2c5abc78 100644 --- a/src/ModuleHomomorphism.jl +++ b/src/ModuleHomomorphism.jl @@ -73,7 +73,7 @@ function kernel(f::Map(FPModuleHomomorphism)) end end # compute the kernel - num_gens, K = left_kernel(N) + num_gens, K = _left_kernel(N) # Construct generators of kernel submodule, reversing rows # and columns so they're correct wrt to original data and # in upper triangular form @@ -151,7 +151,7 @@ function preimage(f::Map(FPModuleHomomorphism), v::FPModuleElem{T}) where end end # Find left inverse of mat - x = solve_left(matr, Generic._matrix(v)) + x = _solve_left(matr, Generic._matrix(v)) if q != 0 x = matrix(R, 1, m, T[x[1, i] for i in 1:m]) end diff --git a/src/Solve.jl b/src/Solve.jl index 73fb8a2fb5..887cd04dd1 100644 --- a/src/Solve.jl +++ b/src/Solve.jl @@ -2,7 +2,7 @@ module Solve using AbstractAlgebra -import AbstractAlgebra: base_ring, nrows, ncols, matrix, rank, Generic +import AbstractAlgebra: base_ring, nrows, ncols, matrix, rank, Generic, kernel ################################################################################ # diff --git a/src/generic/MatRing.jl b/src/generic/MatRing.jl index 61f4bdf866..83712b2067 100644 --- a/src/generic/MatRing.jl +++ b/src/generic/MatRing.jl @@ -43,22 +43,22 @@ end # ############################################################################### -function can_solve_with_solution_lu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} +function _can_solve_with_solution_lu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) MS = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix BS = MatSpaceElem{T}(R, B.entries) - flag, S = can_solve_with_solution_lu(MS, BS) + flag, S = _can_solve_with_solution_lu(MS, BS) SA = MatRingElem{T}(R, S.entries) return flag, SA end -function can_solve_with_solution_fflu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} +function _can_solve_with_solution_fflu(M::MatRingElem{T}, B::MatRingElem{T}) where {T <: RingElement} check_parent(M, B) R = base_ring(M) MS = MatSpaceElem{T}(R, M.entries) # convert to ordinary matrix BS = MatSpaceElem{T}(R, B.entries) - flag, S, d = can_solve_with_solution_fflu(MS, BS) + flag, S, d = _can_solve_with_solution_fflu(MS, BS) SA = MatRingElem{T}(R, S.entries) return flag, SA, d end diff --git a/src/generic/ModuleHomomorphism.jl b/src/generic/ModuleHomomorphism.jl index a226598def..3eea803435 100644 --- a/src/generic/ModuleHomomorphism.jl +++ b/src/generic/ModuleHomomorphism.jl @@ -133,7 +133,7 @@ function ModuleIsomorphism(M1::AbstractAlgebra.FPModule{T}, end # Find left inverse of mat N = identity_matrix(R, n) - X = solve_left(mat, N) + X = _solve_left(mat, N) # Construct matrix of inverse homomorphism from first m columns of X M_inv = X[:, 1:m] end diff --git a/src/generic/Submodule.jl b/src/generic/Submodule.jl index 36b0b7a7ff..911b56579a 100644 --- a/src/generic/Submodule.jl +++ b/src/generic/Submodule.jl @@ -188,7 +188,7 @@ function sub(m::AbstractAlgebra.FPModule{T}, gens::Vector{S}) where {T <: RingEl end end # Rewrite old relations in terms of generators of new submodule - num_rels, K = left_kernel(new_mat) + num_rels, K = _left_kernel(new_mat) new_rels = zero_matrix(base_ring(m), num_rels, num) # we flip rows and columns so that input is in terms of original data and # in upper triangular form, to save time in reduced_form below diff --git a/test/generic/MatRing-test.jl b/test/generic/MatRing-test.jl index 8cde9804e7..05dadcaa38 100644 --- a/test/generic/MatRing-test.jl +++ b/test/generic/MatRing-test.jl @@ -898,7 +898,7 @@ end M = randmat_with_rank(R, dim, -100:100) b = rand(U, -100:100) - flag, x = Generic.can_solve_with_solution_lu(M, b) + flag, x = Generic._can_solve_with_solution_lu(M, b) @test flag && M*x == b end @@ -917,7 +917,7 @@ end MK = T(elem_type(K)[ K(M[i, j]) for i in 1:nrows(M), j in 1:ncols(M) ]) bK = T(elem_type(K)[ K(b[i, j]) for i in 1:nrows(b), j in 1:ncols(b) ]) - flag, x = Generic.can_solve_with_solution_lu(MK, bK) + flag, x = Generic._can_solve_with_solution_lu(MK, bK) @test flag && MK*x == bK end diff --git a/test/generic/Matrix-test.jl b/test/generic/Matrix-test.jl index 65ea3f0570..5c58a8eb35 100644 --- a/test/generic/Matrix-test.jl +++ b/test/generic/Matrix-test.jl @@ -2027,7 +2027,7 @@ end end end -@testset "Generic.Mat.can_solve_with_solution_fflu" begin +@testset "Generic.Mat._can_solve_with_solution_fflu" begin R = ZZ # Test random soluble systems @@ -2053,10 +2053,10 @@ end d2 = rand(R, -20:20) end A *= d2 - flag, X, d = Generic.can_solve_with_solution_fflu(A, B) + flag, X, d = Generic._can_solve_with_solution_fflu(A, B) @test flag && A*X == B*d B = rand(T, -10:10) - flag, X, d = Generic.can_solve_with_solution_fflu(A, B) + flag, X, d = Generic._can_solve_with_solution_fflu(A, B) @test (flag && A*X == B*d) || !flag end @@ -2071,10 +2071,10 @@ end U = matrix_space(R, n, k) A = randmat_with_rank(S, rank, -20:20) B = rand(T, -20:20) - flag1, X, d = Generic.can_solve_with_solution_fflu(A, B) + flag1, X, d = Generic._can_solve_with_solution_fflu(A, B) A2 = change_base_ring(QQ, A) B2 = change_base_ring(QQ, B) - flag2, X2 = can_solve_with_solution(A2, B2) + flag2, X2 = AbstractAlgebra._can_solve_with_solution(A2, B2) @test flag1 == flag2 end end @@ -2100,10 +2100,10 @@ end end X2 = rand(U) B = A*X2 - flag, X = Generic.can_solve_with_solution_lu(A, B) + flag, X = Generic._can_solve_with_solution_lu(A, B) @test flag && A*X == B B = rand(T) - flag, X = Generic.can_solve_with_solution_lu(A, B) + flag, X = Generic._can_solve_with_solution_lu(A, B) @test (flag && A*X == B) || !flag end @@ -2120,7 +2120,7 @@ end A = randmat_with_rank(S, rank, -20:20) X2 = rand(U, -20:20) B = A*X2 - flag, X = Generic.can_solve_with_solution_lu(A, B) + flag, X = Generic._can_solve_with_solution_lu(A, B) @test flag && A*X == B end @@ -2131,7 +2131,7 @@ end M = randmat_with_rank(S, dim, -100:100) b = rand(U, -100:100) - flag, x = Generic.can_solve_with_solution_lu(M, b) + flag, x = Generic._can_solve_with_solution_lu(M, b) @test flag && M*x == b end @@ -2149,13 +2149,13 @@ end MK = matrix(K, elem_type(K)[ K(M[i, j]) for i in 1:nrows(M), j in 1:ncols(M) ]) bK = matrix(K, elem_type(K)[ K(b[i, j]) for i in 1:nrows(b), j in 1:ncols(b) ]) - flag, x = Generic.can_solve_with_solution_lu(MK, bK) + flag, x = Generic._can_solve_with_solution_lu(MK, bK) @test flag && MK*x == bK end end -@testset "Generic.Mat.solve_ff" begin +@testset "Generic.Mat._solve_ff" begin # Exact field R = QQ @@ -2176,7 +2176,7 @@ end end B = divexact(B, d) @test A*divexact(X2, d) == B - X = Generic.solve_ff(A, B) + X = AbstractAlgebra._solve_ff(A, B) @test A*X == B end @@ -2195,12 +2195,12 @@ end A = randmat_with_rank(S, rank, -20:20) X2 = rand(U, -20:20) B = A*X2 - X, d = Generic.solve_ff(A, B) + X, d = AbstractAlgebra._solve_ff(A, B) @test A*X == B*d end end -@testset "Generic.Mat.solve_rational" begin +@testset "Generic.Mat._solve_rational" begin R = ZZ for i = 1:100 @@ -2214,7 +2214,7 @@ end A = randmat_with_rank(S, rank, -20:20) X2 = rand(U, -20:20) B = A*X2 - X, d = solve_rational(A, B) + X, d = AbstractAlgebra._solve_rational(A, B) @test A*X == B*d end @@ -2229,7 +2229,7 @@ end do_test = false try - x, d = solve_rational(M, b) + x, d = AbstractAlgebra._solve_rational(M, b) do_test = true catch e if !(e isa ErrorException) @@ -2259,7 +2259,7 @@ end end M *= d2 b = M*x2 - x, d = solve_rational(M, b) + x, d = AbstractAlgebra._solve_rational(M, b) @test M*x == d*b end @@ -2275,7 +2275,7 @@ end M = randmat_with_rank(S, dim, -100:100) b = rand(U, -100:100) - x = solve(M, b) + x = AbstractAlgebra._solve(M, b) @test M*x == b end @@ -2298,7 +2298,7 @@ end end M *= d2 b = M*x2 - x, d = solve_rational(M, b) + x, d = AbstractAlgebra._solve_rational(M, b) @test M*x == d*b end @@ -2314,7 +2314,7 @@ end M = T([3y*a^2 + (y + 1)*a + 2y (5y+1)*a^2 + 2a + y - 1 a^2 + (-a) + 2y; (y + 1)*a^2 + 2y - 4 3y*a^2 + (2y - 1)*a + y (4y - 1)*a^2 + (y - 1)*a + 5; 2a + y + 1 (2y + 2)*a^2 + 3y*a + 3y a^2 + (-y-1)*a + (-y - 3)]) b = U(permutedims([4y*a^2 + 4y*a + 2y + 1 5y*a^2 + (2y + 1)*a + 6y + 1 (y + 1)*a^2 + 3y*a + 2y + 4], [2, 1])) - x, d = solve_rational(M, b) + x, d = AbstractAlgebra._solve_rational(M, b) @test M*x == d*b end @@ -2334,7 +2334,7 @@ end M = rand(S, -20:20) B = M*X1 - X = solve(M, B) + X = AbstractAlgebra.solve(M, B) @test M*X == B end @@ -2355,14 +2355,14 @@ end M = rand(S, -1:2, -10:10) B = M*X1 - X = solve(M, B) + X = AbstractAlgebra.solve(M, B) @test M*X == B end end end -@testset "Generic.Mat.solve_left" begin +@testset "Generic.Mat._solve_left" begin for R in [ZZ, QQ] for iter = 1:40 for dim = 0:5 @@ -2377,7 +2377,7 @@ end M = rand(U, -20:20) B = X1*M - X = solve_left(M, X1*M) + X = AbstractAlgebra._solve_left(M, X1*M) @test X*M == B end @@ -2399,7 +2399,7 @@ end M = rand(U, -1:2, -10:10) B = X1*M - X = solve_left(M, X1*M) + X = AbstractAlgebra._solve_left(M, X1*M) @test X*M == B end @@ -2420,13 +2420,13 @@ end X2 = rand(U, -1:2, -10:10) b = M*X2 - flag, X = Generic.can_solve_with_solution(M, b) + flag, X = Generic._can_solve_with_solution(M, b) @test flag && M*X == b b = X2*M - flag, X = Generic.can_solve_with_solution(M, b; side=:left) + flag, X = Generic._can_solve_with_solution(M, b; side=:left) @test flag && X*M == b end @@ -2446,7 +2446,7 @@ end M = change_base_ring(S, M) b = change_base_ring(S, b) - flag, X = can_solve_with_solution(M, b) + flag, X = AbstractAlgebra._can_solve_with_solution(M, b) @test flag && M*X == b @@ -2457,7 +2457,7 @@ end M = change_base_ring(S, M) b = change_base_ring(S, b) - flag, X = can_solve_with_solution(M, b; side=:left) + flag, X = AbstractAlgebra._can_solve_with_solution(M, b; side=:left) @test flag && X*M == b end @@ -2477,9 +2477,9 @@ end M = rand(U, -20:20) B = M*X1 - (flag, X) = can_solve_with_solution(M, M*X1) + (flag, X) = AbstractAlgebra._can_solve_with_solution(M, M*X1) - @test can_solve(M, B) + @test AbstractAlgebra._can_solve(M, B) @test flag && M*X == B end @@ -2491,9 +2491,9 @@ end M = rand(U, -20:20) B = X1*M - (flag, X) = can_solve_with_solution(M, X1*M; side = :left) + (flag, X) = AbstractAlgebra._can_solve_with_solution(M, X1*M; side = :left) - @test can_solve(M, B; side = :left) + @test AbstractAlgebra._can_solve(M, B; side = :left) @test flag && X*M == B end end @@ -2516,14 +2516,14 @@ end M = rand(U, -1:2, -10:10) B = M*X1 - (flag, X) = can_solve_with_solution(M, M*X1) + (flag, X) = AbstractAlgebra._can_solve_with_solution(M, M*X1) - @test can_solve(M, B) + @test AbstractAlgebra._can_solve(M, B) @test flag && M*X == B - (flag, X) = can_solve_with_solution(M, M*X1; side = :right) + (flag, X) = AbstractAlgebra._can_solve_with_solution(M, M*X1; side = :right) - @test can_solve(M, B; side = :right) + @test AbstractAlgebra._can_solve(M, B; side = :right) @test flag && M*X == B end @@ -2535,9 +2535,9 @@ end M = rand(U, -1:2, -10:10) B = X1*M - (flag, X) = can_solve_with_solution(M, X1*M; side = :left) + (flag, X) = AbstractAlgebra._can_solve_with_solution(M, X1*M; side = :left) - @test can_solve(M, B; side = :left) + @test AbstractAlgebra._can_solve(M, B; side = :left) @test flag && X*M == B end end @@ -2547,12 +2547,12 @@ end M = matrix(R, 1, 1, [x]) X = matrix(R, 1, 1, [1]) - @assert !can_solve(M, X) - (flag, _) = can_solve_with_solution(M, X) + @assert !AbstractAlgebra._can_solve(M, X) + (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X) @assert !flag - @assert !can_solve(M, X; side = :left) - (flag, _) = can_solve_with_solution(M, X; side = :left) + @assert !AbstractAlgebra._can_solve(M, X; side = :left) + (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X; side = :left) @assert !flag end @@ -2560,12 +2560,12 @@ end M = matrix(ZZ, 2, 2, [1, 1, 1, 1]) X = matrix(ZZ, 2, 1, [1, 0]) - @assert !can_solve(M, X) - (flag, _) = can_solve_with_solution(M, X) + @assert !AbstractAlgebra._can_solve(M, X) + (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X) @assert !flag - @assert !can_solve(M, transpose(X); side = :left) - (flag, _) = can_solve_with_solution(M, transpose(X); side = :left) + @assert !AbstractAlgebra._can_solve(M, transpose(X); side = :left) + (flag, _) = AbstractAlgebra._can_solve_with_solution(M, transpose(X); side = :left) @assert !flag end @@ -2573,18 +2573,18 @@ end M = matrix(ZZ, 3, 3, [2, 0, 0, 0, 1, 0, 0, 0, 1]) X1 = matrix(ZZ, 3, 1, [1, 0, 0]) - @assert !can_solve(M, X1) - (flag, X) = can_solve_with_solution(M, X1) + @assert !AbstractAlgebra._can_solve(M, X1) + (flag, X) = AbstractAlgebra._can_solve_with_solution(M, X1) @assert !flag X2 = matrix(ZZ, 2, 3, [1, 0, 0, 0, 1, 0]) - @assert !can_solve(M, X2; side = :left) - (flag, _) = can_solve_with_solution(M, X2; side = :left) + @assert !AbstractAlgebra._can_solve(M, X2; side = :left) + (flag, _) = AbstractAlgebra._can_solve_with_solution(M, X2; side = :left) @assert !flag end - @test_throws Exception can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = :aaa) - @test_throws TypeError can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = "right") + @test_throws Exception AbstractAlgebra._can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = :aaa) + @test_throws TypeError AbstractAlgebra._can_solve_with_solution(matrix(ZZ, 2, 2, [1, 0, 0, 1]), matrix(ZZ, 2, 1, [2, 3]), side = "right") end @testset "Generic.Mat.can_solve_with_solution_interpolation" begin @@ -2615,13 +2615,13 @@ end end M = M*d2 - flag, X, d = Generic.can_solve_with_solution_interpolation(M, B) + flag, X, d = Generic._can_solve_with_solution_interpolation(M, B) @test flag && M*X == B*d end end -@testset "Generic.Mat.solve_triu" begin +@testset "Generic.Mat._solve_triu" begin R, x = polynomial_ring(QQ, "x") K, = residue_field(R, x^3 + 3x + 1) a = K(x) @@ -2633,13 +2633,13 @@ end M = randmat_triu(S, -100:100) b = rand(U, -100:100) - x = solve_triu(M, b, false) + x = _solve_triu(M, b, false) @test M*x == b end end -@testset "Generic.Mat.solve_left_reduced_triu" begin +@testset "Generic.Mat._solve_left_reduced_triu" begin for iter = 1:40 n = rand(1:6) m = rand(1:n) @@ -2651,7 +2651,7 @@ end M = hnf(M) - flag, x = can_solve_left_reduced_triu(r, M) + flag, x = can__solve_left_reduced_triu(r, M) @test flag == false || x*M == r end @@ -3002,7 +3002,7 @@ end B = M*X - flag, Sol, K = can_solve_with_kernel(M, B) + flag, Sol, K = AbstractAlgebra._can_solve_with_kernel(M, B) @test iszero(M*K) @test flag && M*Sol == B @@ -3010,8 +3010,8 @@ end # fully random B = rand(S, -20:20) - flag, Sol, K = can_solve_with_kernel(M, B) - flag2, Sol2 = can_solve_with_solution(M, B) + flag, Sol, K = AbstractAlgebra._can_solve_with_kernel(M, B) + flag2, Sol2 = AbstractAlgebra._can_solve_with_solution(M, B) @test flag == flag2 if flag @@ -3030,7 +3030,7 @@ end B = X*M - flag, Sol, K = can_solve_with_kernel(M, B; side=:left) + flag, Sol, K = AbstractAlgebra._can_solve_with_kernel(M, B; side=:left) @test iszero(K*M) @test flag && Sol*M == B diff --git a/test/runtests.jl b/test/runtests.jl index 54aeef1f25..f00e12c1ac 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,8 @@ using AbstractAlgebra using Test -if VERSION >= v"1.8.0" +# disable until we encounter GC problems again +if false VERSION >= v"1.8.0" GC.enable_logging(true) # print gc settings