Merge pull request #145 from JuliaAlgebra/bl/MA_div

Exploit mutability in divrem with MutableArithmetics

Author: blegat

Project.toml

@@ -11,7 +11,9 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
 
 [compat]
 DataStructures = "0.17.7"
+DynamicPolynomials = "0.3.12"
 MutableArithmetics = "0.2"
+TypedPolynomials = "0.2.8"
 julia = "1"
 
 [extras]

src/division.jl

@@ -32,33 +32,37 @@ Base.rem(f::APL, g::Union{APL, AbstractVector{<:APL}}; kwargs...) = divrem(f, g;
 
 proddiff(x, y) = x/y - x/y
 function Base.divrem(f::APL{T}, g::APL{S}; kwargs...) where {T, S}
-    rf = convert(polynomialtype(f, Base.promote_op(proddiff, T, S)), f)
-    q = r = zero(rf)
+    # `promote_type(typeof(f), typeof(g))` is needed for TypedPolynomials in case they use different variables
+    rf = convert(polynomialtype(promote_type(typeof(f), typeof(g)), Base.promote_op(proddiff, T, S)), MA.copy_if_mutable(f))
+    q = zero(rf)
+    r = zero(rf)
     lt = leadingterm(g)
     rg = removeleadingterm(g)
     lm = monomial(lt)
     while !iszero(rf)
         ltf = leadingterm(rf)
         if isapproxzero(ltf; kwargs...)
-            rf = removeleadingterm(rf)
+            rf = MA.operate!(removeleadingterm, rf)
         elseif divides(lm, ltf)
             qt = _div(ltf, lt)
-            q += qt
-            rf = removeleadingterm(rf) - qt * rg
+            q = MA.add!(q, qt)
+            rf = MA.operate!(removeleadingterm, rf)
+            rf = MA.operate!(MA.sub_mul, rf, qt, rg)
         elseif lm > monomial(ltf)
             # Since the monomials are sorted in decreasing order,
             # lm is larger than all of them hence it cannot divide any of them
-            r += rf
+            r = MA.add!(r, rf)
             break
         else
-            r += ltf
-            rf = removeleadingterm(rf)
+            r = MA.add!(r, ltf)
+            rf = MA.operate!(removeleadingterm, rf)
         end
     end
     q, r
 end
 function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T, S}
-    rf = convert(polynomialtype(f, Base.promote_op(proddiff, T, S)), f)
+    # `promote_type(typeof(f), eltype(g))` is needed for TypedPolynomials in case they use different variables
+    rf = convert(polynomialtype(promote_type(typeof(f), eltype(g)), Base.promote_op(proddiff, T, S)), MA.copy_if_mutable(f))
     r = zero(rf)
     q = similar(g, typeof(rf))
     for i in eachindex(q)
@@ -71,15 +75,16 @@ function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T
     while !iszero(rf)
         ltf = leadingterm(rf)
         if isapproxzero(ltf; kwargs...)
-            rf = removeleadingterm(rf)
+            rf = MA.operate!(removeleadingterm, rf)
             continue
         end
         divisionoccured = false
         for i in useful
             if divides(lm[i], ltf)
                 qt = _div(ltf, lt[i])
-                q[i] += qt
-                rf = removeleadingterm(rf) - qt * rg[i]
+                q[i] = MA.add!(q[i], qt)
+                rf = MA.operate!(removeleadingterm, rf)
+                rf = MA.operate!(MA.sub_mul, rf, qt, rg[i])
                 divisionoccured = true
                 break
             elseif lm[i] > monomial(ltf)
@@ -90,11 +95,11 @@ function Base.divrem(f::APL{T}, g::AbstractVector{<:APL{S}}; kwargs...) where {T
         end
         if !divisionoccured
             if isempty(useful)
-                r += rf
+                r = MA.add!(r, rf)
                 break
             else
-                r += ltf
-                rf = removeleadingterm(rf)
+                r = MA.add!(r, ltf)
+                rf = MA.operate!(removeleadingterm, rf)
             end
         end
     end

test/division.jl

@@ -69,6 +69,12 @@ const MP = MultivariatePolynomials
         @test q == [x, 0]
         @test r == 0
 
+        @testset "Issue DynamicPolynomials#62" begin
+            p = x^2 + x + 1
+            q = rem(p, [x^2-y])
+            @test q == x + y + 1
+        end
+
         @test (@inferred rem(x^2*y + (1+1e-10)*x*y + 1, [x^2 + x, y + 1])) == 1
         @test (@inferred rem(x^2*y + (1+1e-10)*x*y + 1, [x^2 + x, y + 1]; ztol=1e-11)) == -((1+1e-10)-1)x + 1
     end

test/polynomial.jl

@@ -173,4 +173,11 @@ const MP = MultivariatePolynomials
         p = im * x + 2im * x^2
         @test p * p == -x^2 - 4x^4 - 4x^3
     end
+
+    @testset "$f is a mutable copy, see issue DynamicPolynomials#62" for f in [zero, one]
+        p = 2x + 1
+        q = f(p)
+        q = MA.add!(q, 2y)
+        @test p == 2x + 1
+    end
 end