Unsafe funcs for FreeAssAlgebra, Poly, MPoly (#1847)

* Better conformance tests for FreeAssAlgebra
* Add `zero!`, `one!`, `neg!` for FreeAssAlgElem
* Add `neg!`, `one!` for Poly
* Add `zero!`, `one!`, `neg!` for MPoly
* Add `add!`, `sub!` for FreeAssAlgElem
* Fix `neg!` for `UnivPoly`

src/generic/FreeAssociativeAlgebra.jl

@@ -71,7 +71,7 @@ function isone(a::FreeAssociativeAlgebraElem{T}) where T
     if length(a) < 1
         return isone(zero(base_ring(a)))
     else
-        return a.length == 1 && isone(a.coeffs[1]) && isempty(a.exps[1])
+        return length(a) == 1 && isone(a.coeffs[1]) && isempty(a.exps[1])
     end
 end
 
@@ -646,6 +646,146 @@ function divexact(
     return combine_like_terms!(FreeAssociativeAlgebraElem{T}(R, zcoeffs, copy(a.exps), n))
 end
 
+###############################################################################
+#
+# Unsafe arithmetic functions
+#
+###############################################################################
+
+function zero!(a::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    a.length = 0
+    return a
+end
+
+function one!(a::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    a.length = 1
+    fit!(a, 1)
+    a.coeffs[1] = one(base_ring(parent(a)))
+    a.exps[1] = Int[]
+    return a
+end
+
+function neg!(a::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    for i in 1:length(a)
+        a.coeffs[i] = neg!(a.coeffs[i])
+    end
+    return a
+end
+
+function neg!(z::FreeAssociativeAlgebraElem{T}, a::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    if z === a
+        return neg!(a)
+    end
+    z.length = length(a)
+    fit!(z, length(a))
+    for i in 1:length(a)
+        if isassigned(z.coeffs, i)
+            z.coeffs[i] = neg!(z.coeffs[i], a.coeffs[i])
+        else
+            z.coeffs[i] = -a.coeffs[i]
+        end
+        # mutating z.exps[i] is not allowed since it could be aliased
+        z.exps[i] = a.exps[i]
+    end
+    return z
+end
+
+function add!(a::FreeAssociativeAlgebraElem{T}, b::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    iszero(b) && return a
+    return add!(zero(a), a, b)
+end
+
+function add!(z::FreeAssociativeAlgebraElem{T}, a::FreeAssociativeAlgebraElem{T}, b::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    if z === a
+        return add!(z, b)
+    elseif z === b
+        return add!(z, a)
+    end
+    z.coeffs = empty!(z.coeffs)
+    z.exps = empty!(z.exps)
+    i = j = 1
+    while i <= a.length && j <= b.length
+        c = word_cmp(a.exps[i], b.exps[j])
+        if c < 0
+            push!(z.coeffs, b.coeffs[j])
+            push!(z.exps, b.exps[j])
+            j += 1
+        elseif c > 0
+            push!(z.coeffs, a.coeffs[i])
+            push!(z.exps, a.exps[i])
+            i += 1
+        else
+            s = a.coeffs[i] + b.coeffs[j]
+            if !iszero(s)
+                push!(z.coeffs, s)
+                push!(z.exps, a.exps[i])
+            end
+            i += 1
+            j += 1
+        end
+    end
+    while i <= a.length
+        push!(z.coeffs, a.coeffs[i])
+        push!(z.exps, a.exps[i])
+        i += 1
+    end
+    while j <= b.length
+        push!(z.coeffs, b.coeffs[j])
+        push!(z.exps, b.exps[j])
+        j += 1
+    end
+    z.length = length(z.coeffs)
+    return z
+end
+
+function sub!(a::FreeAssociativeAlgebraElem{T}, b::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    iszero(b) && return a
+    return sub!(zero(a), a, b)
+end
+
+function sub!(z::FreeAssociativeAlgebraElem{T}, a::FreeAssociativeAlgebraElem{T}, b::FreeAssociativeAlgebraElem{T}) where T <: RingElement
+    if z === a
+        return sub!(z, b)
+    elseif z === b
+        return sub!(zero(a), a, b)
+    end
+    z.coeffs = empty!(z.coeffs)
+    z.exps = empty!(z.exps)
+    i = j = 1
+    while i <= a.length && j <= b.length
+        c = word_cmp(a.exps[i], b.exps[j])
+        if c < 0
+            push!(z.coeffs, -b.coeffs[j])
+            push!(z.exps, b.exps[j])
+            j += 1
+        elseif c > 0
+            push!(z.coeffs, a.coeffs[i])
+            push!(z.exps, a.exps[i])
+            i += 1
+        else
+            s = a.coeffs[i] - b.coeffs[j]
+            if !iszero(s)
+                push!(z.coeffs, s)
+                push!(z.exps, a.exps[i])
+            end
+            i += 1
+            j += 1
+        end
+    end
+    while i <= a.length
+        push!(z.coeffs, a.coeffs[i])
+        push!(z.exps, a.exps[i])
+        i += 1
+    end
+    while j <= b.length
+        push!(z.coeffs, -b.coeffs[j])
+        push!(z.exps, b.exps[j])
+        j += 1
+    end
+    z.length = length(z.coeffs)
+    return z
+end
+
 
 ################################################################################
 #

src/generic/MPoly.jl

@@ -3886,6 +3886,43 @@ end
 #
 ###############################################################################
 
+function zero!(a::MPoly{T}) where {T <: RingElement}
+   a.length = 0
+   return a
+end
+
+function one!(a::MPoly{T}) where {T <: RingElement}
+   a.length = 1
+   fit!(a, 1)
+   a.coeffs[1] = one(base_ring(a))
+   a.exps = zero(a.exps)
+   return a
+end
+
+function neg!(a::MPoly{T}) where {T <: RingElement}
+   for i in 1:length(a)
+      a.coeffs[i] = neg!(a.coeffs[i])
+   end
+   return a
+end
+
+function neg!(z::MPoly{T}, a::MPoly{T}) where {T <: RingElement}
+   if z === a
+      return neg!(a)
+   end
+   z.length = length(a)
+   fit!(z, length(a))
+   for i in 1:length(a)
+      if isassigned(z.coeffs, i)
+         z.coeffs[i] = neg!(z.coeffs[i], a.coeffs[i])
+      else
+         z.coeffs[i] = -a.coeffs[i]
+      end
+   end
+   z.exps[:,1:length(a)] .= a.exps[:,1:length(a)]
+   return z
+end
+
 function add!(a::MPoly{T}, b::MPoly{T}, c::MPoly{T}) where {T <: RingElement}
    t = b + c
    a.coeffs = t.coeffs
@@ -3933,11 +3970,6 @@ function fit!(a::MPoly{T}, n::Int) where {T <: RingElement}
    end
    return nothing
 end
-#
-function zero!(a::MPoly{T}) where {T <: RingElement}
-   a.length = 0
-   return a
-end
 
 @doc raw"""
     setcoeff!(a::MPoly{T}, i::Int, c::T) where T <: RingElement

src/generic/Poly.jl

@@ -198,6 +198,29 @@ function zero!(c::Poly{T}) where T <: RingElement
    return c
 end
 
+function one!(c::Poly{T}) where T <: RingElement
+   fit!(c, 1)
+   c = set_length!(c, 1)
+   c.coeffs[1] = one(base_ring(c))
+   return c
+end
+
+function neg!(a::Poly{T}) where T <: RingElement
+   for i in 1:length(a)
+      a.coeffs[i] = neg!(a.coeffs[i])
+   end
+   return a
+end
+
+function neg!(z::Poly{T}, a::Poly{T}) where T <: RingElement
+   fit!(z, length(a))
+   z = set_length!(z, length(a))
+   for i in 1:length(a)
+      z.coeffs[i] = neg!(z.coeffs[i], a.coeffs[i])
+   end
+   return z
+end
+
 function mul!(c::Poly{T}, a::Poly{T}, b::Poly{T}) where T <: RingElement
    lena = length(a)
    lenb = length(b)

src/generic/UnivPoly.jl

@@ -1015,7 +1015,7 @@ function one!(a::UnivPoly{T}) where {T <: RingElement}
 end
 
 function neg!(z::UnivPoly{T}, a::UnivPoly{T}) where {T <: RingElement}
-   z.p = neg!(a.p)
+   z.p = neg!(z.p, a.p)
    return z
 end
 

src/generic/imports.jl

@@ -173,6 +173,7 @@ import ..AbstractAlgebra: mul!
 import ..AbstractAlgebra: mul_classical
 import ..AbstractAlgebra: mul_karatsuba
 import ..AbstractAlgebra: mullow
+import ..AbstractAlgebra: neg!
 import ..AbstractAlgebra: number_of_columns
 import ..AbstractAlgebra: number_of_generators
 import ..AbstractAlgebra: number_of_rows

test/Rings-conformance-tests.jl

@@ -42,6 +42,18 @@ function test_elem(Rx::AbstractAlgebra.PolyRing)
    return Rx(elem_type(R)[test_elem(R) for i in 1:rand(0:6)])
 end
 
+function test_elem(Rx::AbstractAlgebra.MPolyRing)
+   R = base_ring(Rx)
+   num_gens = ngens(Rx)
+   iszero(num_gens) && return Rx(test_elem(R))
+   len_bound = 8
+   exp_bound = rand(1:5)
+   len = rand(0:len_bound)
+   coeffs = [test_elem(R) for _ in 1:len]
+   exps = [[rand(0:exp_bound) for _ in 1:num_gens] for _ in 1:len]
+   return Rx(coeffs, exps)
+end
+
 function test_elem(S::Union{AbstractAlgebra.MatSpace,
                             AbstractAlgebra.MatRing})
    R = base_ring(S)
@@ -72,6 +84,19 @@ function test_elem(Rx::AbstractAlgebra.SeriesRing)
    end
 end
 
+function test_elem(S::AbstractAlgebra.FreeAssociativeAlgebra)
+   f = S()
+   g = gens(S)
+   R = base_ring(S)
+   isempty(g) && return S(test_elem(R))
+   len_bound = 8
+   exp_bound = 6
+   for i in 1:rand(0:len_bound)
+      f += test_elem(R) * prod(rand(g) for _ in 1:rand(0:exp_bound); init = S(1))
+   end
+   return f
+end
+
 
 # helper
 function equality(a::T, b::T) where T <: AbstractAlgebra.NCRingElement

test/generic/FreeAssociativeAlgebra-test.jl

@@ -140,10 +140,6 @@ end
    @test !occursin("\n", sprint(show, R))
 end
 
-function test_elem(R::Generic.FreeAssociativeAlgebra{elem_type(ZZ)})
-   return rand(R, 0:4, 0:5, -10:10)
-end
-
 @testset "Generic.FreeAssociativeAlgebra.change_base_ring" begin
    F5, = residue_ring(ZZ, 5)
    R, varsR = polynomial_ring(F5, ["x"])
@@ -221,6 +217,11 @@ end
 end
 
 @testset "Generic.FreeAssociativeAlgebra.NCRing_interface" begin
-   test_NCRing_interface(free_associative_algebra(ZZ, 3)[1])
+   S, = free_associative_algebra(ZZ, 3)
+   test_NCRing_interface(S)
+
+   R, = QQ[:x, :y]
+   S, = free_associative_algebra(R, :z => 1:3)
+   test_NCRing_interface(S)
 end
 

test/generic/MPoly-test.jl

@@ -1704,3 +1704,18 @@ end
    @test leading_coefficient(f, 1) == coefficients(f, 1)[end]
    @test content(f, 1) == y
 end
+
+@testset "Generic.MPoly.Ring_interface" begin
+  S, = polynomial_ring(QQ, 0)
+  test_Ring_interface_recursive(S)
+
+  S, = polynomial_ring(QQ, 1)
+  test_Ring_interface_recursive(S)
+
+  S, = polynomial_ring(ZZ, 2)
+  test_Ring_interface_recursive(S)
+
+  R, = QQ[:x]
+  S, = polynomial_ring(R, :z => 1:3)
+  test_Ring_interface(S) # _recursive needs too many ressources
+end