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