Enumeration of Misiurewicz points
Define the iterated quadratic polynomial:
f0c(z)=zfn+1c(z)=fc(fnc(z))2+c
The Mandelbrot set is those c for which fnc(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:
fq+pc(0)=fqc(0)p>0q>0
and moreover the period p and the preperiod q of a Misiurewicz point c∈Mq,p are the lowest values that make the equation true. For example, −2∈M2,1 and i∈M2,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 c2+c. Factoring a monic polynomial gives monic polynomials as factors. Factoring over the complex numbers C gives the Mq,p in linear factors, factoring over the integers Z can give irreducible polynomials of degree greater than 1. For example, here's the equation for M2,2:
c3(c+1)2(c+2)(c2+1)
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∈M2,1, so the "real" equation we are interested in is the last term, c2+1, which is irreducible over the integers, but has complex roots ±i. There are two roots, so |M2,2|=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 |Mq,p| (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 |
|M0,p| is known to be A000740. |M2,p| appears to be A038199. |Mq,1| appears to be A000225. |Mq,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 M0,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 fd,c is Mm,n={∑k∣nμ(nk)dk−1m=0(dm−dm−1−d+1)∑k∣nμ(nk)dk−1m≠0 and n∣(m−1)(dm−dm−1)∑k∣nμ(nk)dk−1otherwise
They have fd,c(z)=zd+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.