diff --git a/src/fraction.jl b/src/fraction.jl index 2e7b87f..de54358 100644 --- a/src/fraction.jl +++ b/src/fraction.jl @@ -132,7 +132,7 @@ end function *(x::T, y::T) where T<:Frac a, b, c, d = x.num, x.den, y.num, y.den - g = pgcd(a, d) + g = pgcd(a, d) # TODO call pgcd for GCD domains which are not Eucidean domains a /= g d /= g g = pgcd(b, c) diff --git a/src/linearalgebra.jl b/src/linearalgebra.jl index a07639f..1e01d9c 100644 --- a/src/linearalgebra.jl +++ b/src/linearalgebra.jl @@ -402,16 +402,33 @@ end """ adjugate(A::AbstractMatrix{<:Ring}) -The adjugate of matrix `A`. Invariant is `det(A) * I == adjugate(A) * A == det(A) * I(n)`. +The adjugate of matrix `A`. Invariant is `adjugate(A) * A == det(A) * I`. """ -function adjugate(A::AbstractMatrix{P}) where P +function adjugate(A::AbstractMatrix{P}) where P<:Ring + da = det(A) + if isunit(da) + _adjugate(A, da) + else + _adjugate_fallback(A) + end +end +function _adjugate(A::AbstractMatrix{P}, da) where P<:Union{Polynomial,ZZ} + Q = Frac(P) + B = Q.(A) + C = inv(B) .* Q(da) + numerator.(C) +end +function _adjugate(A::AbstractMatrix{P}, da) where P + inv(A) * da +end +function _adjugate_fallback(A::AbstractMatrix{P}) where P PP = P[:λ] Q = Frac(PP) - B = Q(monom(PP)) + A + B = Q(monom(PP)) + Q.(A) d = det(B) C = inv(B) .* d D = numerator.(C) - evaluate.(D, 0) + convert.(P, evaluate.(D, 0)) end """ diff --git a/src/resultant.jl b/src/resultant.jl index 91a5fa4..598d7b5 100644 --- a/src/resultant.jl +++ b/src/resultant.jl @@ -6,11 +6,11 @@ Pseudo gcd of univariate polynomials `a` and `b`. Uses subresultant sequence to accomplish non-field coeffient types. """ function pgcd(a::T, b::T) where {S,T<:UnivariatePolynomial{S}} -iszero(b) && return a -iszero(a) && return b -a, cc = presultant_seq(a, b, Val(false)) -a = primpart(a) -a / lcunit(a) * cc + iszero(b) && return a + iszero(a) && return b + a, cc = presultant_seq(a, b, Val(false)) + a = primpart(a) + T(a / lcunit(a)) * cc end """ @@ -19,23 +19,23 @@ resultant(a, b) Calculate resultant of two univariate polynomials of general coeffient types. """ function resultant(a::T, b::T) where {S,T<:UnivariatePolynomial{S}} -_resultant(a, b, category_trait(S)) + _resultant(a, b, category_trait(S)) end function _resultant( -a::T, -b::T, -::Type{<:EuclidianDomainTrait}, + a::T, + b::T, + ::Type{<:EuclidianDomainTrait}, ) where {S,T<:UnivariatePolynomial{S}} -_, _, r = presultant_seq(a, b, Val(true)) -r(0) + _, _, r = presultant_seq(a, b, Val(true)) + r(0) end function _resultant( -a::T, -b::T, -::Type{<:CommutativeRingTrait}, + a::T, + b::T, + ::Type{<:CommutativeRingTrait}, ) where {S,T<:UnivariatePolynomial{S}} -resultant_naive(a, b) + resultant_naive(a, b) end resultant(a::T, b::T) where T = iszero(a) || iszero(b) ? zero(T) : oneunit(T) @@ -47,9 +47,9 @@ discriminant(a) Calculate discriminant of a univariate polynomial `a`. """ function discriminant(a::UnivariatePolynomial) -la = LC(a) -s = iseven(deg(a) >> 1) ? la : -la -resultant(a, derive(a)) / s + la = LC(a) + s = iseven(deg(a) >> 1) ? la : -la + resultant(a, derive(a)) / s end """ @@ -62,50 +62,50 @@ See: [`WIKI: Subresultant`](https://en.wikipedia.org/wiki/Polynomial_greatest_co and TAoCP 2nd Ed. 4.6.1. """ function presultant_seq( -a::T, -b::T, -::Val{Usedet}, + a::T, + b::T, + ::Val{Usedet}, ) where {Usedet,S,T<:UnivariatePolynomial{S}} -E = one(S) -da = deg(a) -db = deg(b) -d = da - db -s = E -if d < 0 - a, b, da, db, d = b, a, db, da, -d - if isodd(da) && isodd(db) - s = -s - end -end -cc, a, b = normgcd(a, b) -iszero(b) && return b, cc, b -s *= cc^(da + db) -ψ = -E -β = iseven(d) ? -E : E -det = isodd(da) && isodd(db) ? -E : E -while true - γ = LC(b) - δ = γ^(d + 1) - a = a * δ - c = rem(a, b) - iszero(c) && break - dc = deg(c) - c /= β - # prepare for next turn - if Usedet - det = det * δ^db / γ^(da - dc) / β^db + E = one(S) + da = deg(a) + db = deg(b) + d = da - db + s = E + if d < 0 + a, b, da, db, d = b, a, db, da, -d + if isodd(da) && isodd(db) + s = -s + end end - if isodd(db) && isodd(dc) - det = -det + cc, a, b = normgcd(a, b) + iszero(b) && return b, cc, b + s *= cc^(da + db) + ψ = -E + β = iseven(d) ? -E : E + det = isodd(da) && isodd(db) ? -E : E + while true + γ = LC(b) + δ = γ^(d + 1) + a = a * δ + c = rem(a, b) + iszero(c) && break + dc = deg(c) + c /= β + # prepare for next turn + if Usedet + det = det * δ^db / γ^(da - dc) / β^db + end + if isodd(db) && isodd(dc) + det = -det + end + ψ = iszero(d) ? ψ : (-γ)^d / ψ^(d - 1) + a, b = b, c + da, db = db, dc + d = da - db + β = -γ * ψ^d end - ψ = iszero(d) ? ψ : (-γ)^d / ψ^(d - 1) - a, b = b, c - da, db = db, dc - d = da - db - β = -γ * ψ^d -end -r = db > 0 ? zero(b) : d < 1 ? one(b) : d == 1 ? b : LC(b)^(d - 1) * b -b, cc, r * s / det + r = db > 0 ? zero(b) : d < 1 ? one(b) : d == 1 ? b : LC(b)^(d - 1) * b + b, cc, r * s / det end """ @@ -117,62 +117,62 @@ Its use is restricted to the case of `S` is an integral domain (there is no non-trivial divisor of zero). Effort is `O(deg(p)*deg(q))`. """ function signed_subresultant_polynomials(P::T, Q::T) where {S,T<:UnivariatePolynomial{S}} -# epsi(n) = (-1) ^ (n*(n-1)÷2) -epsi(n::Int) = iseven(n >> 1) ? 1 : -1 -p, q = deg(P), deg(Q) -p > q || throw(ArgumentError("degree of polynomials: $p is not > $q")) -sresp = zeros(T, p + 1) -s = zeros(S, p + 1) -t = zeros(S, p + 1) -sresp[p+1] = P -s[p+1] = t[p+1] = 1 # (sign(LC(P)) -sresp[p] = Q -bq = LC(Q) -sq = bq -t[p] = bq -if p > q - 1 - bqp = bq^(p - q - 1) - sq = bqp * epsi(p - q) - sresp[q+1] = sq * Q - sq *= bq -end -s[q+1] = sq -i = p + 1 -j = p -while j > 0 && !iszero(sresp[j]) - k = deg(sresp[j]) - if k == j - 1 - s[j] = t[j] - if k > 0 - sresp[k] = -rem(s[j]^2 * sresp[i], sresp[j]) / (s[j+1] * t[i]) - end - elseif k < j - 1 - s[j] = 0 - sig = -1 - for d = 1:j-k-1 - t[j-d] = (t[j] * t[j-d+1]) / s[j+1] * sig - sig = -sig - s[k+1] = t[k+1] - sresp[k+1] = s[k+1] * sresp[j] / t[j] - end - for l = j-2:-1:k+1 - sresp[l+1] = 0 - s[l+1] = 0 + # epsi(n) = (-1) ^ (n*(n-1)÷2) + epsi(n::Int) = iseven(n >> 1) ? 1 : -1 + p, q = deg(P), deg(Q) + p > q || throw(ArgumentError("degree of polynomials: $p is not > $q")) + sresp = zeros(T, p + 1) + s = zeros(S, p + 1) + t = zeros(S, p + 1) + sresp[p+1] = P + s[p+1] = t[p+1] = 1 # (sign(LC(P)) + sresp[p] = Q + bq = LC(Q) + sq = bq + t[p] = bq + if p > q - 1 + bqp = bq^(p - q - 1) + sq = bqp * epsi(p - q) + sresp[q+1] = sq * Q + sq *= bq + end + s[q+1] = sq + i = p + 1 + j = p + while j > 0 && !iszero(sresp[j]) + k = deg(sresp[j]) + if k == j - 1 + s[j] = t[j] + if k > 0 + sresp[k] = -rem(s[j]^2 * sresp[i], sresp[j]) / (s[j+1] * t[i]) + end + elseif k < j - 1 + s[j] = 0 + sig = -1 + for d = 1:j-k-1 + t[j-d] = (t[j] * t[j-d+1]) / s[j+1] * sig + sig = -sig + s[k+1] = t[k+1] + sresp[k+1] = s[k+1] * sresp[j] / t[j] + end + for l = j-2:-1:k+1 + sresp[l+1] = 0 + s[l+1] = 0 + end + if k > 0 + sresp[k] = -rem((t[j] * s[k+1]) * sresp[i], sresp[j]) / (s[j+1] * t[i]) + end end if k > 0 - sresp[k] = -rem((t[j] * s[k+1]) * sresp[i], sresp[j]) / (s[j+1] * t[i]) + t[k] = LC(sresp[k]) end + i, j = j, k end - if k > 0 - t[k] = LC(sresp[k]) + for l = 0:j-2 + sresp[l+1] = 0 + s[l+1] = 0 end - i, j = j, k -end -for l = 0:j-2 - sresp[l+1] = 0 - s[l+1] = 0 -end -sresp, s + sresp, s end """ @@ -182,43 +182,43 @@ Extended pseudo GCD algorithm. Return `g == pgcd(a, b)` and `u, v, f` with `a * u + b * v == g * f`. """ function pgcdx(a::T, b::T) where {S,T<:UnivariatePolynomial{S}} -E = one(S) -iszero(b) && return a, E, zero(S), a -iszero(a) && return b, zero(S), E, b -da = deg(a) -db = deg(b) -d = da - db -if d < 0 - g, u, v, f = pgcdx(b, a) - return g, v, u, f -end -cc, a, b = normgcd(a, b) -ψ = -E -β = iseven(d) ? E : -E -EE = one(T) -ZZ = zero(T) -s1, s2, t1, t2 = EE, ZZ, ZZ, EE -while true - γ = LC(b) - γd = γ^(d + 1) - a = a * γd - q, c, h = pdivrem(a, b) - c /= β - a, b = b, c - iszero(b) && break - s1, s2 = s2, (s1 * γd - s2 * q) / β - t1, t2 = t2, (t1 * γd - t2 * q) / β - # prepare for next turn - da = db - db = deg(c) - ψ = d == 0 ? ψ : (-γ)^d / ψ^(d - 1) + E = one(S) + iszero(b) && return a, E, zero(S), a + iszero(a) && return b, zero(S), E, b + da = deg(a) + db = deg(b) d = da - db - β = -γ * ψ^d -end -cs = gcd(content(s2), content(t2)) -a = a / cs -f = content(a) -a / (f / cc), s2 / cs, t2 / cs, f + if d < 0 + g, u, v, f = pgcdx(b, a) + return g, v, u, f + end + cc, a, b = normgcd(a, b) + ψ = -E + β = iseven(d) ? E : -E + EE = one(T) + ZZ = zero(T) + s1, s2, t1, t2 = EE, ZZ, ZZ, EE + while true + γ = LC(b) + γd = γ^(d + 1) + a = a * γd + q, c, h = pdivrem(a, b) + c /= β + a, b = b, c + iszero(b) && break + s1, s2 = s2, (s1 * γd - s2 * q) / β + t1, t2 = t2, (t1 * γd - t2 * q) / β + # prepare for next turn + da = db + db = deg(c) + ψ = d == 0 ? ψ : (-γ)^d / ψ^(d - 1) + d = da - db + β = -γ * ψ^d + end + cs = gcd(content(s2), content(t2)) + a = a / cs + f = content(a) + a / (f / cc), s2 / cs, t2 / cs, f end """ @@ -227,21 +227,21 @@ sylvester_matrix(u::P, v::P) where P<:UnivariatePolynomial Return sylvester matrix of polynomials `u` and `v`. """ function sylvester_matrix(v::P, u::P) where {Z,P<:UnivariatePolynomial{Z}} -nu = deg(u) -nv = deg(v) -n = max(nu, 0) + max(nv, 0) -S = zeros(Z, n, n) -if nv >= 0 - for k = 1:nu - S[k:k+nv, k] .= reverse(v.coeff) + nu = deg(u) + nv = deg(v) + n = max(nu, 0) + max(nv, 0) + S = zeros(Z, n, n) + if nv >= 0 + for k = 1:nu + S[k:k+nv, k] .= reverse(v.coeff) + end end -end -if nu >= 0 - for k = 1:nv - S[k:k+nu, k+nu] .= reverse(u.coeff) + if nu >= 0 + for k = 1:nv + S[k:k+nu, k+nu] .= reverse(u.coeff) + end end -end -S + S end """ @@ -250,8 +250,8 @@ resultant_naive(u, v) Reference implementation of resultant (determinant of sylvester matrix) """ function resultant_naive(u::P, v::P) where {Z,P<:UnivariatePolynomial{Z}} -S = sylvester_matrix(u, v) -det(S) + S = sylvester_matrix(u, v) + det(S) end """ @@ -260,8 +260,8 @@ g, ag, bg = normgcd(a, b) Divided `a` and `b` by the gcd of their contents. """ function normgcd(a, b) -ca = content(a) -cb = content(b) -g = gcd(ca, cb) -isunit(g) ? (one(g), a, b) : (g, a / g, b / g) + ca = content(a) + cb = content(b) + g = gcd(ca, cb) + isunit(g) ? (one(g), a, b) : (g, a / g, b / g) end diff --git a/src/zz.jl b/src/zz.jl index 2dc5cc1..eeb02e1 100644 --- a/src/zz.jl +++ b/src/zz.jl @@ -16,6 +16,13 @@ ZZ(a::T) where T<:Integer = ZZ{T}(a) ZZ{T}(a::ZZ) where T = ZZ{T}(T(a.val)) ZZ(a::ZZ{T}) where T = copy(a) +function ZZ{T}(a::Union{QQ{T},Frac{ZZ{T}}}) where T + a.den != 1 && throw(InexactError(:ZZ, ZZ{T}, a)) + ZZ(a.num) +end +ZZ(a::Union{QQ{T},Frac{ZZ{T}}}) where T = ZZ{T}(a) +ZZ{S}(a::Union{QQ{T},Frac{ZZ{T}}}) where {S,T} = ZZ(promote_type(S,T)(a)) + # promotion and conversion promote_rule(::Type{ZZ{T}}, ::Type{ZZ{S}}) where {S,T} = ZZ{promote_type(S, T)} promote_rule(::Type{ZZ{T}}, ::Type{QQ{S}}) where {S,T} = QQ{promote_type(S, T)} diff --git a/test/linearalgebra.jl b/test/linearalgebra.jl index 5d94c91..9aab6d0 100644 --- a/test/linearalgebra.jl +++ b/test/linearalgebra.jl @@ -79,14 +79,13 @@ end @test length(propertynames(lut)) == 10 end -@testset "determinant" begin - G = ZZ/89 +@testset "determinant and adjugate type $G" for G in (ZZ/89, ZZ{BigInt}, QQ{BigInt}) A = G.([1 2 3; 4 5 6; 7 8 9]) B = G.([1 2 3; 4 5 6; 7 8 8]) @test iszero(det(A)) @test iszero(adjugate(A) * A) @test !iszero(det(B)) - @test adjugate(B) ./ det(B) == inv(B) + @test adjugate(B) * B == det(B) * I(size(B, 1)) p = characteristic_polynomial(B) @test iszero(p(B)) diff --git a/test/zz.jl b/test/zz.jl index 17b1204..8ac34d5 100644 --- a/test/zz.jl +++ b/test/zz.jl @@ -17,7 +17,7 @@ using CommutativeRings @test ZZ(1) + QQ(1, 2) == QQ(3, 2) @test ZZ(Int8(1)) + QQ(1, 2) == QQ(3, 2) @test typeof(ZZ(1) + QQ(1, 2)) == QQ{Int} - @test typeof(ZZ(1) + 1//2) == QQ{Int} + @test typeof(ZZ(1) + 1 // 2) == QQ{Int} @test hash(ZZ(big"123")) == hash(123) end @@ -80,4 +80,12 @@ end @test eltype(factor(z3)) == Pair{Z,Int} end +@testset "conversion from fractions type $T" for T in + (QQ{Int}, Frac{ZZ{Int}}, Rational{Int}) + a = T(1 // 1) + @test_throws InexactError Int(a) + b = T(12) + @test Int(b) == 12 +end + end # module