# Problem with SymPy. Symbolic calculation of Sturm chains (Sturm sequence) for polynomials

I was using `sympy.sturm()` to calculate the Sturm functions of a polynomial and came across an error for a 4th degree polynomial

``````using SymPy
begin
@vars x

# Polynomial coefficients
c4 = Sym[1; -4; 8; -8; 4] # Problem with Sympy sturm() ?
c3 = Sym[1; 2; 1; -2]

p4 = sympy.Poly(c4, x, domain="QQ")
p3 = sympy.Poly(c3, x, domain="QQ")

print("p3.sturm() == my_sturm(p3): ", p3.sturm() == my_sturm(p3),"\n")
print("p4.sturm() == my_sturm(p4): ", p4.sturm() == my_sturm(p4),"\n")
print("\n")
print("p4.sturm()  :", p4.sturm(),"\n\n")
print("my_sturm(p4):", my_sturm(p4),"\n")
end
``````

Output:

p3.sturm() == my_sturm(p3): true
p4.sturm() == my_sturm(p4): false
p4.sturm() :SymPy.Sym[Poly(x^2 - 2x + 2, x, domain=‘QQ’), Poly(2x - 2, x, domain=‘QQ’), Poly(-1, x, domain=‘QQ’)]
my_sturm(p4):SymPy.Sym[Poly(x^4 - 4x^3 + 8x^2 - 8x + 4, x, domain=‘QQ’), Poly(4x^3 - 12x^2 + 16x - 8, x, domain=‘QQ’), Poly(-x^2 + 2*x - 2, x, domain=‘QQ’)]

I believe my implementation

``````function my_sturm(p)
@vars x
n = p.degree()
pₛ = zeros(Sym,n+1)
pₛ = p
pₛ = p.diff()
i = 3
imax = 0

while sympy.degree(pₛ[i-1]) != 0 && i <= n+1
if pₛ[i-1].args != 0
# quotient and remainder
q, r = sympy.div(pₛ[i-2], pₛ[i-1]) #, domain="QQ")
if r.args != 0
pₛ[i] = -(pₛ[i-2] - pₛ[i-1]*q) # -r
imax = i
else
pₛ[i] = sympy.Poly(,x)
end
end
i += 1
end
return pₛ[1:imax]
end
``````

I figured it out myself. In file `rootisolation.py` of `Sympy`, only the square-free part of the polynomial is considered for the calculation of the Sturm chain. The function `dup_sqf_part(f, K)` is the first thing called.