FactoredElem.jl

FactoredElem.jl

using Random
using Test
using Nemo


import Base: +, -, *, /, //, ^, ==, gcd, copy, one, zero, iszero, isone
import Base: checked_add, checked_sub, checked_neg, checked_mul, checked_abs
import Nemo: AbstractAlgebra
import AbstractAlgebra: Field, parent, divexact, expressify, inv, isunit,
                        elem_type, parent_type, base_ring, mul!, addeq!


###############################################################################

mutable struct FactoredField{T <: Union{RingElement}} <: AbstractAlgebra.Field
    base_ring::AbstractAlgebra.Ring

    function FactoredField{T}(R::AbstractAlgebra.Ring) where T <: RingElement
        return new{T}(R)
    end
end

mutable struct FactoredElemTerm{T <: RingElement}
    base::T
    exp::Int
end

mutable struct FactoredElem{T <: RingElement} <: FieldElem
    unit::T
    terms::Vector{FactoredElemTerm{T}}
    parent::FactoredField{T}
end

function FactoredFractionField(R::AbstractAlgebra.Ring)
   return FactoredField{AbstractAlgebra.elem_type(R)}(R)
end

# boilerplate #################################################################

function parent(a::FactoredElem{T}) where T <: RingElement
    return a.parent
end

function parent_type(::Type{FactoredElem{T}}) where {T <: RingElement}
    return FactoredField{T}
end

function elem_type(::Type{FactoredField{T}}) where {T <: RingElement}
    return FactoredElem{T}
end

function base_ring(F::FactoredField{T}) where T <: RingElement
    return F.base_ring::AbstractAlgebra.parent_type(T)
end

function base_ring(a::FactoredElem{T}) where T <: RingElement
    return base_ring(parent(a))
end

function characteristic(F::FactoredField{T}) where T <: RingElement
   return characteristic(base_ring(F))
end

# printing ####################################################################

function Base.show(io::IO, F::FactoredField)
    print(io, "Factored fraction field of ", base_ring(F))
end

function expressify(a::FactoredElem; context = nothing)
    n = Expr(:call, :*, expressify(a.unit, context = context))
    d = Expr(:call, :*)
    for t in a.terms
        b = expressify(t.base; context = context)
        e = checked_abs(t.exp)
        push!((e == t.exp ? n : d).args, isone(e) ? b : Expr(:call, :^, b, e))
    end
    return length(d.args) < 2 ? n : Expr(:call, :/, n, d)
end

function Base.show(io::IO, a::FactoredElem)
    print(io, AbstractAlgebra.obj_to_string(a, context = io))
end

# basic constructors ##########################################################

function (F::FactoredField{T})() where T
    return FactoredElem{T}(zero(base_ring(F)), FactoredElemTerm{T}[], F)
end

function (F::FactoredField{T})(a::T) where T <: RingElem
    parent(a) == base_ring(F) || error("Could not coerce to FactoredElem")
    if iszero(a) || isunit(a)
        return FactoredElem{T}(a, FactoredElemTerm{T}[], F)
    else
        return FactoredElem{T}(one(base_ring(F)), [FactoredElemTerm{T}(a, 1)], F)
    end
end

function (F::FactoredField{T})(a::Integer) where T <: RingElem
    return F(base_ring(F)(a))
end

function (F::FactoredField{T})(a::Rational) where T <: RingElem
    return F(numerator(a))/F(denominator(a))
end

function (F::FactoredField{T})(a::FactoredElem{T}) where T
    F == parent(a) || error("Could not coerce to FactoredElem")
    return a
end

function (F::FactoredField{T})(a::Integer) where T <: Integer
    a = base_ring(F)(a)
    if iszero(a) || isunit(a)
        return FactoredElem{T}(a, FactoredElemTerm{T}[], F)
    else
        return FactoredElem{T}(one(base_ring(F)), [FactoredElemTerm{T}(a, 1)], F)
    end
end

function (F::FactoredField{T})(a::Rational) where T <: Integer
    return F(numerator(a))/F(denominator(a))
end

function Base.deepcopy_internal(a::FactoredElem{T}, dict::IdDict) where T <: RingElement
   return FactoredElem{T}(deepcopy(a.unit), deepcopy(a.terms), a.parent)
end


# implementation ##############################################################

function pow(a::T, e::Int) where T <: RingElement
    if e < 0
        return divexact(one(parent(a)), a)^checked_neg(e)
    else
        return a^e
    end
end

function gcd_cofactors(a::T, b::T) where T
    g = gcd(a, b)
    if iszero(g)
        return (g, a, b)
    else
        return (g, divexact(a, g), divexact(b, g))
    end
end

function copy(a::FactoredElemTerm{T}) where T
    return FactoredElemTerm{T}(a.base, a.exp)
end

function _bases_are_coprime(a::FactoredElem{T}) where T
    a = a.terms
    i = 1
    while i <= length(a)
        j = i + 1
        while j <= length(a)
            if !isunit(gcd(a[i].base, a[j].base))
                return false
            end
            j += 1
        end
        i += 1
    end
    return true
end

function _bases_are_nice(a::FactoredElem{T}) where T
    if !_bases_are_coprime(a)
        return false
    end
    for i in a.terms
        if isunit(i.base)
            return false
        end
    end
    return true
end

# z *= f^e with no normalization
function append_pow!(z::FactoredElem{T}, f::FactoredElem{T}, e::Int) where T
    z.unit *= f.unit^e
    append_pow!(z, f.terms, e)
end

function append_pow!(z::FactoredElem{T}, f::Vector{FactoredElemTerm{T}}, e::Int) where T
    for i in f
        ie = checked_mul(i.exp, e)
        if !iszero(ie)
            push!(z.terms, FactoredElemTerm{T}(i.base, ie))
        end
    end
end

# z *= a^e with normalization
function append_pow_normalize!(z::FactoredElem{T}, a::T, e::Int, i::Int) where T
    iszero(e) && return
    if iszero(a)
        z.unit = a
        empty!(z.terms)
        return
    end
    input_is_good = _bases_are_coprime(z)
    l = z.terms
    while i <= length(l) && !isunit(a)
        (g, lbar, abar) = gcd_cofactors(l[i].base, a)
        # (g*lbar)^l[i].exp * (g*abar)^e
        # (lbar)^l[i].exp * (g)^(l[i].exp+e) * (abar)^e
        if isunit(g)
            i += 1
        elseif isunit(lbar)
            a = abar
            z.unit *= pow(lbar, l[i].exp)
            l[i].base = g
            l[i].exp = checked_add(l[i].exp, e)
            if iszero(l[i].exp)
                l[i] = l[end]
                pop!(l)
            end
        elseif isunit(abar)
            z.unit *= pow(abar, e)
            l[i].base = lbar
            a = g
            e = checked_add(l[i].exp, e)
            if iszero(e)
                return
            end
        else
            a = abar
            append_pow_normalize!(z, g, e, i)
        end
    end
    if isunit(a)
        z.unit *= pow(a, e)
    else
        push!(l, FactoredElemTerm{T}(a, e))
    end
    @assert !input_is_good || _bases_are_coprime(z)
    return
end

# multiply l by p^e and return a such that
# p*a = b*g and gcd(a, g) == 1
function append_coprimefac!(l::FactoredElem{T}, b::T, e::Int, g::T) where T
    @assert !iszero(b)
    append_pow_normalize!(l, g, e, 1)
    a = b
    r = gcd(a, g)
    while !isunit(r)
        @assert !iszero(r)
        append_pow_normalize!(l, r, e, 1)
        a = divexact(a, r)
        r = gcd(r, a)
    end
    return a
end

function gcdhelper(b::FactoredElem{T}, c::FactoredElem{T}) where T
    F = parent(b)
    z = FactoredElem(one(base_ring(F)), FactoredElemTerm{T}[], F)   # gcd(b,c)
    bbar = b.unit                       # expanded form of eventual b/gcd(b,c)
    cbar = c.unit                       # expanded form of eventual c/gcd(b,c)
    b = map(copy, b.terms)
    c = map(copy, c.terms)
    i = 1
    while i <= length(b)
        j = 1
        while j <= length(c)
            (g, b[i].base, c[j].base) = gcd_cofactors(b[i].base, c[j].base)
            if !isunit(g)
                e = Base.checked_sub(b[i].exp, c[j].exp)
                if e >= 0
                    append_pow_normalize!(z, g, c[j].exp, 1)
                    if e > 0
                        push!(b, FactoredElemTerm{T}(g, e))
                    end
                else
                    append_pow_normalize!(z, g, b[i].exp, 1)
                    push!(c, FactoredElemTerm{T}(g, checked_neg(e)))
                end
            end
            if isunit(c[j].base)
                cbar *= pow(c[j].base, c[j].exp)
                c[j] = c[end]
                pop!(c)
            else
                j += 1
            end
        end
        if isunit(b[i].base)
            bbar *= pow(b[i].base, b[i].exp)
            b[i] = b[end]
            pop!(b)
        else
            i += 1
        end
    end
    i = 1
    while i <= length(b)
        if b[i].exp < 0
			append_pow_normalize!(z, b[i].base, b[i].exp, 1)
            cbar *= b[i].base^checked_neg(b[i].exp)
        else
            bbar *= b[i].base^b[i].exp            
        end
        i += 1
    end
    j = 1
    while j <= length(c)
        if c[j].exp < 0
			append_pow_normalize!(z, c[j].base, c[j].exp, 1)
            bbar *= c[j].base^checked_neg(c[j].exp)
        else
            cbar *= c[j].base^c[j].exp            
        end
        j += 1
    end
    return (z, bbar, cbar)
end

function normalize(a::FactoredElem{T}) where T
    z = FactoredElem{T}(a.unit, FactoredElemTerm{T}[], parent(a))
    if !iszero(z.unit)
        for i in a.terms
			append_pow_normalize!(z, a[j].base, a[j].exp, 1)
        end
    end
    return z
end


####################################################################

function one(F::FactoredField{T}) where T
    FactoredElem{T}(one(base_ring(F)), FactoredElemTerm{T}[], F)
end

function zero(F::FactoredField{T}) where T
    FactoredElem{T}(zero(base_ring(F)), FactoredElemTerm{T}[], F)
end

function *(a::FactoredElem{T}, b::T) where T <: RingElem
    F = parent(a)
    parent(b) == base_ring(F) || error("oops")
    if iszero(b)
        return zero(F)
    else
        z = FactoredElem{T}(a.unit, map(copy, a.terms), F)
        append_pow_normalize!(z, b, 1, 1)
        return z
    end
end

function *(a::FactoredElem{T}, b::Union{Integer, Rational}) where T <: RingElem
    F = parent(a)
    b = base_ring(F)(b)
    if iszero(b)
        return zero(F)
    else
        z = FactoredElem{T}(a.unit, map(copy, a.terms), F)
        append_pow_normalize!(z, b, 1, 1)
        return z
    end
end

function *(b::T, a::FactoredElem{T}) where T <: RingElem
   return a*b
end

function *(b::Union{Integer, Rational}, a::FactoredElem{T}) where T <: RingElem
   return a*b
end


function *(b::FactoredElem{T}, c::FactoredElem{T}) where T <: RingElement
    parent(b) == parent(c) || error("oops")
    input_is_good = _bases_are_coprime(b) && _bases_are_coprime(b)
    z = FactoredElem{T}(b.unit*c.unit, FactoredElemTerm{T}[], parent(b))
    if iszero(z.unit)
        return z
    end
    b = map(copy, b.terms)
    c = map(copy, c.terms)
    i = 1
    while i <= length(b)
        j = 1
        while j <= length(c)
            (g, b[i].base, c[j].base) = gcd_cofactors(b[i].base, c[j].base)
            if isunit(g)
                z.unit *= pow(g, checked_add(b[i].exp, c[j].exp))
            else
                b[i].base = append_coprimefac!(z, b[i].base, b[i].exp, g)
                c[j].base = append_coprimefac!(z, c[j].base, c[j].exp, g)
            end
            if isunit(c[j].base)
                z.unit *= pow(c[j].base, c[j].exp)
                c[j] = c[end]
                pop!(c)
            else
                j += 1
            end
        end
        if isunit(b[i].base)
            z.unit *= pow(b[i].base, b[i].exp)
            b[i] = b[end]
            pop!(b)
        else
            i += 1
        end
    end
    append_pow!(z, b, 1)
    append_pow!(z, c, 1)
    @assert !input_is_good || _bases_are_coprime(z)
    return z
end

function inv(a::FactoredElem{T}) where T
    z = FactoredElemTerm{T}[]
    for i in a.terms
        push!(z, FactoredElemTerm{T}(i.base, checked_neg(i.exp)))
    end
    return FactoredElem{T}(inv(a.unit), z, parent(a))
end

function /(a::FactoredElem{T}, b::FactoredElem{T}) where T
    return a*inv(b)
end

function /(a::Union{Integer, Rational}, b::FactoredElem{T}) where T
    return parent(b)(a)/b
end

function /(a::T, b::FactoredElem{T}) where T
    return parent(b)(a)/b
end

function //(a::FactoredElem{T}, b::FactoredElem{T}) where T
    return a*inv(b)
end

function //(a::Integer, b::FactoredElem{T}) where T
    return parent(b)(a)//b
end

function //(a::T, b::FactoredElem{T}) where T
    return parent(b)(a)//b
end


function ^(a::FactoredElem{T}, b::Int) where T
    z = FactoredElemTerm{T}[]
    if !iszero(b)
        for t in a.terms
            push!(z, FactoredElemTerm(t.base, checked_mul(b, t.exp)))
        end
    end
    return FactoredElem{T}(pow(a.unit, b), z, parent(a))
end


function +(a::FactoredElem{T}, b::Union{Integer, Rational}) where T
   return a + parent(a)(b)
end

function +(a::FactoredElem{T}, b::T) where T
   return a + parent(a)(b)
end

function +(b::Union{Integer, Rational}, a::FactoredElem{T}) where T
   return parent(a)(b) + a
end

function +(b::T, a::FactoredElem{T}) where T
   return parent(a)(b) + a
end

function +(a::FactoredElem{T}, b::FactoredElem{T}) where T
    (g, x, y) = gcdhelper(a, b)
    return g*(x + y)
end


function -(a::FactoredElem{T}, b::Union{Integer, Rational}) where T
   return a - parent(a)(b)
end

function -(a::FactoredElem{T}, b::T) where T
   return a - parent(a)(b)
end

function -(b::Union{Integer, Rational}, a::FactoredElem{T}) where T
   return parent(a)(b) - a
end

function -(b::T, a::FactoredElem{T}) where T
   return parent(a)(b) - a
end

function -(a::FactoredElem{T}, b::FactoredElem{T}) where T
    (g, x, y) = gcdhelper(a, b)
    return g*(x - y)
end

function -(a::FactoredElem{T}) where T
    return FactoredElem{T}(-a.unit, a.terms, a.parent)
end

function gcd(a::FactoredElem{T}, b::FactoredElem{T}) where T
    (g, x, y) = gcdhelper(a, b)
    return g*gcd(x, y)
end

function ==(a::FactoredElem{T}, b::FactoredElem{T}) where T
    (g, x, y) = gcdhelper(a, b)
    return x == y    
end

function isunit(a::FactoredElem{T}) where T
    return !iszero(a)
end

function iszero(a::FactoredElem{T}) where T
    return iszero(a.unit)
end

function isone(a::FactoredElem{T}) where T
    if iszero(a.unit)
        return false
    end
    for i in a.terms
        ok = !iszero(i.exp) && !isunit(i.base)
        for j in a.terms
            if j !== i
                ok = ok && isunit(gcd(i.base, j.base))
            end
        end
        ok && return false
    end
    z = normalize(a)
    return isempty(z.terms) && isone(z.unit)
end

function factor(a::FactoredElem{T}) where T
    b = normalize(a)
    z = FactoredElem{T}(b.unit, FactoredElem{T}())
    for i in a.terms
        f = factor(i.base)
        z.unit *= pow(f.unit, i.exp)
        for j in f.fac
            push!(z.terms, FactoredElemTerm{T}(j[1], checked_mul(j[2], i.exp)))
        end
    end
    return z
end


function mul!(a::FactoredElem{T}, b::FactoredElem{T}, c::FactoredElem{T}) where T <: RingElement
    z = b*c
    a.unit = z.unit
    a.terms = z.terms
    a.parent = z.parent
    return a
end

function addeq!(a::FactoredElem{T}, b::FactoredElem{T}) where T <: RingElement
    z = a + b
    a.unit = z.unit
    a.terms = z.terms
    a.parent = z.parent
    return a
end

print("running test code for frac(ZZ) ...")
let test_reps = 100

   Random.seed!(1234)

   FF = FactoredFractionField(ZZ)

   function randomFF(limit::Int)
       t = one(FF)
       for i in 1:abs(rand(Int)%limit)
           s = FF(rand(Int)%(20*limit))
           e = rand(Int)%limit
           t *= iszero(s) ? s^abs(e) : s^e
       end
       return t
   end

   for i in 1:test_reps
      local a, b
      a = randomFF(10)
      b = a - FF(1)
      if isone(a)
        @test iszero(b)
      else
        @test a != FF(1)
      end
      @test _bases_are_nice(a)
      @test _bases_are_nice(b)
   end

   for i in 1:test_reps
      local a, b, c
      a = randomFF(10)
      b = randomFF(10)
      c = a*b
      if iszero(b)
        d = c
        @test iszero(c)
      else
        d = c/b
        @test d == a
        if !iszero(a)
            e = c/a/b
            @test isone(e)
        end
      end
      @test _bases_are_nice(a)
      @test _bases_are_nice(b)
      @test _bases_are_nice(c)
      @test _bases_are_nice(d)
   end


   for i in 1:test_reps
      local a, b, c, d
      a = randomFF(10)
      b = randomFF(10)
      c = a + b
      d = c - b
      @test d == a
      @test _bases_are_nice(a)
      @test _bases_are_nice(b)
      @test _bases_are_nice(c)
      @test _bases_are_nice(d)
   end
end
println("done")



println("trying frac(ZZ[x])")
Zx, x = PolynomialRing(ZZ, "x")
F = FactoredFractionField(Zx)
@show F(x + 1)*F(x + 2)
@show F(x + 1)/F(x + 2)^2 + F(x)

println("trying frac(ZZ[x,y,z])")
Zxyz, (x, y, z) = PolynomialRing(ZZ, ["x", "y", "z"])
F = FactoredFractionField(Zxyz)
(x, y, z) = (F(x), F(y), F(z))
@show 1/(x+y) + 1/(x+y+z)
@show det(matrix(F, [x y z; x^2 y^2 z^2; x^3 y^3 z^3]))


############################# partial fractions ###############################

# prem(A, -B)
function prem(A, B)
   x = gen(parent(A))
   delta_org = degree(A) - degree(B) + 1
   while !iszero(A) && (delta = degree(A) - degree(B)) >= 0
      A = lead(A)*B*x^delta - lead(B)*A
      delta_org -= 1
   end
   return A*(-lead(B))^delta_org
end

# r = n1*A + n2*B
function premx(A, B)
   R = parent(A)
   r = A
   n1 = one(R)
   n2 = zero(R)
   @assert r == n1*A + n2*B
   x = gen(R)
   delta_org = degree(r) - degree(B) + 1
   while !iszero(r) && (delta = degree(r) - degree(B)) >= 0
      n1 =                 - lead(B)*n1
      n2 = lead(r)*x^delta - lead(B)*n2
      r = lead(r)*x^delta*B - lead(B)*r
      @assert r == n1*A + n2*B
      delta_org -= 1
   end
   r *= (-lead(B))^delta_org
   n1 *= (-lead(B))^delta_org
   n2 *= (-lead(B))^delta_org
   @assert r == n1*A + n2*B
   return (r, n1, n2)
end

function pgcd(P, Q)
   @assert degree(P) >= degree(Q) > 0
   s = lead(Q)^(degree(P) - degree(Q))
   A = Q
   B = prem(P, Q)
   Z = A
   while !iszero(B)
      d = degree(A)
      e = degree(B)
      @assert d > e
      Z = C = divexact(B*lead(B)^(d - e - 1), s^(d - e - 1))
      if iszero(e)
         return Z
      end
      B = divexact(prem(A, B), s^(d - e)*lead(A))
      A = C
      s = lead(A)
   end
   return Z
end


function pgcdx(P, Q)
   @assert degree(P) >= degree(Q) > 0
   R = parent(P)
   m11 = one(R)
   m12 = zero(R)
   m21 = zero(R)
   m22 = one(R)
   s = lead(Q)^(degree(P) - degree(Q))
   A = Q
   B, n1, n2 = premx(P, Q)
   m11, m12, m21, m22 = m21, m22, n1*m11 + n2*m21, n1*m12 + n2*m22
   @assert m11*P + m12*Q == A
   @assert m21*P + m22*Q == B
   Z = A
   while !iszero(B)
      d = degree(A)
      e = degree(B)
      @assert d > e
      u = lead(B)^(d - e - 1)
      v = s^(d - e - 1)
      w = s^(d - e)*lead(A)
      Z = C = divexact(B*u, v)
      R, n1, n2 = premx(A, B)
      B = divexact(R, w)
      A = C
      s = lead(A)
      m11, m12, m21, m22 = divexact(m21*u, v), divexact(m22*u, v),
                   divexact(m11*n1 + m21*n2, w), divexact(m12*n1 + m22*n2, w)
      @assert m11*P + m12*Q == A
      @assert m21*P + m22*Q == B
   end
   @assert m11*P + m12*Q == A
   @assert m21*P + m22*Q == B
   return Z, m11, m12, m21, m22
end


# rem = n1*A + n2*B with minimal scaling, requires gcd domain
function pseudoremx(A, B)
   R = parent(A)
   r = A
   n1 = one(base_ring(R))
   n2 = zero(R)
   @assert r == n1*A + n2*B
   x = gen(R)
   while !iszero(r) && (delta = degree(r) - degree(B)) >= 0
      (g, n, d) = gcd_cofactors(lead(r), lead(B))

      n1 =           -d*n1
      n2 = n*x^delta - d*n2
      r = n*x^delta*B - d*r
      @assert r == n1*A + n2*B
   end
   @assert r == n1*A + n2*B
   return (r, n1, n2)
end


# assert that the leading coefficient divisions are exact in divrem
function mydivrem(A, B)
   x = gen(parent(A))
   Q = zero(parent(A))
   R = A
   while !iszero(R) && (delta = degree(R) - degree(B)) >= 0
      q = divexact(lead(R), lead(B))*x^delta
      Q += q
      R -= q*B
   end
   @assert A == Q*B + R
   return Q, R
end

# [a0, a1, a2, ...] -> [a0, a1*s, a2*s^2, ...]
function myrescale!(a::Vector, s)
   t = one(parent(s))
   for i in 2:length(a)
      t *= s
      a[i] *= t
   end
end

# given a0 + a1*(f/z) + a2*(f/z)^2 + ..., reduce the ai mod f. z is a constant
function myreduce!(a::Vector, f, z)
   R = parent(z)
   n = length(a)
   q = zero(R)
   for i in 1:n
      a[i] += z*q
      q, a[i] = mydivrem(a[i], f)
   end
end

# power series multiplication mod (f/z)^n
function mymul(a::Vector, b::Vector, R)
   n = length(a)
   @assert n == length(b)
   r = zeros(R, n)
   for i in 1:n
      for j in 1:n+1-i
         r[i + j - 1] += a[i]*b[j]
      end
   end
   return r
end

mutable struct LazySum
   terms::Vector{Any}
end

function expressify(a::LazySum; context = nothing)
   sum = Expr(:call, :+)
   for t in a.terms
      push!(sum.args, expressify(t, context = context))
   end
   return sum
end

function Base.show(io::IO, a::LazySum)
    print(io, AbstractAlgebra.obj_to_string(a, context = io))
end


function partial_fractions(A::FactoredElem{AbstractAlgebra.Generic.Poly{T}}) where T
   Rx = base_ring(parent(A))
   a = A.unit
   b = elem_type(Rx)[]
   e = Int[]
   for t in A.terms
      if t.exp < 0
         push!(b, t.base)
         push!(e, checked_neg(t.exp))
      else
         a *= t.base^t.exp
      end
   end
   return pfrac(a, b, e, parent(A))
end

function pfrac(aa::E, b::Vector{E}, e::Vector{Int}, F) where E
   answer = LazySum([])

   R = parent(aa)
   n = length(b)
   @assert n == length(e)

   if degree(aa) >= sum(degree(b[i])*e[i] for i in 1:n)
      prodb = prod(b[i]^e[i] for i in 1:n)
      a, n1, n2 = premx(aa, prodb)
      push!(answer.terms, -F(n2)/F(n1))
   else
      n1 = one(R)
      n2 = zero(R)
      a = aa
   end

   # aa/b = -n2/n1 + a/(n1*b)

   for i in 1:n
# set rj = res(bi, bj) = uj*bi + vj*bj
# we have the expansion at bi as (where ck are elements of R)
#
#  r2^e2 * ... * rn^en    ei-1           bi       k-ei
# --------------------- =  sum ck*(--------------)   + ...
#  b1^e1 * ... * bn^en     k=0     prod_{j!=i} rj
#
# substitute bj = (rj - uj*bi)/vj, then the lhs becomes
#
#       prod_{j!=i} vj^ej
#  ----------------------------
#  prod_{j!=i} (1 - uj*bi/rj)^ej
#
# expand this as a power series in bi/prod_{j!=i} rj, being careful to do only
# integral arithmetic at every step
#
# Finally, multiplying by a and reducing mod bi might introduce some power
# of lead(bi) into the denominators of the ck, 0 <= k < ei

      c = [k == 1 ? one(R) : zero(R) for k in 1:e[i]]
      z = one(R)
      r = zeros(R, n)
      for j in 1:n
         if j == i
            continue
         end
         if degree(b[i]) < degree(b[j])
            rj, vj, uj, _, _ = pgcdx(b[j], b[i])
         else
            rj, uj, vj, _, _ = pgcdx(b[i], b[j])
         end
         iszero(degree(rj)) || error("bases not pairwise coprime")
         r[j] = rj
         cj = [vj for k in 1:e[i]]
         myrescale!(cj, uj)
         t = cj
         for k in 2:e[j]
            cj = mymul(cj, t, R)
            myreduce!(cj, b[i], rj)
         end
         myrescale!(c, rj)
         myrescale!(cj, z)
         z *= rj
         c = mymul(c, cj, R)
         myreduce!(c, b[i], z)
      end
      # we now have c[0] + c[1]*b[i]/z + c[2]*(b[i]/z)^2 + ...
      # just need to reduce the c[i] mod b[i], which might introduce more 
      # denominator, hence myreduce! doesn't work
      try
         d = one(base_ring(R))
         v = zero(base_ring(R))
         for k in 1:e[i]
            # (from prev: ck/d - v*z/d*b[i]/z) + c[k]*b[i]/z
            c[k] = a*c[k]*d - v*z
            g = gcd(d, content(c[k]))
            c[k] = divexact(c[k], g)
            d = divexact(d, g)
            ck, u, v = pseudoremx(c[k], b[i])
            # ck = u*c[k] + v*b[i]
            # c[k]/d = ck/u/d - v*z/u/d*b[i]/z
            d *= u
            # write answer terms
            t = F(n1)^-1*F(R(d))^-1
            g = content(ck)
            if iszero(g) || isunit(g)
               t *= F(ck)
            else
               t *= F(R(g))*F(divexact(ck, g))
            end
            for j in 1:n
               if j != i
                  t *= F(r[j])^(1-k-e[j])
               end
            end
            t *= F(b[i])^(-1+k-e[i])
            push!(answer.terms, t)
         end
      catch
         error("TODO: insert non-gcd domain code here")
      end
   end
   return answer
end


println("trying partial fractions")

G, (a, b, c) = PolynomialRing(ZZ, ["a", "b", "c"]);
Zx, x = PolynomialRing(G, "x"); a, b, c = map(Zx, [a,b,c]);
F = FactoredFractionField(Zx); (x, a, b, c) = map(F, [x, a, b, c]);
@show partial_fractions(x/((x-a)^4*(x-b)^3*(x-c)^2))
@show partial_fractions(x^6/((a*x-b)^4*(x-c)^3))