Logo by fractalwizz - Contribute your own Logo!

END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

it was a great time but no longer maintainable by c.Kleinhuis contact him for any data retrieval,
thanks and see you perhaps in 10 years again

this forum will stay online for reference
News: Follow us on Twitter
 
*
Welcome, Guest. Please login or register. September 23, 2019, 10:04:19 AM


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 ... 6 7 [8] 9 10 ... 17   Go Down
  Print  
Share this topic on DiggShare this topic on FacebookShare this topic on GoogleShare this topic on RedditShare this topic on StumbleUponShare this topic on Twitter
Author Topic: SuperFractalThing: Arbitrary precision mandelbrot set rendering in Java.  (Read 66954 times)
0 Members and 1 Guest are viewing this topic.
Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #105 on: July 01, 2013, 04:02:09 AM »

FWIW, Nanoscope's post-e308 scaling is working perfectly on 64-bit Windows. Java, no strictfp. I don't think I've implemented it the same as Kalle though. Kalle seems to be using a pair of doubles whose product is the real value of delta. I'm using one double and an exponent, and most iterations don't need to care about the exponent. The innermost post-e308 loop in Nanoscope is just:

delta_scaled_n+1 = delta_scaled_n * reference_orbit_n + delta_scaled_0.

No delta-squared, because if delta is O(10-308) or smaller, delta squared is smaller than delta by a similar factor and is lost in the noise. No 2* because the copy of the reference orbit used for this loop is premultiplied by 2. If the reference orbit's not about to escape (and if it does, it takes delta with it as they're sure to diverge past the bailout distance), it's confined in a disk of radius 2, so its doubled version is in a disk of radius 4, and the modulus of the first term in the sum above is at most five times larger than the previous, if delta_scaled_0 is still comparable in magnitude to delta_scaled_n, and more usually delta_scaled_0 is smaller. So, basically delta_scaled's modulus grows by a factor of four or less most iterations, and that means it takes delta_scaled around 500 iterations to go from O(1) to O(10308). Only then is delta_scaled (and delta_scaled_0) rescaled and the stored exponent adjusted. This makes for a very fast inner loop below e308.

In principle, this method allows zooming to nearly arbitrary depths; there's no further barrier at e616 or anything like that.

I've also greatly streamlined the memory use, without even having to recalculate the high precision orbit every frame or persist it to disk. High precision points are only needed in three places:

  • The first point, and only the first point, for computing delta_0 per pixel.
  • Each term R_n-1 is needed to compute the series coefficients A_n, B_n, C_n, and those only need wide exponents, not wide mantissas.
  • Each term R_n is needed to compute the premultiplied-by-two, double-precision version of itself used in the iteration loop above and in the above-e308 loop.

So we only need to store:

  • The first point R_0 at full precision.
  • A_i, B_i, C_i up to some point n, as far as will work at the current magnification. We only need wide exponents, not wide mantissas, so these take up much less space.
  • The n we've computed those up to.
  • R_n at full precision, for computing R_n+1, A_n+1, B_n+1, C_n+1 when some deeper image with this reference point may need them.
  • The points 2*R_i, up to R_maxit, at double precision.

And to compute:

  • Per sequence of zooms with a given center point R_0, the double-precision points 2*R_i up to maxit. (Kept.)
  • Per image, any additional coefficients A_i, B_i, C_i that can now be used at the deeper magnification (smaller side length), and R_i along with them, stopping at R_n_new, A_n_new, B_n_new, C_n_new. Only R_n_new and narrow-mantissa, wide-exponent versions of the new A_i need to be kept; the previous R_i can be discarded, except of course for R_0.
  • Per pixel, delta_0 and delta_n (using A_n, B_n, C_n, delta_0).
  • Per iteration, delta_i (or delta_scaled_i, and, every 500 iterations, scale and delta_scaled_0).

The high precision R_i are traversed twice, though, first to build the arrays of doubled doubles and again gradually to build up the sequence of ABC_i values. Small price to pay.

I've also found the need to identify the highest iteration region in an image pre-calculation and sometimes use a second reference point there. The algorithm I'm currently using for this proceeds in two phases. First, it samples 100 points in the image region using perturbation theory (fast!), sorts them by iterations, discards all but the top ten, and then does some graph theory manipulation to identify a big cluster of close-together points, by looking for the closest pair of samples, and then all samples transitively reachable from those two by sample-to-sample hops no longer than 1.5x their separation. The resulting possibly-reduced set of points tends to be a "blob" of high-iter material which, if off-center (i.e. not containing the primary reference point), may not be rendered correctly. This blob's diameter is measured (longest distance between any two of the samples in it) and its barycenter (average of the member samples' coordinates), and if it's significantly bigger than 1 pixel but significantly smaller than the previous examined region, the process is recursively applied to the smaller region to refine the estimate of the high iter region's core's size and center. When the "blob stops shrinking" or reaches only a few pixels in width, the second phase begins. If the "blob" is now tiny (only a few tens of pixels max), they're all iterated with high precision normal Mandelbrot calculations and the highest iteration point is returned; if the blob's still substantial, high precision calculations take over from low but a recursive refinement algorithm that takes logarithmically many points in the blob's diameter is used instead of brute scanning to try to guess the highest-iter point, gradient-climbing to reach it.

Once the high-iter point is identified, its iteration count is checked for being bigger, by more than 10, than the iteration count the PT algorithm assigns. If it is off by ten or more iterations, that means PT is bailing early in that region and we need a second reference point. The highest-iter point is used as that point, which means a second 2*S_i array and D_i, E_i, F_i series coefficients. Also, points in the image are iterated using PT with each set in turn, and the highest iteration wins. This seems to clear up the "blank blobs" in the test images I've thrown it at so far, and even a corrupted minibrot in one test (where the reference point was in an incomparably smaller minibrot some several diameters of the larger minibrot away from it). Sadly, the detection of the need for a second reference point is slowish, taking a minute or so in the e70-e100 few-thousand-iter "shallows" and half an hour at e300+ and 1e5+ iterations (java.math.BigDecimal is slow).

I am, of course, working on figuring out ways to squeeze more speed out of the highest-iter-point-finding algorithm. Even as things are, it seems to be as accurate as UF and dozens to hundreds of times faster, depending on the zoom.

AIUI, Superfractalthing also has some method for determining when a second reference point is needed; I don't know about Kallesfractaler, though, or the methods either use. I've described mine in the hope that it may be helpful to either of their authors, and in the hope that someone might have a useful suggestion for speeding my own up. Similarly, I hope my explanation of Nanoscope's e308+ loop structure can help Kalle to overcome his(?) difficulties with different behavior on different processor architectures, and to be able to go past e616.
« Last Edit: July 01, 2013, 04:47:33 AM by Pauldelbrot » Logged

cKleinhuis
Administrator
Fractal Senior
*******
Posts: 7044


formerly known as 'Trifox'


WWW
« Reply #106 on: July 01, 2013, 07:17:03 AM »

 embarrass great, do you plan to make nanoscope available for others ?
Logged

---

divide and conquer - iterate and rule - chaos is No random!
Kalles Fraktaler
Fractal Senior
******
Posts: 1458



kallesfraktaler
WWW
« Reply #107 on: July 01, 2013, 08:32:23 AM »

Similarly, I hope my explanation of Nanoscope's e308+ loop structure can help Kalle to overcome his(?) difficulties with different behavior on different processor architectures, and to be able to go past e616.
I also made the same kind of custom datatype with a double as the mantissa and an integer as the exponent, to be able to get past e616.
However I did not want to use that until e600 and instead scale ordinary doubles between e300 and e600, which works fine in 32-bit.
The solution is to use the custom datatype earlier in 64-bit, at about e450.

My first tries on series approximation and binary search shows that I should be able to render "tick-tock" in only 1.5 seconds.
But it does not render correct, I probably didn't get the minuses and the pluses correct though... smiley

Edit:
I am beginning to understand that function (2) cannot be used by itself, it can only be used to initiate function (1) so that it does not have to begin from iteration zero. A fractal that only use function (2) will look like http://www.fractalforums.com/images-showcase-(rate-my-fractal)/series-approximation/ which I think is kind of beautiful, however not what is recognized as a Mandelbrot fractal though... So, since function (2) is only accurate to about half of the iterations needed, we will unfortunately not be able to render tick-tock 640x360 in one second, something just below 10 seconds seems more realistic.
« Last Edit: July 01, 2013, 07:40:51 PM by Kalles Fraktaler » Logged

Want to create DEEP Mandelbrot fractals 100 times faster than the commercial programs, for FREE? One hour or one minute? Three months or one day? Try Kalles Fraktaler http://www.chillheimer.de/kallesfraktaler
http://www.facebook.com/kallesfraktaler
Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #108 on: July 01, 2013, 09:47:49 PM »

embarrass great, do you plan to make nanoscope available for others ?

Most likely. Needs more testing and work, and if it's to be usable by non-techies, a real user interface, first. I might be releasing the source code sooner than there's a nontechie-suitable executable artifact, though. Most likely on Github.
Logged

Kalles Fraktaler
Fractal Senior
******
Posts: 1458



kallesfraktaler
WWW
« Reply #109 on: July 02, 2013, 12:18:45 PM »

I have a question to mrflay and Pauldebrot, since you have implemented the "series approximation":
I believe you are validating the approximation in several pixels in the images, in order to find the highest iteration where the approximation is valid? Is that really gaining anything? When I test different pixels I always get the exact same number of valid iterations...

I finally understand mrflays statement
"the time taken rendering Mandelbrot images is largely independent of depth and iteration count, and mainly depends on the complexity of the image being created".
smiley
Logged

Want to create DEEP Mandelbrot fractals 100 times faster than the commercial programs, for FREE? One hour or one minute? Three months or one day? Try Kalles Fraktaler http://www.chillheimer.de/kallesfraktaler
http://www.facebook.com/kallesfraktaler
Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #110 on: July 02, 2013, 12:45:17 PM »

I have a question to mrflay and Pauldebrot, since you have implemented the "series approximation":
I believe you are validating the approximation in several pixels in the images, in order to find the highest iteration where the approximation is valid? Is that really gaining anything? When I test different pixels I always get the exact same number of valid iterations...

Noting that after a short time, A, B, and C are separated by roughly equal numbers of orders of magnitude, I use the current values of B and C to estimate the uncomputed D, coefficient of delta4, and make sure that it times the fourth power of the image side length remains negligibly small. There's no testing against pixels; it just works.

On the other hand, I've overhauled blob-purging (already). Now instead of switching to high precision calculations for whole grids of pixels, it does only a handful at full precision. First, the initial search for a high iter region terminates if either five or more pixels yield identical, non-maxiter smooth iteration values (blob!) or the "net" around it gets tiny or stops shrinking much; then, the center of the "netted" region is iterated at high precision and checked for discrepancy with the PT approximation at the same point. If they differ, the center is used as a reference orbit to continue the process: the net tightens some more, if a blob of constant iter now resolves into structures around a smaller blob (or a point). When the net again stops tightening, the center is again checked at full precision. If the iterations are higher than the current secondary reference point, that becomes the new secondary reference point. In any event, if the net stops shrinking or becomes 1 pixel or so in size, the process terminates and full image rendering can begin. This new method of locating any needed secondary reference point is about twice as fast as the old, and more accurate besides; on that test image where the old took 30 minutes, the new takes 15, despite being single threaded (the old parallelized some things, and I have six-core hardware here). The downside: the slow bits (new reference orbit calculations) aren't amenable to parallelizing. (Keep in mind that that test image is below e308 and has 300,000+ iterations everywhere -- 400,000 in the very centers of the blobs.)

I've also substantially sped up image renders with two reference points. The renderer tracks the lowest iteration count (as returned using the primary reference point) for which a discrepancy was observed between iterations using both reference points, and the highest iteration lower than that for which both were used but no discrepancy was observed. At first, the lowest is set to maxiter and the highest to zero, and points all get both reference points used. Any point that iterates with the primary reference point to a count higher than either "water mark" gets the secondary reference point used, and the max of the iteration values from both iterations used. Any point iterated with both reference points is also used to update the two "water marks", using the iteration count from using the primary reference point to set them. Finally, any point where using the primary reference point gets an iteration count lower than both "water marks" is taken at face value, without using the other reference point. If the blobs don't occupy a large percentage of the image pixels, this results in a near-2x speedup, to near parity with blob-free images.

In the very near future I'll be knocking your socks off in the Images Showcase forum with a result from all of this. wink
Logged

Kalles Fraktaler
Fractal Senior
******
Posts: 1458



kallesfraktaler
WWW
« Reply #111 on: July 02, 2013, 04:16:06 PM »

In the very near future I'll be knocking your socks off in the Images Showcase forum with a result from all of this. wink
Looking forward to that smiley
Logged

Want to create DEEP Mandelbrot fractals 100 times faster than the commercial programs, for FREE? One hour or one minute? Three months or one day? Try Kalles Fraktaler http://www.chillheimer.de/kallesfraktaler
http://www.facebook.com/kallesfraktaler
Roquen
Iterator
*
Posts: 180


« Reply #112 on: July 12, 2013, 11:59:32 AM »

I wondered about that as well, since we use 32-bit and 64-bit Java code at work for scientific calculations.

According to the JLS, if you use the 'strictfp' keyword: "Within an FP-strict expression, all intermediate values must be elements of the float value set or the double value set, implying that the results of all FP-strict expressions must be those predicted by IEEE 754 arithmetic on operands represented using single and double formats. "

As I read this, Java calculations (with strictfp) must be exactly equal on all platforms.

strictfp will disallow higher precision results..aka cannot convert a multiply plus add to fused madd instruction.. all basic operations must be properly rounded results per op.  Java's strictness disallow most FP transforms by default and strictfp only applies to basic ops.
Logged

All code submitted by me is in the public domain. (http://unlicense.org/)
Kalles Fraktaler
Fractal Senior
******
Posts: 1458



kallesfraktaler
WWW
« Reply #113 on: July 15, 2013, 01:56:54 PM »

I am adding this comment to the appropriate thread as well:

Here is a location that might be a challenge for those who are trying to create general fractal explorer programs that are able to take care of all kinds of glitches:
Code:
Re: -1.768534924421349949270981466087798476715853694537361755308932708832468117755950400059653117782799902337478638501072448963628214932441920842176948908374018641322144530898392300564007162704979935800934871014049179313829574155128361677427135196101360338791366078063894990715467902483187609339593823325014485511210158987366776835280397149444804359644147713966607524271943068601381762700362
Im: -0.003238506313847311658390730094059770427468853795857291858449600524173786304597107579640711772650732479651543959102224503193400567765842571403732305799887983178067448569583202109990534768047599808574906753122412433706343662943659711926494919540580535659137842925115520137140901674495412916254795600653116060701559766456079182313472846633378797158212857595572543139303270454793219579761
Magnification: 8.21041081256E354
Since the reference point in the center of these parameters is not the pixel with the highest iterations, it has both blobs and diagonal lines! These diagonal lines start to occur somewhere deeper than e300 and are probably due to the precision limitations in double and, since I did not see this in the previous version of my program, it probably has something to do with series approximation too.

Further, the series approximation function in the linked pdf document contains only the terms A (for delta), B (for delta^2) and C (for delta^3).
Does anybody know what the function would be for calculating a fourth term D for delta^4?


* diagonals.jpg (135.15 KB, 640x360 - viewed 123 times.)
Logged

Want to create DEEP Mandelbrot fractals 100 times faster than the commercial programs, for FREE? One hour or one minute? Three months or one day? Try Kalles Fraktaler http://www.chillheimer.de/kallesfraktaler
http://www.facebook.com/kallesfraktaler
Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #114 on: July 15, 2013, 03:19:32 PM »

I expect that adding more terms gives quickly diminishing returns.
Logged

Roquen
Iterator
*
Posts: 180


« Reply #115 on: July 15, 2013, 07:52:54 PM »

As I've mentioned in a couple of other threads how you compute complex powers is very important to prevent precision loss....up to massive catastrophic cancelling...aka all of the result is garbage although the result can be well represented.
Logged

All code submitted by me is in the public domain. (http://unlicense.org/)
hapf
Fractal Lover
**
Posts: 219


« Reply #116 on: July 31, 2013, 10:49:35 PM »

Can blobs still happen when the series approximation is not used if the reference point is inside or near the image and it does not escape for the valid max iteration for the whole image (which looks correct with that max iteration if sufficient arbitrary precision is used)? If so, why can they still happen?
Logged
Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #117 on: August 01, 2013, 01:49:28 AM »

Blobs have nothing to do with the series approximation. Blobs result when points in the image "resonate differently" from the reference point in substantial ways, traveling far from the reference orbit while still not escaping. Nearby points in a blob end up with deltas with many "static" nonzero significant figures before the first ones that differ, and lack of precision can render them indistinguishable, resulting in the blob. A reference point that travels with them (generally one from the core of the blob itself works best) will avoid producing that blob as those many static figures in the deltas become zeros and thus don't consume any of the available mantissa precision. However, another part of the same image may become a blob with the alternate reference point, so sometimes it is necessary to use two.
Logged

Kalles Fraktaler
Fractal Senior
******
Posts: 1458



kallesfraktaler
WWW
« Reply #118 on: August 01, 2013, 09:02:56 AM »

Yes, if the reference point, even if it eventually ends up in a minibrot, have too few iterations, compared to the blob areas, within depths reachable from the double datatype, there will be blobs.
If a reference point is added in the center of one of the blobs, it usually correct the other blobs of same type/shape.

Series approximation cause other glitches though, the diagonal lines mentioned and also "surfing along curved edges" is not possible at all. I will add a movie later showing what I mean with that.
« Last Edit: August 01, 2013, 09:14:44 AM by Kalles Fraktaler » Logged

Want to create DEEP Mandelbrot fractals 100 times faster than the commercial programs, for FREE? One hour or one minute? Three months or one day? Try Kalles Fraktaler http://www.chillheimer.de/kallesfraktaler
http://www.facebook.com/kallesfraktaler
hapf
Fractal Lover
**
Posts: 219


« Reply #119 on: August 01, 2013, 10:14:50 AM »

Yes, if the reference point, even if it eventually ends up in a minibrot, have too few iterations, compared to the blob areas, within depths reachable from the double datatype, there will be blobs.
If the reference point is in a minibrot it can't have too few iterations compared to other points in the image unless the general maxiter for the image is set too low. That is a separate issue. The problem with the rounding errors I can understand, though. I have not seen such a case yet but I don't do 'very' deep zooms. Does anyone have an example for this that is not zoomed in very deep?
Logged
Pages: 1 ... 6 7 [8] 9 10 ... 17   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
Java applet for exploring the Mandelbrot set Announcements & News Paton 5 3356 Last post March 26, 2007, 06:03:34 PM
by Paton
What range/precision for fractional escape counts for Mandelbrot/Julia sets? Programming Duncan C 7 3985 Last post May 01, 2007, 08:23:13 PM
by Duncan C
Java Mandelbrot segment Help & Support fractalwizz 10 1340 Last post December 29, 2008, 08:01:24 PM
by cKleinhuis
[Java] Double-double library for 128-bit precision. Programming Zom-B 10 12449 Last post December 20, 2010, 04:03:48 AM
by David Makin
*Continued* SuperFractalThing: Arbitrary precision mandelbrot set rendering in Java. Announcements & News « 1 2 ... 23 24 » hapf 347 21649 Last post September 28, 2017, 10:20:19 PM
by claude

Powered by MySQL Powered by PHP Powered by SMF 1.1.21 | SMF © 2015, Simple Machines

Valid XHTML 1.0! Valid CSS! Dilber MC Theme by HarzeM
Page created in 0.184 seconds with 28 queries. (Pretty URLs adds 0.011s, 2q)