Logo by Pauldelbrot - 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: Check out the originating "3d Mandelbulb" thread here
 
*
Welcome, Guest. Please login or register. March 28, 2024, 01:20:39 PM


Login with username, password and session length


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


Pages: 1 [2]   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: Implemented the Mandelbrot loop in Lisp.  (Read 12935 times)
0 Members and 1 Guest are viewing this topic.
Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #15 on: July 08, 2009, 05:31:34 AM »

This code does better than that. Fed a formula in quasi-mathematical notation, it compiles and returns a function that, when called, executes a tight, optimized loop.

I've now got more smarts in it. The optimizer now factors common factors out of sums -- but only if doing so will save a multiply (factoring the x out of (x + 3x2) doesn't, either way there are three multiplies, and factoring it out loses if you want to CSE x2s later. Though not doing so may lose the ability to CSE a (1 + 3x) later, losing the chance to CSE one sum in return for the chance to CSE a product and a different sum is worth it, heuristically.

Furthermore, the grouping now prefers pairs seen elsewhere, prefers to CSE out binary products with large numbers to ones with small numbers, and prefers to CSE out nonnumeric binary products, since other tricks may speed up integer multiplies.

The above test cases, plus two additional ones, now yield:


user=> (optimize-forms-for 55 '[(:= x (* x z))] '[(:= y (+ y x))] false)
[((:= x (* x z))) ((:= y (+ y x)))]
No change.

user=> (optimize-forms-for 55 '[(:= x (* z x x))] '[(:= y (+ y x))] false)
[((:= x (* x (* z x)))) ((:= y (+ y x)))]
No change.

user=> (optimize-forms-for 55 '[(:= x (* z x x))] '[(:= y (+ y x))] true)
[((:= x (* z (* x x)))) ((:= y (+ y x)))]
No change.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b)))] '[(:= y (+ y x))] true)
[((:= x (* (* z (* x x)) (* (/ a) (/ b))))) ((:= y (+ y x)))]
No meaningful change.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b)))] '[(:= y (+ (* y q y) x))] true)
[((:= x (* (* z (* x x)) (* (/ a) (/ b))))) ((:= y (+ (* y (* y q)) x)))]
No meaningful change.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b)))] '[(:= y (+ (* x q x) x))] true)
[((:= x (* (* z (* x x)) (* (/ a) (/ b))))) ((:= y (+ (* q (* x x)) x)))]
No meaningful change. With indiscriminate sum factoring, loses the (* x x) grouping in the second expression, but with sum factoring limited to cases where it saves a multiply, it keeps the squares together for the win.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b)))] '[(:= y (+ (* x q x) (* 7 x)))] true)
[((:= x (* (* z (* x x)) (* (/ a) (/ b))))) ((:= y (* x (+ (* q x) 7))))]
Now factoring the sum does save a multiply. It still keeps the other x^2 in one piece, continuing to save a multiply there if a circular bailout is used.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b)))] '[(:= y (+ (* -1 q x) x (* -1 z) y))] true)
[((:= x (* (* z (* x x)) (* (/ a) (/ b))))) ((:= y (+ (+ (* -1 (* q x)) (* -1 z)) (+ x y))))]
No meaningful change.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b)))] '[(:= y (+ (* -1 q x) x (* -1 z) y k))] true)
[((:= x (* (* z (* x x)) (* (/ a) (/ b))))) ((:= y (+ (+ x y) (+ k (+ (* -1 (* q x)) (* -1 z))))))]
No meaningful change.

user=> (optimize-forms-for 55 '[(:= x (* z (/ a) x x (/ b) (+ x k)))] '[(:= y (+ (* -1 q x) x (* -1 z) y k))] true)
[((:= x (* (* (/ a) (/ b)) (* (* x x) (* z (+ x k)))))) ((:= y (+ (+ (* -1 (* q x)) (* -1 z)) (+ y (+ x k)))))]
No meaningful change.

user=> (optimize-forms-for 55 '[(:= x (* x y (/ (+ (* 3 y x z) 1)) z))] [] true)
[((:= x (* (/ (+ (* 3 (* x (* y z))) 1)) (* x (* y z))))) ()]
No meaningful change.

user=> (optimize-forms-for 55 '[(:= x (* x y (/ (+ (* 3 y w z) 1)) z))] [] true)
[((:= x (* (/ (+ (* 3 (* w (* y z))) 1)) (* x (* y z))))) ()]
Now groups the yz products.

user=> (optimize-forms-for 55 '[(:= x (+ x -1 (* x y (/ (+ a 1))))) (:= y (* 2 y (/ b) y (/ (+ a 1)) x))] [] true)
[((:= x (+ (* x (* y (/ (+ a 1)))) (+ x -1))) (:= y (* 2 (* (/ b) (* (* y y) (* (/ (+ a 1)) x)))))) ()]
Now groups y2 instead of xy/(a + 1) unfortunately. But ...

user=> (optimize-forms-for 55 '[(:= x (+ x -1 (* x y (/ (+ a 1))))) (:= y (* 2 y (/ b) y (/ (+ a 1)) x))] [] false)
[((:= x (+ (* x (* y (/ (+ a 1)))) (+ x -1))) (:= y (* 2 (* (/ b) (* y (* x (* y (/ (+ a 1))))))))) ()]
Square bailout version keeps that fraction together, and would probably be faster anyway with this system. And the formula compiler will test both and keep the faster one, automatically. Heh, heh, heh.

user=> (optimize-forms-for 55 '[(:= x (* 3 x z))] '[(:= y (* 3 y z))] true)
[((:= x (* x (* 3 z)))) ((:= y (* y (* 3 z))))]
Now groups the 3z products.

user=> (optimize-forms-for 55 '[(:= x (+ a (* 3 x)))] '[(:= y (+ b (* 3 x y)))] true)
[((:= x (+ a (* 3 x)))) ((:= y (+ b (* y (* 3 x)))))]
No change.

user=> (optimize-forms-for 55 '[(:= x (+ a (* 3 x)))] '[(:= y (+ b (* 6 x y)))] true)
[((:= x (+ a (* 3 x)))) ((:= y (+ b (* 2 (* y (* 3 x))))))]
No change.

user=> (optimize-forms-for 55 '[(:= x (+ a (* 3 x) (* 2 y)))] '[(:= y (+ b (* 6 x y)))] true)
[((:= x (+ (* 2 y) (+ a (* 3 x))))) ((:= y (+ b (* (* 3 x) (* 2 y)))))]
No change.

user=> (optimize-forms-for 55 (simplify-forms (real-forms *n-complex* '[(:= z (+ (** z 2) c))])) [] true)
[((:= z-r (+ (* -1 (* z-i z-i)) (+ c-r (* z-r z-r)))) (:= z-i (+ c-i (* 2 (* z-i z-r))))) ()]
Standard M-set iteration still gives usual result.

user=> (optimize-forms-for 55 (simplify-forms (real-forms *n-complex* '[(:= z (+ (** z 3) c))])) [] true)
[((:= z-r (+ (* -1 (* (* z-i z-i) (* 3 z-r))) (+ c-r (* z-r (* z-r z-r))))) (:= z-i (+ (* (* z-i z-i) (* -1 z-i)) (+ c-i (* (* z-r z-r) (* 3 z-i)))))) ()]
Cubic M-set iteration still gives useful grouping of the squares. Different grouping of numeric terms, because it now sees opportunities to CSE two occurrences each of (* -1 z-i) and (* 3 z-r) but ends up choosing to keep the squares together instead. The (* -1 (* (* z-i z-i) (* 3 z-r))) is inefficient, but will be cleaned up by the not-yet-implemented fourth pass, which will move numeric factors out of non-CSE-able terms like that non-replicated (* 3 z-r).

user=> (optimize-forms-for 55
          (simplify-forms (real-forms *n-complex* '[(:= z (+ (** z 2) c))]))
          (simplify-forms (real-forms *n-complex* '[(:= dz (+ (* 2 z dz) 1))])) true)
[((:= z-r (+ (* -1 (* z-i z-i)) (+ c-r (* z-r z-r)))) (:= z-i (+ c-i (* z-i (* 2 z-r)))))
 ((:= dz-r (+ (* -1 (* dz-i (* 2 z-i))) (+ 1 (* dz-r (* 2 z-r))))) (:= dz-i (* 2 (+ (* dz-r z-i) (* dz-i z-r)))))]
Now groups (* 2 z-r) in two places, but at the third possible site prefers to factor the 2 out of the sum instead of factor (* 2 z-i) and (*2 z-r). This saves only one multiply by 2 instead of two, unfortunately.
Logged

Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #16 on: July 08, 2009, 07:38:15 AM »


user=> (optimize-forms-for 55
          (simplify-forms (real-forms *n-complex* '[(:= z (+ (** z 2) c))]))
          (simplify-forms (real-forms *n-complex* '[(:= dz (+ (* 2 z dz) 1))])) true)
[((:= z-r (+ (* -1 (* z-i z-i)) (+ c-r (* z-r z-r)))) (:= z-i (+ c-i (* z-i (* 2 z-r)))))
 ((:= dz-r (+ (* -1 (* dz-i (* 2 z-i))) (+ 1 (* dz-r (* 2 z-r))))) (:= dz-i (* 2 (+ (* dz-r z-i) (* dz-i z-r)))))]
Now groups (* 2 z-r) in two places, but at the third possible site prefers to factor the 2 out of the sum instead of factor (* 2 z-i) and (*2 z-r). This saves only one multiply by 2 instead of two, unfortunately.


Fixed. I made it smarter. Pass 1 does sum common-factor extraction. Pass 2 gathers repeated possible binary subexpressions. Pass 3 is the one I improved further: when it sees a product, it sees if it simplifies into a sum via the distributive rule. If so, then it tries both versions, counting how many binary product subexpressions it can eliminate because they occur elsewhere in each case. It weighs the change in this number from undoing the common-factor extraction (extra multiplies) with the multiplies saved by the extraction (one fewer than the number of terms in the post-extraction sum that are not 1 -- one fewer for the single multiply by the whole sum, and ignoring sum terms equal to 1 because 1*common-factor is not a real multiply). Then it chooses. (This glosses over some major subtleties I have spent the past two hours working out.)

All of the tests except the last continue to give the same results as before, while the last yields:

[((:= z-r (+ (* -1 (* z-i z-i)) (+ c-r (* z-r z-r)))) (:= z-i (+ c-i (* z-i (* 2 z-r)))))
 ((:= dz-r (+ (* -1 (* dz-i (* 2 z-i))) (+ 1 (* dz-r (* 2 z-r))))) (:= dz-i (+ (* dz-i (* 2 z-r)) (* dz-r (* 2 z-i)))))]


Saved: two squarings (assuming circular bailout), one doubling of z-i, two doublings of z-r, and a partridge in a pear tree. Distance estimator can blaze now.

Next up: make it preferentially group expressions involving no variables appearing on the LHS of assignments, overriding all else. These can be factored out of the iteration loop entirely. Then I'll implement and test the fourth pass and run integration tests.
Logged

Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #17 on: July 08, 2009, 08:13:57 AM »

Next up: make it preferentially group expressions involving no variables appearing on the LHS of assignments, overriding all else. These can be factored out of the iteration loop entirely.

Done:

user=> (optimize-forms-for 55 '[(:= x (* z x w))] '[(:= y (* z y w))] true)
[((:= x (* x (* z w)))) ((:= y (* y (* z w))))]
user=> (optimize-forms-for 55 '[(:= x (* z x w)) (:= w z) (:= z w)] '[(:= y (+ z y w))] true)
[((:= x (* w (* z x))) (:= w z) (:= z w)) ((:= y (+ w (+ z y))))]

Note how it preferentially clumps the uniterated, per-pixel constants z and w in the first instance (green); in the second, when those have been made no longer constant (red, added assignment expressions) the preference disappears (blue and purple).

This means mappings expressed with lots of (* 2 pi) or similarly won't have that being calculated every iteration. It will get, or stay, grouped, and as a big common subexpression with only constant terms it will get factored out of the inner-loop code.

This will be particularly important for Antimatter, which has an expensive (exp (* 2 pi i a)) in its expression.
Logged

Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #18 on: July 08, 2009, 12:55:20 PM »


user=> (optimize-forms-for 55 (simplify-forms (real-forms *n-complex* '[(:= z (+ (** z 3) c))])) [] true)
[((:= z-r (+ (* -1 (* (* z-i z-i) (* 3 z-r))) (+ c-r (* z-r (* z-r z-r))))) (:= z-i (+ (* (* z-i z-i) (* -1 z-i)) (+ c-i (* (* z-r z-r) (* 3 z-i)))))) ()]
Cubic M-set iteration still gives useful grouping of the squares. Different grouping of numeric terms, because it now sees opportunities to CSE two occurrences each of (* -1 z-i) and (* 3 z-r) but ends up choosing to keep the squares together instead. The (* -1 (* (* z-i z-i) (* 3 z-r))) is inefficient, but will be cleaned up by the not-yet-implemented fourth pass, which will move numeric factors out of non-CSE-able terms like that non-replicated (* 3 z-r).


That fourth pass has now been implemented and tested fairly thoroughly, and it works. In particular, it changes the above to:

[((:= z-r (+ (* -3 (* (* z-i z-i) z-r)) (+ c-r (* z-r (* z-r z-r))))) (:= z-i (- (+ c-i (* 3 (* (* z-r z-r) z-i))) (* (* z-i z-i) z-i)))) ()]


Note that the 3 has been pulled up and combined with the -1. The 3 would have stayed put if it had been part of a repeated subexpression, but it was not.

The distance estimator becomes:

[((:= z-r (- (+ c-r (* z-r z-r)) (* z-i z-i))) (:= z-i (+ c-i (* z-i (* 2 z-r))))) ((:= dz-r (- (+ 1 (* dz-r (* 2 z-r))) (* dz-i (* 2 z-i)))) (:= dz-i (+ (* dz-i (* 2 z-r)) (* dz-r (* 2 z-i)))))]

Note that here and above some expressions like (* -1 foo) have been changed into subtraction operations in their enclosing additive expressions. Additionally, note that (* dz-i (* 2 z-i)) did not have its 2 brought out and combined with the -1 to produce a -2; instead, the factor of -1 became a subtraction. This is in contrast to the 3 and -1 combining in the previous example, and is a result of (* 2 z-i) being a recurrent subexpression. As a result of this behavior, a multiply is saved in preference to saving a subtraction. (This looks slightly bad with bignums, since the (* 2 z-i) would become a (+ z-i z-i) and it therefore appears to be saving a mere add at the cost of a more expensive subtraction in this case. But if it did hoist the 2 and make a -2 factor, this would either be a very expensive multiply, or else a subtraction and add. So the subtraction is actually unavoidable, except at the cost of a more expensive multiply, so by keeping the 2 where it is it saves a bignum add.)

The upshot of these last few posts is this: the optimizer is essentially complete. Oh, it will probably get the odd tweak now and again but the bulk of the job is done there and it's time to move on to integration testing, and then start layering functionality above the iteration loops.
Logged

Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #19 on: July 08, 2009, 03:51:01 PM »

Results of the integration tests are in:

.............................................................
...........................*#*...............................
............................*................................
.......................**#######**#..........................
.....................**############..........................
.............*.#....*##############*.........................
............*#####***##############*.........................
#*********#########*##############...........................
............*#####***##############*.........................
.............*.#....*##############*.........................
.....................**############..........................
.......................**#######**#..........................
............................*................................
...........................*#*...............................
.............................................................

That's with an iteration max of 1000 and 72-bit mantissas. Yeah, it's a very shallow deep zoom, and it's using the new expression optimizer. The optimizer's only effect on this particular mapping being to change 2xy to (xy + xy) in the next-y part of the calculations (which it does because bignum multiplies are so much slower than bignum adds). (*s indicate points whose periods were not detected.)

Still, everything seems to hang together and work together at this point. Now work can begin on building the higher-level stuff.

The code I executed at my REPL to generate that, by the way, being

(def *l* (make-iteration-loops 72 *n-complex* 2 2 '[(:= z (+ (** z 2) c))] '[] nil '(0 0) '(-1.9913844699254 0)))
(asciibrot 30 15 *l* 1000)

Yes, it really does transform '[(:= z (+ (** z 2) c))] on the fly into an optimized, for non-bignums only 7 nanosecond, loop that can be called upon at any future time, all by itself without further human intervention.

Well, almost.

The quoted lists at the end of the first line are a z and c pair that gives rise to an ergodic orbit on the spike. The benchmarking make-iteration-loops uses to choose alternative loop generation parameters didn't work too well with dummy loops -- either for code cache or branch prediction reasons, I'm not sure which. Using a live loop presents a problem if it bails out too early. So an interior point needs to be supplied, and for it to decide between circular and square trap around finite attractors/period checks, it needs not to converge too fast.

I'll experiment to see how sensitive the timing data are to the speed of escape/convergence of the supplied point. If worst comes to worst, for complex single-variable maps and quats, hypers, and any other algebraic rings that embed the complex plane, solving for the fixed points and then for parameter points such that the derivative at some fixed point is e^{\pi i (\sqrt{5} - 1)}, which is irrationally indifferent. (For the quadratic Mandelbrot, there are two such parameter points, belonging to the M-set and lying between seahorse valley and the period-3 bulbs; zooming there shows cauliflower shapes and gets into very high iterations very quickly, and the most common images of Siegel disk Julia sets correspond to these points. The orbit of zero is ergodic for all z2 + c Siegel disk parameters.) But I'm not sure what I'll do about systems of more than one variable or of one real variable in that eventuality.
Logged

Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #20 on: July 10, 2009, 03:41:58 PM »

Added another loop-generation option: fixed attractor search. The loop in this case is generated to check every iteration for approach to a particular point that is baked in, rather than a parameter to the loop function. The advantage is that the tests can be CSE'd with the rest of the calculations. If the point is the origin, and both tests are for circular rather than square targets, the same squares and sum used for bailout get recycled, for instance.

I also have derivatives and code to take an input expression and spit out an expression that, when iterated, performs Newton's method on the original one. I've used this to make a working root-finder. Combining these all yields the critical points for Mandelbrot calculations. Oh, and the derivative code can output Jacobian matrix calculations for multi-variable systems, too.

The bad news: tests have determined that the loop timings code really does need a trapped point to do its job at all well. Even adding code to calculate the function call overhead (for the call itself, and the loop function's preamble code that executes before entering the loop proper) and account for it, the results were all over the board for z = 2.0001 as a test point (worst case -- 1 iteration and it escapes).

I think I shall have to code it to just generate a normal escaping-points loop with a circular bailout, start rendering, and if it hits a trapped point, do the timings and everything and switch to the original "plan" then. With a modification in that it should identify basins whose attractors are unvarying with pixel coordinates and generate loops for those attractors to use whenever it hits those basins again.
Logged

Nahee_Enterprises
World Renowned
Fractal Senior
******
Posts: 2250


use email to contact


nahee_enterprises Nahee.Enterprises NaheeEnterprise
WWW
« Reply #21 on: January 02, 2010, 11:02:23 PM »

It's shocking that nobody much else is apparently developing
fractal software in Lisp, given what it's capable of.   smiley

The first part of the optimizer is done.  It behaves more or less as described....

Whatever happened to this ??  Was another thread started somewhere else, or was this project dropped to the wayside ??
Logged

Pauldelbrot
Fractal Senior
******
Posts: 2592



pderbyshire2
« Reply #22 on: January 03, 2010, 04:24:40 AM »

Got back-burnered as I ran into a thorny problem to solve (unrelated to iteration loop generation) and got busy with other things.

Don't worry, more will be done, and perhaps soon.
Logged

Nahee_Enterprises
World Renowned
Fractal Senior
******
Posts: 2250


use email to contact


nahee_enterprises Nahee.Enterprises NaheeEnterprise
WWW
« Reply #23 on: January 03, 2010, 04:47:16 AM »

Got back-burnered... and got busy with other things.
Don't worry, more will be done, and perhaps soon.

Thanks for the quick reply.  No rush on my part, was just curious after going through some of the postings over the past few months.  (Still catching up on several hundred postings.)     wink

Will keep an eye open for when you do get back around to it.
Logged

erstwhile
Alien
***
Posts: 21


@backtracemusic
WWW
« Reply #24 on: September 03, 2014, 08:44:22 PM »

I'm probably going to transfer my Supercollider Mandelbrot Music generator to Clojure / Overtone since Supercollider doesn't support arbitrary-precision numbers. This could be a good starting point for implementing m-set calculations in Clojure.
Logged
Pages: 1 [2]   Go Down
  Print  
 
Jump to:  

Related Topics
Subject Started by Replies Views Last post
Help with Fractal zoom loop Movies Showcase (Rate My Movie) OdedS 0 428 Last post October 01, 2008, 02:32:20 PM
by OdedS
little loop Animation taurus 0 4564 Last post June 01, 2012, 10:19:05 AM
by taurus
My old animation loop Animations Showcase (Rate My short Animation) Alef 1 1143 Last post June 07, 2015, 08:54:14 PM
by 3dickulus
One line and a loop Images Showcase (Rate My Fractal) 4rsenic 0 822 Last post November 01, 2013, 01:05:00 PM
by 4rsenic
Mandelbrot phase angle GIF loop Animations Showcase (Rate My short Animation) billtavis 4 9238 Last post February 24, 2015, 06:29:58 AM
by Chris Thomasson

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.201 seconds with 24 queries. (Pretty URLs adds 0.011s, 2q)