Title: The Mandelbrot Polynomial Roots Challenge
Post by: bugman on January 02, 2010, 06:15:03 PM
The quadratic Mandelbrot set is defined by iterating the following function:
f(z) = z² + c
For example, the 2nd level Mandelbrot polynomial is given by:
F(z) = f(f(z)) = (z² + c)² + c
Now, suppose we arbitrarily choose a "seed" value of z = 0 and solve for the roots of c such that:
F(0) = 0
Here is a plot of the roots of the 12th level Mandelbrot polynomial. Altogether, there are 2048 roots. This may not look like much, but it was fairly difficult to calculate. I found these roots using the Durand-Kerner method. Total runtime was 4 hours, 45 minutes, although I think it could be made much faster:
http://en.wikipedia.org/wiki/Durand-Kerner_method
For those of you computer math wizards out there, I challenge you to find more roots than this. The first one to find 1 million roots wins! 8)
Of course, this is still a far cry from the inverse Julia set pruning methods which can probe hundreds or even thousands of levels deep. Unfortunately, I don't know of any similar way to prune Mandelbrot set roots.
f(z) = z² + c
For example, the 2nd level Mandelbrot polynomial is given by:
F(z) = f(f(z)) = (z² + c)² + c
Now, suppose we arbitrarily choose a "seed" value of z = 0 and solve for the roots of c such that:
F(0) = 0
Here is a plot of the roots of the 12th level Mandelbrot polynomial. Altogether, there are 2048 roots. This may not look like much, but it was fairly difficult to calculate. I found these roots using the Durand-Kerner method. Total runtime was 4 hours, 45 minutes, although I think it could be made much faster:
http://en.wikipedia.org/wiki/Durand-Kerner_method
For those of you computer math wizards out there, I challenge you to find more roots than this. The first one to find 1 million roots wins! 8)
Of course, this is still a far cry from the inverse Julia set pruning methods which can probe hundreds or even thousands of levels deep. Unfortunately, I don't know of any similar way to prune Mandelbrot set roots.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: bugman on January 02, 2010, 06:15:36 PM
The way I found the roots was by first solving the roots of a lower level Mandelbrot polynomial, and then making an array containing 2 duplicates of each of those roots and plugging this new array back in as my initial guess before trying to solve the higher degree polynomial (for stability, I had to slightly perturb the duplicated roots so there would be no duplicates in my intial guess). The convergence time varied depending on how much I perturbed the duplicated roots, but in the better scenarios, I found that this accelerated my overall convergence time by a factor of 40. I stopped calculating when the solution stopped changing by any more than 10^-7.
I understand Steven Fortune has a program called Eigensolve that can find the roots of high order polynomials very quickly:
http://ect.bell-labs.com/who/sjf/eigenvalueiterexamples.html
He told me his program wasn't too much faster, because convergence time increases by roughly 10x for each level. However, I still suspect that in the specific case of solving Mandelbrot polynomials, it might be 30 to 120 times faster than my program if the same approach is used.
I understand Steven Fortune has a program called Eigensolve that can find the roots of high order polynomials very quickly:
http://ect.bell-labs.com/who/sjf/eigenvalueiterexamples.html
He told me his program wasn't too much faster, because convergence time increases by roughly 10x for each level. However, I still suspect that in the specific case of solving Mandelbrot polynomials, it might be 30 to 120 times faster than my program if the same approach is used.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: bugman on January 02, 2010, 06:24:32 PM
What about the triplex roots of the Mandelbulb polynomial? Well, I tried that but convergence was way too slow. I couldn't even solve the roots of the 2nd level Mandelbulb polynomial.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: BradC on January 06, 2010, 04:16:59 PM
I tried to find the roots of the 2nd level Mandelbulb (degree 8) function by numerically finding roots starting from a whole mess of random seed points and then removing duplicate results. I got 64 unique answers, and they appear to be situated in a configuration that looks kinda like the Mandelbulb (see attached image below).
The coordinates of these 64 roots to double precision (rounded from higher precision) are:
The coordinates of these 64 roots to double precision (rounded from higher precision) are:
Code:
{{-1., 0., 0.},
{-0.8117449009293668, -0.3909157412340149, -0.4338837391175581},
{-0.8117449009293668, -0.3909157412340149, 0.4338837391175581},
{-0.8117449009293668, 0.3909157412340149, -0.4338837391175581},
{-0.8117449009293668, 0.3909157412340149, 0.4338837391175581},
{-0.766044443118978, 0., -0.6427876096865394},
{-0.766044443118978, 0., 0.6427876096865394},
{-0.6234898018587335, -0.7818314824680298, 0.},
{-0.6234898018587335, 0.7818314824680298, 0.},
{-0.4776208980552355, -0.5989176626001069, -0.6427876096865394},
{-0.4776208980552355, -0.5989176626001069, 0.6427876096865394},
{-0.4776208980552355, 0.5989176626001069, -0.6427876096865394},
{-0.4776208980552355, 0.5989176626001069, 0.6427876096865394},
{-0.20048443395120957, -0.8783796973249267, -0.4338837391175581},
{-0.20048443395120957, -0.8783796973249267, 0.4338837391175581},
{-0.20048443395120957, -0.0965482148568969, -0.9749279121818236},
{-0.20048443395120957, -0.0965482148568969, 0.9749279121818236},
{-0.20048443395120957, 0.0965482148568969, -0.9749279121818236},
{-0.20048443395120957, 0.0965482148568969, 0.9749279121818236},
{-0.20048443395120957, 0.8783796973249267, -0.4338837391175581},
{-0.20048443395120957, 0.8783796973249267, 0.4338837391175581},
{-0.17364817766693036, 0., -0.984807753012208},
{-0.17364817766693036, 0., 0.984807753012208},
{-0.10826786788668456, -0.13576361217320798, -0.984807753012208},
{-0.10826786788668456, -0.13576361217320798, 0.984807753012208},
{-0.10826786788668456, 0.13576361217320798, -0.984807753012208},
{-0.10826786788668456, 0.13576361217320798, 0.984807753012208},
{-0.049515566048790434, -0.21694186955877906, -0.9749279121818236},
{-0.049515566048790434, -0.21694186955877906, 0.9749279121818236},
{-0.049515566048790434, 0.21694186955877906, -0.9749279121818236},
{-0.049515566048790434, 0.21694186955877906, 0.9749279121818236},
{0., 0., 0.},
{0.03864035467425736, -0.16929445530699877, -0.984807753012208},
{0.03864035467425736, -0.16929445530699877, 0.984807753012208},
{0.03864035467425736, 0.16929445530699877, -0.984807753012208},
{0.03864035467425736, 0.16929445530699877, 0.984807753012208},
{0.1387395330218428, -0.17397387167523584, -0.9749279121818236},
{0.1387395330218428, -0.17397387167523584, 0.9749279121818236},
{0.1387395330218428, 0.17397387167523584, -0.9749279121818236},
{0.1387395330218428, 0.17397387167523584, 0.9749279121818236},
{0.15645160204589237, -0.07534312061707779, -0.984807753012208},
{0.15645160204589237, -0.07534312061707779, 0.984807753012208},
{0.15645160204589237, 0.07534312061707779, -0.984807753012208},
{0.15645160204589237, 0.07534312061707779, 0.984807753012208},
{0.17046092493487977, -0.746838109568473, -0.6427876096865394},
{0.17046092493487977, -0.746838109568473, 0.6427876096865394},
{0.17046092493487977, 0.746838109568473, -0.6427876096865394},
{0.17046092493487977, 0.746838109568473, 0.6427876096865394},
{0.2225209339563144, -0.9749279121818236, 0.},
{0.2225209339563144, 0., -0.9749279121818236},
{0.2225209339563144, 0., 0.9749279121818236},
{0.2225209339563144, 0.9749279121818236, 0.},
{0.5617449009293668, -0.7044058256496909, -0.4338837391175581},
{0.5617449009293668, -0.7044058256496909, 0.4338837391175581},
{0.5617449009293668, 0.7044058256496909, -0.4338837391175581},
{0.5617449009293668, 0.7044058256496909, 0.4338837391175581},
{0.6901821946798448, -0.3323742273106898, -0.6427876096865394},
{0.6901821946798448, -0.3323742273106898, 0.6427876096865394},
{0.6901821946798448, 0.3323742273106898, -0.6427876096865394},
{0.6901821946798448, 0.3323742273106898, 0.6427876096865394},
{0.9009688679024191, -0.4338837391175581, 0.},
{0.9009688679024191, 0., -0.4338837391175581},
{0.9009688679024191, 0., 0.4338837391175581},
{0.9009688679024191, 0.4338837391175581, 0.}}Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: kram1032 on January 06, 2010, 06:22:21 PM
could you look at those positions? Just thinking about the Mandelbrot set and looking at the two points in the center of the set lets me guess that those are exactly positions of minibulbs. :)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: bugman on January 06, 2010, 07:10:51 PM
Hey, that is really interesting BradC! This might prove useful for finding more Mandelbulb 3D orbits (so far I have only been able to find the main 3D cardioid):
http://www.fractalforums.com/theory/mandelbrot-pearls/msg11043/#msg11043
I am curious what method did you use to converge the random points and how long did it take to solve?
http://www.fractalforums.com/theory/mandelbrot-pearls/msg11043/#msg11043
I am curious what method did you use to converge the random points and how long did it take to solve?
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: BradC on January 06, 2010, 09:22:49 PM
I used FindRoot with two starting values for each variable, so Mathematica wouldn't try to take derivatives. It took about 10 minutes on my old P4. Two starting values for each of the three variables = 6 values total, I chose them all independently, uniformly distributed in (-1.5, 1.5). This strategy could probably be improved. It found most of the points fairly quickly, but then took a long time to find the last few values near the poles. Better starting values would probably help that.
Some other implementation details:
With a real or complex polynomial of one variable, I think you can divide out zeros as you find them to avoid finding the same zero twice. I'm not sure how to do this with a function of more than one variable.
This approach probably wouldn't be able to realistically find all solutions for the 3rd level Mandelbulb or higher, but it should be able to find lots of solutions farily easily.
Some other implementation details:
- WorkingPrecision everywhere is 4 $MachinePrecision; then I used N[] at the very end to get doubles.
- I mapped all failed convergences to a special value so I could recognize and remove them at the end.
- I removed duplicates using Gather, considering two points to be equal if they were less than 10^-8 apart. This was smaller than necessary.
- FindRoot was sometimes calling my function with x=y=0, so I had to handle that case specially to avoid messages and Indeterminate values.
- I ran this procedure 1000 points at a time until the resulting solution set looked symmetrical and seemed to stop growing. It took 5 runs (5000 FindRoot calls).
With a real or complex polynomial of one variable, I think you can divide out zeros as you find them to avoid finding the same zero twice. I'm not sure how to do this with a function of more than one variable.
This approach probably wouldn't be able to realistically find all solutions for the 3rd level Mandelbulb or higher, but it should be able to find lots of solutions farily easily.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: bugman on January 07, 2010, 03:58:43 AM
Triplex roots can also be found analytically:
http://www.fractalforums.com/theory/triplex-algebra/msg9477/#msg9477
There are 64 roots for the 8th order power formula. I wonder if this can be applied to exactly solve for the 64 points that you found.
http://www.fractalforums.com/theory/triplex-algebra/msg9477/#msg9477
There are 64 roots for the 8th order power formula. I wonder if this can be applied to exactly solve for the 64 points that you found.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: kram1032 on January 07, 2010, 09:50:47 PM
only 64 roots?
I'd expect a triplex to have power³ roots which would be 512 for 8...
The number of roots would also increase per iterations, I guess?
And as the set fully includes the Mset of power, it should be fully including all the roots for power 8 too...
Or how's that 64 meant?
I'd expect a triplex to have power³ roots which would be 512 for 8...
The number of roots would also increase per iterations, I guess?
And as the set fully includes the Mset of power, it should be fully including all the roots for power 8 too...
Or how's that 64 meant?
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on May 26, 2013, 10:09:36 PM
Sorry for bumping this old thread but:
It uses some kind of interval arithmetics and a divide and conquere strategy.
BTW! The bottleneck with Durand-Kerner method comes the product in the denominators which makes it O(N²), the evaluation of the Mandelbrot polynomial itself is O(log(N)). It should be possible to transform the Durand-Kerner method to make it O(N E(N)) (where E(N) is the complexity of the polynomial evaluation) by using a Fast multipole method to compute the denominators. This is because the "potential" from a cluster of approximated roots that are far from the curent approximation root have a small magnitude and doesn't change too much in the vicinity of the current apprximation.
Because I'm not a math wizard, I challenge those of you computer math wizards to write such algorithm or prove me wrong (which is equally probable). ;D
For those of you computer math wizards out there, I challenge you to find more roots than this. The first one to find 1 million roots wins! 8)
Ok! here is attached an evaldraw script that can compute 1 million roots (and maybe 2). :tool: ;DIt uses some kind of interval arithmetics and a divide and conquere strategy.
BTW! The bottleneck with Durand-Kerner method comes the product in the denominators which makes it O(N²), the evaluation of the Mandelbrot polynomial itself is O(log(N)). It should be possible to transform the Durand-Kerner method to make it O(N E(N)) (where E(N) is the complexity of the polynomial evaluation) by using a Fast multipole method to compute the denominators. This is because the "potential" from a cluster of approximated roots that are far from the curent approximation root have a small magnitude and doesn't change too much in the vicinity of the current apprximation.
Because I'm not a math wizard, I challenge those of you computer math wizards to write such algorithm or prove me wrong (which is equally probable). ;D
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: M Benesi on May 27, 2013, 01:34:44 AM
Are there roots at the center of each p/q bulb? Then you could generate roots with (p/q is rational) -
})
From the wiki. (http://en.wikipedia.org/wiki/Mandelbrot_set#Main_cardioid_and_period_bulbs)
From the wiki. (http://en.wikipedia.org/wiki/Mandelbrot_set#Main_cardioid_and_period_bulbs)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on May 27, 2013, 01:57:04 PM
AFAIU, that formula gives the points of tangency of the disks with the main cardioid not roots.
Of course, each root of the N-th Mandelbrot polynomial is the center of an hyperbolic componenent. Because N is finite, There is a finite number of roots (2^(N-1) roots) so it's impossible to have the roots for all p/q which is infinite.
Of course, each root of the N-th Mandelbrot polynomial is the center of an hyperbolic componenent. Because N is finite, There is a finite number of roots (2^(N-1) roots) so it's impossible to have the roots for all p/q which is infinite.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: M Benesi on May 27, 2013, 07:45:02 PM
</snip>
We could have
. Seems like the roots would tend to 
as well. </snip>
Could we follow the tangent out from the main Cardioid, determining the size of the bulb by the p/q ratio to find the approximate center.... and then do Newton's approx (or some other approximation function).
</snip>
This is sortof neat:
http://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set#color_is_proportional_to_magnitude_of_zn
Above it is other info about the roots.
We could have
Could we follow the tangent out from the main Cardioid, determining the size of the bulb by the p/q ratio to find the approximate center.... and then do Newton's approx (or some other approximation function).
</snip>
This is sortof neat:
http://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set#color_is_proportional_to_magnitude_of_zn
Above it is other info about the roots.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on May 29, 2013, 05:13:56 PM
Could we follow the tangent out from the main Cardioid, determining the size of the bulb by the p/q ratio to find the approximate center.... and then do Newton's approx (or some other approximation function).
Very nice idea and I think it's possible: There is a formula (http://www.mrob.com/pub/muency/muatomsizes.html) giving an estimation of the size of a mu-atom attached to a cardioid.One have to be carful about the accuracy of floating number format used when solving for very small mu-atoms.
In the script I posted above, using double floats will probably miss some roots when the order of Mandelbrot polynomial is above 24 or so. Anyway, the script is limited to the 22th order because of the limitations to the size of the array storing the roots. Notice that the algorithm I used is naturally parallelizable so one can write a multithreaded application to speed up things. I'm wondering if it's possible to use cuda or opencl (with say, an implementation of double-double or quad-floats) in order to compute even more roots (maybe on billion ;D).
It looks like Piers Lawrence (http://publish.uwo.ca/~plawren/) have already computed more than 1 million roots. There is in particular, a nice animation of a kind of buddhabrot obtained from them.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: M Benesi on May 30, 2013, 07:30:02 AM
At least it looks like he's trying for it:
:D
You know, that root animation on his main page begins to look like a 2d Mandelbulb towards the end.
And here is a bit of planetary buddhabrot nebulousness f (http://infinity-imagined.tumblr.com/post/51437486231/a-collection-of-planetary-nebulae-clouds-of)or you.
Quote from: PiersLawrence
The goal of this project is ultimately compute over one million roots of the Mandelbrot polynomials (i.e. compute the roots of $p_{20}(z)$). This will be achieved through the use of multiple CPU and GPU computations using sparse iterative matrix algorithms. (http://publish.uwo.ca/~plawren/coding/mandel.html)
:D
You know, that root animation on his main page begins to look like a 2d Mandelbulb towards the end.
And here is a bit of planetary buddhabrot nebulousness f (http://infinity-imagined.tumblr.com/post/51437486231/a-collection-of-planetary-nebulae-clouds-of)or you.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on June 01, 2013, 03:37:07 PM
Thank you M Benesi. :)
This poster (http://publish.uwo.ca/~plawren/papers/eccad2012_poster_Piers_Lawrence.pdf) says:
The "was" and "were" suggest that the computations are over, thought it is not very clear to me.
This poster (http://publish.uwo.ca/~plawren/papers/eccad2012_poster_Piers_Lawrence.pdf) 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.
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.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: M Benesi 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. :D
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
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty 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. :D
Humm... yes! Piers Lawrence is definitly the winner. :'( :D 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
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
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:
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: M Benesi 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.
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.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty 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 :sad1:). Now it computes 1M roots in a little more than 3 mn on my computer instead of 20 mn for the previous version.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: plawrence 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 (http://people.cs.kuleuven.be/~piers.lawrence/posters/IAP_study_day_2013_poster.pdf) or here: IAP Poster 2013 (alternative) (http://research.lawrences.info/posters/IAP_study_day_2013_poster.pdf).
Let me say something about the progression of this work:
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)
cost rather than }-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.
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)
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 }-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.
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,
say, such that as we vary from zero to one, the roots trace a path from the (double) roots of )
(plus one extra point at zero) to those of )
.
From this new equation, we can form a differential equation describing the path of the roots (see the poster (http://research.lawrences.info/posters/IAP_study_day_2013_poster.pdf)). 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-domain, but what was surprising was the speed: to compute all roots from 
to 
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
, since we have to worry about the singularities due to the higher derivatives of .
I attach the code to compute the roots via the homotopy continuation method, it seems clear that one could get over 1 billion roots (
) 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 , 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
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 (http://people.cs.kuleuven.be/~piers.lawrence/posters/IAP_study_day_2013_poster.pdf) or here: IAP Poster 2013 (alternative) (http://research.lawrences.info/posters/IAP_study_day_2013_poster.pdf).
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
- 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
- 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,
From this new equation, we can form a differential equation describing the path of the roots (see the poster (http://research.lawrences.info/posters/IAP_study_day_2013_poster.pdf)). 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
The complication is the singularities, these only get more dense as we increase
I attach the code to compute the roots via the homotopy continuation method, it seems clear that one could get over 1 billion roots (
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on November 11, 2015, 03:54:03 PM
Thx. It's interesting topic.
Please add informastions about progress.
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
Adam
See also :
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
Please add informastions about progress.
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
Adam
See also :
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: plawrence 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
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski 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
Adam
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
Adam
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: plawrence on November 11, 2015, 07:18:28 PM
Hi Adam,
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
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski 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
Adam
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
Adam
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: plawrence on November 12, 2015, 01:05:27 PM
There is no removal of roots of lower periods.
For period 14 there should be}-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
For period 14 there should be
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: JArtherton 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 (http://www.scideveloper.com/Fractals/fractals.html)
These are found using 320 bit math for higher precision.
Data available at www.scideveloper.com/Fractals/fractals.html (http://www.scideveloper.com/Fractals/fractals.html)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on November 24, 2015, 04:11:12 PM
:joy:
Hello all!
@ plawrence: Congratulations for winning the contest! :D
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).
Hello all!
@ plawrence: Congratulations for winning the contest! :D
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).
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: JArtherton 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
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
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on November 24, 2015, 08:24:23 PM
Perhaps this (excellent) paper is of interest for this programme: http://webserver2.tecgraf.puc-rio.br/~lhf/ftp/doc/ij.pdf
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on November 24, 2015, 08:26:17 PM
@ JArtherton
What do you think aboout atom domain :
http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
https://commons.wikimedia.org/wiki/File:Mandelbrot_Atom_Domains_Animation.gif
?
See also :
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers
What do you think aboout atom domain :
http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
https://commons.wikimedia.org/wiki/File:Mandelbrot_Atom_Domains_Animation.gif
?
See also :
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: cKleinhuis on November 24, 2015, 10:06:50 PM
people, i would like to collect the roots and store them in a public useable database provided here through a forum api, the mandelbrot roots serve very interesting aspects, for once the roots make nice julia sets ;)
regarding webspace and api programming i am totally available in doing so (at least for a while - 1gig mysql database space available here), a realtime rendering for julia seeds of that root locations could provide interesting explorations, or importing the roots into other programs like ultrafractal or kalles fractaller (remember to go deep into the 0,0 location of that julia sets ;) )
regarding webspace and api programming i am totally available in doing so (at least for a while - 1gig mysql database space available here), a realtime rendering for julia seeds of that root locations could provide interesting explorations, or importing the roots into other programs like ultrafractal or kalles fractaller (remember to go deep into the 0,0 location of that julia sets ;) )
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on November 26, 2015, 06:07:48 PM
@ cKleinhuis
Great idea.
The mininal set of data is : period, cx, cy
but it could be extended to :
* radius see : https://en.wikibooks.org/wiki/Fractals/Computer_graphic_techniques/2D/plane#radius
* angled internal addresses
* angles of rays landing on the root points of such component
* ....
Great idea.
The mininal set of data is : period, cx, cy
but it could be extended to :
* radius see : https://en.wikibooks.org/wiki/Fractals/Computer_graphic_techniques/2D/plane#radius
* angled internal addresses
* angles of rays landing on the root points of such component
* ....
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: cKleinhuis on November 26, 2015, 08:19:11 PM
correct, the roots database could contain pretty much information, grouped by location, depth, i hope this does not explode too early ;) we can chat in private about setting up a minimal api, where you could send solutions to, each solution is then marked as new, some ways to determine the correctness could be to save versions corresponding to program names, versions and dates that suggested them
the thing is going to be handling of the huge arbitrary precision logic in that database
the thing is going to be handling of the huge arbitrary precision logic in that database
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on November 26, 2015, 08:50:55 PM
cKleinhuis
Posts: 6666
:evil1:
Posts: 6666
:evil1:
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on November 29, 2015, 12:45:15 PM
Hi,
I've implemented my mandelbrot polynomial roots finder in c++ (attached) but there is someting I really don't undestand:
- when compilling with -O3 it takes 6.8s to compute the 8192 roots of the 14th mandelbrot polynomial on i7 processor.
- When compiling with -Ofast it takes only 0.3s seconds.
Why such a difference in speed?
When using float or long double instead of double, -Ofast doesn't improve computation speed. :/
I'm using GCC 5.2.0
It looks like the compiler doesn't use hardware fp unless one uses doubles and -Ofast switch... strange!
I did once something similar. See attached evaldraw scripts.
I've implemented my mandelbrot polynomial roots finder in c++ (attached) but there is someting I really don't undestand:
- when compilling with -O3 it takes 6.8s to compute the 8192 roots of the 14th mandelbrot polynomial on i7 processor.
- When compiling with -Ofast it takes only 0.3s seconds.
Why such a difference in speed?
When using float or long double instead of double, -Ofast doesn't improve computation speed. :/
I'm using GCC 5.2.0
It looks like the compiler doesn't use hardware fp unless one uses doubles and -Ofast switch... strange!
@ JArtherton
What do you think aboout atom domain :
http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
https://commons.wikimedia.org/wiki/File:Mandelbrot_Atom_Domains_Animation.gif
?
See also :
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers
What do you think aboout atom domain :
http://mathr.co.uk/blog/2014-11-02_practical_interior_distance_rendering.html
https://commons.wikimedia.org/wiki/File:Mandelbrot_Atom_Domains_Animation.gif
?
See also :
https://en.wikibooks.org/wiki/Fractals/Iterations_in_the_complex_plane/Mandelbrot_set/centers
I did once something similar. See attached evaldraw scripts.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on November 29, 2015, 01:36:54 PM
Perhaps this (excellent) paper is of interest for this programme: http://webserver2.tecgraf.puc-rio.br/~lhf/ftp/doc/ij.pdf
Hahaha, and of course the mighty Knighty jumps in with not just interval arithmetic, but affine arithmetic! What a legend O0
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on November 29, 2015, 03:39:26 PM
O0 :embarrass:
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on November 29, 2015, 06:00:13 PM
Quote
When using float or long double instead of double, -Ofast doesn't improve computation speed.
one thing is that long double is going to use the 80bit fpu and as such cant be vectorized. not sure why there would be no difference with single precision float though. do you use any other compiler options?Quote
Why such a difference in speed?
ffast-math (brought in by Ofast) is where it's at for floating point performance ;D if you want speed and don't need strict IEEE compliance definitely use fast relaxed math options whenever you can get away with it.i ran a few tests of your code on my haswell system, gcc 5.2:
Code:
O3 0.02602
O3 std=gnu++14 0.02560
O3 march=haswell 0.02420
O3 march=haswell std=gnu++14 0.02462
Ofast 0.00838
Ofast std=gnu++14 0.00854
Ofast march=haswell 0.00805
Ofast march=haswell std=gnu++14 0.00803
of course running them 300 times and averaging the results would be a more useful benchmark if one was so inclined :)
PS. theres some really good stuff happening in c++ land in recent years. c++98 is nearing 20 years old now. let's start writing some modern c++ :D
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on November 29, 2015, 07:19:49 PM
g++ -Ofast -std=c++11 m.cpp -Wall // time taken : 0.007739 s
g++ -Ofast -march=haswell -std=gnu++14 m.cpp -Wall // m.cpp:1:0: error: bad value (haswell) for -march= switch
g++ -Ofast -march=native -std=c++1y m.cpp -Wall // time taken : 0.026236 s
g++ -Ofast -march=haswell -std=gnu++14 m.cpp -Wall // m.cpp:1:0: error: bad value (haswell) for -march= switch
g++ -Ofast -march=native -std=c++1y m.cpp -Wall // time taken : 0.026236 s
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on November 29, 2015, 07:52:25 PM
adam, gross! old version of gcc i guess :'(
as a followup to my last post, ffast-math also (SHOULD) greatly affect the performance of std::complex. in fact this is where the majority of the overall performance difference is coming from here. unfortunately this is not currently the case for non-gcc users. there was a rather eye-opening presentation on this issue at CppCon 2015:
https://youtu.be/he-XVt1xIE0
to illustrate this on knighty's code:
clang you turd, get your shit together! and micros~1, well, just die in a fire. :angel1:
as a followup to my last post, ffast-math also (SHOULD) greatly affect the performance of std::complex. in fact this is where the majority of the overall performance difference is coming from here. unfortunately this is not currently the case for non-gcc users. there was a rather eye-opening presentation on this issue at CppCon 2015:
https://youtu.be/he-XVt1xIE0
to illustrate this on knighty's code:
Code:
gcc-5.2 Ofast march=haswell std=gnu++14 0.00803
clang-3.7 Ofast march=haswell std=gnu++14 0.02023
clang you turd, get your shit together! and micros~1, well, just die in a fire. :angel1:
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on November 29, 2015, 09:37:46 PM
What about small modification :
Code:
void saveThem(int period){
std::string name= std::to_string(period) + ".txt"; // C++11 for std::to_string
ofstream SaveFile(name);
switch (sizeof(T)){
case 4 : {SaveFile.precision(8); break;}
case 8 : {SaveFile.precision(16); break;}
case 10:
case 12: {SaveFile.precision(20); break;}
case 16: {SaveFile.precision(25); break;} //__float128. doesn't work ... need more info
default : SaveFile.precision(16);// should not happen :/
}
SaveFile.setf(ios::fixed);
for (unsigned int i = 0; i < m_tab.size() /*m_cIndex*/; i++){
SaveFile << m_tab[i].real() << "\t" << m_tab[i].imag() << endl;
}
}
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on November 29, 2015, 09:42:25 PM
adam, gross! old version of gcc i guess :'(
gcc version 4.8.4 (Ubuntu 4.8.4-2ubuntu1~14.04
not the latest gcc 5.2 version
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: zebastian on November 29, 2015, 10:43:18 PM
when we are getting down on performance, i thought i could just throw in some openmp pragmas to make use of multicore.
binary is by factor (~2 / 5 * core count) faster, but calculates wrong and crashes from time to time :D
probably some thread unsafety...
the main part is:
if(depth <= 2)
{
#pragma omp parallel for
for(int i = 0; i < 4; i++){
findRootsRecurse(cc + complex<T>(-hw + (i / 2) * hw,-hw + (i % 2) * hw), hw, depth + 1);
}
}
else{
findRootsRecurse(cc + complex<T>(-hw,-hw),hw, depth + 1);
findRootsRecurse(cc + complex<T>( hw,-hw),hw, depth + 1);
findRootsRecurse(cc + complex<T>( hw, hw),hw, depth + 1);
findRootsRecurse(cc + complex<T>(-hw, hw),hw, depth + 1);
}
if someone is interested, full source code in attachment.
dont forget to put -fopenmp as compiler flag!
binary is by factor (~2 / 5 * core count) faster, but calculates wrong and crashes from time to time :D
probably some thread unsafety...
the main part is:
if(depth <= 2)
{
#pragma omp parallel for
for(int i = 0; i < 4; i++){
findRootsRecurse(cc + complex<T>(-hw + (i / 2) * hw,-hw + (i % 2) * hw), hw, depth + 1);
}
}
else{
findRootsRecurse(cc + complex<T>(-hw,-hw),hw, depth + 1);
findRootsRecurse(cc + complex<T>( hw,-hw),hw, depth + 1);
findRootsRecurse(cc + complex<T>( hw, hw),hw, depth + 1);
findRootsRecurse(cc + complex<T>(-hw, hw),hw, depth + 1);
}
if someone is interested, full source code in attachment.
dont forget to put -fopenmp as compiler flag!
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 01, 2015, 02:09:34 PM
Thank you for the tests (and the switches.)
I gess you tested the program with default value of period. You can change it by passing the period as a command line argument (for a period of 25 : MandelRoots 25 ... beweare that's more than 16 millions roots. Using double float you can go for 31 period which means more than 1 billion roots... if it works correctly)
The numbers I gave were with a 32bit build. -Ofast gives 25 times faster results than -O3. for a 64bit buil its "only" about 3 times faster.
Using long double or floats there is no big difference between -O3 and -Ofast and they are much slower than double.
I don't think vectorization have something to do with this discreapancy because the evaldraw script is almost as fast as -Ofast (with doubles). Evaldraw's compiler, which is a JIT, does'nt do any vectorization optimization. It uses SSE though.
@ Adam Majewski: Thank you for the tip. I didn't know about the -std=c++11 switch. :D :embarrass:
@ zebastian: Thank you! I gess the problem comes from store() which is called from findRootsRecurse(). store() stores the root in a vector array. It is not thread safe at all.
Maybe using multiple instances of MRootsCalc class ?
I gess you tested the program with default value of period. You can change it by passing the period as a command line argument (for a period of 25 : MandelRoots 25 ... beweare that's more than 16 millions roots. Using double float you can go for 31 period which means more than 1 billion roots... if it works correctly)
The numbers I gave were with a 32bit build. -Ofast gives 25 times faster results than -O3. for a 64bit buil its "only" about 3 times faster.
Using long double or floats there is no big difference between -O3 and -Ofast and they are much slower than double.
I don't think vectorization have something to do with this discreapancy because the evaldraw script is almost as fast as -Ofast (with doubles). Evaldraw's compiler, which is a JIT, does'nt do any vectorization optimization. It uses SSE though.
@ Adam Majewski: Thank you for the tip. I didn't know about the -std=c++11 switch. :D :embarrass:
@ zebastian: Thank you! I gess the problem comes from store() which is called from findRootsRecurse(). store() stores the root in a vector array. It is not thread safe at all.
Maybe using multiple instances of MRootsCalc class ?
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on December 01, 2015, 04:20:01 PM
g++ -Ofast -std=c++11 m.cpp -fopenmp -Wall
a@zalman:~/cpp/mroots/m2$ time ./a.out
tab size: 1024
nbrRoots (expected) = 512
Searchig for roots...
*** Error in `./a.out': double free or corruption (!prev): 0x0000000000dd5080 ***
Przerwane (core dumped)
real 0m0.096s
user 0m0.025s
sys 0m0.073s
a@zalman:~/cpp/mroots/m2$ g++ m.cpp -std=c++11 -fopenmp -Wall
a@zalman:~/cpp/mroots/m2$ time ./a.out
tab size: 1024
nbrRoots (expected) = 512
Searchig for roots...
roots found (with duplicates): 1017
Removing duplicates...
roots found: 243
243
time taken : 0.069925 s
Saving results to save.txt
a@zalman:~/cpp/mroots/m2$ time ./a.out
tab size: 1024
nbrRoots (expected) = 512
Searchig for roots...
*** Error in `./a.out': double free or corruption (!prev): 0x0000000000dd5080 ***
Przerwane (core dumped)
real 0m0.096s
user 0m0.025s
sys 0m0.073s
a@zalman:~/cpp/mroots/m2$ g++ m.cpp -std=c++11 -fopenmp -Wall
a@zalman:~/cpp/mroots/m2$ time ./a.out
tab size: 1024
nbrRoots (expected) = 512
Searchig for roots...
roots found (with duplicates): 1017
Removing duplicates...
roots found: 243
243
time taken : 0.069925 s
Saving results to save.txt
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 08:05:32 PM
Quote
Evaldraw's compiler, which is a JIT, does'nt do any vectorization optimization. It uses SSE though.
not sure what you mean? sse is a vector unitTitle: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 01, 2015, 08:59:51 PM
not sure what you mean? sse is a vector unit
Sorry... It only uses scalar sse instructions.@zebastian:
see: http://stackoverflow.com/questions/9042571/is-stdvector-or-boostvector-thread-safe (http://stackoverflow.com/questions/9042571/is-stdvector-or-boostvector-thread-safe)
It is not straightforward to make a multithreaded version. Maybe one way is to use err... a critical section or a mutex around store().
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 09:07:34 PM
vectors are most definitely not thread safe. if one was to write proper threaded code you would definitely use a thread synchronization utility like a semaphore or mutex for such things. if you are just going to throw an openmp pragma around a block of code, there is also openmp stuff for specifying how to synchronize things.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 01, 2015, 10:02:28 PM
A simple way around this is to make N vectors (for N threads) and then just write into the vector indexed by omp_get_thread_num(). Combine after processing is done.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 10:05:31 PM
my own preference is to just use a normal C array and write into that. it is thread safe as long as your threads stick to their own designated ranges of the array.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 01, 2015, 10:09:58 PM
it is thread safe as long as your threads stick to their own designated ranges of the array.
Somewhat circular argument :tongue1:
Besides, the idea is to not have mutex locking, atomics, etc.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 01, 2015, 10:26:04 PM
Somewhat circular argument :tongue1:
Besides, the idea is to not have mutex locking, atomics, etc.
Simple! Obvious! Genius! ;DBesides, the idea is to not have mutex locking, atomics, etc.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 10:29:12 PM
Somewhat circular argument :tongue1:
Besides, the idea is to not have mutex locking, atomics, etc.
Besides, the idea is to not have mutex locking, atomics, etc.
sorry, what? i was simply pointing out that multiple threads writing to different sections of one array is thread safe and does not require any thread synchronization techniques, like mutexes or atomics as you said. simple, fast, and does not require any other screwing around like creating multiple vectors and then combining them later.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 01, 2015, 10:33:24 PM
sorry, what?
The statement is roughly, "C arrays are thread safe as long as you use them in a thread-safe way" :)
So my suggestion is basically a way to do that, with whatever structure you like (C array, vector, whatever), without any contention or locking. Note also that we don't know the number of roots ahead of time, so the simple C array would turn into something like std::vector anyway.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 10:46:42 PM
Quote
Note also that we don't know the number of roots ahead of time, so the simple C array would turn into something like std::vector anyway.
ah, indeed that puts a damper on using a C array :)Quote
The statement is roughly, "C arrays are thread safe as long as you use them in a thread-safe way"
i guess im just confused as to why you feel the need to be hostile about this...lol...it wasnt a "circular argument." i was simply pointing out that writing to C arrays is not thread UNsafe in the way that writing to std containers like std::vector is. did you already know that? thats wonderful. O0 i thought this was a discussion board where people discuss things, my bad. besides, a lot of people here tend to be mathematicians first, with the code writing often happening in a more by-gee and by-golly sort of manner, so fostering code-related discussion seems like a reasonable thing to me :)Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 01, 2015, 10:50:23 PM
Awwww, I wasn't being hostile I think? No hostility intended - I hope you see what I mean about circular argument (argument not as in "having a heated argument", rather "point being made"). Actually I would say calling "my" approach "screwing around" is more hostile than suggesting a circular argument is circular...
PS. God knows I prefer low-level arrays too, where possible ;) But there's nothing to fear in std::vector, it's indistinguishable (compiles to just a MOV) from an array if you use resize() and only index in valid range, with range checks turned off (in MSVC).
PS. God knows I prefer low-level arrays too, where possible ;) But there's nothing to fear in std::vector, it's indistinguishable (compiles to just a MOV) from an array if you use resize() and only index in valid range, with range checks turned off (in MSVC).
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 01, 2015, 10:59:27 PM
oops!
wrong quote. I was actually referring to quazor's idea: getting rid of the syncronization. Doing it using multiple arrays also solves the problem. The problem is not that we don't the number of roots (because we do: 2^(m_maxiter) but it is that the number of duplicates is not obvious to estimate (potentially a root can be reported up to 8 times because of the affine arithmetics and rounding errors).
wrong quote. I was actually referring to quazor's idea: getting rid of the syncronization. Doing it using multiple arrays also solves the problem. The problem is not that we don't the number of roots (because we do: 2^(m_maxiter) but it is that the number of duplicates is not obvious to estimate (potentially a root can be reported up to 8 times because of the affine arithmetics and rounding errors).
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 01, 2015, 11:01:42 PM
*sigh* I'd better leave this discussion 88)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 11:08:54 PM
yes i know the word argument has more meanings than "a heated confrontation." jeez, you just cant stop talking down to me can you ;D
and yes, i am also familiar with the idea of a circular argument. an example on the wikipedia page for circular reasoning is "Whatever is less dense than water will float, because whatever is less dense than water will float." my statement was not really intended (nor written?) as "writing to C arrays is thread safe because writing to C arrays is thread safe," the point was simply "writing to C arrays is thread safe whereas writing to std containers in the same manner is not thread safe." to use the stuff that floats analogy, i said "boats float, whereas rocks sink," not "boats float because boats float." if you wanted to respond to my statement in a weird, unprovoked, passive aggressive way, it would have made more sense to just say "that is obvious. i already knew that." calling it a circular argument is neither accurate nor...called for? as again, i thought this was a discussion board, not a debate club or...whatever...
and yes, i am also familiar with the idea of a circular argument. an example on the wikipedia page for circular reasoning is "Whatever is less dense than water will float, because whatever is less dense than water will float." my statement was not really intended (nor written?) as "writing to C arrays is thread safe because writing to C arrays is thread safe," the point was simply "writing to C arrays is thread safe whereas writing to std containers in the same manner is not thread safe." to use the stuff that floats analogy, i said "boats float, whereas rocks sink," not "boats float because boats float." if you wanted to respond to my statement in a weird, unprovoked, passive aggressive way, it would have made more sense to just say "that is obvious. i already knew that." calling it a circular argument is neither accurate nor...called for? as again, i thought this was a discussion board, not a debate club or...whatever...
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 01, 2015, 11:21:50 PM
:headbatting:
:D
:D
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 01, 2015, 11:27:06 PM
lol :)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 02, 2015, 02:19:54 PM
Ok!
Yesterday, it was late, I was tired...
Returning back and reading more carefully:
I did a fast test with <array> and found that the overhead of not using reserve() is not terrible. Therefore and because using c arrays raises the problem of correct amount of memory to allocate. I will use lycium's idea.
Finally, big apologies to both of you and thank you very much for your awesome ideas. :)
Yesterday, it was late, I was tired...
Returning back and reading more carefully:
A simple way around this is to make N vectors (for N threads) and then just write into the vector indexed by omp_get_thread_num(). Combine after processing is done.
And:(...) the idea is to not have mutex locking, atomics, etc.
It was these ideas That I found genius.my own preference is to just use a normal C array and write into that. it is thread safe as long as your threads stick to their own designated ranges of the array.
I also found this idea awesome.I did a fast test with <array> and found that the overhead of not using reserve() is not terrible. Therefore and because using c arrays raises the problem of correct amount of memory to allocate. I will use lycium's idea.
Finally, big apologies to both of you and thank you very much for your awesome ideas. :)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 02, 2015, 02:49:02 PM
You're too kind, and please don't apologise on my behalf :beer:
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 02, 2015, 03:57:16 PM
I didn't. I was talking for myself.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: lycium on December 02, 2015, 04:04:38 PM
Oh, ehm, I meant, not need to apologise to me / because of anything I said etc. Anyhow, glad we're cool :beer: (I don't actually drink beer either, but the symbolism is there)
Edit: one last thing about the std::vector, on Windows/MSVC you'll want to build with _SECURE_SCL=0 defined in the preprocessor settings. This is what I was referring to earlier about making them same perf to access as arrays; afaik GCC etc doesn't do range checking.
Edit: one last thing about the std::vector, on Windows/MSVC you'll want to build with _SECURE_SCL=0 defined in the preprocessor settings. This is what I was referring to earlier about making them same perf to access as arrays; afaik GCC etc doesn't do range checking.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: quaz0r on December 02, 2015, 04:57:18 PM
Quote from: http://en.cppreference.com/w/cpp/container/vector
at access specified element with bounds checking
(public member function)
(public member function)
Quote from: http://en.cppreference.com/w/cpp/container/vector
operator[] access specified element
(public member function)
(public member function)
is micros~1 really so nonstandard that you have to turn off bounds checking for operator[] ? i guess it wouldnt surprise me :sad1:
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: hgjf2 on December 05, 2015, 09:23:49 AM
is micros~1 really so nonstandard that you have to turn off bounds checking for operator[] ? i guess it wouldnt surprise me :sad1:
Wow! This topic have many replies. COOL!:peacock: :star:
Reminder: This chapter isn't competition. The real competition will begin on 1 may at chapter "Competitions and contests" when we will say: "Go faster else will beating our the horses!" or "Fly like the wind Artax!"
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: cyseal on December 05, 2015, 12:21:34 PM
What tools and algortihms are you using ?
Can it be done in Mathematica ?
Can it be done in Mathematica ?
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on December 05, 2015, 04:17:24 PM
IMHO results are more important that tools.
Questions about algorithms, precision of computations are very important.
BTW :
http://arxiv.org/pdf/1508.02935.pdf
HTH
Questions about algorithms, precision of computations are very important.
BTW :
http://arxiv.org/pdf/1508.02935.pdf
HTH
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 11, 2015, 05:19:23 PM
Finally got multithreading to work using c++11, so not with pthreads. I'll try to post the code next sunday after some cleaning.
What tools and algortihms are you using ?
Can it be done in Mathematica ?
It should be possible.Can it be done in Mathematica ?
BTW :
http://arxiv.org/pdf/1508.02935.pdf
Interesting.http://arxiv.org/pdf/1508.02935.pdf
Quote
The total number of iterations required to find all roots was 3 056 825 939 654 ≈ 2.78d² (...) The total number of Newton iterations required, while theoretically bounded above by O(d2log4d)
Lawrence's method is probably O(d log d). Mine is maybe between O(d log d) and O(d log2d). Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 13, 2015, 12:07:58 PM
Here it is.
Computes one million roots in 22.5 seconds on a i7. O0
Computes one million roots in 22.5 seconds on a i7. O0
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: cKleinhuis on December 13, 2015, 02:01:04 PM
hey people, regarding a database, if you still think it is useful, i prepared a project able to store roots for various formulas, it is not finished yet
and user accounts need to be created, but here you have a swagger2 api definition from the live backend:
http://qa.feuerware.com:8100/ufp-api/v2/api-docs
the idea is:
- create formulas
- store polynomal roots with formula reference, rootindex, iteration depth
- numbers are submitted as arbitrary long strings, and get stored as 18 digits of long for searching through the results
as soon as i am happy with the setup i will provide a tutorial on how to use the api, in general every user would need an account,
but perhaps we can keep the reading of results in a public manner as long as server says "yes" ;)
and user accounts need to be created, but here you have a swagger2 api definition from the live backend:
http://qa.feuerware.com:8100/ufp-api/v2/api-docs
the idea is:
- create formulas
- store polynomal roots with formula reference, rootindex, iteration depth
- numbers are submitted as arbitrary long strings, and get stored as 18 digits of long for searching through the results
as soon as i am happy with the setup i will provide a tutorial on how to use the api, in general every user would need an account,
but perhaps we can keep the reading of results in a public manner as long as server says "yes" ;)
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: Adam Majewski on December 13, 2015, 03:30:42 PM
BTW
There is also MPFR program by Claude :
https://en.wikibooks.org/wiki/Fractals/Mathematics/Newton_method#center
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: knighty on December 13, 2015, 04:48:25 PM
@ cKleinhuis:
To save storage space, I think it would be better to compute the root when the user asks for it. Given a region, the computations, to get the smallest period root inside it, would take some milliseconds, even with big precision and a java implementation. Claude's program is a good candidate for that purpose but there are cases where it fails (encoutering a zero derivative).
Here is attached a new version. I used my own implementation of complex class. Now using long double is only 3 times slower than using doubles. order 14 is computed in about 0.6s instead of 2.62s on an i7 (both multithreaded).
Moreover the difference in speed between using -O3 and -Ofast is small now.
To save storage space, I think it would be better to compute the root when the user asks for it. Given a region, the computations, to get the smallest period root inside it, would take some milliseconds, even with big precision and a java implementation. Claude's program is a good candidate for that purpose but there are cases where it fails (encoutering a zero derivative).
Here is attached a new version. I used my own implementation of complex class. Now using long double is only 3 times slower than using doubles. order 14 is computed in about 0.6s instead of 2.62s on an i7 (both multithreaded).
Moreover the difference in speed between using -O3 and -Ofast is small now.
Title: Re: The Mandelbrot Polynomial Roots Challenge
Post by: cKleinhuis on December 13, 2015, 04:57:14 PM
the idea would be to actually store roots and provide a fast reusable database for those numbers