In the old forum, there was an exciting challenge 10 years ago (when I was only computing Lyapunov diagrams, so I missed it) about one million roots in the Mandelbrot set (http://www.fractalforums.com/theory/the-mandelbrot-polynomial-roots-challenge/) and a few years back there was a similar article by Schleicher et al. Newton's method in practice: Finding all roots of polynomials of degree one million efficiently, 2017.

I have yet to read those articles completely and in depth, but I thought, I try it with my interval arithmetics based subdivision implementation to see how far/fast I can get there - and what types of optimizations need to be implemented to reach a million (if at all possible).

Below are the first levels to detect the hyperbolic centers of period 7 and 1 using the z-constant term of the 7-fold composition f[7] of f(z)=z^2+c, i.e. g[7] for g[n]=g[n-1]^2+c, g[0]=0 with degree 64 (so quite a long distance to a million to go ).

Upper row shows the development of some levels. The gray region shrinks in an interesting fashion from outwards in streaks and at early levels follow roughly the Mset shape. It looks a bit like a crowded game of life.

At deepest level L16 so far it found 40 definite root regions (lower left, red rectangles around black pixels; probably 40 roots then from the image, but see details' last item). Some regions are yet to be judged (lower middle, gray circumferenced gray regions) which interestingly cover a lot of the components sitting on the real axis (including the main cardioid).

Computationally the number of boxes per level increases exponentially as expected (lower right) - about 2^1.8 fold at start, so not quite the maximal 2^2-fold increase due to subdivision into 4 squares - but at some point the complexity drops rapidly and box numbers go down, so interestingly no plateau there. It seems to be a "climb-the-hill-and-then-roll-down-fast" type of complexity. There are only 248 boxes currently left in total at end of L16, so the present 64 root regions might be in reach soon.

Has anyone attempted this (with or w/o IA) recently? Probably one can get to more than a million roots with nowadays computers?

Technical details

• the image shows the complex 2-square. It is for visualization and not intended to provide a reliable filter (but region coordinates are available).
• regions fully in the positive imaginary half are not analyzed due to symmetry.
• I use my custom-made fixed-point complex interval type with arbitrary (at compile time) precision of about 270 decimal digits for now with simulated outward rounding when needed.
• subdivision uses the interval Newton operator as described in forum link.
• the polynomial g[n] is implemented in iterative form, the inverse Jacobian is computed using elementary relations between complex and partial derivatives.
• computation was performed on-disk in a single-thread manner and took 4 hrs in total.
• as a plausibility check, the midpoint (double precision calculated) of a definite root region was analyzed by the TSA cell mapping. In 2*13/40, a cycle of length 7 could be verified at level 18 (no other cycles were found).
• root regions have not yet been checked whether they are pairwise disjoint (they could touch if a root lies on a rationally complex border of two adjacent regions).

Linkback: https://fractalforums.org/programming/11/the-mandelbrot-set-root-challenge-10-years-later/3952/

"A thousand roots"

All the 1024 roots of the hyperbolic centers of periods 11 and 1 have been reliably enclosed. The automatically created output of the verification reads:

--- BEGIN root status -----------

  559 DEFINITE root regions in CURRENT level L62
    of which 94 touch/straddle the real axis
    and 465 purely in the negative imag complex part

  CONCLUSION: The expected number of 1024 (=2*465 + 94) roots have been found (including symmetry).

  ALL root regions are pairwise disjoint.

  Each root in its root region is at most 10^-95 away from any point therein

--- END -----------

Attached is a text-file containing all the root regions. Black pixels in the image denote the root regions.

Algorithmic details

• a quick interval-evaluated escape test was performed to discard boxes before computing the expensive IA operator (identical to the exteriorIA method in forum link). Escape max it was only slightly increased from 20-30 as it will only be used in lower levels.
• precision used was about 410 decimal digits (although not tested, how much is sufficient to find the roots).
• computation took about 3 days in total with ongoing optimizations.
• boxes were analyzed until only regions touching the real axis were left undecided. Then one wobbly expansion event was performed.
• overlappung regions after expansion were fused - they were removed from the analysis queue and a rectangular superset covering both was newly introduced.
• subdivision for 2 more levels sufficed to get to the above mentioned accuracy.
• accuracy was conservatively determined by calculating the maximum of the sum of the width and height (Manhattan distance) over all root containing regions and then rounded to the next larger pure negative power of 10. Individual boxes can provide a higher accuracy.
• the process was terminated when the expected number of disjoint definitely root-containing regions was reached - factoring in symmetry of boxes that lie properly in the negative imaginary complex hemisphere.

@pauldelbrot: Interesting classification result. I would have expected more bulbs (as there are - in a sense - infinitely times more bulbs than cardioids), but hyperbolic centers then seem to be distributed with some regularity. I'm currently computing centers for periods 14/7/2/1 (the hill of increasing numbers of box regions in the queue has almost been reached) - what ratio between bulbs and cardioids occurs there?


-----------------
"Heuristics"

Analyzing the data for last post's period 11, it turned out that the most powerful method to discard a box as not containing a root (and hence emptying the queue) was: simply checking, whether the polynomial evaluation overlaps with the origin (light blue bars in image below). Two more methods are currently employed to discard a box: the interval escape test (red), which has an effect in the early stages, but later diminishes because I only increase the max it slightly to save time there.

Interestingly, the 3rd method, the Newton IA subdivision operator (green) - which works both as an inclusion and exclusion predicate - only shows an effect once the box number has reached vicinity of its peak - and then the "discarding power" is quite slim.

For now, I use this as a heuristics: Omitting the calculation of the IAoperator until I see a clear drop in box numbers. This is the type of optimizations I like most, as they do not cut off from the search space and do not lose a root, and hopefully help finding them faster.

Technical details

• this has only been analyzed with one period, so results might not be applicable in general.
• besides being a powerful predicate, the Newton operator also serves to shrink the region to subdivide: If the box B contains a root, so does the result of op(B). Hence it suffices to subdivide the intersection of both. This could then lead to faster convergence and an overall lower number of boxes (not analyzed).
• calculations were performed with ongoing optimizations. Peak numbers therefore can now be lower with all optimizations on in every level.