diff --git a/src/DynamicPolynomials.jl b/src/DynamicPolynomials.jl index f4eab05..5df9ab0 100644 --- a/src/DynamicPolynomials.jl +++ b/src/DynamicPolynomials.jl @@ -84,6 +84,7 @@ include("promote.jl") include("operators.jl") include("comp.jl") +include("anti_diff.jl") include("diff.jl") include("subs.jl") diff --git a/src/anti_diff.jl b/src/anti_diff.jl new file mode 100644 index 0000000..46ce047 --- /dev/null +++ b/src/anti_diff.jl @@ -0,0 +1,34 @@ +function _div_by_power(x::T, y::Int) where {T} + x / y +end + +function _div_by_power(x::T, y::Int)::Rational{T} where {T<:Integer} + x // y +end + +function MP.antidifferentiate(m::Monomial{V,M}, x::Variable{V,M}) where {V,M} + z = copy(m.z) + i = findfirst(isequal(x), MP.variables(m)) + if (i === nothing || i == 0) || m.z[i] == 0 + Monomial(MP.variables(m), z) * x + else + z[i] += 1 + (1 // (m.z[i] + 1)) * Monomial(MP.variables(m), z) + end +end + +function MP.antidifferentiate(p::Polynomial{V,M,T}, x::Variable{V,M}) where {V,M,T} + i = something(findfirst(isequal(x), MP.variables(p)), 0) + S = typeof(_div_by_power(zero(T), Int(1))) + if iszero(i) + x * p + else + Z = copy.(p.x.Z) + a = Vector{S}(undef, length(p.a)) + for j in 1:length(Z) + a[j] = _div_by_power(p.a[j], (Z[j][i] + 1)) + Z[j][i] += 1 + end + Polynomial(a, MonomialVector(MP.variables(p), Z)) + end +end diff --git a/test/mono.jl b/test/mono.jl index 25b2361..17dd672 100644 --- a/test/mono.jl +++ b/test/mono.jl @@ -141,6 +141,28 @@ import MultivariatePolynomials as MP @test m.z == [2, 1, 1, 0, 0] end + @testset "Antidifferentiation" begin + @ncpolyvar x y z + + m = x + mi = DynamicPolynomials.MP.antidifferentiate(m, y) + @test mi == x * y + + # Antidifferentiation is product => Integral coefficients + @test MP.coefficient_type(mi) == Int + + # General antidifferentiation => Rational coefficients + m = x^3 + mi = DynamicPolynomials.MP.antidifferentiate(m, x) + @test mi == (x^4 / 4) + @test MP.coefficient_type(mi) == Rational{Int} + + m = Monomial([x, y, z], [1, 2, 3]) + mi = DynamicPolynomials.MP.antidifferentiate(m, z) + @test mi == (x*y^2*z^4) / 4 + @test MP.coefficient_type(mi) == Rational{Int} + end + @testset "Evaluation" begin @polyvar x y @test (x^2 * y)(3, 2) == 18 diff --git a/test/poly.jl b/test/poly.jl index 9cd4f28..80c55a3 100644 --- a/test/poly.jl +++ b/test/poly.jl @@ -81,6 +81,34 @@ @inferred polynomial(2.0u, Int) end + @testset "Antiderivative" begin + @polyvar x y + + p = (x^2 + 4 * y^3) + @test MP.coefficient_type(p) == Int + + pi = DynamicPolynomials.antidifferentiate(p, y) + @test pi == (x^2 * y + y^4) + + pi = DynamicPolynomials.antidifferentiate(p, x) + @test MP.coefficient_type(pi) == Rational{Int} + + p = (1.0 * x^2 + 2.0 * y^2) + @test MP.coefficient_type(p) == Float64 + + pi = DynamicPolynomials.antidifferentiate(p, x) + @test MP.coefficient_type(pi) == Float64 + + p = 2 * y + pi = DynamicPolynomials.antidifferentiate(p, y) + @test pi == y^2 + + p = x^2 + pi = DynamicPolynomials.antidifferentiate(p, y) + @test pi == x^2 * y + + end + @testset "Evaluation" begin @polyvar x y @test (x^2 + y^3)(2, 3) == 31