Loading [MathJax]/jax/output/HTML-CSS/jax.js

mathr / blog / #

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 cMq,p are the lowest values that make the equation true. For example, 2M2,1 and iM2,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 2M2,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:

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

|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={knμ(nk)dk1m=0(dmdm1d+1)knμ(nk)dk1m0 and n(m1)(dmdm1)knμ(nk)dk1otherwise

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.