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

Specific Domains Numerics
, ,
1

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

is giving the correct answer.

Does SymPy use a different definition of Sturm functions or am I missing something here?

2

I don’t know, but did verify that both give the same answer to the example problem here: Polynomials Manipulation Module Reference - SymPy 1.11 documentation

1 Like
3

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.

1 Like