fixed SubResultant, now works with integral domains

Author: HaraldHofstaetter

src/general.jl

@@ -28,7 +28,7 @@ end
 Base.show(io::IO, t::IdTerm) = show(io, t.arg)
 
 struct FunctionTerm <: Term
-    op::Function
+    op # ::Function
     coeff::RingElement
     arg::RingElement
 end
@@ -252,12 +252,9 @@ function PartialFraction(a::T, d::Vector{T}, e::Vector{Int}) where T <: RingElem
     a0, A
 end
 
-function SubResultant(A::PolyElem{T}, B::PolyElem{T}) where T <: FieldElement
+function SubResultant(A::PolyElem{T}, B::PolyElem{T}) where T <: RingElement
     # See Bronstein's book, Section 1.5, p. 24 
-    # Note: This implementation requires that A, B are polynomials over a field
-    # (and not mereley over an integral domain), because some intermediate
-    # calculations use divisons. For example, negative exponents of γ[i-1] 
-    # can occur, or the division in the definition of q.
+    degree(A) >= degree(B) || error("degree(A) must be >= degree(B)")
     T_one = one(leading_coefficient(A)) # 1 of type T
     Rs = [A, B] # indices shifted! 
     i = 1
@@ -271,7 +268,8 @@ function SubResultant(A::PolyElem{T}, B::PolyElem{T}) where T <: FieldElement
         q = divexact(R, β[i])        
         push!(Rs, q)
         i += 1
-        push!(γ, (-r[i-1])^δ[i-1]*γ[i-1]^(1-δ[i-1]) )
+        h = δ[i-1]<=1 ? (-r[i-1])^δ[i-1]*γ[i-1]^(1-δ[i-1]) : divexact((-r[i-1])^δ[i-1], γ[i-1]^(δ[i-1]-1))
+        push!(γ, h)
         push!(δ, degree(Rs[i-1+1]) - degree(Rs[i+1]) )
         push!(β, -r[i-1]*γ[i]^δ[i])       
     end
@@ -295,22 +293,6 @@ function SubResultant(A::PolyElem{T}, B::PolyElem{T}) where T <: FieldElement
     constant_coefficient((s*c)*Rs[k+1]^degree(Rs[k-1+1])), vcat(Rs[1:k+1], zero_poly)
 end
 
-function SubResultant(A::PolyElem{T}, B::PolyElem{T}) where T <: RingElement
-    # If we have polynomials over integral domain we upgrade them to polynomials
-    # over the fraction field. Then we call the above implementation of SubResultant 
-    # for polynomials over fields. Finally we downgrade the results to elements
-    # over the base integral domain.
-    Z = zero(leading_coefficient(A))//one(leading_coefficient(A))
-    # upgrade A, B \in R[t] to elements a, b of K[t] where K is the fraction ring of R
-    a = map_coefficients(c->(c+Z), A)
-    b = map_coefficients(c->(c+Z), B)
-    r0, Rs0 = SubResultant(a, b) 
-    # downgrade r, RS from elements of K[t] to elements  of R[t]
-    r = numerator(r0)
-    Rs = [map_coefficients(c->numerator(c), R) for R in Rs0]
-    r, Rs
-end
-
 function Squarefree_Musser(A::PolyElem{T}) where T <: RingElement
     # See Bronstein's book, Section 1.7, p. 29 
     c = content(A)