using SymPy

# Override the Base.show definition of SymPy.jl:
# https://github.com/JuliaPy/SymPy.jl/blob/29c5bfd1d10ac53014fa7fef468bc8deccadc2fc/src/types.jl#L87-L105

@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::SymbolicObject)
    print(io, as_markdown("\\displaystyle " * sympy.latex(x, mode="plain", fold_short_frac=false)))
end
@eval SymPy function Base.show(io::IO, ::MIME"text/latex", x::AbstractArray{Sym})
    function toeqnarray(x::Vector{Sym})
        a = join(["\\displaystyle " * sympy.latex(x[i]) for i in 1:length(x)], "\\\\")
        """\\left[ \\begin{array}{r}$a\\end{array} \\right]"""
    end
    function toeqnarray(x::AbstractArray{Sym,2})
        sz = size(x)
        a = join([join("\\displaystyle " .* map(sympy.latex, x[i,:]), "&") for i in 1:sz[1]], "\\\\")
        "\\left[ \\begin{array}{" * repeat("r",sz[2]) * "}" * a * "\\end{array}\\right]"
    end
    print(io, as_markdown(toeqnarray(x)))
end

all_coeffs(f, x) = sympy.poly(f, x).all_coeffs()

function all_divisors(x)
    @assert x != 0
    A = typeof(x)[]
    for i in 1:abs(x)
        if mod(x, i) == 0
            push!(A, i, -i)
        end
    end
    A
end

function factor_for_all_divisors(f, var, constant_term)
    A = all_divisors(constant_term)
    stack([k, f(var=>k).factor()] for k in A; dims=1)
end

@vars a b c d e p q r s t u v x

f = x^5 + p*x^4 + q*x^3 + r*x^2 + s*x + t

g = x^2 + a*x + b

h = x^3 + c*x^2 + d*x + e

quotient, remainder = sympy.div(f, g, x)
sympy.poly.([quotient, remainder], x)

_, C, D, E = all_coeffs(quotient, x)

eqrem0 = remainder.coeff(x, 0)/b |> expand

eqrem1 = remainder.coeff(x, 1)

sympy.div(f, x^3 + a*x^2 + b*x + c, x) |> collect

sympy.div(x^6 + p*x^5 + q*x^4 + r*x^3 + s*x^2 + t*x + u, g, x) |> collect

Q63, R63 = sympy.div(x^6 + p*x^5 + q*x^4 + r*x^3 + s*x^2 + t*x + u, x^3 + a*x^2 + b*x + c, x) |> collect

C63 = all_coeffs(R63, x)

B63 = sympy.solve(C63[end], b)[end]

@show eqrem63_1 = -(C63[1](b=>B63)(u=>c*v) * (-2*a + p)^2).expand().factor();

all_coeffs(eqrem63_1, a)

Problems = [
    x^5 + x + 1
    x^5 - 5x + 3
    x^5 + x + 6
    x^5 - x + 15
    x^5 - 9x - 27
    x^5 - 9x + 27
    x^5 - 5x + 28
    x^5 + 16x + 32
    x^5 + 16x - 32
    x^5 - 5x^2 + 2
    x^5 - x^2 + 4
    x^5 + 2x^2 + 4
    x^5 - 3x^2 + 18
    x^5 + x^2 + 18
]

stack([f, f.factor()] for f in Problems; dims=1)

GG = x^2 + 3x - 2

HH = (x+1)^3-2 |> expand

FF = GG*HH

FF = expand(FF)

factor(FF)

_, P, Q, R, S, T = all_coeffs(FF, x)

Eqrem0 = eqrem0(p=>P, q=>Q, r=>R, s=>S, t=>T)

Eqrem1 = eqrem1(p=>P, q=>Q, r=>R, s=>S, t=>T)

factor_for_all_divisors(Eqrem0, b, N(T))

BB, AA = -2, 3

C, D, E

CC, DD, EE = C(a=>AA, b=>BB, p=>P), D(a=>AA, b=>BB, p=>P, q=>Q), E(a=>AA, b=>BB, p=>P, q=>Q, r=>R, t=>T)

Eqrem1(a=>AA, b=>BB)

sol = (x^2 + AA*x + BB)*(x^3 + CC*x^2 + DD*x + EE)

sol - FF |> simplify

FF = x^5 - x^4 - 1

factor(FF)

_, P, Q, R, S, T = all_coeffs(FF, x)

Eqrem0 = eqrem0(p=>P, q=>Q, r=>R, s=>S, t=>T)

Eqrem1 = eqrem1(p=>P, q=>Q, r=>R, s=>S, t=>T)

factor_for_all_divisors(Eqrem0, b, N(T))

BB, AA = 1, -1
Eqrem1(a=>AA, b=>BB)

BB, AA = -1, -1
Eqrem1(a=>AA, b=>BB)

AA, BB = -1, 1
CC, DD, EE = C(a=>AA, b=>BB, p=>P), D(a=>AA, b=>BB, p=>P, q=>Q), E(a=>AA, b=>BB, p=>P, q=>Q, r=>R, t=>T)

sol = (x^2 + AA*x + BB)*(x^3 + CC*x^2 + DD*x + EE)

sol - FF |> simplify

FF = x^5 + 16x + 32

factor(FF)

_, P, Q, R, S, T = all_coeffs(FF, x)

Eqrem0 = eqrem0(p=>P, q=>Q, r=>R, s=>S, t=>T)

Eqrem1 = eqrem1(p=>P, q=>Q, r=>R, s=>S, t=>T)

factor_for_all_divisors(Eqrem0, b, N(T))

BB, AA = 4, 2
Eqrem1(a=>AA, b=>BB)

BB, AA = -2, 2
Eqrem1(a=>AA, b=>BB)

AA, BB = 2, 4
CC, DD, EE = C(a=>AA, b=>BB, p=>P), D(a=>AA, b=>BB, p=>P, q=>Q), E(a=>AA, b=>BB, p=>P, q=>Q, r=>R, t=>T)

sol = (x^2 + AA*x + BB)*(x^3 + CC*x^2 + DD*x + EE)

sol - FF |> simplify

FF = x^5 + 3x^4 - 2x^3 - 2x^2 - 6x + 4

factor(FF)

factor(FF)

_, P, Q, R, S, T = all_coeffs(FF, x)

eqrem0

Eqrem0 = eqrem0(p=>P, q=>Q, r=>R, s=>S, t=>T)

eqrem1

Eqrem1 = eqrem1(p=>P, q=>Q, r=>R, s=>S, t=>T)

factor_for_all_divisors(Eqrem0, b, N(T))

BB, AA = -2, 3
Eqrem1(a=>AA, b=>BB)

AA, BB = 3, -2
CC, DD, EE = C(a=>AA, b=>BB, p=>P), D(a=>AA, b=>BB, p=>P, q=>Q), E(a=>AA, b=>BB, p=>P, q=>Q, r=>R, t=>T)

sol = (x^2 + AA*x + BB)*(x^3 + CC*x^2 + DD*x + EE)

sol - FF |> simplify