Mandlebrot set

[Robert J Woodhead]

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 at biar.UUCP
  ``I can read your mind - right now, you're thinking I'm full of it...''