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

**URL:** https://discourse.julialang.org/t/problem-with-sympy-symbolic-calculation-of-sturm-chains-sturm-sequence-for-polynomials/84921
**Category:** Numerics
**Tags:** question, polynomials, sympy
**Created:** [July 28, 2022, 3:16pm UTC](https://discourse.julialang.org/t/problem-with-sympy-symbolic-calculation-of-sturm-chains-sturm-sequence-for-polynomials/84921 "2022-07-28T15:16:13Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)
#### Post date: [July 28, 2022, 3:16pm UTC](https://discourse.julialang.org/t/problem-with-sympy-symbolic-calculation-of-sturm-chains-sturm-sequence-for-polynomials/84921/1 "2022-07-28T15:16:13Z")

</div>

I was using `sympy.sturm()` to calculate the [Sturm functions](https://mathworld.wolfram.com/SturmFunction.html) of a polynomial and came across an error for a 4th degree polynomial

```julia
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 - 2_x + 2, x, domain=‘QQ’), Poly(2_x - 2, x, domain=‘QQ’), Poly(-1, x, domain=‘QQ’)]  
> my\_sturm(p4):SymPy.Sym[Poly(x^4 - 4_x^3 + 8_x^2 - 8_x + 4, x, domain=‘QQ’), Poly(4_x^3 - 12_x^2 + 16_x - 8, x, domain=‘QQ’), Poly(-x^2 + 2\*x - 2, x, domain=‘QQ’)]

I believe my implementation

```julia
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?

---

<div class="post-metadata">

### Author: ![j\_verzani](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/j_verzani/32/8551_2.png) [@j\_verzani](https://discourse.julialang.org/u/j_verzani)
#### Post date: [July 28, 2022, 4:05pm UTC](https://discourse.julialang.org/t/problem-with-sympy-symbolic-calculation-of-sturm-chains-sturm-sequence-for-polynomials/84921/2 "2022-07-28T16:05:13Z")

</div>

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](https://docs.sympy.org/latest/modules/polys/reference.html)

---

<div class="post-metadata">

### Author: ![mendeli](https://avatars.discourse-cdn.com/v4/letter/m/6de8d8/32.png) [@mendeli](https://discourse.julialang.org/u/mendeli)
#### Post date: [August 3, 2022, 4:44pm UTC](https://discourse.julialang.org/t/problem-with-sympy-symbolic-calculation-of-sturm-chains-sturm-sequence-for-polynomials/84921/3 "2022-08-03T16:44:54Z")

</div>

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.