Critical polynomial roots bigger than 2 In the question Why study critical polynomials?, it was said that each hyperbolic component of the Mandelbrot set contains a root for a critical polynomial. That is, the sequence $$\{0, 0^2 + c, c^2 + c, (c^2 + c)^2 + c, \dots, p_n, p_n^2 + c, \dots\}$$
the sequence of critical polynomials has roots corresponding to one of the blobs/mini-mandelbrots in the Mandelbrot set.
However, I also know that the mandelbrot set is contained inside $\{z \in \mathbb{C} : |z|<2\}$. 
So we should have that the roots of each polynomial is bounded inside this set.
In Sage I typed 
sage: R.<c> = PolynomialRing(ZZ)
sage: R
Univariate Polynomial Ring in c over Integer Ring
sage: p = 0
sage: P = []
sage: for i in xrange(10):
....:         p = p^2 + c
....:         P.append( points( (p.complex_roots()), hue=i/10, size = 20 ))
sage: sum(P)


I get roots with size greater than 2. What do these roots mean? Have I
  messed up somewhere in the code? Are the degrees so big that sage can't
  handle it (unlikely).

This is for the first 7 polynomials 
This is all 10 
 A: This is just numerical error - a standard thing that you've got
to get used to dealing with in this context. On my computer, the
errors become apparent around 7 or 8 iterations. If you expand the
polynomial at 8 iterations and check the largest coefficient, you
should get 2676118542978972739644, which is larger than $2^{63}$ or the
maximum machine integer. Mathematica and Sage can both deal with this
but it's not automatic. You'll have to trigger extra precision
somehow.
In Mathematica, you could do something like so:
f[c_][z_] = z^2 + c;
F[c][z_] = Nest[f[c], z, 8];
pts = c /. NSolve[F[c][0] == 0, c,
    WorkingPrecision -> 20];
ListPlot[{Re[#], Im[#]} & /@ pts,
 PlotRange -> All]

Notice the WorkingPrecision option. If you remove it, you'll find you
get roots outside the disk of radius two. It quickly gets worse as the
number of iterates increases.
In sage there is a way to do this too:
sage: R.<c> = PolynomialRing(ZZ)
sage: p = 0
sage: P = []
sage: for i in xrange(10):
....:         p = p^2 + c
....:         P.append( points((p.roots(ring=ComplexField(400), multiplicities=False)), hue=i/10, size = 20 ))
sage: sum(P)

Where we change ComplexField(400) to be how precise we want it. In this case, it would be 400 bits of precision.
Doing this gives us the great pictures like

And here's a live version with the iteration count reduced by one to get it to run a little faster.
