# mathr

## Enumeration of Misiurewicz points

$f_c^0(z) = 0 \\ 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 ]

f :: P -> P
f z = z^2 + c

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)).

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):
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
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.