search

JuliaAlgebra / DynamicPolynomials.jl

Multivariate polynomials implementation of commutative and non-commutative variables
Other
60 stars 21 forks source link

Matrix inverse #49

Closed ahumenberger closed 5 years ago

ahumenberger commented 5 years ago 
using DynamicPolynomials, MultivariatePolynomials, LinearAlgebra
@polyvar n
M = [n 0; 0 n]
inv(M)

The above piece of code causes two errors which can be fixed by adding the following two lines. I assume that these additions belong to MultivariatePolynomials.jl, right?

RationalPoly{S,T}(x::Bool) where {S,T} = ifelse(x, one(RationalPoly{S,T}), zero(RationalPoly{S,T}))
LinearAlgebra.adjoint(r::RationalPoly) = adjoint(numerator(r))/adjoint(denominator(r))

However, after adding the above two methods I get another error which is actually coming from DynamicPolynomials:

ERROR: UndefVarError: C not defined
Stacktrace:
 [1] Polynomial(::Array{Union{},1}, ::Array{Union{},1}) at /Users/ahumenberger/.julia/packages/DynamicPolynomials/ISFei/src/poly.jl:42
 [2] polynomial(::Array{Union{},1}, ::Array{Union{},1}, ::MultivariatePolynomials.SortedUniqState) at /Users/ahumenberger/.julia/packages/DynamicPolynomials/ISFei/src/poly.jl:174
 [3] polynomial at /Users/ahumenberger/.julia/packages/MultivariatePolynomials/BxDgV/src/polynomial.jl:53 [inlined]
 [4] polynomial(::Array{Union{},1}, ::MultivariatePolynomials.SortedState) at /Users/ahumenberger/.julia/packages/MultivariatePolynomials/BxDgV/src/polynomial.jl:107
 [5] polynomial at /Users/ahumenberger/.julia/packages/MultivariatePolynomials/BxDgV/src/polynomial.jl:108 [inlined] (repeats 2 times)
 [6] adjoint(::Polynomial{true,Union{}}) at /Users/ahumenberger/.julia/packages/MultivariatePolynomials/BxDgV/src/operators.jl:117
 [7] \(::RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}}, ::RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}}) at ./operators.jl:536
 [8] naivesub!(::UpperTriangular{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},Array{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},2}}, ::Array{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},1}, ::Array{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},1}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:1148
 [9] naivesub! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:1141 [inlined]
 [10] ldiv! at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/bidiag.jl:511 [inlined]
 [11] ldiv!(::UpperTriangular{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},Array{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},2}}, ::Array{RationalPoly{Polynomial{true,Int64},Polynomial{true,Int64}},2}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/bidiag.jl:524
 [12] inv(::UpperTriangular{RationalPoly{Polynomial{true,Int64},Term{true,Int64}},Array{RationalPoly{Polynomial{true,Int64},Term{true,Int64}},2}}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/triangular.jl:630
 [13] inv(::Array{Term{true,Int64},2}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.1/LinearAlgebra/src/dense.jl:728
 [14] top-level scope at none:0
blegat commented 5 years ago

I assume that these additions belong to MultivariatePolynomials.jl, right?

Indeed!

However, after adding the above two methods I get another error which is actually coming from DynamicPolynomials:

The issue is with [6] adjoint(::Polynomial{true,Union{}}), it should be [6] adjoint(::Polynomial{true,Int64}) so it seems numerator is returning this which is weird

ahumenberger commented 5 years ago

I tracked down the issue. A method for dividing two RationalPoly is missing. That is the following piece of code does not work yet.

r = (n+1)/(n^2+1)
s = (n+2)/(n^2+2)
r / s

And the absence of this method also causes the problems above. Therefore, after adding something like

Base.:/(r::RationalPoly, s::RationalPoly) = (r.num * s.den) / (s.num * r.den)

everything works fine.

So, there are only changes to be made in MultivariatePolynomials.