Add antidifferentiation (#158)

* Add antidifferentiation for monomials

* Add antidifferentiation for polynomials

* Add support for rational coefficients in antidifferentiation

* Update src/anti_diff.jl

Co-authored-by: Benoît Legat <[email protected]>

* Remove superfluous print

* Use coefficient_type in polynomial antidifferentiaton tests

* Make monomial antidifferentiation coefficient types rationals

---------

Co-authored-by: Benoît Legat <[email protected]>

Author: chrhansk

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")
 

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

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

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