Number of real zeroes of iterated polynomial: x^3-2x+1

math2962644.sage

from sage.rings.polynomial.real_roots import *
x = polygen(ZZ)
def conv_p(p, n):
    if n == 0:
        return x
    else:
        return p(conv_p(p, n-1))

for n in range(1, 21): print add(map(lambda (a,b): b, real_roots(conv_p(x^3-2*x-1, n))))

math2962644b.sage

# This doesn't isolate the roots, but it actually runs slower
x = polygen(QQ)

# Implemented by Polynomial.compose_power if you have a recent enough version
def conv_p(p, n):
    if n == 0:
        return x
    else:
        return p(conv_p(p, n-1))

def count_real_roots(p):
    if p == 0:
        return infinity
    if p.degree() == 0:
        return 0
    count = 0
    for (factor, multiplicity) in p.squarefree_decomposition():
        sturm_seq = [factor, derivative(factor)]
        while (sturm_seq[-1].degree() > 0):
            sturm_seq.append(-(sturm_seq[-2].quo_rem(sturm_seq[-1])[1]))
        for i in range(len(sturm_seq) - 1):
            a = sturm_seq[i]
            b = sturm_seq[i+1]
            if (a.degree() + b.degree()) % 2 == 1:
                count += multiplicity * sign(a.leading_coefficient() * b.leading_coefficient())
    return count

for n in range(1, 21): print count_real_roots(conv_p(x^3-2*x-1, n))