Okay, technical details.
There are multiple numerical methods to find complex roots of polynomials. I'm aware about Durand-Kerner method (https://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method) and Dandelin-Graeffe method (https://en.wikipedia.org/wiki/Graeffe's_method), where "aware" means "I know that they exist". I did not implement root finding myself; instead I've used GNU scientific library (http://www.gnu.org/software/gsl/) and built-in function of it called gsl_poly_complex_solve (http://www.gnu.org/software/gsl/manual/html_node/General-Polynomial-Equations.html#General-Polynomial-Equations). The documentation claims that they use "balanced-QR reduction of the companion matrix" and I have no idea what is it and was not aware about it existence before. Since I wrote my code in OCaml, I've used OCaml bindings for GSL (from https://github.com/mmottl/gsl-ocaml).
After the roots are calculated, I transform them to pixel coordinates, making sure that the whole thing scaled to the given image size (I must admit that I crop out some purely real roots, nothing interesting there). Then, according to pixel coordinates, I write 'FFFFFF' (i.e. "white") bytes to the corresponding positions in an array. This array backs up Cairo (http://cairographics.org/) surface of the 20000x20000 size, and Cairo PNG module is used to write the result to a file. Again, I needed OCaml bindings (https://github.com/Chris00/ocaml-cairo).
This approach limits the resolution of the image; for example, Cairo fails to create surface of 40000x40000 size. See below for future plans.
Most of the time is taken by calculating the roots and scaling them to fill the image of the given size; contribution of drawing and writing png file is noticeable but insignificant. On Macbook Pro with 3GHz Intel Core i7 the whole process takes about 40 minutes; drawing the picture and writing it to file takes about 2-3 minutes.
- I preferred convenience in writing to optimization, so there's a lot of work for garbage collector when it kicks in -- the pauses are quite noticeable.
- The whole thing is written in a pretty dumb way, it runs in one thread only, so no multi-CPU multi-core benefits.
- Everything is kept in memory until image rendered. As as result I cannot do degree 25 or more, as I'm running out of memory (I have 16 Gb of RAM). Memory usage definitely can be optimized. See "convenience over optimization" point above.
I have started to rewrite everything in a parallel fashion. I also plan to do rendering in chunks (like "top left quarter image, top right quarter image, bottom left quarter image, bottom right quarter image") -- this would allow me to have much higher resolution. When this is done, I'll rent an hour from some Amazon high-memory high-CPU server and will do much better resolution and higher degree.
I have expected that drawing the picture of degree 23 on top of the picture of degree 24 would be interesting. Not so. Here is the degree 23 (in green) on top of degree 24 (in cyan):
Green is basically distributed everywhere -- more in the middle of the ring, but also in all interesting fractal-like parts.
Here is degree 22 (in orange) on top of degree 23 (in green) on top of degree 24 (in cyan):
Same story. Adding more degrees creates the mess of colored pixels.
But what's clear (from bigger pictures which I do not share) is that my resolution is too rough to see how exactly degree 23 roots overlap degree 24 roots in interesting areas. Stay tuned :)