Skip to content

Instantly share code, notes, and snippets.

@tthsqe12
Last active Feb 16, 2021
Embed
What would you like to do?
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))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment