Use shared MP monomials and SA terms and polynomials

Project.toml

@@ -17,9 +17,13 @@ StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 LinearAlgebra = "1.10"
 MacroTools = "0.5"
 MultivariatePolynomials = "0.5.15"
-MutableArithmetics = "1"
+MutableArithmetics = "1.9"
 Pkg = "1.10"
 Random = "1.10"
 Reexport = "1"
-StarAlgebras = "0.3"
+StarAlgebras = "0.4"
 julia = "1.10"
+
+[sources]
+MultivariatePolynomials = {url = "https://github.com/JuliaAlgebra/MultivariatePolynomials.jl", rev = "bl/sa_poly"}
+StarAlgebras = {url = "https://github.com/JuliaAlgebra/StarAlgebras.jl", rev = "bl/term"}

README.md

@@ -30,30 +30,32 @@ julia> typeof(y)
 TypedPolynomials.Variable{:y}
 ```
 
-Multiplying variables creates a `Monomial{V}` where `V` is the vector of variables contained in the monomial:
+Multiplying variables creates a `Monomial{V}` where `V` is the tuple type of the
+variables contained in the monomial. Monomials share MP's implementation, while
+terms and polynomials use StarAlgebras:
 
 ```julia
 julia> x * y
 xy
 
-julia> typeof(x * y)
-TypedPolynomials.Monomial{(x, y),2}
+julia> typeof(x * y) === Monomial{typeof((x, y)),2}
+true
 
-julia> typeof(x^2)
-TypedPolynomials.Monomial{(x,),1}
+julia> typeof(x^2) === Monomial{typeof((x,)),1}
+true
 ```
 
-Multiplying a monomial (or variable) by anything other than another `Variable` or `Monomial` creates a `Term`:
+Multiplying a monomial (or variable) by a number creates a `Term`:
 
 ```julia
 julia> 3 * x
 3x
 
-julia> typeof(3 * x)
-TypedPolynomials.Term{Int64,TypedPolynomials.Monomial{(x,),1}}
+julia> 3 * x isa Term{Int}
+true
 
-julia> typeof(3.0 * x^2 * y)
-TypedPolynomials.Term{Float64,TypedPolynomials.Monomial{(x, y),2}}
+julia> 3.0 * x^2 * y isa Term{Float64}
+true
 ```
 
 Addition or subtraction of terms, monomials, and/or variables creates a `Polynomial`:
@@ -73,6 +75,9 @@ x^3 + 3.0y^2 + x + 0.5z
 
 ### Differentiation and Substitution
 
+Use `p(...)` to provide values for every variable, and `subs(p, ...)` for
+partial substitution.
+
 ```julia
 using TypedPolynomials
 using Test
@@ -81,7 +86,7 @@ p = 2x + 3.0x*y^2 + y
 @test differentiate(p, x) == 3y^2 + 2 # compute the derivative of p with respect to x
 @test differentiate.(p, (x, y)) == (3y^2 + 2, 6*x*y + 1) # compute the gradient of p
 @test p((x, y)=>(y, x)) == 2y + 3y*x^2 + x  # replace any x by y and y by x
-@test p(y=>x^2) == 2x + 3x*(x^4) + x^2 # replace any occurence of y by x^2
+@test subs(p, y=>x^2) == 2x + 3x*(x^4) + x^2 # replace any occurrence of y by x^2
 @test p(x=>1, y=>2) == 2 * 1 + 3 * 1 * 2^2 + 2 # evaluate p at [1, 2]
 ```
 
@@ -94,6 +99,6 @@ using TypedPolynomials
 A = rand(3, 3)
 @polyvar x[1:3] # assign x to a tuple of variables x1, x2, x3
 p = sum(x .* x) # x_1^2 + x_2^2 + x_3^2
-p(x[1]=>2, x[3]=>3) # x_2^2 + 13
+subs(p, x[1]=>2, x[3]=>3) # x_2^2 + 13
 p(x=>A*vec(x)) # corresponds to dot(A*x, A*x), need vec to convert the tuple to a vector
 ```

src/TypedPolynomials.jl

@@ -6,36 +6,47 @@ const MA = MutableArithmetics
 using Reexport
 @reexport using MultivariatePolynomials
 const MP = MultivariatePolynomials
-const Term = MP.Term
-const Polynomial = MP.Polynomial
+import StarAlgebras as SA
 export Term, Polynomial
 
 using MacroTools
-import Base: *, +, -, /, ^, ==,
-    promote_rule, show, isless, size, getindex,
-    one, zero, iszero, @pure, copy, exponent, vec
+import Base:
+    *,
+    +,
+    -,
+    /,
+    ^,
+    ==,
+    promote_rule,
+    show,
+    isless,
+    size,
+    getindex,
+    one,
+    zero,
+    iszero,
+    @pure,
+    copy,
+    exponent,
+    vec
 
 using LinearAlgebra
 import LinearAlgebra: dot, adjoint
 export @polyvar,
-       Variable,
-       Monomial,
-       Term,
-       Polynomial,
-       variables,
-       terms,
-       differentiate,
-       antidifferentiate,
-       subs
+    Variable,
+    Monomial,
+    Term,
+    Polynomial,
+    variables,
+    terms,
+    differentiate,
+    antidifferentiate,
+    subs
 
 include("types.jl")
-include("operators.jl")
 include("substitution.jl")
-include("calculus.jl")
 include("conversion.jl")
 include("promotion.jl")
-include("call.jl")
 include("macros.jl")
-include("sa.jl")
 
 end # module

src/conversion.jl

@@ -1,61 +1,23 @@
-Base.convert(::Type{M}, v::Variable) where {Vars, M <: Monomial{Vars}} = convert(M, Monomial(v))
-Base.convert(::Type{Monomial}, v::Variable) = Monomial(v)
-
-Base.convert(T::Type{Monomial{V}}, m::Monomial) where {V} = convert(Monomial{V, nvariables(T)}, m)
-
-@pure function matchindices(::Type{Monomial{V1, N1}}, ::Type{Monomial{V2, N2}}) where {V1, N1, V2, N2}
-    i2 = 1
-    inds = ntuple(i -> begin
-        if i2 > N2
-            0
-        elseif V1[i] === V2[i2]
-            i2 += 1
-            i2 - 1
-        else
-            0
-        end
-    end, Val{N1}())
-    if i2 <= N2
-        throw(InexactError(:matchindices, Monomial{V1, N1}, Monomial{V2, N2}))
-    end
-    inds
+function Base.convert(::Type{M}, v::Variable) where {M<:Monomial}
+    return convert(M, monomial(v))
 end
 
-function Base.convert(::Type{Monomial{V1, N1}}, m::Monomial) where {V1, N1}
-    inds = matchindices(Monomial{V1, N1}, typeof(m))
-    exps = ntuple(i -> begin
-        @inbounds ii = inds[i]
-        if ii == 0
-            0
-        else
-            m.exponents[ii]
-        end
-    end, Val{N1}())
-    Monomial{V1, N1}(exps)
+function Base.convert(::Type{A}, a::SA.StarAlgebra) where {A<:_Algebra}
+    return MP.algebra(MP.FullBasis{MP.Monomial}(variables(monomial_type(A))))
 end
+Base.convert(::Type{A}, a::A) where {A<:_Algebra} = a
+Base.convert(::Type{Monomial}, v::Variable) = monomial(v)
+Base.convert(::Type{Monomial}, m::Monomial) = m
+Base.convert(::Type{M}, m::M) where {M<:Monomial} = m
 
-function Base.convert(::Type{Variable}, mono::Monomial)
-    var = _mono2var(mono, powers(mono)...)
-    var === nothing && throw(InexactError(:convert, Variable, mono))
-    return var
-end
-
-function _checknovar(mono) end
-function _checknovar(mono, ve, ves...)
-    if iszero(ve[2])
-        _checknovar(mono, ves...)
-    else
-        throw(InexactError(:convert, Variable, mono))
+function Base.convert(::Type{M}, m::Monomial) where {M<:Monomial}
+    vars = variables(M)
+    for v in variables(m)
+        v in vars || throw(InexactError(:convert, M, m))
     end
-end
-function _mono2var(mono, ve, ves...)
-    if iszero(ve[2])
-        _mono2var(mono, ves...)
-    elseif isone(ve[2])
-        _checknovar(mono, ves...)
-        ve[1]
-    else
-        throw(InexactError(:convert, Variable, mono))
+    exps = map(vars) do v
+        i = findfirst(isequal(v), variables(m))
+        return isnothing(i) ? 0 : exponents(m)[i]
     end
+    return M(exps)
 end
-_mono2var(mono) = throw(InexactError(:convert, Variable, mono))

src/macros.jl

@@ -2,13 +2,14 @@ _new_variable(name::Symbol) = :(Variable{$(esc(Expr(:quote, name)))}())
 _assign_new_variable(name::Symbol) = :($(esc(name)) = $(_new_variable(name)))
 
 function _split(expr::Expr)
-    @capture(expr, var_Symbol[lb_Int:ub_Int]) || error("Unsupported syntax. Expected x[a:b] for literal integers a and b")
+    @capture(expr, var_Symbol[lb_Int:ub_Int]) ||
+        error("Unsupported syntax. Expected x[a:b] for literal integers a and b")
     var, lb, ub
 end
 
 function _assign_new_variable(expr::Expr)
     var, lb, ub = _split(expr)
-    :($(esc(var)) = $(Expr(:tuple, [_new_variable(Symbol("$var[$i]")) for i in lb:ub]...)))
+    :($(esc(var)) = $(Expr(:tuple, [_new_variable(Symbol("$var[$i]")) for i = lb:ub]...)))
 end
 
 _return_name(s::Symbol) = esc(s)