Define the iterated quadratic polynomial:

\[ f_c^0(z) = z \\ f_c^{n+1}(z) = f_c(f_c^n(z))^2 + c \]

The Mandelbrot set is those \(c\) for which \(f_c^n(0)\) remains bounded for all \(n\). Misiurewicz points are dense in the boundary of the Mandelbrot set. They are strictly preperiodic, which means they satisfy this polynomial equation:

\[ f_c^{q+p}(0) = f_c^{q}(0) \\ p > 0 \\ q > 0\]

and moreover the period \(p\) and the preperiod \(q\) of a Misiurewicz point \( c \in M_{q,p} \) are the lowest values that make the equation true. For example, \(-2 \in M_{2,1}\) and \(i \in M_{2,2}\), which can be verified by iterating the polynomial (exercise: do that).

Misiurewicz points are algebraic integers (a subset of the algebraic numbers), which means they are the roots of a monic polynomial with integer coefficients. A monic polynomial is one with leading coefficient \(1\), for example \(c^2+c\). Factoring a monic polynomial gives monic polynomials as factors. Factoring over the complex numbers \(\mathbb{C}\) gives the \(M_{q,p}\) in linear factors, factoring over the integers \(\mathbb{Z}\) can give irreducible polynomials of degree greater than \(1\). For example, here's the equation for \(M_{2,2}\):

\[c^3\,\left(c+1\right)^2\,\left(c+2\right)\,\left(c^2+1\right)\]

Note that the repeated root \(0\) corresponds to a hyperbolic component of period \(1\) (the nucleus of the top level cardioid of the Mandelbrot set), and the repeated root \(-1\) corresponds to the period \(2\) circle to the left. And \(-2 \in M_{2,1}\), so the "real" equation we are interested in is the last term, \(c^2+1\), which is irreducible over the integers, but has complex roots \(\pm i\). There are two roots, so \(\left|M_{2,2}\right| = 2\).

So, a first **attempt** at enumerating Misiurewicz points works like this:

-- using numeric-prelude and MemoTrie from Hackage type P = MathObj.Polynomial.T Integer -- h with all factors g removed divideAll :: P -> P -> P divideAll h g | isZero h = h | isOne g = h | isZero g = error "/0" | otherwise = case h `divMod` g of (di, mo) | isZero mo -> di `divideAll` g | otherwise -> h -- h with all factors in the list removed divideAlls :: P -> [P] -> P divideAlls h [] = h divideAlls h (g:gs) = divideAlls (h `divideAll` g) gs -- the variable for the polynomials c :: P c = fromCoeffs [ 0, 1 ] -- the base quadratic polynomial f :: P -> P f z = z^2 + c -- the iterated quadratic polynomial fn :: Int -> P fn = memo fn_ where fn_ 0 = 0 ; fn_ n = f (fn (n - 1)) -- the raw M_{q,p} polynomial m_raw :: Int -> Int -> P m_raw = memo2 m_raw_ where m_raw_ q p = fn (q + p) - fn q -- the M_{q,p} polynomial with lower (pre)periods removed m :: Int -> Int -> P m = memo2 m_ where m_ q p = m_raw q p `divideAlls` [ mqp | q' <- [ 0 .. q ] , p' <- [ 1 .. p ] , q' + p' < q + p , p `mod` p' == 0 , let mqp = m q' p' , not (isZero mqp) ] -- |M_{q,p}| d :: Int -> Int -> Int d q p = case degree (m q p) of Just k -> k ; Nothing -> -1

This is using numeric-prelude and MemoTrie from Hackage, but with a reimplemented divMod for monic polynomials that doesn't try to divide by an Integer (which will always be \(1\) for monic polynomials). The core polynomial divMod from numeric-prelude needs a Field for division, and the integers don't form a field.

Tabulating this **attempt** at \(\left|M_{q,p}\right|\) (`d q p`

)
for various small \(q,p\) gives:

q | p | |||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | |

0 | 1 | 1 | 3 | 6 | 15 | 27 | 63 | 120 | 252 | 495 | 1023 | 2010 | 4095 | 8127 | 16365 | 32640 |

1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | |

2 | 1 | 2 | 6 | 12 | 30 | 54 | 126 | 240 | 504 | 990 | 2046 | 4020 | 8190 | 16254 | ||

3 | 3 | 3 | 12 | 24 | 60 | 108 | 252 | 480 | 1008 | 1980 | 4092 | 8040 | 16380 | |||

4 | 7 | 8 | 21 | 48 | 120 | 216 | 504 | 960 | 2016 | 3960 | 8184 | 16080 | ||||

5 | 15 | 15 | 48 | 90 | 240 | 432 | 1008 | 1920 | 4032 | 7920 | 16368 | |||||

6 | 31 | 32 | 96 | 192 | 465 | 864 | 2016 | 3840 | 8064 | 15840 | ||||||

7 | 63 | 63 | 189 | 384 | 960 | 1701 | 4032 | 7680 | 16128 | |||||||

8 | 127 | 128 | 384 | 768 | 1920 | 3456 | 8001 | 15360 | ||||||||

9 | 255 | 255 | 768 | 1530 | 3840 | 6912 | 16128 | |||||||||

10 | 511 | 512 | 1533 | 3072 | 7680 | 13824 | ||||||||||

11 | 1023 | 1023 | 3072 | 6144 | 15345 | |||||||||||

12 | 2047 | 2048 | 6144 | 12288 | ||||||||||||

13 | 4095 | 4095 | 12285 | |||||||||||||

14 | 8191 | 8192 | ||||||||||||||

15 | 16383 |

\(|M_{0,p}|\) is known to be A000740. \(|M_{2,p}|\) appears to be A038199. \(|M_{q,1}|\) appears to be A000225. \(|M_{q,2}|\) appears to be A166920.

**HOWEVER there is a fatal flaw**. The polynomials might
not be irreducible, which means that `divideAlls`

might not be
removing all of the lower (pre)period roots! A proper solution would be
to port the code to a computer algebra system that can factor polynomials
into irreducible polynomials. Or alternatively, mathematically prove that
the polynomials in question will always be irreducible (as far as I know
this is an open question, verified only for \(M_{0,p}\) up to \(p = 10\),
according to
Corollary 5.6 (Centers of Components as Algebraic Numbers)).

You can download my full Haskell code.

**UPDATE** I wrote some Sage code (Python-based) with an
improved algorithm (I think it's perfect now). The values all matched the
original table, and I extended it with further values and links to OEIS.
All the polynomials in question are irreducible, up to the \(p + q < 16\)
limit. No multiplicities greater than one were reported. Code:

@parallel(16) def core(q, p, allroots): mpq = 0 roots = set() R.<x> = ZZ[] w = 0*x for i in range(q): w = w^2 + x wq = w for i in range(p): w = w^2 + x wqp = w f = wqp - wq r = f.factor() for i in r: m = i[0] k = i[1] if not (m in allroots) and not (m in roots): roots.add(m) mpq += m.degree() if k > 1: print(("multiplicity > 1", k, "q", q, "p", p, "degree", m.degree())) return (q, p, mpq, roots) allroots = set() for n in range(16): print(n) res = sorted(list(core([(q, n - q, allroots) for q in range(n)]))) for r in res: t = r[1] q = t[0] p = t[1] mpq = t[2] roots = t[3] print((q, p, mpq, len(roots), [root.degree() for root in roots])) allroots |= roots

**UPDATE2** I bumped the table to \(q + p < 17\). I
ran into some OOM-kills, so I had to run it with less parallelism to get
it to finish.

**UPDATE3** I found a simple function that fits all the
data in the table, but I don't know if it is correct or will break for
larger values. Code (the function is called `f`

):

import Math.NumberTheory.ArithmeticFunctions (divisors, moebius, runMoebius) -- arithmoi import Data.Set (toList) -- containers mu :: Integer -> Integer mu = runMoebius . moebius mqps :: [[Integer]] mqps = [[1,1,3,6,15,27,63,120,252,495,1023,2010,4095,8127,16365,32640] ,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0] ,[1,2,6,12,30,54,126,240,504,990,2046,4020,8190,16254] ,[3,3,12,24,60,108,252,480,1008,1980,4092,8040,16380] ,[7,8,21,48,120,216,504,960,2016,3960,8184,16080] ,[15,15,48,90,240,432,1008,1920,4032,7920,16368] ,[31,32,96,192,465,864,2016,3840,8064,15840] ,[63,63,189,384,960,1701,4032,7680,16128] ,[127,128,384,768,1920,3456,8001,15360] ,[255,255,768,1530,3840,6912,16128] ,[511,512,1533,3072,7680,13824] ,[1023,1023,3072,6144,15345] ,[2047,2048,6144,12288] ,[4095,4095,12285] ,[8191,8192] ,[16383] ] m :: Integer -> Integer -> Integer m q p = mqps !! fromInteger q !! fromInteger (p - 1) f :: Integer -> Integer -> Integer f 0 p = sum [ mu (p `div` d) * 2 ^ (d - 1) | d <- toList (divisors p) ] f 1 _ = 0 f q 1 = 2 ^ (q - 1) - 1 f q p = (2 ^ (q - 1) - if q `mod` p == 1 then 1 else 0) * f 0 p check :: Bool check = and [ f q p == m q p | n <- [1 .. 16], p <- [1 .. n], let q = n - p ] main :: IO () main = print check

**UPDATE4** Progress! I found a paper with the answer:

Misiurewicz Points for Polynomial Maps and Transversality

Benjamin Hutz, Adam Towsley

Corollary 3.3.The number of \((m,n)\) Misiurewicz points for \(f_{d,c}\) is \[ M_{m,n} = \begin{cases} \sum_{k \mid n} \mu\left(n \over k \right) d^{k-1} & m = 0 \\ (d^m - d^{m-1} - d + 1) \sum_{k \mid n} \mu\left(n \over k \right) d^{k-1} & m \ne 0 \text{ and } n \mid (m - 1) \\ (d^m - d^{m-1}) \sum_{k \mid n} \mu\left(n \over k \right) d^{k-1} & \text{otherwise} \end{cases} \]

They have \(f_{d,c}(z) = z^d + c\), so this result is more general than the case \(d = 2\) I was researching in this post. The formula I came up with is the same, with minor notational differences.