News: Visit us on facebook
 Welcome, Guest. Please login or register. September 19, 2019, 09:26:08 AM 1 Hour 1 Day 1 Week 1 Month Forever Login with username, password and session length

## The All New FractalForums is now in Public Beta Testing! Visit FractalForums.org and check it out!

 Pages: 1 [2] 3 4 ... 6   Go Down
 Author Topic: The Mandelbrot Polynomial Roots Challenge  (Read 11112 times) Description: 0 Members and 1 Guest are viewing this topic.
knighty
Fractal Iambus

Posts: 819

 « Reply #15 on: June 01, 2013, 03:37:07 PM »

Thank you M Benesi.

This poster says:
Quote
The aim of this work was to compute the roots
of the polynomials pk(zeta), with the ultimate goal
of computing all 2^20-1 = 1,048,575 roots of
p21(zeta) shown in Figure 1. Parallel eigenvalue com-
putations were run on Sharcnet taking approxi-
mately 31 serial computing years.

The "was" and "were" suggest that the computations are over, thought it is not very clear to me.
 Logged
M Benesi
Fractal Schemer

Posts: 1075

 « Reply #16 on: June 01, 2013, 09:47:48 PM »

Yeah.  Maybe the project was completed and the web page (of P.W.) I quoted is not updated.

I wonder if the following would be faster than Piers Lawrence's method (not that I really understand the method that he's using):

a) locate tangent line from cardioid
b) calculate "circle" size and approximate center of each bud
......  thinking "circle" size should be approximated by a function related to p/q???
......  could absolute center be on tangent line, or at least calculated with p/q???
c) Newtonian root approximation (or other approximation method) started from approximate center

 « Last Edit: June 01, 2013, 09:49:41 PM by M Benesi » Logged

knighty
Fractal Iambus

Posts: 819

 « Reply #17 on: June 03, 2013, 02:42:36 PM »

Yeah.  Maybe the project was completed and the web page (of P.W.) I quoted is not updated.
Humm... yes! Piers Lawrence is definitly the winner.

I wonder if the following would be faster than Piers Lawrence's method (not that I really understand the method that he's using):

a) locate tangent line from cardioid
b) calculate "circle" size and approximate center of each bud
......  thinking "circle" size should be approximated by a function related to p/q???
......  could absolute center be on tangent line, or at least calculated with p/q???
c) Newtonian root approximation (or other approximation method) started from approximate center

This is valid only for (potentially all) the discs that are attached to the main cardioid. all other hyberbolic components, especially the minibrots are not attainable this way.

The script I've presented is already much faster than P. Lawrence's algorithm. The difference, at least, is that his method could be used for other sparse polynomials. It's an academic work not a hobbyist hack .

some pictures:
 EVAL0000.png (6.85 KB, 1152x845 - viewed 295 times.)  EVAL0001.png (7.12 KB, 1152x845 - viewed 504 times.)  EVAL0002.png (4.56 KB, 1152x845 - viewed 309 times.) « Last Edit: December 01, 2015, 09:26:36 PM by knighty, Reason: Wrong statment. » Logged
M Benesi
Fractal Schemer

Posts: 1075

 « Reply #18 on: June 03, 2013, 09:29:45 PM »

I'm assuming the 'millionth' root around the cardioid would require a bit more precision than the millionth root around the whole Mset.  Could get costly.

Have to admit, a few days ago I downloaded, but took no gander at, your root algorithm.  Will check it out after a much needed nap and see if I understand it.

I do like the root images.  Something nice about the stark B&W images, and the lack of sharp delineation due to the distribution of roots.
 Logged

knighty
Fractal Iambus

Posts: 819

 « Reply #19 on: May 25, 2014, 10:08:51 PM »

Sorry for reviving again this thread, but there is an update that is rotting on my hard drive for six months just because I wanted to optimize it further, implement the technique suggested by Benessi (thank you BTW) and do a multithreaded standalone applications (things that may never happen  ). Now it computes 1M roots in a little more than 3 mn on my computer instead of 20 mn for the previous version.
 MBPolyroots.zip (13.13 KB - downloaded 136 times.) Logged
plawrence
Forums Newbie

Posts: 4

 « Reply #20 on: November 11, 2015, 03:21:26 PM »

Hello,

I have no idea if anybody is still following this thread any more. I am Piers Lawrence, and I can add some of our new advancements since our first poster.

You can find some of the methodology here: IAP Poster 2013 or here: IAP Poster 2013 (alternative).

Let me say something about the progression of this work:

• Lagrange interpolation based approach:

We found the challenge a nice example for testing our Lagrange interpolation based rootfinder since the recurrence definition allows one to evaluate the polynomial in $O(k)$ cost rather than $O(degree=2^{(k-1)}-1)$. This approach worked up to about k=12, the problem being that the roots get too close together, making the computation of some auxiliary values (the barycentric weights) inaccurate.

• Mandelbrot matrices:

Since the polynomials have such a sparse recurrence, I started looking to see if we could form a similarly sparse matrix whose characteristic polynomial is exactly the $p_k(\zeta)$ whose roots we are after. This lead to the the construction you can see in both posters on the topic. Because the size of the matrices increase at a rate of $2^{(k-1)}-1$, it was clear that we needed some tricks. Essentially, this boils down to using Krylov based eigenvalue solvers (i.e. ARPACK) in shift and invert mode to find a small number of roots near some shift. It was particularly nice that the matrices have a very simple and sparse LU decomposition, which made it possible to get over 1 million roots. This solution, while the structure of the matrices is quite beautiful, is not really a good one. Basically since we ran the code on a parallel cluster (it is trivially parallelizable over the shifts chosen) taking around 31 serial years of computing resources, I think it took a little under a week on the cluster.

• Homotopy based approach

We started to investigate other ways of computing roots of polynomials, and we discovered that a homotopy based approach is the way of the future! The basic idea is to rewrite the polynomial with a new parameter, $\epsilon$ say, such that as we vary $\epsilon$  from zero to one, the roots trace a path from the (double) roots of $p_{k}(\zeta)$ (plus one extra point at zero) to those of $p_{k+1}(\zeta)$.

From this new equation, we can form a differential equation describing the path of the roots (see the poster). The homotopy continuation method is not without other problems (namely singularities in the path), some of which can be alleviated by integrating in the complex $\epsilon$-domain, but what was surprising was the speed: to compute all roots from $k=2$ to $k=14$ takes now only around 3 seconds on a desktop! Something that would take on the order of hours via the Mandelbrot matrices!

The complication is the singularities, these only get more dense as we increase $k$, since we have to worry about the singularities due to the higher derivatives of $p_k(\zeta)$.

I attach the code to compute the roots via the homotopy continuation method, it seems clear that one could get over 1 billion roots ($k\approx 31$) amd perhaps even further (although I guess higher precision will be needed). There is a shell script run_mandel.sh that will compute all the roots up to $k=14$, the rest is written in C, although not well documented (sorry for that), it also will require the installation of the Gnu scientific library.

Anyway, I hope you find this interesting, and perhaps it will be useful for computing other fractals, too, who really knows.

Kind regards,

Dr. Piers Lawrence
 mandel.zip (10.13 KB - downloaded 98 times.) « Last Edit: November 11, 2015, 04:19:31 PM by plawrence, Reason: Missing license file » Logged
Fractal Lover

Posts: 221

 « Reply #21 on: November 11, 2015, 03:54:03 PM »

Thx. It's interesting topic.
Is it possible ( do you plan it ) to make a procedure ( or library ) with free licence ( code ) which use your method ?
Is it posible to see the results of your computations ? ( like list )

TIA

https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers
http://fraktal.republika.pl/mset_centers.html
http://fraktal.republika.pl/eigensolve.html
 « Last Edit: November 11, 2015, 04:10:38 PM by Adam Majewski, Reason: new important link » Logged
plawrence
Forums Newbie

Posts: 4

 « Reply #22 on: November 11, 2015, 04:26:24 PM »

Sorry I prematurely pressed post, but now everything seems to be there.

The code is attached to the post, with a permissive license.

The results of the computation can be seen on the poster, basically it outputs the paths that have been computed. You should have all the points in the data directory within a few seconds, and there is a small gnuplot script to plot the paths and roots.

Kind regards,

Piers
 Logged
Fractal Lover

Posts: 221

 « Reply #23 on: November 11, 2015, 05:48:30 PM »

Hi,

Thx for the code.

I have compiled it

make
cc    -c -o mandelhom.o mandelhom.c
cc    -c -o mandel.o mandel.c
cc    -c -o mandelutil.o mandelutil.c
cc  -o mandelhom mandelhom.o mandel.o mandelutil.o -lgsl -lgslcblas -lm
cc    -c -o mandeleuler.o mandeleuler.c
cc  -o mandeleuler mandeleuler.o mandel.o -lgsl -lgslcblas -lm
cc    -c -o mandelconvert.o mandelconvert.c
cc  -o mandelconvert mandelconvert.o

and run

./run_mandel.sh
Read 1 initial values from data/fin2.bin
Read 3 initial values from data/fin3.bin
Read 7 initial values from data/fin4.bin
Read 15 initial values from data/fin5.bin
Read 31 initial values from data/fin6.bin
Read 63 initial values from data/fin7.bin
Read 127 initial values from data/fin8.bin
Read 255 initial values from data/fin9.bin
Read 511 initial values from data/fin10.bin
Read 1023 initial values from data/fin11.bin
Read 2047 initial values from data/fin12.bin
Read 4095 initial values from data/fin13.bin
singular minus: -1.6223428021976933e+00 -2.3975009778627543e-02 (-nan -nan)
Read 8191 initial values from data/fin14.bin
singular plus: -nan -nan (-nan -nan)
singular minus: -nan -nan (-nan -nan)

then I wanted to make image :

cd data
a@zalman:~/c/mandel/centers/plawrence/mandelhom/data$ls fin10.bin fin12.bin fin14.bin fin2.bin fin3.bin fin5.bin fin7.bin fin9.bin path11.dat path13.dat path15.dat path4.dat path6.dat path8.dat plothom.pg fin11.bin fin13.bin fin15.bin fin2.dat fin4.bin fin6.bin fin8.bin path10.dat path12.dat path14.dat path3.dat path5.dat path7.dat path9.dat zero.dat a@zalman:~/c/mandel/centers/plawrence/mandelhom/data$ ./plothom.pg
line 0: warning: Skipping unreadable file "./fin.bin"
line 0: warning: Skipping unreadable file "./fin1.bin"
line 0: warning: Skipping unreadable file "./path1.dat"
line 0: warning: Skipping unreadable file "./cont1.dat"
Warning: empty x range [0:0], adjusting to [-1:1]
Warning: empty y range [0:0], adjusting to [-1:1]

Help is welcome .

I suppose centers are in a path files .
What is in these files.
How cen I extract only centers ( complex number) from path file

TIA

 Logged
plawrence
Forums Newbie

Posts: 4

 « Reply #24 on: November 11, 2015, 07:18:28 PM »

Sorry that was not so clear. The plotting script takes an argument being the level that you want to plot.

The data is written in binary format, with the real and imaginary parts.

There is a conversion program mandelconvert which allows you to print out the ascii real and imaginary parts, one on each line. It can be used like this:

./mandelconvert -mb < data/fin3.bin

The nan values indicate that a singularity has been encountered in the path. And the plus and minus indicate if it is on the positive or negative square root used in the formula to move off the starting value. (This is because you start at something like a saddle point).

Regards,

Piers

 Logged
Fractal Lover

Posts: 221

 « Reply #25 on: November 11, 2015, 08:46:05 PM »

OK

I have :
./mandelconvert -mb < data/fin14.bin> 14.txt

removed nan lines

Now file 14.txt has :

wc -l 14.txt
8191 14.txt

8191 lines.

If I'm not wrong for period 14 there is a 8127 centers.

Do you remove centers of lower periods ?
http://fraktal.republika.pl/mset_components.html

 Logged
plawrence
Forums Newbie

Posts: 4

 « Reply #26 on: November 12, 2015, 01:05:27 PM »

There is no removal of roots of lower periods.

For period 14 there should be $2^{(14-1)}-1=8191$ roots.

I think what you are getting at is if we count the period 7 roots also, indeed these are included in the 8191 roots.

Piers
 Logged
JArtherton
Forums Newbie

Posts: 2

 « Reply #27 on: November 23, 2015, 10:45:35 PM »

Hello, I have been working on finding roots for the Mandelbrot set for some time and have developed a search program that has now found 522,921 roots to period 19.
These are found using 320 bit math for higher precision.
Data available at www.scideveloper.com/Fractals/fractals.html

 Logged
knighty
Fractal Iambus

Posts: 819

 « Reply #28 on: November 24, 2015, 04:11:12 PM »

Hello all!
@ plawrence: Congratulations for winning the contest!
Thank you very much for the informations about your work and the source code.

@ JArtherton: Could you give us more information about the method used.
BTW! Assuming that the smallest distance between two roots occure near -2, Dr. Lawrence's formula (Largest root of P_k) can be used to get an idea about the needed accuracy. For one, I used 2^(-(2*k+2)). So using doubles, maximal k would be 26. With 80-bit extended precision max k would be 30 or 31. The roots can then be refined to the needed accuracy (using Newton method for example).
 Logged
JArtherton
Forums Newbie

Posts: 2

 « Reply #29 on: November 24, 2015, 06:44:15 PM »

My newest update has exceeded 1 million roots.
I use a quadtree subdivision of the grid to increase the search resolution in areas with higher root density. This adapts to very fine resolution in the tips of the antennae.

James
 Logged
 Pages: 1 [2] 3 4 ... 6   Go Down