mathr / blog / #

Enumeration of Misiurewicz points

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:

qp
12345678910111213141516
0113615276312025249510232010409581271636532640
1000000000000000
212612305412624050499020464020819016254
333122460108252480100819804092804016380
478214812021650496020163960818416080
515154890240432100819204032792016368
631329619246586420163840806415840
7636318938496017014032768016128
812712838476819203456800115360
925525576815303840691216128
1051151215333072768013824
11102310233072614415345
1220472048614412288
134095409512285
1481918192
1516383

\(|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.