Speedup for content (#200)

* Speedup for content

Co-authored-by: Yingbo Ma <[email protected]>

* Add gcd for terms

* Fixes

* Fix

* Add credits to SIMDPolynomials

* Fix test

* Exclude test where it was overflowing anyway

Co-authored-by: Yingbo Ma <[email protected]>

Author: blegat

src/gcd.jl

@@ -129,6 +129,10 @@ function Base.gcd(p1::APL{T}, p2::APL{S}, algo::AbstractUnivariateGCDAlgorithm=G
     return inflate(g, shift, defl)::MA.promote_operation(gcd, typeof(p1), typeof(p2))
 end
 
+function Base.gcd(t1::AbstractTermLike{T}, t2::AbstractTermLike{S}, algo::AbstractUnivariateGCDAlgorithm=GeneralizedEuclideanAlgorithm()) where {T, S}
+    return term(_coefficient_gcd(coefficient(t1), coefficient(t2)), gcd(monomial(t1), monomial(t2)))
+end
+
 # Inspired from to `AbstractAlgebra.deflation`
 function deflation(p::AbstractPolynomialLike)
     if iszero(p)
@@ -508,6 +512,18 @@ _simplifier(a, b, algo) = _gcd(a, b, algo)
 # which makes the size of the `BigInt`s grow significantly which slows things down.
 _simplifier(a::Rational, b::Rational, algo) = gcd(a.num, b.num) // gcd(a.den, b.den)
 
+# Largely inspired from from `YingboMa/SIMDPolynomials.jl`.
+function termwise_content(p::APL)
+    ts = terms(p)
+    length(ts) == 1 && return first(ts)
+    g = gcd(ts[1], ts[2])
+    isone(g) || for i in 3:length(ts)
+        g = gcd(g, ts[i])
+        isone(g) && break
+    end
+    return g
+end
+
 """
     content(poly::AbstractPolynomialLike{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}
 
@@ -524,11 +540,28 @@ function content(poly::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {T}
     P = MA.promote_operation(gcd, T, T)
     # This is tricky to infer a `content` calls `gcd` which calls `content`, etc...
     # To help Julia break the loop, we annotate the result here.
-    return reduce(
-        (a, b) -> _simplifier(a, b, algo),
-        coefficients(poly),
-        init = zero(P),
-    )::P
+    coefs = coefficients(poly)
+    length(coefs) == 0 && return zero(T)
+    length(coefs) == 1 && return first(coefs)
+    # Largely inspired from from `YingboMa/SIMDPolynomials.jl`.
+    if T <: APL
+        for i in eachindex(coefs)
+            if nterms(coefs[i]) == 1
+                g = gcd(termwise_content(coefs[1]), termwise_content(coefs[2]))
+                isone(g) || for i in 3:length(coefs)
+                    g = gcd(g, termwise_content(coefs[i]))
+                    isone(g) && break
+                end
+                return convert(P, g)
+            end
+        end
+    end
+    g = gcd(coefs[1], coefs[2])::P
+    isone(g) || for i in 3:length(coefs)
+        g = _simplifier(g, coefs[i], algo)::P
+        isone(g) && break
+    end
+    return g::P
 end
 function content(poly::APL{T}, algo::AbstractUnivariateGCDAlgorithm) where {T<:AbstractFloat}
     return one(T)

test/division.jl

@@ -151,7 +151,7 @@ function _mult_test(a, b)
 end
 function mult_test(expected, a, b, algo)
     g = @inferred MP._simplifier(a, b, algo)
-    @test g isa promote_type(polynomialtype(a), polynomialtype(b))
+    @test g isa Base.promote_typeof(a, b)
     _mult_test(expected, g)
 end
 function mult_test(expected, a::Number, b, algo)
@@ -213,7 +213,7 @@ function multivariate_gcd_test(::Type{T}, algo=GeneralizedEuclideanAlgorithm())
         y^2*z^3 + y*z^4 + y^3 + y^3*z - y - z,
         algo,
     )
-    if T != Int || (algo != GeneralizedEuclideanAlgorithm(false, false) && algo != GeneralizedEuclideanAlgorithm(true, false))
+    if T != Int
         test_relatively_prime(
             -3o*y^2*z^3 - 3*y^4 + y*z^3 + z^4 + 2*y^3 + 2*y^2*z - y,
             3o*y^3*z^3 - 2*y^5 + y^2*z^3 + y*z^4 + y^3 + y^3*z - y - z,