; 18-bit math
;
; The 18-bit evaluation routines represent numbers from -2.0 to +2.0 using
; a fixed-point format containing one sign bit, 14 "0" bits, one bit to the
; left of the binary point, and 16 bits to the right of the binary point.
; This format takes 32 bits to store in memory or a register, but because
; 14 bits are zero, there are actually only 18 bits of precision.
;
; This format is equivalent to Apple's Fixed format when the absolute value
; of the number being represented is less than 2.  Because the integer part
; of the number is always either 0 or 1, most of the partial products can
; be computed without using multiply instructions.  However, the speed
; advantage over the 30-bit routines is minimal.
;
; See the file ÒMZ-Math16Ó for an explanation of the other number formats.

; Macro to compute the square of an 18-bit quantity.
;
macro	square_18 s,ans,tmp =
	tst.l {s}
	bpl @1
	neg.l {s}

@1	move {s},{tmp}
	mulu {s},{tmp}
	swap {tmp}

	move.l {s},{ans}
	move {tmp},{ans}

	move.l {s},{tmp}
	clr {tmp}
	swap {tmp}
	neg {tmp}
	and {s},{tmp}
	add.l {tmp},{ans}
	add.l {tmp},{ans}
	|

; Subroutine to multiply one 18-bit value by another.
;
; factor 1		a	b
; times factor 2     x	c	d
;		     ------------
; b * d			p	q
; a * d			s
; b * c			u
; a * c	     +	w
;	     --------------------
; product	w	p+s+u	q
;
_long_mult_18
	tst.l d0	; Make both factors positive and remember sign
	smi d4		; of product.
	bpl.s @1
	neg.l d0
@1	tst.l d1
	bpl.s @2
	not.b d4
	neg.l d1
			; D0	D1	D2	D3	D4
			;------	-------	-------	-------	-------
@2	move d0,d2	; a:b	c:d	*:b
	mulu d1,d2	; a:b	c:d	p:q
	swap d2		; a:b	c:d	q:p

	move.l d0,d3	; a:b	c:d	q:p	a:b
	and.l d1,d3	; a:b	c:d	q:p	w:*
	move d2,d3	; a:b	c:d	q:p	w:p

	move.l d0,d2	; a:b	c:d	a:b	w:p
	clr d2		; a:b	c:d	a:0	w:p
	swap d2		; a:b	c:d	0:a	w:p
	neg d2		; a:b	c:d	0:[a]	w:p
	and d1,d2	; a:b	c:d	0:s	w:p
	add.l d2,d3	; a:b	c:d	0:s	w:ps

	clr d1		; a:b	c:0	0:s	w:ps
	swap d1		; a:b	0:c	0:s	w:ps
	neg d1		; a:b	0:[c]	0:s	w:ps
	and d0,d1	; a:b	0:u	0:s	w:ps
	add.l d1,d3	; a:b	0:u	0:s	w:psu

	tst.b d4	; Set sign of product correctly.
	beq.s _lm1_exit
	neg.l d3
_lm1_exit
	rts

; Subroutine to compute Z', given Z and C   (18 bit version) 
;
;               Z = a + b i
;
;                     2         2    2
;               Z' = Z + C = ( a  - b  + Cr ) + ( 2ab + Ci ) i
;
;			; D0	D1	D2	D3	D4	D5	D6	A0	A1	A2	A3
; At startup		;		a	b				Cr	Ci
; At finish		; *	*	a'	b'	*	a'		Cr	Ci	a	b
;			;------	-------	-------	-------	-------	-------	-------	-------	-------	-------	-------
; At startup		;		a	b				Cr	Ci
_z_prime_18 
	move.l d2,a2	;		a	b				Cr	Ci	a
	move.l d3,a3	;		a	b				Cr	Ci	a	b
_z1_1	;
     square_18 d3,d0,d1	; b^2	*	a	*				Cr	Ci	a	b
_z1_2	;
     square_18 d2,d5,d1	; b^2	*	*			a^2		Cr	Ci	a	b
_z1_3	;
	sub.l d0,d5	; b^2					a^-b^		Cr	Ci	a	b
			
	add.l a0,d5	; b^2					a'	 	Cr	Ci	a	b
	compare_l d5,#$20000
	bge.s _z1_err
	compare_l d5,#$-20000
	ble.s _z1_err
	
			; D0	D1	D2	D3	D4	D5	D6	A0	A1	A2	A3
			;------	-------	-------	-------	-------	-------	-------	-------	-------	-------	-------
	move.l a2,d0	; a					a'		Cr	Ci	a	b
	move.l a3,d1	; a	b				a'	 	Cr	Ci	a	b
      jsr _long_mult_18	; *	*	*	ab	*	a'		Cr	Ci	a	b

	add.l d3,d3	;  			2ab		a'		Cr	Ci	a	b

	add.l a1,d3	;			b'		a'		Cr	Ci	a	b
	compare_l d3,#$20000
	bge.s _z1_err
	compare_l d3,#$-20000
	ble.s _z1_err

	move.l d5,d2	;		a'	b'		a'		Cr	Ci	a	b
_z1_exit moveq #0,d0
	rts

_z1_err	moveq #1,d0	; Set overflow condition
	rts

; 18-bit version of _in_set.
;
; Reg.	Input				Output
; ----	-------------------------	------
;  D0	Real component of point		Unchanged.
;	   (Ls2.29)
;  D1	Imag. component of point	Unchanged.
;	   (Ls2.29)
;  D2	   used as longint real component of Z
;  D3	   used as longint imag. component of Z
;  D4	   used as longint real component of Z'
;  D5	   used as longint imag. component of Z'
;  D6	   used as loop counter
;  D7	Ignored.			Number of iterations it took.
;					(If this is -1, assume that the
;					point is in the Mandelbrot Set.)
_in_set_18 movem.l d0-d6/a0-a4,-(sp)
	moveq #0,d2	; Initial Z is 0.
	moveq #0,d3
	move n_max,d7
	asr.l #8,d0	; Shift right 13 bits to put C in desired format
	asr.l #5,d0	; (which is 16 frac. bits)
	move.l d0,a0
	asr.l #8,d1
	asr.l #5,d1
	move.l d1,a1

_is18_loop 
	jsr _z_prime_18
	tst d0
	bne _is18_exit

_is18_next dbf d7,_is18_loop
	move #-1,d7		; -1 means it's in the set.
	bra _is18_rts

_is18_exit sub n_max,d7		; Compute n = n_max - d7
	neg d7

	move.l a2,d2		; Get old a and b
	move.l a3,d3
_is18_z1				; Compute a^2 + b^2 for old Z
	square_18 d3,d0,d1	; d0 = b^2 (28 frac. bits)
_is18_z2
	square_18 d2,d3,d1	; d3 = a^2 (28 frac. bits)
_is18_z3
	add.l d3,d0		; d0 = MS word of a^2 + b^2, 28 frac. bits
	compare_l d0,#$40000	; Is it less than 4.0000 ?
	blt.s _is18_rts		; If so, it was still inside the circle.
	subq #1,d7	; Otherwise Z was outside the circle; adjust count
			; to make countour lines smooth.
_is18_rts movem.l (sp)+,d0-d6/a0-a4
	rts