Maybe I should just be working through the excercises in Graham's ACL
book. However, I wanted to play around with something that I know
how to do by brute force in C. I expect that my C-think shows in the
code below. I haven't tried to optimize away some of the redundant
calculations yet. I did try the idea of memoization in
mandelbrot-escape-fn, but it so far has only served well for a few
specific input cases. I shall have to google around a bit to see if
there are decent ways of detecting when I'm orbiting some attractor
or set of attractors.
Oh, I'm talking about calculating whether or not a complex number is
in the Mandelbrot set. When I did this in C, I was using IEEE-754
doubles. I even made a movie where I zoom in just short of where I
can go with 64 bit doubles. I thought perhaps using CL's rational
numbers, I could go a lot further.
The first problem I ran into was that squaring fractions generates a
lot of digits very quickly. If you want a way to consume CPU and
memory resources, I recomend it. I expect there are fixed precision
libraries available that I should look at. In the meantime, I hacked
together a function that does the rational equivalent of rounding
numbers to some settable precision.
Anyway, even though this is code in progress, or to be more accurate,
I am exploring code options, I am surprised at how much time is spent
in GC. It looks to be about 10% of the total time. How typical is
that for a Lisp program? I would expect that good Lisp programmers
know how to keep GC time down to a minimum.
As for the huge amount of memory I allocate, I think that is mostly
going into the hash table. I'm going to be backing that thing out
and seeing what my numbers are without it. I don't think I have the
memory or address space to handle millions of iterations in the
mandelbrot-escape-fn function.
Another question on technique: When do you decide whether to use a
macro rather than a function call? I'm already thinking I might
benefit from some optimization using macros instead of functions with
what I have so far, but I just don't know yet. I imagine macros are
also something that can be abused fairly easily.
I also would hope the compiler warns me about variable capture. So
far, every warning I've seen from OpenMCL has been because I made a
coding error. Even with show-paren-mode turned on in Emacs, my lack
of familiarity with Lisp still has me put symbols in the wrong
place. It is also easy to forget when I need a list of lists rather
than just a list for a function.
One thing I am happy about is that even though I find I have to test
my code little bits at a time, Lisp makes it much easier for me to do
that than even Perl. Having the top level right there is very handy
as is the ability to do C-x C-e in SLIME. I did manage to crash
Emacs though when I did C-g a couple times to stop a lengthy
computation. This prompted me to update SLIME from CVS. I haven't
updated Emacs again though. Maybe I should.
OK, I'm rambling on a bit too much now. Here is the code I am
talking about. Perhaps it will show where my mind is right now.
;;; clamp-rational -- reduce the number of digits ("precision") of a rational number by dividing by one
(defun clamp-rational (r &key (clamp most-positive-fixnum))
(if (rationalp r)
(if (> (denominator r) clamp)
(let ((scale-divisor (/ (denominator r) clamp)))
(/ (round (/ (numerator r) scale-divisor)) (round (/ (denominator r) scale-divisor))))
r)
r))
;;; clamp-complex -- clamp-rational the parts of a complex number
(defun clamp-complex (c &key (clamp most-positive-fixnum))
(complex (clamp-rational (realpart c) :clamp clamp)
(clamp-rational (imagpart c) :clamp clamp)))
;;; cplx-num-escaped -- test that z has not escaped
(defun cplx-num-escaped (z)
(>= (+ (* (realpart z) (realpart z)) (* (imagpart z) (imagpart z))) 4))
;;; mandelbrot-escape-fn -- (values escaped iterations)
(defun mandelbrot-escape-fn (c &key (max-iterations (- most-positive-fixnum 1)))
(do ((z #C(0 0) (clamp-complex (+ (* z z) c)))
(iterations 0 (+ iterations 1))
(escaped nil (cplx-num-escaped z))
(ht (make-hash-table :test #'equal)))
((or (= iterations max-iterations)
escaped
(gethash z ht))
(values escaped iterations))
(setf (gethash z ht) t)))
(time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
(MANDELBROT-ESCAPE-FN #c(-1 1/4) :MAX-ITERATIONS 100000) took 23,912 milliseconds (23.912 seconds) to run.
Of that, 17,770 milliseconds (17.770 seconds) were spent in user mode
5,170 milliseconds (5.170 seconds) were spent in system mode
972 milliseconds (0.972 seconds) were spent executing other OS processes.
2,345 milliseconds (2.345 seconds) was spent in GC.
106,854,528 bytes of memory allocated.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> (defun clamp-complex (c &key (clamp most-positive-fixnum))
> (complex (clamp-rational (realpart c) :clamp clamp)
> (clamp-rational (imagpart c) :clamp clamp)))
This can be a handy way to pass unchanged keyword arguments through to
another function:
(defun clamp-complex (c &rest keys)
(complex (apply #'clamp-rational (realpart c) keys)
(apply #'clamp-rational (imagpart c) keys)))
--
Lars Brinkhoff, Services for Unix, Linux, GCC, HTTP
Brinkhoff Consulting http://www.brinkhoff.se/
Lars Brinkhoff <·········@nocrew.org> writes:
> David Steuber <·············@verizon.net> writes:
> > (defun clamp-complex (c &key (clamp most-positive-fixnum))
> > (complex (clamp-rational (realpart c) :clamp clamp)
> > (clamp-rational (imagpart c) :clamp clamp)))
>
> This can be a handy way to pass unchanged keyword arguments through to
> another function:
>
> (defun clamp-complex (c &rest keys)
> (complex (apply #'clamp-rational (realpart c) keys)
> (apply #'clamp-rational (imagpart c) keys)))
Thanks! I snagged this and substituted it in. This seems like a
tremendous way to improve code maintainability.
I think Frode is also correct about my use of the hash table. I
thought eql was just a pointer comparator which is why I chose equal
instead of the default eql. I'll give that a shot. I want only the
strictest comparison of equality that I can get away with. I might
also want to create my own function that basicly says, "close enough".
Although I may not need to do that. I had a minor burst of
inspiration and decided to actually *use* the :clamp key:
(defun mandelbrot-escape-fn (c &key (max-iterations (- most-positive-fixnum 1)))
(do ((z #C(0 0) (clamp-complex (+ (* z z) c) :clamp 10000))
(iterations 0 (+ iterations 1))
(escaped nil (cplx-num-escaped z))
(ht (make-hash-table :test #'equal)))
((or (= iterations max-iterations)
escaped
(gethash z ht))
(values escaped iterations))
; (print z)
(setf (gethash z ht) t)))
This made a huge difference:
(time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
(MANDELBROT-ESCAPE-FN #c(-1 1/4) :MAX-ITERATIONS 100000) took 46 milliseconds (0.046 seconds) to run.
Of that, 20 milliseconds (0.020 seconds) were spent in user mode
0 milliseconds (0.000 seconds) were spent in system mode
26 milliseconds (0.026 seconds) were spent executing other OS processes.
181,232 bytes of memory allocated.
NIL
499
Look, ma! No GC!
The NIL return value means it did not escape, that is correct. 499
is the number of iterations it took to figure that out. That sure
beats going the full 100000. So I just have to figure out how small
I can make the clamp without getting a false positive.
I should point out that I learned two things doing this with Lisp
that I did not learn doing this in C:
* Going with strictly arbitrary precision arithmetic can get me into
real trouble fast.
* I can detect an attractor providing I can use a 'near enough'
method to detect it in the first place.
The actual edge of the Mandelbrot Set is supposed to be incalculable
because the number of attractors rises to infinity.
I am liking Lisp more and more. I really am beginning to wonder why
the heck it is not a so called mainstream language. It seems foolish
to ignore Lisp as a tool for exploration.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber wrote:
> I am liking Lisp more and more. I really am beginning to wonder why
> the heck it is not a so called mainstream language. It seems foolish
> to ignore Lisp as a tool for exploration.
After you fill in your newby survey (start search at second link in
sig), it sounds as if you should add yourself to the "What Do I Know
That They Do Not Know?" subcategory of the highlight film.
:)
kt
--
http://tilton-technology.com
Why Lisp? http://alu.cliki.net/RtL%20Highlight%20Film
Your Project Here! http://alu.cliki.net/Industry%20Application
David Steuber <·············@verizon.net> wrote in message news:<··············@david-steuber.com>...
> Lars Brinkhoff <·········@nocrew.org> writes:
> > This can be a handy way to pass unchanged keyword arguments through to
> > another function:
> >
> > (defun clamp-complex (c &rest keys)
> > (complex (apply #'clamp-rational (realpart c) keys)
> > (apply #'clamp-rational (imagpart c) keys)))
>
> Thanks! I snagged this and substituted it in. This seems like a
> tremendous way to improve code maintainability.
>
A slightly better (I think) way of doing this is the following
slightly obscure idiom:
(defun clamp-complex (c &rest keys &key &allow-other-keys)
... (apply ... keys) ...)
This allows CLAMP-COMPLEX to validate that all its arguments are kw
args (so they're in pairs) without caring what they are. If you want,
you can now actually specify what args are allowed, using something
like:
(defun clamp-complex (x &rest args &key a b c)
(declare (ignore a b c))
...)
Unfortunately if you do this you don't get to provide default values.
Defaulting and passing down keyword args correctly is fairly fiddly to
get right while remaining maintainable.
··········@tfeb.org (Tim Bradshaw) writes:
> [..]
> (defun clamp-complex (x &rest args &key a b c)
> (declare (ignore a b c))
> ...)
>
> Unfortunately if you do this you don't get to provide default
> values. Defaulting and passing down keyword args correctly is
> fairly fiddly to get right while remaining maintainable.
Another thing to consider is that when you find you have variables
that you need to pass around verbatim like this, a special variable
might be more appropriate. That way, the clamping value (or whatever)
is implicit, and won't bother code that doesn't really need to care
about it.
--
Frode Vatvedt Fjeld
David Steuber <·············@verizon.net> writes:
>> (defun clamp-complex (c &key (clamp most-positive-fixnum))
>> (complex (clamp-rational (realpart c) :clamp clamp)
>> (clamp-rational (imagpart c) :clamp clamp)))
Lars Brinkhoff <·········@nocrew.org> writes:
> (defun clamp-complex (c &rest keys)
> (complex (apply #'clamp-rational (realpart c) keys)
> (apply #'clamp-rational (imagpart c) keys)))
This technique can sometimes be appropriate, but I think in this case
the true argument-list to clamp-complex is obscured, something I'd
consider a serious problem. I like to be able to tell a function's
argument-list by a key-press, and for this list to provide concise
information that is semantically meaningful. Also, I'd not be
surprised if there's a performance penalty to be paid with this
technique.
Another thing is the rather peculiar use of most-positive-fixnum,
which means that the behavior of this calculation will change with
implementations for no (to me) apparent reason.
--
Frode Vatvedt Fjeld
Frode Vatvedt Fjeld <······@cs.uit.no> writes:
> Another thing is the rather peculiar use of most-positive-fixnum,
> which means that the behavior of this calculation will change with
> implementations for no (to me) apparent reason.
The fixnum type is to be used when efficiency is desired. Of course
the user of anything that is dependent on the value of
most-positive-fixnum must ensure that the size of the limitation is
not critical to the output of their algorithm, but that doesn't I
think make it peculiar. When would you use most-positive-fixnum?
Christophe
--
http://www-jcsu.jesus.cam.ac.uk/~csr21/ +44 1223 510 299/+44 7729 383 757
(set-pprint-dispatch 'number (lambda (s o) (declare (special b)) (format s b)))
(defvar b "~&Just another Lisp hacker~%") (pprint #36rJesusCollegeCambridge)
David Steuber <·············@verizon.net> writes:
> (defun clamp-rational (r &key (clamp most-positive-fixnum))
> (if (rationalp r)
> (if (> (denominator r) clamp)
> (let ((scale-divisor (/ (denominator r) clamp)))
> (/ (round (/ (numerator r) scale-divisor)) (round (/ (denominator r) scale-divisor))))
> r)
> r))
Frode Vatvedt Fjeld <······@cs.uit.no> writes:
>> Another thing is the rather peculiar use of most-positive-fixnum,
>> which means that the behavior of this calculation will change with
>> implementations for no (to me) apparent reason.
Christophe Rhodes <·····@cam.ac.uk> writes:
> The fixnum type is to be used when efficiency is desired. Of course
> the user of anything that is dependent on the value of
> most-positive-fixnum must ensure that the size of the limitation is
> not critical to the output of their algorithm, but that doesn't I
> think make it peculiar.
What I initially reacted to was the use of most-positive-fixnum as one
would use MAXINT etc. in C. Looking a bit closer at the code above, I
agree that this use of most-positive-fixnum makes some sense, although
I'd put the clamp somewhat lower than than the fixnum maximum (are
there no negatives?) on the grounds that when you've created your
first non-fixnum ratio it's already too late (you've had to muck about
with bignums), and clamping it to just within the fixnum range again
is only likely to yield another bignum-ratio after the next
iteration. As the case often is, the profiler is your friend when it
comes to knowing whether bignums are generated.
Personally, I would have chosen a true constant as the default value
in this case, leaving the accuracy/efficiency trade-off to the caller
rather than the callee. But this consideration only applies to
general-purpose library API design, not my-quick-mandelbrot-generator,
of course.
Btw (for David Steuber), if speed really is of great interest, I would
also avoid using &key in such (presumably) often-called functions. And
the round function takes a divisor argument, too. (And, while I'm not
that good with number-crunching, isn't that newly-generated
denominator always going to end up the same as clamp?)
> When would you use most-positive-fixnum?
When writing system code (i.e. parts of an implementation), or when
optimizing something specifically for a certain implementation. Or as
a user-level descision. For example in this instance, I'd consider
this appropriate (if it was a truly general library):
(defvar *default-clamp* 10000)
(defun clamp-rational (r &key (clamp *default-clamp*))
..)
And then, in my ~/.mandel-clamper.lisp or something like that:
;; Try to balance efficiency and accuracy
(setf *default-clamp* (ash most-positive-fixnum -8))
--
Frode Vatvedt Fjeld
Frode Vatvedt Fjeld <······@cs.uit.no> writes:
> David Steuber <·············@verizon.net> writes:
>
> > (defun clamp-rational (r &key (clamp most-positive-fixnum))
> > (if (rationalp r)
> > (if (> (denominator r) clamp)
> > (let ((scale-divisor (/ (denominator r) clamp)))
> > (/ (round (/ (numerator r) scale-divisor)) (round (/ (denominator r) scale-divisor))))
> > r)
> > r))
>
> Btw (for David Steuber), if speed really is of great interest, I would
> also avoid using &key in such (presumably) often-called functions. And
> the round function takes a divisor argument, too. (And, while I'm not
> that good with number-crunching, isn't that newly-generated
> denominator always going to end up the same as clamp?)
I'll have to check, but I think you're right on all counts. I expect
I shall always pass in the :clamp parameter so I might as well make it
required and do away with the default entirely.
Actually, after that last division, the rational may be reduced.
Unless it is an implementation feature, the denominator is always a
positive number and the rational is returned in reduced form. I
should probably look at the CLHS about that.
One thing does seem to be true. A newbie at Lisp can easily write
inefficient code.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Frode Vatvedt Fjeld <······@cs.uit.no> writes:
>
> > David Steuber <·············@verizon.net> writes:
> >
> > > (defun clamp-rational (r &key (clamp most-positive-fixnum))
> > > (if (rationalp r)
> > > (if (> (denominator r) clamp)
> > > (let ((scale-divisor (/ (denominator r) clamp)))
> > > (/ (round (/ (numerator r) scale-divisor)) (round (/ (denominator r) scale-divisor))))
> > > r)
> > > r))
> >
> > Btw (for David Steuber), if speed really is of great interest, I would
> > also avoid using &key in such (presumably) often-called functions. And
> > the round function takes a divisor argument, too. (And, while I'm not
> > that good with number-crunching, isn't that newly-generated
> > denominator always going to end up the same as clamp?)
>
> I'll have to check, but I think you're right on all counts. I expect
> I shall always pass in the :clamp parameter so I might as well make it
> required and do away with the default entirely.
Trust me to forget how to do simple arithmetic. This new
clamp-rational function returns the same result:
(defun clamp-rational (r clamp)
(if (rationalp r)
(if (> (denominator r) clamp)
(let ((scale-divisor (/ (denominator r) clamp)))
(/ (round (numerator r) scale-divisor) clamp))
r)
r))
I've also gone ahead and made the clamp required.
So the returned rational should be the closest rational with clamp as
the denominator. There may be cases where the rational is further
reduced but I have no idea either way.
I wonder if clamp is the correct verb to use. It makes sense to me,
but I may be misunderstanding the meaning of clamp.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
Frode Vatvedt Fjeld <······@cs.uit.no> writes:
> David Steuber <·············@verizon.net> writes:
>
> >> (defun clamp-complex (c &key (clamp most-positive-fixnum))
> >> (complex (clamp-rational (realpart c) :clamp clamp)
> >> (clamp-rational (imagpart c) :clamp clamp)))
>
> Lars Brinkhoff <·········@nocrew.org> writes:
>
> > (defun clamp-complex (c &rest keys)
> > (complex (apply #'clamp-rational (realpart c) keys)
> > (apply #'clamp-rational (imagpart c) keys)))
>
>
> This technique can sometimes be appropriate, but I think in this case
> the true argument-list to clamp-complex is obscured, something I'd
> consider a serious problem. I like to be able to tell a function's
> argument-list by a key-press, and for this list to provide concise
> information that is semantically meaningful. Also, I'd not be
> surprised if there's a performance penalty to be paid with this
> technique.
I see your point about the argument list. I didn't think of that. I
do like the feature I have seen in SLIME of showing function/macro
arguments in the minibuffer. I also don't want to incur a
performance penalty to boot, but I don't know if that would really be
the case. I guess that depends on how clever the compiler is about
code generation.
> Another thing is the rather peculiar use of most-positive-fixnum,
> which means that the behavior of this calculation will change with
> implementations for no (to me) apparent reason.
I wanted to have a default value for the key and that seemed like the
way to go. Normally I don't like to have "magic" numbers burried in
my code so I picked this ANSI defined symbol as what seemed like a
reasonable default. It didn't actually occur to me that the value
would change with different implementations on a given platform.
I guess it is something I need to be aware of when I pick default
values. In this case, I'm not sure that a guaranteed value makes
much difference. Some decisions are arbitrary.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Another question on technique: When do you decide whether to use a
> macro rather than a function call? I'm already thinking I might
> benefit from some optimization using macros instead of functions
> with what I have so far, but I just don't know yet. I imagine
> macros are also something that can be abused fairly easily.
Only use macros when there is some substantial _syntactic_ advantages
to be gained. Using macros as "a speedy function" is only for C
programming etc. If you are in cycle-counting mode, you can declare
functions to be inlinable. Or, in extra-special cases, you can write
compiler-macros for the function, which is a way to customize how the
compiler inlines the function.
> ;;; mandelbrot-escape-fn -- (values escaped iterations)
> (defun mandelbrot-escape-fn (c &key (max-iterations
> (- most-positive-fixnum 1)))
> (do ((z #C(0 0) (clamp-complex (+ (* z z) c)))
> (iterations 0 (+ iterations 1))
> (escaped nil (cplx-num-escaped z))
> (ht (make-hash-table :test #'equal)))
> ((or (= iterations max-iterations)
> escaped
> (gethash z ht))
> (values escaped iterations))
> (setf (gethash z ht) t)))
I don't understand too much of this code, but I suspect you're not
using the hash-table optimally here. You're using numbers as keys,
suggesting to me you should use eql as the test (for numbers, eql is
equivalent to =). And you might want to experiment with the size
parameters to make-hash-table, because having to do many re-hashes can
(very likely) be expensive both in time and space. For example, you
might want to derive the hash-table's initial size from the value of
iterations.
And you might want to know about the functions 1+ and 1-.
--
Frode Vatvedt Fjeld
David Steuber <·············@verizon.net> writes:
> The first problem I ran into was that squaring fractions generates a
> lot of digits very quickly. If you want a way to consume CPU and
> memory resources, I recomend it. I expect there are fixed precision
> libraries available that I should look at. In the meantime, I hacked
> together a function that does the rational equivalent of rounding
> numbers to some settable precision.
> (time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
If you don't want rational, don't use rationals!
(time (mandelbrot-escape-fn #C(-1.0 0.25) :max-iterations 100000))
[8]> (defun mandelbrot-escape-fn (c &key (max-iterations (- most-positive-fixnum 1)))
(do ((z #C(0 0) (clamp-complex (+ (* z z) c)))
(iterations 0 (+ iterations 1))
(escaped nil (cplx-num-escaped z))
(ht (make-hash-table :test #'equal)))
((or (= iterations max-iterations)
escaped
(gethash z ht))
(values escaped iterations))
(setf (gethash z ht) t)))
MANDELBROT-ESCAPE-FN
[9]> (time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
Real time: 45.138042 sec.
Run time: 15.74 sec.
Space: 102872896 Bytes
GC: 46, GC time: 6.06 sec.
NIL ;
100000
[10]> (defun mandelbrot-escape-fn (c &key (max-iterations (- most-positive-fixnum 1)))
"Returns: (values escaped iterations)"
(do ((z #C(0.0 0.0) (+ (* z z) c))
(iterations 0 (+ iterations 1))
(escaped nil (cplx-num-escaped z))
(ht (make-hash-table :test #'equal)))
((or (= iterations max-iterations)
escaped
(gethash z ht))
(values escaped iterations))
(setf (gethash z ht) t)))
MANDELBROT-ESCAPE-FN
[11]> (time (mandelbrot-escape-fn #C(-1.0 0.25) :max-iterations 100000))
Real time: 8.486377 sec.
Run time: 3.17 sec.
Space: 20148724 Bytes
GC: 12, GC time: 1.35 sec.
NIL ;
89135
[12]>
--
__Pascal_Bourguignon__ http://www.informatimago.com/
There is no worse tyranny than to force a man to pay for what he doesn't
want merely because you think it would be good for him.--Robert Heinlein
http://www.theadvocates.org/
Pascal Bourguignon wrote:
> David Steuber <·············@verizon.net> writes:
>
>>The first problem I ran into was that squaring fractions generates a
>>lot of digits very quickly. If you want a way to consume CPU and
>>memory resources, I recomend it. I expect there are fixed precision
>>libraries available that I should look at. In the meantime, I hacked
>>together a function that does the rational equivalent of rounding
>>numbers to some settable precision.
>
>
>>(time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
>
>
> If you don't want rational, don't use rationals!
I think he does want rationals: "I even made a movie where I zoom in
just short of where I can go with 64 bit doubles. I thought perhaps
using CL's rational numbers, I could go a lot further."
kt
--
http://tilton-technology.com
Why Lisp? http://alu.cliki.net/RtL%20Highlight%20Film
Your Project Here! http://alu.cliki.net/Industry%20Application
Kenny Tilton <·······@nyc.rr.com> writes:
> Pascal Bourguignon wrote:
>
> > David Steuber <·············@verizon.net> writes:
> >
> >>The first problem I ran into was that squaring fractions generates a
> >>lot of digits very quickly. If you want a way to consume CPU and
> >>memory resources, I recomend it. I expect there are fixed precision
> >>libraries available that I should look at. In the meantime, I hacked
> >>together a function that does the rational equivalent of rounding
> >>numbers to some settable precision.
> >
> >>(time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
> > If you don't want rational, don't use rationals!
>
> I think he does want rationals: "I even made a movie where I zoom in
> just short of where I can go with 64 bit doubles. I thought perhaps
> using CL's rational numbers, I could go a lot further."
Oops, sorry, I overlooked that last sentence.
Clamping only makes things worse: it creates two more rationals in
addition to the two rationals created in each iteration.
In addition, I don't think that clamping to fixnums will improve the
precision vs. double-floats. It's ok for 1/3, but n/1000000000 will
not have more precision than n.0e9.
The best would be to use higher precision floating
points. Unfortunately, they are available only in clisp:
#+clisp
(progn
(setf CUSTOM:*FLOATING-POINT-CONTAGION-ANSI* t)
(setf CUSTOM:*WARN-ON-FLOATING-POINT-CONTAGION* nil)
(SETF (EXT:LONG-FLOAT-DIGITS) 512) ;; 512-bit long-float
(defparameter c-zero #c(0L0 0L0))
(defmacro clamp (x) x))
#-clisp
(progn
(defparameter c-zero #c(0 0))
(defmacro clamp (x) `(clamp-complex x)))
(defun cplx-num-escaped (z)
"test that z has not escaped"
(>= (+ (* (realpart z) (realpart z)) (* (imagpart z) (imagpart z))) 4))
(defun mandelbrot-escape-fn (c &key (max-iterations (- most-positive-fixnum 1)))
"Returns: (values escaped iterations)"
(do ((z c-zero (clamp (+ (* z z) c)))
(iterations 0 (+ iterations 1))
(escaped nil (cplx-num-escaped z))
(ht (make-hash-table :test #'eql)))
((or (= iterations max-iterations)
escaped
(gethash z ht))
(values escaped iterations))
(setf (gethash z ht) t)));;mandelbrot-escape-fn
[22]> (time (mandelbrot-escape-fn #C(-1.0 0.25) :max-iterations 100000))
Real time: 10.429873 sec.
Run time: 7.98 sec.
Space: 117894952 Bytes
GC: 32, GC time: 4.35 sec.
NIL ;
100000
There, you get a run time intermediary between double-float (64-bit)
and rationals, but we still cons about the same as with rationals.
--
__Pascal_Bourguignon__ http://www.informatimago.com/
There is no worse tyranny than to force a man to pay for what he doesn't
want merely because you think it would be good for him.--Robert Heinlein
http://www.theadvocates.org/
Pascal Bourguignon wrote:
> Kenny Tilton <·······@nyc.rr.com> writes:
>
>
>>Pascal Bourguignon wrote:
>>
>>>If you don't want rational, don't use rationals!
>>
>>I think he does want rationals: "I even made a movie where I zoom in
>>just short of where I can go with 64 bit doubles. I thought perhaps
>>using CL's rational numbers, I could go a lot further."
>
>
> Oops, sorry, I overlooked that last sentence.
No problemo, it was easy to miss. But now we have to help with the
rationals optimization. I am just a humble applications guy who is busy
having his ass kicked by the continuing "Mystery of GL-Ortho and the
Z-Axis", but I have done my part by changing the subject such that
everyone will dive in to optimize his code beyond belief.
The crew was all over a C++ astronomy benchmark like ugly on an ape when
it was offered as an insult, but now that a Lisp newbie enthusiast has a
genuine, interesting problem ("unlimited fractal zoom") you can hear a
pin drop.
Hmmm, should I have made it "A challenge to Pascal Bourguignon!"?
:)
kenneth
--
http://tilton-technology.com
Why Lisp? http://alu.cliki.net/RtL%20Highlight%20Film
Your Project Here! http://alu.cliki.net/Industry%20Application
Kenny Tilton <·······@nyc.rr.com> writes:
> Pascal Bourguignon wrote:
> > Kenny Tilton <·······@nyc.rr.com> writes:
> >
> >>Pascal Bourguignon wrote:
> >>
> >>>If you don't want rational, don't use rationals!
> >>
> >>I think he does want rationals: "I even made a movie where I zoom in
> >>just short of where I can go with 64 bit doubles. I thought perhaps
> >>using CL's rational numbers, I could go a lot further."
> > Oops, sorry, I overlooked that last sentence.
>
> No problemo, it was easy to miss. But now we have to help with the
> rationals optimization. I am just a humble applications guy who is
> busy having his ass kicked by the continuing "Mystery of GL-Ortho and
> the Z-Axis", but I have done my part by changing the subject such that
> everyone will dive in to optimize his code beyond belief.
>
> The crew was all over a C++ astronomy benchmark like ugly on an ape
> when it was offered as an insult, but now that a Lisp newbie
> enthusiast has a genuine, interesting problem ("unlimited fractal
> zoom") you can hear a pin drop.
>
> Hmmm, should I have made it "A challenge to Pascal Bourguignon!"?
LOL
There's a thought.
There are basicly two problems to sort out in my code so far. One is
how to reduce the number of actual numerical operations for a single
iteration of the do loop. That includes the comparison to see if Z
has escaped. The other is to find a way to short circuit the loop
(this is what the hash table is for) without generating a false
positive. That is, I don't want to say Z never escapes when it does.
As for the precision, I chose rationals because they are there and
offer unlimited precision. 1/3 does not suffer from rounding error.
Nor does 1/10 which also has problems in binary representation. A
good alternative would be a fixed radix point representation that
lets me choose how big the number can get (values over 4 are out of
the set) and how many digits are in the mantissa. I once did this in
80386 assembler using 32bit ints. The trick was simply to keep the
high 32 bits of an int multiply rather than the low 32 bits that C
gives you. However, I didn't take the extra leap to also do numbers
that consisted of a vector of ints. I doubt I would have coded
anything as good as what Lisp has for bignums.
The clamp that I am using is intended to be expanded with the zoom.
Pascal is correct that I have limited the effective precision. That
was the entire point. The problem I had was that doing Z = Z^2 + C
with complex rational numbers gave me an explosion of digits. The
rationals in reduced form were consuming all my memory in fewer than
100 iterations. The clamp does introduce inexactitude but it also
makes the numbers computable in finite time and space. The clamp can
be a bignum though. I can make it 10^23 or 10^100. That leads to
another problem of course.
The hash table takes up space. At some point, it may take up too
much space and I will have to do without it. So there is some
experimentation to be done there.
I'm sure I'm not doing anything that hasn't already been done at
IBM's Watson Research Center (where Mandelbrot work[s|ed]). For me
this is an excercise in learning how to use Lisp to solve a real
problem within the limits of finite computing power. I think it is
educational and giving me some appreciation for the language beyond
what Graham and other Lisp evangelists say. You can read about how
wonderful a Ferrari is to drive. But until you get behind the wheel
of one, it's just another car.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
The problem is that the lisp compilers we use don't "put the
intermediary results on the stack". All big number is on the heap and
left for the garbage collector to deal with.
The more operations you compute on these rational or big numbers, the
more garbage you generate. The big numbers are not "interned". For
example, running this on clisp:
(dotimes (i 10 ) (time (prog1 (values) (ext:! 20000))))
;; ext:! == factorial
shows you that you regenerate the 20000-1 intermediary numbers each time.
And long-floats in clisp are managed in the same way.
One problem may be that you don't have _mutable_ numbers. More
exactly, numbers can't be mutable, but you cannot store them in a data
block of a special type that you could overwrite with a new number
value. Well, you could do something like:
(defun store (bv num)
(declare (type (array bit (*)) bv)
(type (integer 0) num))
(dotimes (i (integer-length num))
(setf (aref bv i) (logbitp num i))))
and vice-versa, but you cannot do: (add-into bv1 bv2)
instead of: (setf n1 (+ n1 n2))
Well, yes, you can do it in Lisp:
(defun add-into (bv1 bv2)
(do ((carry 0)
(sum 0)
(i 0 (1+ i)))
((>= i (length bv1)) bv1)
(setf sum (+ carry (aref bv1 i) (aref bv2 i)))
(if (< sum 2)
(setf carry 0
(aref bv1 i) sum)
(setf carry 1
(aref bv1 i) (- sum 2)))))
but it's not as fast as a primitive +. The best you could do is
working with arrays of (unsigned-byte +word-size+).
Conclusion:
You'd have to write your own big number primitives to avoid consing.
David Steuber <·············@verizon.net> writes:
> The hash table takes up space. At some point, it may take up too
> much space and I will have to do without it. So there is some
> experimentation to be done there.
In the case of long-floats of 512 bits, the hash table represents only
about 7% of the space used. The question here is that I'm not sure
that this function loops so often that it's worthwhile to check
it. Does it really? In any case, you have already limited the number
of iterations...
[Note how sbcl is slower and conses more than clisp:
* (time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
Evaluation took:
26.095 seconds of real time
20.49 seconds of user run time
0.39 seconds of system run time
1 page fault and
311022904 bytes consed.
NIL
100000
[9]> (time (mandelbrot-escape-fn #C(-1 1/4) :max-iterations 100000))
Real time: 10.850633 sec.
Run time: 8.37 sec.
Space: 124799164 Bytes
GC: 146, GC time: 5.85 sec.
NIL ;
100000
]
--
__Pascal_Bourguignon__ http://www.informatimago.com/
There is no worse tyranny than to force a man to pay for what he doesn't
want merely because you think it would be good for him.--Robert Heinlein
http://www.theadvocates.org/
Pascal Bourguignon <····@thalassa.informatimago.com> writes:
> The problem is that the lisp compilers we use don't "put the
> intermediary results on the stack". All big number is on the heap and
> left for the garbage collector to deal with.
>
> The more operations you compute on these rational or big numbers, the
> more garbage you generate.
Yes, this is a problem. The way the code is now, the GC will be kept
rather busy when mandelbrot-escape-fn is called for every pixel in a
bitmap. Ouch.
> Conclusion:
> You'd have to write your own big number primitives to avoid consing.
If this is what I have to do, then I will do it. I would have to do
the same thing in C (or assembler). Hopefully it is easier to do in
Lisp.
> David Steuber <·············@verizon.net> writes:
> > The hash table takes up space. At some point, it may take up too
> > much space and I will have to do without it. So there is some
> > experimentation to be done there.
>
>
> In the case of long-floats of 512 bits, the hash table represents only
> about 7% of the space used. The question here is that I'm not sure
> that this function loops so often that it's worthwhile to check
> it. Does it really? In any case, you have already limited the number
> of iterations...
It does loop often. If a number is in the set, the loop will go
until it hits max-iterations. With a suitable value of :clamp (I
used 10000 in a test case), the loop can bail earlier because the
same number came up again and I've detected a cycle. In the test
that I posted earlier, I got out in 499 iterations instead of the
full 100000. That makes a big difference in time. What I don't know
yet is how the integrity of the function is affected. If I make the
clamp too small, then I could bail early when I shouldn't. Also
there may be no cycle for a given point.
I may need to have code that abandons the hash table after it reaches
a certain size.
The function
Z = Z^2 + C
has rather wild behavior. The closer to the edge of the fractal you
get, the more iterations are required to tell if you are in or out.
When I made this movie:
http://www.david-steuber.com/~david/fractal-movie/
I didn't bother with a bail out. I did the computations with a C
program and chose the max-iterations value based on the zoom factor.
As a point in the set will loop infinitely if there is no
max-iterations value, the challenge becomes one of where to draw the
line. Generating the frame data for that movie took about three days
of computation with both CPUs in a dual Athlon box dedicated to the
task. The last frame had the max-iteration value set to 10,000,000.
The last frame is magnified by 40,000,000,000,000 times. A lot of
numbers were crunched.
I want to go much, much further.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber wrote:
> Pascal Bourguignon <····@thalassa.informatimago.com> writes:
>
>
>>The problem is that the lisp compilers we use don't "put the
>>intermediary results on the stack". All big number is on the heap and
>>left for the garbage collector to deal with.
>>
>>The more operations you compute on these rational or big numbers, the
>>more garbage you generate.
>
>
> Yes, this is a problem. The way the code is now, the GC will be kept
> rather busy when mandelbrot-escape-fn is called for every pixel in a
> bitmap. Ouch.
If your ephemeral GC is working right it shouldn't cost much at all.
paul
Paul Wallich <··@panix.com> writes:
> David Steuber wrote:
>
> > Pascal Bourguignon <····@thalassa.informatimago.com> writes:
> >
> >
> >>The problem is that the lisp compilers we use don't "put the
> >>intermediary results on the stack". All big number is on the heap and
> >>left for the garbage collector to deal with.
> >>
> >>The more operations you compute on these rational or big numbers, the
> >>more garbage you generate.
> >
> >
> > Yes, this is a problem. The way the code is now, the GC will be kept
> > rather busy when mandelbrot-escape-fn is called for every pixel in a
> > bitmap. Ouch.
>
> If your ephemeral GC is working right it shouldn't cost much at all.
More specifically, the cost of an ephemeral GC run is the cost of
scanning the root set, tracing the live objects in the nursery
generation, and copying those live objects to the next generation.
The amount of garbage only triggers the GC, it doesn't cost any time
once you've done that. It's not garbage that's expensive, it's live
objects.
So, in your case, I'd try tuning the GC. What you want is to have a
longish period of time between runs, and to have the GC run when
you're not referencing much data that you're going to throw away. Try
setting the egc threshold to something like 20 MB.
--
/|_ .-----------------------.
,' .\ / | No to Imperialist war |
,--' _,' | Wage class war! |
/ / `-----------------------'
( -. |
| ) |
(`-. '--.)
`. )----'
In article <···············@famine.OCF.Berkeley.EDU>,
Thomas F. Burdick <···@famine.OCF.Berkeley.EDU> wrote:
> More specifically, the cost of an ephemeral GC run is the cost of
> scanning the root set, tracing the live objects in the nursery
> generation, and copying those live objects to the next generation.
> The amount of garbage only triggers the GC, it doesn't cost any time
> once you've done that. It's not garbage that's expensive, it's live
> objects.
>
> So, in your case, I'd try tuning the GC. What you want is to have a
> longish period of time between runs, and to have the GC run when
> you're not referencing much data that you're going to throw away. Try
> setting the egc threshold to something like 20 MB.
If you're using CMUCL's or SBCL's gencgc, it also scavenges the entire
static space for each collection - there is no write barrier for static
space. At least for SBCL, that is somewhat large (~5 megs?), though it
doesn't grow.
I think there was some testing done for not running purify (which makes
the static and read-only space) and instead using a tenured 6th
generation (which will use the normal write barrier), but I'm not sure
what the results were.
-bcd
--
*** Brian Downing <bdowning at lavos dot net>
David Steuber <·············@verizon.net> writes:
> The function
>
> Z = Z^2 + C
>
> has rather wild behavior. The closer to the edge of the fractal you
> get, the more iterations are required to tell if you are in or out.
> When I made this movie:
>
> http://www.david-steuber.com/~david/fractal-movie/
I always was amazed that the round-off errors does apparently not kill you
during those calculations. Does anyone know of a reference where this
phenomenon is explained?
Yours, Nicolas
Nicolas Neuss <·······@iwr.uni-heidelberg.de> writes:
> I always was amazed that the round-off errors does apparently not kill you
> during those calculations. Does anyone know of a reference where this
> phenomenon is explained?
I've wondered about that myself. One would expect the round-off
errors to infect the significant digits. I guess that doesn't
happen. I don't know why. I do know that with binary numbers,
rounding and truncating are really the same thing. Perhaps that has
something to do with it.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
On 09 Mar 2004 19:01:05 +0100, Nicolas Neuss
<·······@iwr.uni-heidelberg.de> wrote:
>I always was amazed that the round-off errors does apparently not kill you
>during those calculations. Does anyone know of a reference where this
>phenomenon is explained?
I don't think I understand the question; the original poster mentioned
100,000 iterations on a point, I think; that should only lose you
about 8 bits of precision (out of a 53 bit mantissa); but I may be
misunderstanding the issues involved?
--
"Sore wa himitsu desu."
To reply by email, remove
the small snack from address.
http://www.esatclear.ie/~rwallace
················@eircom.net (Russell Wallace) writes:
> On 09 Mar 2004 19:01:05 +0100, Nicolas Neuss
> <·······@iwr.uni-heidelberg.de> wrote:
>
> >I always was amazed that the round-off errors does apparently not kill you
> >during those calculations. Does anyone know of a reference where this
> >phenomenon is explained?
>
> I don't think I understand the question; the original poster mentioned
> 100,000 iterations on a point, I think; that should only lose you
> about 8 bits of precision (out of a 53 bit mantissa); but I may be
> misunderstanding the issues involved?
The problem is that rounding (clamping) consists in mapping all the
points in some area to one center and go on from there. For the
mandelbrot set, this should be catastrophic, because of it's fractal
nature: whatever the scale, whatever the precision you have, the set
boundary is always as complicated and you must have a lot of these
rounding area across the boundary. For a circle or any regular
figure, you'd just get a good approximation. For the mandelbrot set,
unless we get a mathematical proof validating the transformation, I'd
say that it would lead to meaningless results. (The drawn set will not
be an approximation of the manderbrot set at the suggested scale).
--
__Pascal_Bourguignon__ http://www.informatimago.com/
There is no worse tyranny than to force a man to pay for what he doesn't
want merely because you think it would be good for him.--Robert Heinlein
http://www.theadvocates.org/
On 10 Mar 2004 19:19:12 +0100, Pascal Bourguignon
<····@thalassa.informatimago.com> wrote:
>The problem is that rounding (clamping) consists in mapping all the
>points in some area to one center and go on from there. For the
>mandelbrot set, this should be catastrophic, because of it's fractal
>nature: whatever the scale, whatever the precision you have, the set
>boundary is always as complicated and you must have a lot of these
>rounding area across the boundary.
But wouldn't that apply only to those points close to the boundary?
i.e. be a limiting factor on how close to the boundary you can zoom in
and still get meaningful results?
--
"Sore wa himitsu desu."
To reply by email, remove
the small snack from address.
http://www.esatclear.ie/~rwallace
················@eircom.net (Russell Wallace) writes:
> On 10 Mar 2004 19:19:12 +0100, Pascal Bourguignon
> <····@thalassa.informatimago.com> wrote:
>
> >The problem is that rounding (clamping) consists in mapping all the
> >points in some area to one center and go on from there. For the
> >mandelbrot set, this should be catastrophic, because of it's fractal
> >nature: whatever the scale, whatever the precision you have, the set
> >boundary is always as complicated and you must have a lot of these
> >rounding area across the boundary.
>
> But wouldn't that apply only to those points close to the boundary?
> i.e. be a limiting factor on how close to the boundary you can zoom in
> and still get meaningful results?
When I made a fractal zoom movie, what I did was increase the
max-iterations for the escape function as a function of the zoom
factor. So at lower level zooms, I had fewer iterations that would
get chewed up inside the set.
As far as actual rounding error, that is inevitable with floats
anyway. It doesn't seem to be a real issue until adjacent pixels no
longer have a unique coordinate because the precision ran out. In my
movie, I stoped the zoom short of that happening. But not by much.
The numbers I was using were good for about 15 decimal digits.
So far as I can tell, to determin the boundry, you need to increase
the iterations exponentially with proximity. I expect precision
needs to increase in the same fashion in terms of number of unique
numbers not number of digits. That is, the number of significant
digits should increase linearly with zoom to get a well defined
boundry.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> I may need to have code that abandons the hash table after it reaches
> a certain size.
I wonder how efficient it would be to detect periodicity with
fast and slow pointers, as in LIST-LENGTH.
Kalle Olavi Niemitalo <···@iki.fi> writes:
> David Steuber <·············@verizon.net> writes:
>
> > I may need to have code that abandons the hash table after it reaches
> > a certain size.
>
> I wonder how efficient it would be to detect periodicity with
> fast and slow pointers, as in LIST-LENGTH.
I'm afraid I don't follow.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Kalle Olavi Niemitalo <···@iki.fi> writes:
>
>> I wonder how efficient it would be to detect periodicity with
>> fast and slow pointers, as in LIST-LENGTH.
>
> I'm afraid I don't follow.
I'm sorry. The standard function LIST-LENGTH detects whether a
linked list is circular. The Common Lisp HyperSpec contains a
sample implementation that scans the list which two pointers, one
of which moves faster than the other; if the pointers ever meet
again, the list is circular. If I remember correctly, Knuth
described a similar algorithm for testing the periodicity of
random-number generators.
You could also try letting the numbers stabilize for a few
hundred iterations before you start checking for periodicity.
The length of this initial delay would affect only the speed of
the program and not the final result, so you might be able to
adjust it automatically, based on how fast other nearby pixels
have lead to periodic values.
Kalle Olavi Niemitalo <···@iki.fi> writes:
> David Steuber <·············@verizon.net> writes:
>
> > Kalle Olavi Niemitalo <···@iki.fi> writes:
> >
> >> I wonder how efficient it would be to detect periodicity with
> >> fast and slow pointers, as in LIST-LENGTH.
> >
> > I'm afraid I don't follow.
>
> I'm sorry. The standard function LIST-LENGTH detects whether a
> linked list is circular. The Common Lisp HyperSpec contains a
> sample implementation that scans the list which two pointers, one
> of which moves faster than the other; if the pointers ever meet
> again, the list is circular. If I remember correctly, Knuth
> described a similar algorithm for testing the periodicity of
> random-number generators.
Yes, I looked up LIST-LENGTH in the CLHS when you mentioned it. I
see how it deals with a circular list to return the number of
elements in the list rather than running in circles forever. What I
don't see is the jump to detecting periodicity of numbers. I don't
have Knuth on hand, so I can't check that algorithm for random
numbers. However bad the hash table is, I know that as soon as I see
a number that is already in there, I have detected a period.
> You could also try letting the numbers stabilize for a few
> hundred iterations before you start checking for periodicity.
> The length of this initial delay would affect only the speed of
> the program and not the final result, so you might be able to
> adjust it automatically, based on how fast other nearby pixels
> have lead to periodic values.
I was wondering if there is a way of having a hash table just clear
itself rather than resizing when N number of slots have been filled.
Waiting to create the hash table (or other structure for detecting
periodicity) is a reasonable idea. It certainly takes time for Z to
get "captured" by a period/orbit. The closer to the boundry C is, the
longer it takes. Some special values get captured straight away, but
they are exceptions like #C(-1 0).
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Yes, I looked up LIST-LENGTH in the CLHS when you mentioned it.
> I see how it deals with a circular list to return the number of
> elements in the list rather than running in circles forever.
> What I don't see is the jump to detecting periodicity of numbers.
Just substitute the operations:
(cdr z) -> (clamp (+ (* z z) c))
(endp z) -> (cplx-num-escaped z)
(eq fast slow) -> (= fast slow)
> However bad the hash table is, I know that as soon as I see
> a number that is already in there, I have detected a period.
Yes, the hash table is likely to detect the period in fewer
iterations, and with less arithmetic per iteration. But it uses
more memory.
On 2004-03-13, David Steuber <·············@verizon.net> wrote:
|
| Yes, I looked up LIST-LENGTH in the CLHS when you mentioned it. I
| see how it deals with a circular list to return the number of
| elements in the list rather than running in circles forever. What I
| don't see is the jump to detecting periodicity of numbers. I don't
| have Knuth on hand, so I can't check that algorithm for random
| numbers. However bad the hash table is, I know that as soon as I see
| a number that is already in there, I have detected a period.
Below is some Mandelbrot code I wrote a while ago which does
periodicity checking. It's based on the algorithm used in Fractint.
Basically, ax and ay hold old values of zx and zy, and are updated
periodically (and the period increases over time). If zx=ax and
zy=ay, then we've detected a cycle and can bail out. If you're using
rationals in your code, then you can probably just use
(and (= zx ax) (= zy zy)) instead of (< (abs (- zx ax)) dx) etc.
Also, because it can be relatively time-consuming to do periodicity
checking for points outside the set, it's only done if the last point
calculated was inside the set (that's what the basin variable is for).
(defparameter *max-dwell* 512) ; max number of iterations
(defparameter *bailout* 4.0d0)
(defparameter basin nil) ; was last point calculated inside set?
(defun mandel (x y)
(let ((period 1)
(check 3)
(zx 0.0)
(zy 0.0)
(ax 0.0)
(ay 0.0)
;; dx and dy are small thresholds for
;; checking float equality
(dx (float (/ 2 (* *max-x* *mag*))))
(dy (float (/ 2 (* *max-y* *mag*)))))
(loop for n below *max-dwell*
do (let ((zx2 (* zx zx))
(zy2 (* zy zy)))
(when (> (+ zx2 zy2) *bailout*)
(setq basin nil)
(return n))
(setq zy (+ y (* 2 zx zy))
zx (+ x (- zx2 zy2)))
(when basin
(when (and (< (abs (- zx ax)) dx)
(< (abs (- zy ay)) dy))
(return *max-dwell*))
(incf period)
(when (> period check)
(setq period 1
check (* 2 check)
ax zx
ay zy))))
finally (setq basin t)
(return *max-dwell*))))
--
Cliff Crawford *** ·····@cornell.edu
"The perfection of art is to conceal art." -- Quintilian
David Steuber <·············@verizon.net> writes:
> http://www.david-steuber.com/~david/fractal-movie/
Have you seen the XaoS program by Jan Hubicka? It cheats when
zooming, by reusing pixels from the previous frames. The result
looks still quite good. If I remember correctly, it can also
output motion vectors to help MPEG encoders.
Kalle Olavi Niemitalo <···@iki.fi> writes:
> David Steuber <·············@verizon.net> writes:
>
> > http://www.david-steuber.com/~david/fractal-movie/
>
> Have you seen the XaoS program by Jan Hubicka? It cheats when
> zooming, by reusing pixels from the previous frames. The result
> looks still quite good. If I remember correctly, it can also
> output motion vectors to help MPEG encoders.
I have not seen or even heard of it. You'll probably get a kick out
of this, but I actually calculated every single pixel in that movie.
I split the work in such a way that the even frames were handled by
one CPU and the odd frames by the other in a dual Athlon system.
Seeing what the MPEG-4 compression did to it almost made me cry. I
was thinking that for my next attempt, I would just generate every
other frame like Hanna-Barbara. Another short cut would be to
calculate bitmaps that are lower resolution than the final output. I
was producing 480x360 pixel bitmaps. I could probably get away with
something smaller. Image scaling is cheap.
My data files were basicly raw iteration counts that I post processed
into image files that I then processed with QuickTime and iMovie.
I'm not entirely sure how I can take advantage of previous frames,
but it should be doable.
My movie has a fixed center, so all the pixels appear to fly away
radialy from that point. I don't really know how MPEG encoders
work. QuickTime is able to work with frame sequences, so that is how
I did it.
I've been meaning to post the code that did my end of the work but I
haven't gotten around to it yet. It's mainly a C file for the main
number crunching, a Perl file that drives it for each frame, and
another Perl file that assigned colors to the iterations.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber wrote:
> Oh, I'm talking about calculating whether or not a complex number is
> in the Mandelbrot set. When I did this in C, I was using IEEE-754
> doubles. I even made a movie where I zoom in just short of where I
> can go with 64 bit doubles. I thought perhaps using CL's rational
> numbers, I could go a lot further.
>
> The first problem I ran into was that squaring fractions generates a
> lot of digits very quickly. If you want a way to consume CPU and
> memory resources, I recomend it. I expect there are fixed precision
> libraries available that I should look at. In the meantime, I hacked
> together a function that does the rational equivalent of rounding
> numbers to some settable precision.
Someone's already suggested that you consider using CLISP,
whose LONG-FLOAT type is of user-settable size. That will
probably be faster than using "clamped rationals".
You might want to consider making the threshold at which
you start clamping bigger than the size you clamp to.
That way, you can avoid making *every* operation after
clamping starts take the extra time needed for clamping.
The downside is that some of those operations will be
done on larger numbers. I suspect it's a win overall,
especially if you go with my proposal for more accurate
clamping below.
Your clamping code is unnecessarily expensive; instead
of
(defun clamp-rational (r &key (clamp most-positive-fixnum))
(if (rationalp r)
(if (> (denominator r) clamp)
(let ((scale-divisor (/ (denominator r) clamp)))
(/ (round (/ (numerator r) scale-divisor))
(round (/ (denominator r) scale-divisor))))
r)
r))
try
(defun clamp-rational (r &key (clamp most-positive-fixnum))
(if (and (ratiop r)
(> (denominator r) clamp))
(let ((scale-divisor (round (denominator r) clamp)))
(/ (round (numerator r) scale-divisor)
(round (denominator r) scale-divisor)))
r))
or perhaps
(defun clamp-rational (r &key (clamp most-positive-fixnum))
(if (and (ratiop r)
(> (denominator r) clamp))
(let ((scale-divisor (round (denominator r) clamp)))
(/ (round (* (numerator r) clamp) (denominator r))
clamp))
r))
which is in fact equivalent to what you wrote, though probably
not to what you intended to write :-).
You can make clamping less accuracy-destroying, at some
cost in time, by using continued fractions to find the
smallest "good enough" approximation. I suspect this is
a loss overall, but it might be worth a go.
I should explain that last paragraph a bit better. Here's
some code.
(defun rationalize (x max-denominator)
"Return the closest approximation to X whose denominator
is <= MAX-DENOMINATOR."
(let ((p 0) (q 1) (r 1) (s 0))
;; p/q, r/s are successive c.f. convergents to x
(loop while (<= s max-denominator) do
(multiple-value-bind (intpart fracpart)
(floor x)
(shiftf p r (+ p (* intpart r)))
(shiftf q s (+ q (* intpart s)))
(when (zerop fracpart)
(return-from rationalize (/ r s)))
(setq x (/ fracpart))))
(/ p q)))
Now you could do your clamping by calling RATIONALIZE,
getting better approximations and smaller denominators
than you do at present. But it does quite a lot more work
per clamping than you do. What do you get for that extra
work? The main thing is that if you clamp to denominators
no bigger than N, the error with your method can be on
the order of 1/N whereas the error with mine is O(1/N^2).
So for a given level of accuracy I can afford to have a
clamping limit with half as many bits (ish), so all my
arithmetic operations are (once we're into bignum-land,
and assuming the numbers aren't big enough for sophisticated
arithmetic algorithms to be being used) 2-4 times faster
than yours. Whether that's worth it or not depends on
how often we clamp.
Cool though having native complex numbers is, you *may* find
that you actually get better performance by separating out
the real and imaginary parts. Especially as then you can
avoid squaring the real and imaginary parts of z twice
(once for the "attracted to infinity?" check, once for
the calculation of the real part of z^2). In your application,
especially once you start using bigna, the latter is likely
to matter more than the former.
I suspect you may get better performance without the hash
table than with it, because surely most points -- even
rational ones of small denominator -- in the Mandelbrot
set aren't periodic, and none of the points outside it.
For plots that include substantial lumps of the set itself,
there are ways to prove that (e.g.) a whole disc lies within
the set; it may be possible to use this to do the calculations
more efficiently than you could do by processing each pixel
separately.
Consider using LOOP instead of DO, since it's more readable
(except to LOOP-phobics):
(let ((seen (make-hash-table :test #'eql)))
(loop for z = #C(0 0) then (clamp-complex (+ (* z z) c))
for iterations upfrom 0 below max-iterations
for escaped = (cplx-num-escaped z)
do (when escaped (return (values t iterations)))
(when (gethash z seen) (return (values nil iterations)))
(setf (gethash z seen) t)
finally (return (values nil iterations))))
I agree with whoever suggested that the clamping parameter(s)
should be special variables rather than being passed around
everywhere explicitly.
--
Gareth McCaughan
.sig under construc
Gareth McCaughan <················@pobox.com> writes:
> Someone's already suggested that you consider using CLISP,
> whose LONG-FLOAT type is of user-settable size. That will
> probably be faster than using "clamped rationals".
That may well be the most appropriate thing for this problem. Except
I am also trying to learn Lisp and a CLISP only feature might add to
much to the mix.
> You might want to consider making the threshold at which
> you start clamping bigger than the size you clamp to.
> That way, you can avoid making *every* operation after
> clamping starts take the extra time needed for clamping.
> The downside is that some of those operations will be
> done on larger numbers. I suspect it's a win overall,
> especially if you go with my proposal for more accurate
> clamping below.
I've thought about that and not yet decided about it. More
experimentation is required.
> Your clamping code is unnecessarily expensive;
Yes. I've fixed it a bit before seeing this post.
> You can make clamping less accuracy-destroying, at some
> cost in time, by using continued fractions to find the
> smallest "good enough" approximation. I suspect this is
> a loss overall, but it might be worth a go.
>
> I should explain that last paragraph a bit better. Here's
> some code.
>
> (defun rationalize (x max-denominator)
> "Return the closest approximation to X whose denominator
> is <= MAX-DENOMINATOR."
> (let ((p 0) (q 1) (r 1) (s 0))
> ;; p/q, r/s are successive c.f. convergents to x
> (loop while (<= s max-denominator) do
> (multiple-value-bind (intpart fracpart)
> (floor x)
> (shiftf p r (+ p (* intpart r)))
> (shiftf q s (+ q (* intpart s)))
> (when (zerop fracpart)
> (return-from rationalize (/ r s)))
> (setq x (/ fracpart))))
> (/ p q)))
I'll have to look at this. I don't understand it on the face of it.
> Cool though having native complex numbers is, you *may* find
> that you actually get better performance by separating out
> the real and imaginary parts. Especially as then you can
> avoid squaring the real and imaginary parts of z twice
> (once for the "attracted to infinity?" check, once for
> the calculation of the real part of z^2). In your application,
> especially once you start using bigna, the latter is likely
> to matter more than the former.
This is what I did in my C code. I had a bunch of temp variables
though. I don't know if that was a good thing overall or bad.
The C version of the main loop looked like this:
// Z = Z^2 + C
int calcPoint(int maxItt, double c_real, double c_imag)
{
int count;
double temp, temp_2real, temp_2imag, z_real, z_imag;
temp = temp_2real = temp_2imag = z_real = z_imag = 0.0L;
for (count = 0; count < maxItt && temp_2real + temp_2imag < 4.0L; count++) {
temp = temp_2real - temp_2imag + c_real;
z_imag = z_real * z_imag * 2.0L + c_imag;
z_real = temp;
temp_2real = z_real * z_real;
temp_2imag = z_imag * z_imag;
}
return count;
}
> I suspect you may get better performance without the hash
> table than with it, because surely most points -- even
> rational ones of small denominator -- in the Mandelbrot
> set aren't periodic, and none of the points outside it.
I'm not sure either way. The main reason I return two values instead
of just one is to see how often I bail because I detected a cycle
rather than running out of iterations. Otherwise, I would just
return NIL instead of a number when I was in the set.
> For plots that include substantial lumps of the set itself,
> there are ways to prove that (e.g.) a whole disc lies within
> the set; it may be possible to use this to do the calculations
> more efficiently than you could do by processing each pixel
> separately.
The problem is I don't know how to find the edges of a disc or
cardioid and then fill it in. Then there are all those discs that
are on the edge of any disc or cardioid.
> Consider using LOOP instead of DO, since it's more readable
> (except to LOOP-phobics):
>
> (let ((seen (make-hash-table :test #'eql)))
> (loop for z = #C(0 0) then (clamp-complex (+ (* z z) c))
> for iterations upfrom 0 below max-iterations
> for escaped = (cplx-num-escaped z)
> do (when escaped (return (values t iterations)))
> (when (gethash z seen) (return (values nil iterations)))
> (setf (gethash z seen) t)
> finally (return (values nil iterations))))
I haven't inerned loop yet as well as do. But I will try it.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber wrote:
[I said:]
>> I should explain that last paragraph a bit better. Here's
>> some code.
>>
>> (defun rationalize (x max-denominator)
>> "Return the closest approximation to X whose denominator
>> is <= MAX-DENOMINATOR."
>> (let ((p 0) (q 1) (r 1) (s 0))
>> ;; p/q, r/s are successive c.f. convergents to x
>> (loop while (<= s max-denominator) do
>> (multiple-value-bind (intpart fracpart)
>> (floor x)
>> (shiftf p r (+ p (* intpart r)))
>> (shiftf q s (+ q (* intpart s)))
>> (when (zerop fracpart)
>> (return-from rationalize (/ r s)))
>> (setq x (/ fracpart))))
>> (/ p q)))
>
> I'll have to look at this. I don't understand it on the face of it.
Um, it's not meant to be obvious how it works :-). It was
meant to be a black box that has the effect of returning
the rational number closest to X with denominator <= MAX-DENOMINATOR.
But, since you ask ...
---------- begin digression ----------
Every non-negative number can be written as a "continued fraction",
in the form
n0 + 1/(n1 + 1/(n2 + 1/(...)))
where n0,n1,n2,... are non-negative integers called its "coefficients",
only the first of which is allowed to be 0. It's easy to find
these integers; split the number into its integer part n0 and
its fractional part x0, and then observe that the continued fraction
for 1/x0 has coefficients n1,n2,... so that we can (according to
taste) iterate or recurse.
(defun cf (x)
"Return the continued fraction for X, which should be
a non-negative rational number."
(multiple-value-bind (intpart fracpart)
(floor x)
(cons x (if (zerop fracpart) nil (cf (/ fracpart))))))
If the number is irrational then the continued fraction never
terminates; there are infinitely many coefficients. However,
you can truncate it at any point, just throwing away all the
subsequent coefficients, and then the continued fraction of
what remains is a rational approximation to x. It's called a
"convergent" to x; thus the first convergent is just n0, the
integer part of x, the next takes account of n1 too, and so
on.
Here's an inverse to CF that indicates one way of calculating
successive convergents. You should recognize about half of the
RATIONALIZE function above in here, the other half being the
CF function I gave a moment ago.
(defun df (coeffs)
(let ((p 0) (q 1) (r 1) (s 0))
(loop for n in coeffs do
(shiftf p r (+ p (* n r)))
(shiftf q s (+ q (* n s))))
(/ r s)))
Why does this work? Sketchy explanation: suppose df(coeffs) = p/q;
then df(n . coeffs) = n + q/p. Represent rational numbers as 2-element
vectors; then the operation that maps p/q to n+q/p corresponds to
premultiplication by the 2x2 matrix (n 1; 1 0). The DF function
is accumulating the product of those matrices, from left to right.
The reason why all this is interesting (well, one reason, anyway)
is the following fact. Suppose p/q and r/s are successive convergents
to some number x. Then p/q is the nearest rational number to x whose
denominator is < s. Hence you can find the nearest rational number
to x whose denominator is <= any particular bound by calculating
its continued fraction coefficients (as with CF above) and, as you
go, calculating the convergents (as with DF above), and stopping
when a convergent's denominator exceeds the bound. That's exactly
what RATIONALIZE does.
There's another nice theorem, which says that if p/q and r/s are
successive convergents to x, then x lies between p/q and r/s and
|p/q-x| is at most 1/qs. (I think there's an embarrassingly short
proof of this, but I forget.) Hence the convergents are not only
*best* approximations to x, but *very good* approximations. And
when one of the coefficients in a continued fraction is large,
the preceding convergent is very close; that's why 355/113 is
such a good approximation to pi. (Work out the beginning of the
continued fraction for pi and you'll see.) Conversely, if all
the coefficients are small, then that limits how closely x can
be approximated by rational numbers; for instance, the continued
fraction for sqrt(2) has coefficients 1,2,2,2,2,2,..., so sqrt(2)
isn't very well approximated by rational numbers ... which (kinda)
is why the musical interval of an augmented fourth sounds so
nasty. :-) Anyway, we're now well and truly off topic, so I'll
shut up about continued fractions.
---------- end digression ----------
> The C version of the main loop looked like this:
>
> // Z = Z^2 + C
> int calcPoint(int maxItt, double c_real, double c_imag)
> {
> int count;
> double temp, temp_2real, temp_2imag, z_real, z_imag;
>
> temp = temp_2real = temp_2imag = z_real = z_imag = 0.0L;
>
> for (count = 0; count < maxItt && temp_2real + temp_2imag < 4.0L; count++) {
> temp = temp_2real - temp_2imag + c_real;
> z_imag = z_real * z_imag * 2.0L + c_imag;
> z_real = temp;
> temp_2real = z_real * z_real;
> temp_2imag = z_imag * z_imag;
> }
> return count;
> }
Yes, I'm afraid that's the sort of thing I'm suggesting :-).
>> For plots that include substantial lumps of the set itself,
>> there are ways to prove that (e.g.) a whole disc lies within
>> the set; it may be possible to use this to do the calculations
>> more efficiently than you could do by processing each pixel
>> separately.
>
> The problem is I don't know how to find the edges of a disc or
> cardioid and then fill it in. Then there are all those discs that
> are on the edge of any disc or cardioid.
There are, as I say, definitely algorithms for this sort of thing.
I don't know how practical they are for your purpose. My hazy memory
says that there was a book many years ago by someone with a name
like Hans-Otto Peitgen (I'm more sure of the surname than the first
names), which might have been called "The beauty of fractals" (there
certainly was such a book, but I'm not sure whether it's the one
I'm thinking of), which discussed some algorithms along these lines.
But I could be spectacularly wrong :-). (Which is why I'm giving you
this brain dump rather than checking with Google and giving a precise
reference complete with the misleading impression of authority.)
--
Gareth McCaughan
.sig under construc
Gareth McCaughan <················@pobox.com> writes:
> David Steuber wrote:
>
> [I said:]
> >> I should explain that last paragraph a bit better. Here's
> >> some code.
> >>
> >> (defun rationalize (x max-denominator)
> >> "Return the closest approximation to X whose denominator
> >> is <= MAX-DENOMINATOR."
> >> (let ((p 0) (q 1) (r 1) (s 0))
> >> ;; p/q, r/s are successive c.f. convergents to x
> >> (loop while (<= s max-denominator) do
> >> (multiple-value-bind (intpart fracpart)
> >> (floor x)
> >> (shiftf p r (+ p (* intpart r)))
> >> (shiftf q s (+ q (* intpart s)))
> >> (when (zerop fracpart)
> >> (return-from rationalize (/ r s)))
> >> (setq x (/ fracpart))))
> >> (/ p q)))
> >
> > I'll have to look at this. I don't understand it on the face of it.
>
> Um, it's not meant to be obvious how it works :-). It was
> meant to be a black box that has the effect of returning
> the rational number closest to X with denominator <= MAX-DENOMINATOR.
> But, since you ask ...
I really do wish my math (such that I had) did not decide to go
away. I'm probably reading something wrong, but it seems like x is a
real number rather than a rational. I think CL has a function for
turning floats into the closest approximate ratio.
> > The C version of the main loop looked like this:
> >
> > // Z = Z^2 + C
> > int calcPoint(int maxItt, double c_real, double c_imag)
> > {
> > int count;
> > double temp, temp_2real, temp_2imag, z_real, z_imag;
> >
> > temp = temp_2real = temp_2imag = z_real = z_imag = 0.0L;
> >
> > for (count = 0; count < maxItt && temp_2real + temp_2imag < 4.0L; count++) {
> > temp = temp_2real - temp_2imag + c_real;
> > z_imag = z_real * z_imag * 2.0L + c_imag;
> > z_real = temp;
> > temp_2real = z_real * z_real;
> > temp_2imag = z_imag * z_imag;
> > }
> > return count;
> > }
>
> Yes, I'm afraid that's the sort of thing I'm suggesting :-).
I am looking at LOOP in the CLHS. I also need to check out the
modified BNF notation section because I'm having a heck of a time
with it. Recasting the C function above into Lisp is trivial
enough. It is only more work for me because I am not yet thinking in
Lisp.
> >> For plots that include substantial lumps of the set itself,
> >> there are ways to prove that (e.g.) a whole disc lies within
> >> the set; it may be possible to use this to do the calculations
> >> more efficiently than you could do by processing each pixel
> >> separately.
> >
> > The problem is I don't know how to find the edges of a disc or
> > cardioid and then fill it in. Then there are all those discs that
> > are on the edge of any disc or cardioid.
>
> There are, as I say, definitely algorithms for this sort of thing.
> I don't know how practical they are for your purpose. My hazy memory
> says that there was a book many years ago by someone with a name
> like Hans-Otto Peitgen (I'm more sure of the surname than the first
> names), which might have been called "The beauty of fractals" (there
> certainly was such a book, but I'm not sure whether it's the one
> I'm thinking of), which discussed some algorithms along these lines.
> But I could be spectacularly wrong :-). (Which is why I'm giving you
> this brain dump rather than checking with Google and giving a precise
> reference complete with the misleading impression of authority.)
I'll certainly be Googling around for such info. I think I have
heard of that book, but I do not have a copy.
It is certainly desirable to know that a point is in the set quickly
regardless of the zoom factor. Consuming all the iterations does
take quite a bit of time. At least outside the set I get a number
that I can then use to map to a color to make a pretty picture. I've
had some fun with logarithms doing that.
I feel like I've regressed back to high school level programming
somehow. Perhaps I never really advanced beyond that point. I
wonder what my percentile ranking is in the set of people who program
for money (or try to).
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber wrote:
[I said:]
> > Um, it's not meant to be obvious how it works :-). It was
> > meant to be a black box that has the effect of returning
> > the rational number closest to X with denominator <= MAX-DENOMINATOR.
> > But, since you ask ...
>
> I really do wish my math (such that I had) did not decide to go
> away. I'm probably reading something wrong, but it seems like x is a
> real number rather than a rational. I think CL has a function for
> turning floats into the closest approximate ratio.
I explained in terms of rationals, but when you don't need
to either produce or consume the entire potentially-infinite
continued fraction you can do all the same stuff for any
positive real number.
> I feel like I've regressed back to high school level programming
> somehow. Perhaps I never really advanced beyond that point. I
> wonder what my percentile ranking is in the set of people who program
> for money (or try to).
The fact that you're capable of wondering that suggests to me
that your ranking is probably pretty good :-).
--
Gareth McCaughan
.sig under construc
Gareth McCaughan <················@pobox.com> writes:
> I don't know how practical they are for your purpose. My hazy memory
> says that there was a book many years ago by someone with a name
> like Hans-Otto Peitgen (I'm more sure of the surname than the first
> names), which might have been called "The beauty of fractals" (there
> certainly was such a book, but I'm not sure whether it's the one
> I'm thinking of), which discussed some algorithms along these lines.
I don't know the book you mention, but there's on of Heinz-Otto
Peitgen and Dietmar Saupe with the title "The Science of Fractal
Images". Maybe you are thinking of that one? It's a beautiful book and
in chapter 2 it indeed discusses algorithms for fractals. (But it talks
more about brownian motion than generating nice pictures.)
Andreas
--
Wherever I lay my .emacs, there's my $HOME.
Andreas Eder <············@t-online.de> writes:
> Gareth McCaughan <················@pobox.com> writes:
> > My hazy memory says that there was a book many years ago by
> > someone with a name like Hans-Otto Peitgen (I'm more sure of the
> > surname than the first names), which might have been called "The
> > beauty of fractals" (there certainly was such a book, but I'm not
> > sure whether it's the one I'm thinking of), which discussed some
> > algorithms along these lines.
> I don't know the book you mention, but there's on of Heinz-Otto
> Peitgen and Dietmar Saupe with the title "The Science of Fractal
> Images". Maybe you are thinking of that one? It's a beautiful book and
> in chapter 2 it indeed discusses algorithms for fractals. (But it talks
> more about brownian motion than generating nice pictures.)
I think Gareth refers to sections D.1 Bounding the distance to M ("In
this section we compute a lower bound on the distance from a point c
to M..."), and D.2 Finding disks in the interior of M ("In this
section we derive a lower bound on the distance from a point c \in M
to the boundary of M.")
--
Lars Brinkhoff, Services for Unix, Linux, GCC, HTTP
Brinkhoff Consulting http://www.brinkhoff.se/
Lars Brinkhoff <·········@nocrew.org> writes:
> Andreas Eder <············@t-online.de> writes:
> > Gareth McCaughan <················@pobox.com> writes:
> > > My hazy memory says that there was a book many years ago by
> > > someone with a name like Hans-Otto Peitgen (I'm more sure of the
> > > surname than the first names), which might have been called "The
> > > beauty of fractals" (there certainly was such a book, but I'm not
> > > sure whether it's the one I'm thinking of), which discussed some
> > > algorithms along these lines.
> >
> > I don't know the book you mention, but there's on of Heinz-Otto
> > Peitgen and Dietmar Saupe with the title "The Science of Fractal
> > Images". Maybe you are thinking of that one? It's a beautiful book and
> > in chapter 2 it indeed discusses algorithms for fractals. (But it talks
> > more about brownian motion than generating nice pictures.)
>
> I think Gareth refers to sections D.1 Bounding the distance to M ("In
> this section we compute a lower bound on the distance from a point c
> to M..."), and D.2 Finding disks in the interior of M ("In this
> section we derive a lower bound on the distance from a point c \in M
> to the boundary of M.")
Aha, that sounds quite familiar. I should mention that it's at least
12 years, probably more, since I saw the book :-).
--
Gareth McCaughan
.sig under construc
Lars Brinkhoff <·········@nocrew.org> writes:
> I think Gareth refers to sections D.1 Bounding the distance to M ("In
> this section we compute a lower bound on the distance from a point c
> to M..."), and D.2 Finding disks in the interior of M ("In this
> section we derive a lower bound on the distance from a point c \in M
> to the boundary of M.")
I entered "finding disks in the interior of m" and this is one of the
hits I got:
http://www.fractalus.com/kerry/articles/area/mandelbrot-area.html
I don't know if it will help me or not, particularly on zooms > 100,000,000 or whatever.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> I don't know if it will help me or not, particularly on zooms >
> 100,000,000 or whatever.
The algorithm in the book involves finding periodic cycles in the
iterative computation. As far as I can see, you usually don't have a
lot of "interior M" area in your images, so I don't think you should
bother.
--
Lars Brinkhoff, Services for Unix, Linux, GCC, HTTP
Brinkhoff Consulting http://www.brinkhoff.se/
David Steuber wrote:
>
> The problem is I don't know how to find the edges of a disc or
> cardioid and then fill it in. Then there are all those discs that
> are on the edge of any disc or cardioid.
>
One of the mandelbrot's nifty properties is that all these regions
are contiguous; if all the pixels surrounding any region are in the
set (or for that matter in any "iterates N times before escaping"
group) then all the pixels in the region are in the set (or in that
group).
Bear
Ray Dillinger <····@sonic.net> writes:
> David Steuber wrote:
> >
> > The problem is I don't know how to find the edges of a disc or
> > cardioid and then fill it in. Then there are all those discs that
> > are on the edge of any disc or cardioid.
>
> One of the mandelbrot's nifty properties is that all these regions
> are contiguous; if all the pixels surrounding any region are in the
> set (or for that matter in any "iterates N times before escaping"
> group) then all the pixels in the region are in the set (or in that
> group).
This is true. While the Gausian like field lines outside the set are
rather turbulant, perhaps not making it worth doing edge detection,
your mention of this property reminds me of convex hulls in graphics
and the algorithm for filling them in. If I work around the boundry
of an in M artifact, I then have a convex hull I can fill in.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber wrote:
> The first problem I ran into was that squaring fractions generates a
> lot of digits very quickly. If you want a way to consume CPU and
> memory resources, I recomend it. I expect there are fixed precision
> libraries available that I should look at. In the meantime, I hacked
> together a function that does the rational equivalent of rounding
> numbers to some settable precision.
Here's another simple approach which I think will probably be
faster than anything based on rationals. Roll your own fixed-precision
arithmetic.
(defvar *bits* 64
"Number of fractional bits to compute with.")
(defvar *threshold* (ash 4 (* 2 *bits*))
"The value beyond which iteration is deemed to have escaped.")
(defun shift-down (x)
(ash x (- *bits*)))
;; or, more accurately, (round x (ash 1 *bits*))
;; or (if (logbitp x (1- *bits*)) (1+ (ash x (- *bits*))) (ash x (- *bits*)))
;; or any of a number of other equivalent formulations whose speed and
;; accuracy will vary.
(defun next (x y cx cy)
;; x,y are integers, representing the number (x+i.y)/2^n.
;; cx,cy are integers, representing the number (cx+i.cy)/4^n.
;; Either return (values new-x new-y) where (newx+i.newy)/2^n
;; ~= [(x+iy)/2^n]^2 + c/4^n, or (if |x+iy|/2^n >= 2)
;; return NIL instead.
(let ((x2 (* x x)) (y2 (* y y)))
(when (>= (+ x2 y2) *threshold*)
(return-from next nil))
(values (shift-down (+ (- x2 y2) cx))
(shift-down (+ (* 2 x y) cy)))))
And so on. (All the code above is untested and may be horribly broken,
but I hope the idea's reasonably clear.)
--
Gareth McCaughan
.sig under construc
Gareth McCaughan <················@pobox.com> writes:
> Here's another simple approach which I think will probably be
> faster than anything based on rationals. Roll your own fixed-precision
> arithmetic.
Rolling my own fixed-precision arithmetic has always been on my
mind. Figuring out how to do it in an efficient way has always been
a problem (not that I've worked hard on it). Thanks for the code.
ASH is one of those Lisp words I do not know yet. I have the CLHS
locally, so I can look it up no problem.
I think I mentioned this before. I did 32 bit fixed-precision
arithmetic in 80386 assembler. The result of MUL on a 32 bit number
saved the entire 64 bit result. All I had to do was grab the high
word. It was the next step that I never tried; working with numbers
that were vectors of words.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Gareth McCaughan <················@pobox.com> writes:
>
> > Here's another simple approach which I think will probably be
> > faster than anything based on rationals. Roll your own fixed-precision
> > arithmetic.
>
> Rolling my own fixed-precision arithmetic has always been on my
> mind. Figuring out how to do it in an efficient way has always been
> a problem (not that I've worked hard on it). Thanks for the code.
> ASH is one of those Lisp words I do not know yet. I have the CLHS
> locally, so I can look it up no problem.
Well, I just spent a bunch of hours hacking at this. It's amazing
how few lines of code I have to show for my time. I don't know if
that's a good thing or a bad thing at this stage. My fractal code is
now completely rewritten. It looks a lot like the C code! Well, the
loop part does anyway. I'm sure there is still room for improvement
with both the algorithm and the style. Periodicity checking is gone
for now. I may bring it back. I don't know. I'm thinking more
along the lines of doing edge detection for Cs that are in M.
Here is what I have now with some test runs.
;;;; fractal.lisp -- an exploration of the Mandelbrot Set using
;;;; arbitrary sized integer math.
(defvar *fractional-bits* 16
"Number of fractional bits in our number")
(defvar *escape-value* (ash 4 *fractional-bits*)
"Numbers larger than this have escaped and are not part of the Mandelbrot Set")
(defvar *my-two* (ash 2 *fractional-bits*))
(defun set-fractional-bits (fb)
"A nice way to set both *fractional-bits* and *escape-value* at the same time"
(setf *escape-value* (ash 4 fb)
*my-two* (ash 2 fb)
*fractional-bits* fb))
(defun shifted-multiply (&rest args)
"Multiplies the args as integers left-shifted by *fractional-bits* and right-shifts accordingly"
(let ((n (apply #'* args)))
(ash n (- (* (max 0 (- (list-length args) 1)) *fractional-bits*)))))
;;; mandelbrot-escape-fn -- (values escaped iterations)
(defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
"Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
(loop for iterations below max-iterations
with z-real = 0
and z-imag = 0
and z-real-next = 0
and z-real-squared = 0
and z-imag-squared = 0
unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
do (setf z-real-next (+ (- z-real-squared z-imag-squared) c-real))
(setf z-imag (+ (shifted-multiply z-real z-imag *my-two*) c-imag))
(setf z-real z-real-next)
(setf z-real-squared (shifted-multiply z-real z-real))
(setf z-imag-squared (shifted-multiply z-imag z-imag))
finally (return (values nil iterations))))
; (time (mandelbrot-escape-fn (ash -1 *fractional-bits*) (ash 1 (- *fractional-bits* 2)) :max-iterations 1000000))
CL-USER> (set-fractional-bits 32)
32
CL-USER> *my-two*
8589934592
CL-USER> *escape-value*
17179869184
CL-USER> (shifted-multiply *my-two* *my-two*)
17179869184
CL-USER> (set-fractional-bits 8)
8
;;;; (time (mandelbrot-escape-fn (ash -1 *fractional-bits*) (ash ...
(MANDELBROT-ESCAPE-FN (ASH -1 *FRACTIONAL-BITS*) (ASH 1 (- *FRACTIONAL-BITS* 2)) :MAX-ITERATIONS 1000000) took 2,162 milliseconds (2.162 seconds) to run.
Of that, 1,450 milliseconds (1.450 seconds) were spent in user mode
580 milliseconds (0.580 seconds) were spent in system mode
132 milliseconds (0.132 seconds) were spent executing other OS processes.
351 milliseconds (0.351 seconds) was spent in GC.
56,000,000 bytes of memory allocated.
CL-USER> (set-fractional-bits 16)
16
;;;; (time (mandelbrot-escape-fn (ash -1 *fractional-bits*) (ash ...
(MANDELBROT-ESCAPE-FN (ASH -1 *FRACTIONAL-BITS*) (ASH 1 (- *FRACTIONAL-BITS* 2)) :MAX-ITERATIONS 1000000) took 4,203 milliseconds (4.203 seconds) to run.
Of that, 2,980 milliseconds (2.980 seconds) were spent in user mode
1,040 milliseconds (1.040 seconds) were spent in system mode
183 milliseconds (0.183 seconds) were spent executing other OS processes.
631 milliseconds (0.631 seconds) was spent in GC.
99,999,976 bytes of memory allocated.
CL-USER> (set-fractional-bits 32)
32
;;;; (time (mandelbrot-escape-fn (ash -1 *fractional-bits*) (ash ...
(MANDELBROT-ESCAPE-FN (ASH -1 *FRACTIONAL-BITS*) (ASH 1 (- *FRACTIONAL-BITS* 2)) :MAX-ITERATIONS 1000000) took 12,101 milliseconds (12.101 seconds) to run.
Of that, 9,130 milliseconds (9.130 seconds) were spent in user mode
2,440 milliseconds (2.440 seconds) were spent in system mode
531 milliseconds (0.531 seconds) were spent executing other OS processes.
1,347 milliseconds (1.347 seconds) was spent in GC.
208,027,400 bytes of memory allocated.
CL-USER> (set-fractional-bits 64)
64
;;;; (time (mandelbrot-escape-fn (ash -1 *fractional-bits*) (ash ...
(MANDELBROT-ESCAPE-FN (ASH -1 *FRACTIONAL-BITS*) (ASH 1 (- *FRACTIONAL-BITS* 2)) :MAX-ITERATIONS 1000000) took 14,884 milliseconds (14.884 seconds) to run.
Of that, 11,030 milliseconds (11.030 seconds) were spent in user mode
3,280 milliseconds (3.280 seconds) were spent in system mode
574 milliseconds (0.574 seconds) were spent executing other OS processes.
1,953 milliseconds (1.953 seconds) was spent in GC.
292,027,360 bytes of memory allocated.
CL-USER> (set-fractional-bits 128)
128
;;;; (time (mandelbrot-escape-fn (ash -1 *fractional-bits*) (ash ...
(MANDELBROT-ESCAPE-FN (ASH -1 *FRACTIONAL-BITS*) (ASH 1 (- *FRACTIONAL-BITS* 2)) :MAX-ITERATIONS 1000000) took 18,909 milliseconds (18.909 seconds) to run.
Of that, 13,210 milliseconds (13.210 seconds) were spent in user mode
4,850 milliseconds (4.850 seconds) were spent in system mode
849 milliseconds (0.849 seconds) were spent executing other OS processes.
2,876 milliseconds (2.876 seconds) was spent in GC.
420,027,320 bytes of memory allocated.
I must confess that I am surprised at the way the time increased with
*fractional-bits*. It's not linear. It doesn't look like log(n)
either. Oh well. I'm still not sure if the code will be fast enough
to do what I want to do.
SLIME doesn't seem to give me the ability to interrupt the process.
Maybe it's a simple case of RTFM if there was one. I found myself
killing the Lisp (OpenMCL) more than once and restarting SLIME.
I tried C-g.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> SLIME doesn't seem to give me the ability to interrupt the process.
> Maybe it's a simple case of RTFM if there was one. I found myself
> killing the Lisp (OpenMCL) more than once and restarting SLIME.
>
> I tried C-g.
Try C-c C-c instead. C-g is for interrupting Emacs; C-c C-c
tells Emacs to send a C-c (SIGINT) to the Lisp process. Note: this key
mapping works in (some of) the SLIME buffers.
C-h b is highly recommended to find out what keybindings are
active in the current buffers - I discovered something about the SLIME
debugger window a couple of days back this way :-)
SLIME is wonderful, SBCL is wonderful, MacOSX is wonderful...
and I have to work on porting !"#$%& C++ code from Windows to *nix.
Life is hard.
--
Raymond Wiker Mail: ·············@fast.no
Senior Software Engineer Web: http://www.fast.no/
Fast Search & Transfer ASA Phone: +47 23 01 11 60
P.O. Box 1677 Vika Fax: +47 35 54 87 99
NO-0120 Oslo, NORWAY Mob: +47 48 01 11 60
Try FAST Search: http://alltheweb.com/
Raymond Wiker <·············@fast.no> writes:
> SLIME is wonderful, SBCL is wonderful, MacOSX is wonderful...
> and I have to work on porting !"#$%& C++ code from Windows to *nix.
> Life is hard.
MFC code? If so, I don't envy you.
I took a stab at Objective-C + Cocoa for OSX. There are some things
I don't get about who is responsible for deallocating things (gotta
love Lisp for taking that worry away). I decided to step back from
Cocoa to learn Lisp first. Once I have some comfort with that, then
it will be back to Cocoa. The end goal is to be able to write
commercial applications in Lisp that play nice with the OSX user
interface guidelines.
I would soooooo love to be able to use Interface Builder for a Lisp
app. That would really rock. One method I heard someone propose (I
don't remember who) was to write a Cocoa GUI using Xcode and have it
talk to a Lisp listener for doing the major heavy lifting. I'm not
sure how that would work in the context of an app that needs to use
OpenGL though.
Maybe Kenny's Cello project will provide an answer there. I haven't
looked at it yet. I'm still in the Hello World stage.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Raymond Wiker <·············@fast.no> writes:
>
> > SLIME is wonderful, SBCL is wonderful, MacOSX is wonderful...
> > and I have to work on porting !"#$%& C++ code from Windows to *nix.
> > Life is hard.
>
> MFC code? If so, I don't envy you.
>
> I took a stab at Objective-C + Cocoa for OSX. There are some things
> I don't get about who is responsible for deallocating things (gotta
> love Lisp for taking that worry away).
In OpenStep, it's reference counting, and there is no designated
"owner" of an object. Objects interested in an object just increase
the reference count sending it the retain message, and when they're
done, indicate it decreasing the reference count sending a release
message. When the reference count reaches 0, the object itself commit
suicide (which may lead to the releasing of the objects it was
interested in and further deallocating). The interesting bits is how
an object can return a new objects it's not interested in. This is
done with the autorelease mechanism, where the object is placed in an
autorelease pool (which is interested in the objects, therefore which
retains them!). The autorelease pool is released and deallocated in
the event loop, so this leave time to the caller to retain the
autoreleased object if it wants to keep it. When the autorelease pool
is deallocated, it releases all the objects it contains, and those
only whose reference count reaches 0 are deallocated (the temporary
objects, the newest generation). The remaining item is that some
constructors return retained objects, assuming the caller wants to
keep them, while other return autoreleased objects, assuming temporary
objects. [[Class alloc] initXxxx] return retained objects, while
[Class classXxxx] return autoreleased objects.
That said, in your programs you want to avoid circular chains of
retain, so you must either structure your data in an acyclic graph, or
you must take note that some objects are "owned" by others and then
have back references not counted in the reference count (where the
"owned" object does not retain its "owner").
But OpenStep contains no cycle.
> I decided to step back from Cocoa to learn Lisp first. Once I have
> some comfort with that, then it will be back to Cocoa. The end goal
> is to be able to write commercial applications in Lisp that play
> nice with the OSX user interface guidelines.
>
> I would soooooo love to be able to use Interface Builder for a Lisp
> app. That would really rock.
What's a shame is that the original Interface Builder was a Lisp program!
> One method I heard someone propose (I don't remember who) was to
> write a Cocoa GUI using Xcode and have it talk to a Lisp listener
> for doing the major heavy lifting. I'm not sure how that would work
> in the context of an app that needs to use OpenGL though.
openmcl has Objective-C bindings allowing to mix Objective-C and CLOS
classes and objects. You should be able to write OpenStep GUI
directly in Common-Lisp, including loading and using .nib files
designed in Interface Builder.
--
__Pascal_Bourguignon__ http://www.informatimago.com/
There is no worse tyranny than to force a man to pay for what he doesn't
want merely because you think it would be good for him.--Robert Heinlein
http://www.theadvocates.org/
Pascal Bourguignon <····@thalassa.informatimago.com> writes:
> openmcl has Objective-C bindings allowing to mix Objective-C and CLOS
> classes and objects. You should be able to write OpenStep GUI
> directly in Common-Lisp, including loading and using .nib files
> designed in Interface Builder.
This would indeed be a good solution.
Of course, I'm still climbing the learning curve. The OS X
documentation already has more info in it than I can cache. The same
thing goes for the CLHS. I'm looking for commonalities (be they
synonyms or whatever) that I can use to factor the information down
into a smaller number of chunks.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
Peter Seibel <·····@javamonkey.com> writes:
> > David Steuber <·············@verizon.net> writes:
>
> > ;;; mandelbrot-escape-fn -- (values escaped iterations)
> > (defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
> > "Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
> > (loop for iterations below max-iterations
> > with z-real = 0
> > and z-imag = 0
> > and z-real-next = 0
> > and z-real-squared = 0
> > and z-imag-squared = 0
> > unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
> > do (setf z-real-next (+ (- z-real-squared z-imag-squared) c-real))
> > (setf z-imag (+ (shifted-multiply z-real z-imag *my-two*) c-imag))
> > (setf z-real z-real-next)
> > (setf z-real-squared (shifted-multiply z-real z-real))
> > (setf z-imag-squared (shifted-multiply z-imag z-imag))
> > finally (return (values nil iterations))))
>
> I think you can right this (IMO) more clearly like this (not tested):
>
> (defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
> "Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
> (loop for iterations below max-iterations
> for z-imag = 0 then (+ (shifted-multiply z-real z-imag *my-two*) c-imag)
> for z-real = 0 then (+ (- z-real-squared z-imag-squared) c-real)
> for z-real-squared = 0 then (shifted-multiply z-real z-real)
> for z-imag-squared = 0 then (shifted-multiply z-imag z-imag)
> unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
> finally (return (values nil iterations)))
You are shy one closing paren (copy paste error?), but it looks
equivalent providing those steppings are parallel. I always hated
that z-real-next variable. This also gets back to the brevity of the
C version. I tried it out and it looks like it works.
Thanks.
Since Lisp is written in Lisp data structures, are there any Lisp
CASE tools that can analyse Lisp code and offer ways to refactor it
to perform fewer operations, clean up the style, and perhaps even
offer better algorithms (that last sounds like AI)?
Such a tool would totally rock. Imagine being able to type in what
you mean and the code analyser being able to make the code much
better for you (probably with guidance from you but perhaps also with
some knowledge of your Lisp implimentation as well). I imagine such
a tool could even help with writing the Lisp portions of Lisp
compilers.
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> writes:
> Peter Seibel <·····@javamonkey.com> writes:
>
>> > David Steuber <·············@verizon.net> writes:
>>
>> > ;;; mandelbrot-escape-fn -- (values escaped iterations)
>> > (defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
>> > "Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
>> > (loop for iterations below max-iterations
>> > with z-real = 0
>> > and z-imag = 0
>> > and z-real-next = 0
>> > and z-real-squared = 0
>> > and z-imag-squared = 0
>> > unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
>> > do (setf z-real-next (+ (- z-real-squared z-imag-squared) c-real))
>> > (setf z-imag (+ (shifted-multiply z-real z-imag *my-two*) c-imag))
>> > (setf z-real z-real-next)
>> > (setf z-real-squared (shifted-multiply z-real z-real))
>> > (setf z-imag-squared (shifted-multiply z-imag z-imag))
>> > finally (return (values nil iterations))))
>>
>> I think you can right this (IMO) more clearly like this (not tested):
>>
>> (defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
>> "Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
>> (loop for iterations below max-iterations
>> for z-imag = 0 then (+ (shifted-multiply z-real z-imag *my-two*) c-imag)
>> for z-real = 0 then (+ (- z-real-squared z-imag-squared) c-real)
>> for z-real-squared = 0 then (shifted-multiply z-real z-real)
>> for z-imag-squared = 0 then (shifted-multiply z-imag z-imag)
>> unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
>> finally (return (values nil iterations)))
>
> You are shy one closing paren (copy paste error?), but it looks
> equivalent providing those steppings are parallel. I always hated
> that z-real-next variable. This also gets back to the brevity of the
> C version. I tried it out and it looks like it works.
Hmmm. I'm not sure what you mean by parallel. I'd call them
"sequential" in the sense that each time through the LOOP each the
assignments are made in the order the "for" clauses appear, and if one
of the later clauses refers to one of the variables in an earlier
clause it will see the new value. Which is the way your SETF version
was used. If you wanted parallel assignment you'd replace all but one
of the "for"'s with "and" but in this case I don't think it matters
because of the way the clauses are ordered.
-Peter
--
Peter Seibel ·····@javamonkey.com
Lisp is the red pill. -- John Fraser, comp.lang.lisp
Peter Seibel <·····@javamonkey.com> writes:
> David Steuber <·············@verizon.net> writes:
>
> > Peter Seibel <·····@javamonkey.com> writes:
> >
> >> > David Steuber <·············@verizon.net> writes:
> >>
> >> > ;;; mandelbrot-escape-fn -- (values escaped iterations)
> >> > (defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
> >> > "Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
> >> > (loop for iterations below max-iterations
> >> > with z-real = 0
> >> > and z-imag = 0
> >> > and z-real-next = 0
> >> > and z-real-squared = 0
> >> > and z-imag-squared = 0
> >> > unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
> >> > do (setf z-real-next (+ (- z-real-squared z-imag-squared) c-real))
> >> > (setf z-imag (+ (shifted-multiply z-real z-imag *my-two*) c-imag))
> >> > (setf z-real z-real-next)
> >> > (setf z-real-squared (shifted-multiply z-real z-real))
> >> > (setf z-imag-squared (shifted-multiply z-imag z-imag))
> >> > finally (return (values nil iterations))))
> >>
> >> I think you can right this (IMO) more clearly like this (not tested):
> >>
> >> (defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations (1- most-positive-fixnum)))
> >> "Return (values escaped iterations) for c-real + c-imag i as integers left shifted by *fractional-bits*"
> >> (loop for iterations below max-iterations
> >> for z-imag = 0 then (+ (shifted-multiply z-real z-imag *my-two*) c-imag)
> >> for z-real = 0 then (+ (- z-real-squared z-imag-squared) c-real)
> >> for z-real-squared = 0 then (shifted-multiply z-real z-real)
> >> for z-imag-squared = 0 then (shifted-multiply z-imag z-imag)
> >> unless (< (+ z-real-squared z-imag-squared) *escape-value*) return (values t iterations)
> >> finally (return (values nil iterations)))
> >
> > You are shy one closing paren (copy paste error?), but it looks
> > equivalent providing those steppings are parallel. I always hated
> > that z-real-next variable. This also gets back to the brevity of the
> > C version. I tried it out and it looks like it works.
>
> Hmmm. I'm not sure what you mean by parallel. I'd call them
> "sequential" in the sense that each time through the LOOP each the
> assignments are made in the order the "for" clauses appear, and if one
> of the later clauses refers to one of the variables in an earlier
> clause it will see the new value. Which is the way your SETF version
> was used. If you wanted parallel assignment you'd replace all but one
> of the "for"'s with "and" but in this case I don't think it matters
> because of the way the clauses are ordered.
Looking at the code again, my original version seems to use
z-real-next as a totally superfluous variable. It looks like I made
the same mistake in my C code as well.
What's important is that Z^2 looks like this:
(a+bi)(a+bi) = (aa - bb) + (2ab)i
--
Those who do not remember the history of Lisp are doomed to repeat it,
badly.
> (dwim x)
NIL
David Steuber <·············@verizon.net> wrote in message news:<··············@david-steuber.com>...
>
> What's important is that Z^2 looks like this:
>
> (a+bi)(a+bi) = (aa - bb) + (2ab)i
For high-precision calculations, it may be useful to organize this as
(a+bi)(a+b1)= (a+b)(a-b) + (2ab)i
and trade off one multiplication for one addition.
Brad
······@math.purdue.edu (Brad Lucier) writes:
> David Steuber <·············@verizon.net> wrote in message news:<··············@david-steuber.com>...
> >
> > What's important is that Z^2 looks like this:
> >
> > (a+bi)(a+bi) = (aa - bb) + (2ab)i
>
> For high-precision calculations, it may be useful to organize this as
>
> (a+bi)(a+b1)= (a+b)(a-b) + (2ab)i
>
> and trade off one multiplication for one addition.
That looks like it might be a good idea. I'm just not sure how to put
it into this context:
(defun mandelbrot-escape-fn (c-real c-imag &key (max-iterations 1000))
"Return iterations for c-real + c-imag i as integers left shifted by *fractional-bits*"
(loop for iterations below max-iterations
for z-imag = 0 then (+ (shifted-multiply-and-double z-real z-imag) c-imag)
for z-real = 0 then (+ (- z-real-squared z-imag-squared) c-real)
for z-real-squared = 0 then (shifted-multiply z-real z-real)
for z-imag-squared = 0 then (shifted-multiply z-imag z-imag)
unless (< (+ z-real-squared z-imag-squared) *escape-value*) return iterations
finally (return iterations)))
I am using aa and bb to decide whether to bail out or not.
The complete code in its current state is here (for now):
http://www.david-steuber.com/~david/fractal-movie/fractal.lisp
I'm not spending a lot of time hacking at it right now. I've been
busy with less fun things. I am also thinking about another Lisp
project which will probably be more useful and also help me learn
other aspects of the language.
--
It would not be too unfair to any language to refer to Java as a
stripped down Lisp or Smalltalk with a C syntax.
--- Ken Anderson
http://openmap.bbn.com/~kanderso/performance/java/index.html
In article <··············@david-steuber.com>,
David Steuber <·····@david-steuber.com> wrote:
>······@math.purdue.edu (Brad Lucier) writes:
>
>> David Steuber <·············@verizon.net> wrote in message
>news:<··············@david-steuber.com>...
>> >
>> > What's important is that Z^2 looks like this:
>> >
>> > (a+bi)(a+bi) = (aa - bb) + (2ab)i
>>
>> For high-precision calculations, it may be useful to organize this as
>>
>> (a+bi)(a+b1)= (a+b)(a-b) + (2ab)i
>>
>> and trade off one multiplication for one addition.
>
>I am using aa and bb to decide whether to bail out or not.
Yes, I realized this later; but you don't really need a*a+b*b to high
accuracy to know when to bail, you can shift a and b right a lot (maybe
so that a*a+b*b is fixnum size) and compare *this* to 4 (shifted appropriately).
There will be a few times that you do an extra iteration because you're
underestimating a^2+b^2, but otherwise you should go faster.
Brad
David Steuber <·············@verizon.net> writes:
[snip]
> SLIME doesn't seem to give me the ability to interrupt the process.
> Maybe it's a simple case of RTFM if there was one. I found myself
> killing the Lisp (OpenMCL) more than once and restarting SLIME.
>
> I tried C-g.
The key is C-c C-g (or C-c C-b in XEmacs). If you can reproduce
problems with interrupting, please send a bug report to the
slime-devel list. Working interrupts are very high on our priority
list and we'll try to fix them.
Helmut.