; 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