Division with remainder for power series (#1545)

* `divrem` for power series

* Add test for HNF/SNF of power series

src/AbsSeries.jl

@@ -738,6 +738,25 @@ function Base.inv(a::AbsPowerSeriesRingElem{T}) where T <: FieldElement
     return x
 end
 
+###############################################################################
+#
+#   Division with remainder
+#
+###############################################################################
+
+function Base.divrem(a::AbsPowerSeriesRingElem{T}, b::AbsPowerSeriesRingElem{T}) where {T <: FieldElement}
+   check_parent(a, b)
+   if length(b) == 0
+      throw(DivideError())
+   end
+   if valuation(a) < valuation(b)
+      return zero(parent(a)), a
+   end
+   # valuation(a) >= valuation(b), so the exact division works
+   q = divexact(a, b)
+   return q, a - q*b
+end
+
 ###############################################################################
 #
 #   Composition

src/RelSeries.jl

@@ -965,6 +965,25 @@ function Base.inv(a::RelPowerSeriesRingElem{T}) where T <: FieldElement
     return x
 end
 
+###############################################################################
+#
+#   Division with remainder
+#
+###############################################################################
+
+function Base.divrem(a::RelPowerSeriesRingElem{T}, b::RelPowerSeriesRingElem{T}) where {T <: FieldElement}
+   check_parent(a, b)
+   if pol_length(b) == 0
+      throw(DivideError())
+   end
+   if valuation(a) < valuation(b)
+      return zero(parent(a)), a
+   end
+   # valuation(a) >= valuation(b), so the exact division works
+   q = divexact(a, b)
+   return q, a - q*b
+end
+
 ###############################################################################
 #
 #   Composition

test/generic/AbsSeries-test.jl

@@ -1247,3 +1247,37 @@ end
       @test isequal(h, R())
    end
 end
+
+@testset "Generic.AbsSeries.euclidean" begin
+   R, x = power_series_ring(QQ, 20, "x", model = :capped_absolute)
+   S, y = power_series_ring(GF(5), 20, "y", model = :capped_absolute)
+
+   T, z = power_series_ring(GF(7), 20, "z", model = :capped_absolute)
+   @test_throws ErrorException divrem(y, z)
+
+   for (T, t) in [ (R, x), (S, y) ]
+      @test divrem(t^2, t) == (t, zero(T))
+      @test divrem(t + 1, t) == (zero(T), t + 1)
+      @test_throws DivideError divrem(t, zero(T))
+
+      @test gcd(t^2*(t - 1), t^2) == t^2
+      @test gcd((t + 1)*(t - 1), t + 1) == one(T)
+      @test gcdx(t^2*(t - 1), t^2) == (t^2, zero(T), one(T))
+      g, u, v = gcdx((t + 1)*(t - 1), t + 1)
+      @test u*(t + 1)*(t - 1) + v*(t + 1) == g
+
+      M = matrix(T, 4, 3, [ 0 0 0; t^3 + 1 t^2 0; 0 t^2 t^5; t^4 + 1 t^2 t^5 + t^3 ])
+      @test is_hnf(hnf(M))
+      H, U = hnf_with_transform(M)
+      @test is_hnf(H)
+      @test U*M == H
+      @test is_unit(det(U))
+
+      @test is_snf(snf(M))
+      S, U, V = snf_with_transform(M)
+      @test is_snf(S)
+      @test U*M*V == S
+      @test is_unit(det(U))
+      @test is_unit(det(V))
+   end
+ end

test/generic/RelSeries-test.jl

@@ -1283,3 +1283,37 @@ end
       @test isequal(h, R())
    end
 end
+
+@testset "Generic.RelSeries.euclidean" begin
+   R, x = power_series_ring(QQ, 20, "x", model = :capped_relative)
+   S, y = power_series_ring(GF(5), 20, "y", model = :capped_relative)
+
+   T, z = power_series_ring(GF(7), 20, "z", model = :capped_relative)
+   @test_throws ErrorException divrem(y, z)
+
+   for (T, t) in [ (R, x), (S, y) ]
+      @test divrem(t^2, t) == (t, zero(T))
+      @test divrem(t + 1, t) == (zero(T), t + 1)
+      @test_throws DivideError divrem(t, zero(T))
+
+      @test gcd(t^2*(t - 1), t^2) == t^2
+      @test gcd((t + 1)*(t - 1), t + 1) == one(T)
+      @test gcdx(t^2*(t - 1), t^2) == (t^2, zero(T), one(T))
+      g, u, v = gcdx((t + 1)*(t - 1), t + 1)
+      @test u*(t + 1)*(t - 1) + v*(t + 1) == g
+
+      M = matrix(T, 4, 3, [ 0 0 0; t^3 + 1 t^2 0; 0 t^2 t^5; t^4 + 1 t^2 t^5 + t^3 ])
+      @test is_hnf(hnf(M))
+      H, U = hnf_with_transform(M)
+      @test is_hnf(H)
+      @test U*M == H
+      @test is_unit(det(U))
+
+      @test is_snf(snf(M))
+      S, U, V = snf_with_transform(M)
+      @test is_snf(S)
+      @test U*M*V == S
+      @test is_unit(det(U))
+      @test is_unit(det(V))
+   end
+ end