maclab@reed.UUCP (S.Gillespie-Mac DLab) (05/26/87)
[MandelBrot Routine]
For anyone interested in creating their own Mandelbrot set generating
program, the following is a short assembly routine which uses the Mac
II's hardware floating point processor to determine whether a given
point on the complex plane is in the Mandelbrot set. It is quite fast
(we have found that a full Mac II screen with 150 iterations, and with,
say 20% black pixels, takes about a half hour). You'll need an MPW to
assemble this. You'll also need to figure out how to write a user
interface to wrap around this routine...
Scott Gillespie
Reed College
{decvax,ihnp4, etc}!tektronix!reed!maclab
; ---------------------------------------------------------------
;
; MPW Asm Mandelbrot function using 68881 processor
; by S. Gillespie @ Reed College
;
MC68881
MandelBrot Proc Export
;
;
;
; Function MandelBrot(LoopMax: Integer; c: Float[2]): Integer;
;
; Note: Float = Byte[12] -- a 96-bit Extended 68881 float.
; All float operations assume Extended floating point
; variables.
;
; This function does the necessary calculation to determine
; if a given complex number (x + iy) is in the Mandelbrot set,
; depending on the number of iterations requested.
;
; Returns 0 if complex number c is IN the set at iteration
; LoopMax. Otherwise, returns the number of the iteration
; at which the point blows up.
;
; Given:
; LoopMax: Maximum number of iterations to perform
; c: The complex number (c[0] + i*c[1]) to be tested
; x: Temporary complex number (a,i*b)
; Modsquare(x) = a^2 + b^2
; j: Loop variable
;
; Here is a pseudo-Pascal translation of the routine (note,
; however, that the asm routine handles the loop counter
; a little differently):
;
; Begin
; x := c;
; j := 0;
;
; Repeat
; x := x^2 + c;
; j := j+1;
; Until (j=LoopMax) or (Modsqare(x)>4);
;
; If Modsquare(x)>4 Then
; Mandelbrot := j
; Else
; Mandelbrot := 0;
; End;
;
; --------------------------------------------
;
imoff equ 12 ; Offset to imaginary part of the number
Creal equ fp0 ; This reg holds the real part of c
Cimag equ fp1 ; " " " " imaginary part of c
Xreal equ fp2 ; Real part of x
Ximag equ fp3 ; Imaginary part of x
Move.l (sp)+,a1 ; Return address
Move.l (sp)+,a0 ; Pointer to c
Clr.l d0
Move.w (sp)+,d0 ; LoopMax
Subq.w #1,d0 ; LoopCounter in d0
Clr.w (sp) ; Clear function result location
Suba.w #12,sp ; Allocate space for temporary
; float on the stack. This float
; will be used in the Modsquare test,
; so it is set equal to 4 in the next two
; lines.
FMove.w #4,fp0
FMove.x fp0,(sp)
FMove.x (a0),Creal ; Load up the fp registers with c
FMove.x imoff(a0),Cimag
FMove.x Creal,Xreal ; x := c
FMove.x Cimag,Ximag
Clr.l d1 ; keeps track of no. of iterations
Loop
Addq.w #1,d1 ; bump iteration count
; Now, must do the x := x^2 + c calculation,
; so load up some temporary registers for the operations
FMove.x Xreal,fp4
FMove.x Ximag,fp5
FMove.x Xreal,fp6
FMove.x Ximag,fp7
; for complex number x = (a,i*b), x^2 = (a^2-b^2,i*2ab)
FMul.x fp4,fp4 ; a^2
FMul.x fp5,fp5 ; b^2
FMove.x fp4,Xreal
FSub.x fp5,Xreal ; a^2-b^2
FMul.x fp6,fp7 ; a*b
FMul.w #2,fp7 ; 2*ab
FMove.x fp7,Ximag
FAdd.x Creal,Xreal ; Add c to x^2
FAdd.x Cimag,Ximag
; Now, do the Modsquare test...
FMove.x Xreal,fp4
FMove.x Ximag,fp5
FMul.x fp4,fp4 ; a^2
FMul.x fp5,fp5 ; b^2
Fadd.x fp4,fp5 ; a^2 + b^2
; The stack pointer points to temporary float which
; equals 4.
FCmp.x (sp),fp5
Fdbgt d0,Loop ; Break if loop counter goes negative
; or if modsquare(x)>4
Fbgt NotInSet ; If not > 4, return blowup value
Adda.w #12,sp
Bra.s GetOut
NotInSet
Adda.w #12,sp
Move.w d1,(sp) ; Return value...
GetOut
Jmp (a1)
endp
End