; 				MZ-Math16
;
;    The files ÒMZ-Math16Ó, ÒMZ-Math18Ó, ÒMZ-Math30Ó, and ÒMZ-Math32Ó
; contain all of the evaluation routines.  Each file contains the routines
; for computing the Mandelbrot Set calculations for one of the four
; levels of precision offered to the user by the program.  The 16-bit
; version is the simplest, and only uses one subroutine.  Each of the
; other three uses a set of subroutines and macros.
;
;
;			     NUMBER FORMATS
;
; There are several different integer and fixed-point formats used in
; this program.  The formats are described here.
;
; In the diagrams, the following symbols are used:
;
;	s	Sign bit
;	i	Integer bit (to the left of the binary point)
;	f	Fraction bit (to the right of the binary point)
;	0	this bit must be zero
;	x	this bit is a "don't-care"
;
; Note also that the binary point is always assumed to be immediately
; to the right of the rightmost "i" bit.
;
; Ws15 :
; +-+-----------------------------+
; |s|i i i i i i i i i i i i i i i|
; +-+-----------------------------+
;    This is an ordinary 16-bit signed integer.  It is used for
; various counters, QuickDraw coordinates, etc.
;
;
; Ls31 :
; +-+-------------------------------------------------------------+
; |s|i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i i|
; +-+-------------------------------------------------------------+
;    32-bit signed integer.
;
; Ws1.14 :
; +-+-+---------------------------+
; |s|i|f f f f f f f f f f f f f f|
; +-+-+---------------------------+
;    This is the format used internally by the 16-bit evaluation
; routines.
;
; Ls3.28 :
; +-+-----+-------------------------------------------------------+
; |s|i i i|f f f f f f f f f f f f f f f f f f f f f f f f f f f f|
; +-+-----+-------------------------------------------------------+
;    This format is used
;
; Ls2.29 :
; +-+---+---------------------------------------------------------+
; |s|i i|f f f f f f f f f f f f f f f f f f f f f f f f f f f f f|
; +-+---+---------------------------------------------------------+
;    This is the input format for all of the evaluation routines.
;
; Ls1.30 :
; +-+-+-----------------------------------------------------------+
; |s|i|f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f|
; +-+-+-----------------------------------------------------------+
;

; 16-bit math
;
; The 16-bit evaluation routine represents numbers from -2.0 to +2.0 using
; a fixed-point format containing one sign bit, one bit to the left of
; the binary point, and 14 bits to the right of the binary point.

; Subroutine to determine if a given point is in the Mandelbrot Set, or
; if it isn't, how many interations of  Z'=Z^2+p  it took to get a value
; more than 2 units from the origin.  NOTE -- the test for |Z|>2 is not
; performed every time; instead, in the interest of speed, |Zr| and |Zi|
; are compared to 2 in the loop, causing the loop to exit when Z is outside
; the square.  Then Z(n-1) is checked retrospectively to see if it was
; inside or outside the circle.  This gives us a reasonable but imperfect
; approximation, and it means that although the contours are smooth, it is
; possible to find small angular protrusions in some places.  These are
; found in areas which map onto the left-hand edge of the square (those
; points with real coordinates less than -1 .)
;
;               Z = a + b i
;
;                     2         2    2
;               Z' = Z + C = ( a  - b  + Cr ) + ( 2ab + Ci ) i
;
; Reg.	Input					Output
; ----	-------------------------		------
; D0.l	Real component of C (Ls2.29)		Unchanged.
; D1.l	Imag. component of C (Ls2.29)		Unchanged.
; D2	   used as real component of Z		trashed
; D3	   used as imag. component of Z		trashed
; D4	   used as real component of Z'		trashed
; D5	   used as imag. component of Z'	trashed
; D6	Unused.					Unchanged.
; D7.w	Ignored.			Number of iterations it took.
;					(If this is -1, assume that the
;					point is in the Mandelbrot Set.)
; A0.l	   used as ptr to array
;
_in_set moveq #0,d7	; Count will be zero if one of these first two bvs's
			; is taken.
	move.l d0,d2	; Initial Z is C.
	add.l d2,d2	; Ls1.30
	bvs.s _is_rts
	swap d2		; Ws1.14
	move.l d1,d3
	add.l d3,d3	; Ls1.30
	bvs.s _is_rts
	swap d3		; Ws1.14
	
	move n_max,d7	; Set up counter
	subq #1,d7	; Decrement one because of initial Z=C instead of
			; Z=0
	
_is_loop move d2,d4	; Calculate  a^2 - b^2 + Cr  term.
	muls d4,d4	; D4 = a^2

	move d3,d5
	muls d5,d5	; D5 = b^2

	sub.l d5,d4	; D4 = a^2 - b^2 (Ls3.28)
	add.l d4,d4	; (Ls2.29, standard format for C)
	add.l d0,d4	; D4 = a^2 - b^2 + Cr
	bvs.s _is_exit
	add.l d4,d4	; Ls1.30
	bvs.s _is_exit	; check for overflow
	swap d4		; Ws1.14

			; Calculate  2*a*b + Ci  term.
	move d3,d5	; D5 = b
	muls d2,d5	; D5 = a*b (Ls3.28)
	asl.l #2,d5	; D5 = 2*a*b (Ls2.29)
	bvs.s _is_exit	; [ Leave out this bvs to get neat patterns near
			; boundaries... ]
	add.l d1,d5	; Add Ci
	bvs.s _is_exit
	add.l d5,d5	; Ls1.30
	bvs.s _is_exit
	swap d5		; Ws1.14

	move d4,d2	; Set  Z = Z'
	move d5,d3

_is_next dbf d7,_is_loop

_is_inset move #-1,d7
	bra _is_rts

_is_exit sub n_max,d7	; Compute D7 = n_max - d7
	neg d7

	move d2,d4	; Compute   a^2 + b^2  (for previous Z)
	muls d4,d4
	move d3,d5
	muls d5,d5
	add.l d5,d4	; Ls3.28
	cmp.l #$40000000,d4	; Compare to +4.0 in Ls3.28 format.
	blt @1		; If less, it was still inside the circle.
	subq #1,d7	; Otherwise Z was outside the circle; adjust count
			; to make countour lines smooth.
@1	
_is_rts	rts

; Subroutine to convert an Ls3.28 number to a string.
;
; Reg.	Input				Output
; ----	------------------------------	------
; D0.L	The number (Ls3.28)		Trashed.
; D1	  Used as temp. register.	Unchanged.
; D2	  Used as temp counter		Unchanged.
; D3	  Used to hold a constant (28)	Unchanged.
; D4	Number of characters that the	Unchanged.
;         answer will have.  This is 7 
;	  for 16-bit input or 11 for
;	  32-bit input.
; A0	Pointer to buffer for string.	Unchanged.
; A1.L	  Used as temp pointer		Unchanged.
;
dc.w 0
_fixed_to_str 
	movem.l d1-d4/a0-a1,-(sp)
	move.l a0,a1		; remember pointer to length byte.
	moveq #3,d2		; init every-third-digit counter.

	moveq #28,d3
	move.b d4,(a0)+		; Store length
	tst.l d0
	bpl _fts_1
	move.b #'-',(a0)+
	neg.l d0
	bra _fts_2

_fts_1	move.b #'+',(a0)+

_fts_2	move.l d0,d1
	lsr.l d3,d1
	add #'0',d1
	move.b d1,(a0)+
	and.l #$0fffffff,d0

	move.b #'.',(a0)+

	subq #4,d4	; Set up d4 for loop.
_fts_3	add.l d0,d0	; d0=2*d0
	move.l d0,d1
	lsl.l #2,d1	; d1=8*d0
	add.l d1,d0	; d0=10*d0
	move.l d0,d1
	lsr.l d3,d1
	add #'0',d1
	move.b d1,(a0)+
	
	subq #1,d2	; Every third time, write out an extra
	bne @1		; space.
	move.b #' ',(a0)+
	addq.b #1,(a1)
	moveq #3,d2
@1
	and.l #$0fffffff,d0
	dbf d4,_fts_3

	movem.l (sp)+,d1-d4/a0-a1
	rts

temp_string	dc.l 0,0,0,0,0,0,0,0,0,0
; Subroutine to convert number to string and draw on the screen.
;
; Reg.	Input				Output
; ----	------------------------------	------
;  D0	the number (16-bit signed)	unchanged.
;
_draw_number 
	movem.l d0/a0,-(sp)
	lea temp_string, a0
	jsr _numtostr
	DrawString temp_string
	movem.l (sp)+,d0/a0
	rts

; Subroutine to convert a 16-bit signed integer to a string.
;
; Reg.	Input				Output
; ----	------------------------------	------
;  D0	the number (16-bit signed)	unchanged.
;  D1      (Used as temp)		unchanged.
;  A0	pointer to buffer for string.	unchanged.
;  A1      (Used as temp)		unchanged.
_numtostr 
	movem.l d0/d1/a0/a1,-(sp)
	clr.w d1	; D1 will be nonzero if number is negative.
	ext.l d0
	tst.l d0
	bpl _nts1
	neg.l d0
	addq.w #1,d1	; Set flag.
_nts1	clr -(sp)	; Push a zero - this will be the end of the string.
_nts2	divs #10,d0	; Divide by 10 and get remainder.
	swap d0
	add.w #'0',d0	; Convert to ASCII and push on stack.
	move.w d0,-(sp)
	clr.w d0
	swap d0
	tst d0
	bne _nts2
	tst d1		; If it was negative, put a '-' on the stack.
	beq _nts3
	move.w #'-',-(sp)
_nts3	clr.w d0	; Use D0 to count number of characters.
	move.l a0,a1	; Remember starting address of destination buffer.
	addq #1,a0	; Skip place where length field will go.
_nts5	move.w (sp)+,d1
	beq _nts4	; Zero means end of string.
	addq.w #1,d0
	move.b d1,(a0)+
	bra _nts5
_nts4	move.l a1,a0	; Restore A0
	move.b d0,(a0)	; Set length field of string.
	movem.l (sp)+,d0/d1/a0/a1
	rts

; Subroutine to convert a string to a number.  The string is a Pascal
; string, starting with an optional sign bit, one digit (optional),
; a period, and an arbitrary number of fraction digits.  Spaces are
; ignored; other illegal characters result in an error return from this
; routine.
;
; Reg.	Input					Output
; ----	------------------------------		------
; D0.L	ignored.				The number (Ls2.29)
; D1.W	ignored.				Zero unless format error.
; D2	  (Used for current char)		Unchanged.
; D3	  (Used or negative flag)		Unchanged.
; D4	  (Used for integer part)		Unchanged.
; D5	  (Used for accumulating fraction)	Unchanged.
; D6	  (Used for frac. digits counter)	Unchanged.
; D7	  (Used for chars-left flag)		Unchanged.
; A0.L	pointer to the Pascal string.		Unchanged.
; A1.L	  (Used for fractional part pointer)	Unchanged.
;
_string_to_fixed 
	movem.l d2-d7/a0-a4,-(sp)
	
	moveq #0,d4
	
	moveq #0,d7
	move.b (a0)+,d7		; Get length byte
	clr d3			; clear minus flag
	lea stf_table,a1
	moveq #8,d6		; Allow at most 9 frac. digits
	moveq #0,d5		; Init frac. part to zero
	
	jsr _stf_get		; Try to get '+','-', digit, or '.'
	beq _stf_error
	compare d2,#'+'
	beq _stf_1		; plus found - try to get first digit
	compare d2,#'-'
	bne _stf_2		; minus not found - check for digit
	moveq #1,d3		; set minus flag
	
_stf_1	jsr _stf_get		; Try to get digit or '.'
	beq _stf_error
_stf_2	compare d2,#'.'		; Dot found, look for fractional digits
	beq _stf_5
	sub #'0',d2		; Check for less than zero
	blt _stf_error
	compare d2,#4		; Integer part may not exceed four
	bgt _stf_error
	move.b d2,d4		; Remember integer part.

_stf_4
	jsr _stf_get		; Try to get '.'
	beq _stf_okay
_stf_45	compare d2,#'.'
	bne _stf_65		; No dot found, check for digit
	
_stf_5				; '.' found - look for digits

_stf_6	jsr _stf_get		; Get a char
	beq _stf_okay		; Last one - return OK.
_stf_65
	sub #'0',d2		; Check for less than zero
	blt _stf_error
	compare d2,#9		; Check for greater than nine
	bgt _stf_error
	
_stf_7	add.l (a1),d5		; Add 1/(10^n) to accumulating fraction.
	dbf d2, _stf_7
	sub.l (a1)+,d5		; Subtract to make up for extra add due to
				; DBF loop; bump pointer to components table.
	
	dbf d6,_stf_6

_stf_okay			; We're done -- create the final answer.
	move.l d5,d0		; This is most of the answer.
	moveq #29,d1		; Shift d4 left 29 bits
	asl.l d1,d4
	add.l d4,d0		; Add in integer part.
	tst d3			; Check for minus
	beq @1
	neg.l d0
@1
	moveq #0,d1		; No error return
	movem.l (sp)+,d2-d7/a0-a4
	rts

_stf_error
	moveq.l #0,d0
	moveq #1,d1
	movem.l (sp)+,d2-d7/a0-a4
	rts

stf_table
	dc.l	53687091, 5368709, 536871, 53687, 5369, 537, 54, 5, 1, 0

; Subroutine of _string.to_num - gets next character, filters out
; spaces.
;
_stf_get
	moveq #0,d2	; Clear high bytes of d2
	move.b (a0)+,d2	; Get char.
	subq #1,d7	; Decrement chars-left count.
	bge @1
	moveq #0,d2	; No chars left.
@1	compare d2,#32	; If space, get another.
	beq _stf_get
	tst d2		; So caller can branch on end-of-string.
	rts

; Subroutine to convert a string to a 32-bit signed integer.  The string
; is a Pascal string, starting with an optional sign bit and up to ten
; digits.  Spaces are ignored; other illegal characters result in an error
; return from this routine.
;
; Reg.	Input					Output
; ----	------------------------------		------
; D0.L	ignored.				The number.
; D1.W	ignored.				Zero unless format error.
; D2	  (Used for current char)		Unchanged.
; D3	  (Used or negative flag)		Unchanged.
; D4.L	  (Used for temp.)			Unchanged.
; D7	  (Used for chars-left flag)		Unchanged.
; A0.L	pointer to the Pascal string.		Unchanged.
;
_string_to_long 
	movem.l d2-d7/a0-a4,-(sp)
	
	moveq #0,d0
	
	moveq #0,d7
	move.b (a0)+,d7		; Get length byte
	clr d3			; clear minus flag
	
	jsr _stf_get		; Try to get '+','-', or digit.
	beq _stl_error
	compare d2,#'+'
	beq _stl_1		; plus found - try to get first digit
	compare d2,#'-'
	bne _stl_2		; minus not found - check for digit
	moveq #1,d3		; set minus flag
	
_stl_1	jsr _stf_get		; Try to get digit.
	beq _stl_error
	
_stl_2	sub #'0',d2		; Check for less than zero
	blt _stl_error
	compare d2,#9		; Check for greater than nine.
	bgt _stl_error

_stl_3
	add.l d0,d0		; Multiply d0 by ten
	bvs _stl_error
	move.l d0,d4
	add.l d0,d0
	bvs _stl_error
	add.l d0,d0
	bvs _stl_error
	add.l d4,d0
	bvs _stl_error
	ext.w d2		; Add latest digit.
	ext.l d2
	add.l d2,d0
	
	jsr _stf_get		; Get next char
	beq _stl_okay		; Last one - return OK.
	sub #'0',d2		; Check for less than zero
	blt _stl_error
	compare d2,#9		; Check for greater than nine
	bgt _stl_error
	bra _stl_3		; Add it in to the number and get another char.
	
_stl_okay			; We're done -- create the final answer.
	tst d3			; Check for minus
	beq @1
	neg.l d0
@1
	moveq #0,d1		; No error return
	movem.l (sp)+,d2-d7/a0-a4
	rts

_stl_error
	moveq.l #0,d0
	moveq #1,d1
	movem.l (sp)+,d2-d7/a0-a4
	rts