themod@blake.acs.washington.edu (Chris Hinton) (06/22/89)
Hello to all you netters out there, maybe one of you can help me. I am in search of a function that, given a point in the complex plane, will return the number of iterations it takes to determine if the given point is in the Mandlebrot (sp?) set. Put more simply, I would call the function in this way : s = function(real,imag); where s is the number of iterations real is the real part of the complex number imag is the imaginary part of the complex number. Any help is greatly appreciated. +---------------------------------+------------------------------------------+ + Chris Hinton (aka The Mod) | The veiws presented in the above posting | + Mod Software Systems | are my own, and should not be take as | + (206)547-1342 | those of the University of Wasington. | + | These ideas are available for sale or | + themod@blake.acs.washington.edu | rent at a reasonable price. | +---------------------------------+------------------------------------------+
gwyn@smoke.BRL.MIL (Doug Gwyn) (06/22/89)
In article <2492@blake.acs.washington.edu> themod@blake.acs.washington.edu (Chris Hinton) writes: >I am in search of a function that, given a point in the complex plane, will >return the number of iterations it takes to determine if the given point is >in the Mandlebrot (sp?) set. Not asking much, are you? A solution to this would probably qualify you for the Fields medal. A more practical requirement would be: Given coordinates of a point in the complex plane and an iteration count limit, determine whether the iterative Mandelbrot computation (z <- z^2+c starting with z=0) produces a z with norm > 2 within the iteration limit.
mpl@cbnewsl.ATT.COM (michael.p.lindner) (06/23/89)
In article <2492@blake.acs.washington.edu>, themod@blake.acs.washington.edu (Chris Hinton) writes: > Hello to all you netters out there, maybe one of you can help me. I am > in search of a function that, given a point in the complex plane, will return > the number of iterations it takes to determine if the given point is in the > Mandlebrot (sp?) set. > Put more simply, I would call the function in this way : > > s = function(real,imag); . . > + Chris Hinton (aka The Mod) | The veiws presented in the above posting | Well, the number of iterations may be infinite, so if you wrote a function like the one you describe, it might never terminate. In general, one specifies a maximum number of iterations to try, after which one assumes the point is in the Mandelbrot (sp!) set. In addition, your function as defined does not say whether or not the point is actually in the set. So, I give you the following, which returns the number of iterations after which it is known the point is NOT in the set (if the value ever goes outside a circle of radius 2 it is not in the set). If the function returns MAX_ITER, the maximum number of iterations to try, assume the point is in the set. double func(real, imag) double real, imag; { int n; double newreal, newimag, re2, im2; for (n = 0; n < MAX_ITER; n++) { if ((re2 = real * real) + (im2 = imag * imag) > 4.0) break; /* point is not in Mandelbrot set */ /* Xn+1 = Xn^2 + Xn */ newreal = re2 - im2 + real; newimag = 2.0 * real * imag + imag; if (newreal == real && newimag == imag) { n = MAX_ITER; break; } } return n; } For a good discussion of the Mandelbrot set, see the January 1989 Scientific American, or read "Chaos: Making a New Science" by James Gleick, published by Penguin Books. Mike Lindner attunix!mpl AT&T Bell Laboratories 190 River Rd. Summit, NJ 07901
burkie@uxe.cso.uiuc.edu (06/26/89)
>double >func(real, imag) >double real, imag; >{ > int n; > double newreal, newimag, re2, im2; > > for (n = 0; n < MAX_ITER; n++) { > if ((re2 = real * real) + (im2 = imag * imag) > 4.0) > break; /* point is not in Mandelbrot set */ > > /* Xn+1 = Xn^2 + Xn */ > newreal = re2 - im2 + real; > newimag = 2.0 * real * imag + imag; > if (newreal == real && newimag == imag) { > n = MAX_ITER; > break; > } > } > return n; >} I think there is something missing in the above code. ----------------------- The general algorithm is: given complex number: c z = 0 { do (z = z ^ 2 + c) until |z| > 2 } <--- don't check beyond MAX_ITER iterations; assume beyond MAX_ITER is for all practical purposes 'stuck forever' return the number of iterations it took for |z| to exceed 2 and if never exceeded 2 during the MAX_ITER iterations that we tried then return MAX_ITER ------------------------ double func(c_real, c_imag) double c_real, c_imag; /* c = (c_real + i*c_imag) */ { int n; /* iteration count */ double z_real, z_imag, /* z = (z_real + i*z_imag) */ re2, im2; /* z_real^2 , z_imag^2 */ z_real = 0; z_real = 0; n = 0; while ( (n < MAX_ITER) && ((re2 = z_real * z_real) + (im2 = z_imag * z_imag) <= 4.0) ) { /* iterate z = z^2 + c */ /* = (z_real + i*z_imag)^2 + (c_real + i*c_imag) */ z_imag = 2.0 * z_real * z_imag + c_imag; z_real = re2 - im2 + c_real; n++; } /* if n < MAX_ITER then |z| exceeded 2 after iteration n * (and 'c' is definitely not in mandelbrot set) * if n = MAX_ITER then we were not able to see it diverge * (so we may assume it is part of mandelbrot set) * (for lack of knowing any better) */ return(n); } ---------------------------- I realize most people already know this but here is a short explanation: the mandelbrot set is the set of c's in the complex plane which never diverge when iterated in the manner described above. (i.e. the z remains finite, oscillating between values, converging to a value etc. destined never to escape beyond the circle) in most pictures of the mandelbrot set - the set is colored black and other areas are colored according to the number of iterations it took for the condition |z| > 2 to be reached using the 'c' value at that point. the condition |z| > 2 guarantees that future iterations of z = z^2 + c will diverge more and more away and so we use it as a convenient condition to test for divergence. however for certain points 'c' it may take a very long time to verify that |z| eventually exceeds 2 and for some this condition may never become true (those points which are truly within the mandelbrot set) So, to avoid waiting forever for this condition, we only test upto MAX_ITER iterations and then if |z| still < 2 , we assume for our purposes that the z never diverges and give up testing further for divergence. (and include 'c' in the mandelbrot set) the more detailed your examination of the boundary of the mandelbrot set, the higher MAX_ITER will have to be to give a good representation. For a given pixel resolution, MAX_ITER beyond a certain value will not increase the apparent accuracy of the picture and will only slow down the program. -Salman burkie@uxe.cso.uiuc.edu
trebor@biar.UUCP (Robert J Woodhead) (06/26/89)
When I wrote the Mac ``MandelColor'' program, I had a little fun writing some assembly language code for the 68881 FPU chip to compute the Mandelbrot set. Here it is. It has the nice feature that it manages to fit all the variables into the 68881 registers, so it is reasonably fast. ; ; function mCalc(p,q,myM:Extended; maxK:integer):integer; ; ; returns mandelbrot value for point p,q with escaped value myM and ; maximum value maxK. Value is number of loops through the computation ; needed before end result is >myM, or maxK, whichever comes first. ; ; maxK should be the # of colors you want to use in your plot, and ; myM is some arbitrary cutoff. I always use 100, as that is the ; number found in Mandelbrot's book. ; mCalc func export ; ; register definitions ; retAddr equ A0 ; pascal return address param equ A1 ; parameter pointer k equ D0 ; loop index maxK equ D1 ; maximum K value p equ FP0 ; mandelbrot p value q equ FP1 ; mandelbrot q value myM equ FP2 ; myM loop top value xv equ FP3 ; x value yv equ FP4 ; y value xv2 equ FP5 ; x value squared yv2 equ FP6 ; y value squared scratch equ FP7 ; temp results register ; ; start of code ; ; stack contents are as follows: ; ; top of stack > return address long pointer ; maxK word ; myM long pointer ; q long pointer ; p long pointer ; function result word ; move.l (sp)+,retAddr ; get return address move.w (sp)+,maxK ; get maxK from the stack move.l (sp)+,param ; get pointer to myM fmove.x (param),myM ; and save it move.l (sp)+,param ; get pointer to q fmove.x (param),q ; and save it move.l (sp)+,param ; get pointer to p fmove.x (param),p ; and save it fmovecr #$0f,xv ; xv := 0 fmovecr #$0f,yv ; yv := 0 fmovecr #$0f,xv2 ; xv2 := 0 fmovecr #$0f,yv2 ; yv2 := 0 moveq.l #0,k ; k := 0 fmove k,fpcr ; set IEEE standards for the math ; ; loop performing the mandelbrot calculations. loop ends when ; either maxK iterations have been performed or the result ; exceeds myM ; mLoop fmul xv,yv ; compute yv:=2*xv*yv+q addq.w #1,k ; while this is going on, increment k and cmp.w maxK,k ; see if we are done beq.s eLoop ; if so, exit fadd yv,yv ; yv:=2*xv*yv fadd q,yv ; yv:=2*xv*yv+q fmove xv2,xv ; compute xv:=xv2-yv2+p fsub yv2,xv fadd p,xv fmove xv,xv2 ; xv2:=xv*xv fmul xv,xv2 fmove yv,yv2 ; yv2:=yv*yv fmul yv,yv2 fmove yv2,scratch ; compute scratch:=yv2+xv2 fadd xv2,scratch fcmp scratch,myM ; if myM>=scratch, reloop fbge mLoop ; ; at this point, k contains our return result ; eLoop move.w k,(sp) ; store function result jmp (retAddr) ; and return to Pascal ENDF -- (^;-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-;^) Robert J Woodhead, Biar Games, Inc. !uunet!biar!trebor | trebor@biar.UUCP ``I can read your mind - right now, you're thinking I'm full of it...''
mpl@cbnewsl.ATT.COM (michael.p.lindner) (06/26/89)
In article <225800193@uxe.cso.uiuc.edu>, burkie@uxe.cso.uiuc.edu writes: > > >double > >func(real, imag) > >double real, imag; > >{ erroneous code deleted > >} > I think there is something missing in the above code. correct code deleted > -Salman > burkie@uxe.cso.uiuc.edu Oops! I seem to have posted the wrong code! Also, my reference to the January 1989 Scientific American should be to the Februaury 1989 Scientific American. Mike Lindner attunix!mpl AT&T Bell Laboratories 190 River Rd. Summit, NJ 07901
chris@ssbell.UUCP (Chris Olson) (06/27/89)
In article <690@biar.UUCP> trebor@biar.UUCP (Robert J Woodhead) writes: >When I wrote the Mac ``MandelColor'' program, I had a little fun writing >some assembly language code for the 68881 FPU chip to compute the [lot's of bandwidth wasted here with assembly code] >-- >(^;-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-;^) >Robert J Woodhead, Biar Games, Inc. !uunet!biar!trebor | trebor@biar.UUCP > ``I can read your mind - right now, you're thinking I'm full of it...'' ^ | No lie. #ifdef FLAME_ON! There must be some mistake here. For a moment there I thought I was reading comp.lang.c. Obviously, it's comp.lang.68881. Come on folks! Let's at least keep the language this group is based on as 'C'. What's next? Mandelbrot's in COBOL? #endif FLAME_ON! Ahem. At least now I feel better... -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- Chris Olson / (aka, Lord Valentine, O Electric Samurai, OXXXXXXI>>>>>>>>>>>>>>>>>>>>>>>>>>>>> & other foolishness) O [uunet!ssbell!chris] \
mfinegan@uceng.UC.EDU (michael k finegan) (07/01/89)
While the calculation of members of the Mandelbrot set may be well known, do you have pseudo-code for calculation of members of a Julia set? I implemented the algorithms outlined in Scientific American (forget issue #) for Mandelbrot and Julia sets; the Mandelbrot set was correct, but the Julia set results were always wrong. I realize this is not comp.graphics, but since you are discussing it here as C code ... Mike Finegan mfinegan@uceng.UC.EDU
gwyn@smoke.BRL.MIL (Doug Gwyn) (07/02/89)
In article <1354@uceng.UC.EDU> mfinegan@uceng.UC.EDU (michael k finegan) writes: > While the calculation of members of the Mandelbrot set may be well known, do >you have pseudo-code for calculation of members of a Julia set? There are all sorts of clever schemes, but the following extract from a program that creates Mandelbrot & Julia convergence displays in the "obvious" way shows that the calculations are quite similar: if ( julia ) { zx = cx; zy = cy; for ( k = 0; k < iters; ++k ) { register double x = zx * zx - zy * zy + jx; zy = 2.0 * zx * zy + jy; zx = x; if ( zx * zx + zy * zy >= thresh ) break; } } else { /* mandelbrot */ zx = zy = 0.0; for ( k = 0; k < iters; ++k ) { register double x = zx * zx - zy * zy + cx; zy = 2.0 * zx * zy + cy; zx = x; if ( zx * zx + zy * zy >= thresh ) break; } }