dos33fsprogs/asm_routines/multiply_fast.s

352 lines
6.7 KiB
ArmAsm

; Fast mutiply
; Note for our purposes we only care about 8.8 x 8.8 fixed point
; with 8.8 result, which means we only care about the middle two bytes
; of the 32 bit result. So we disable generation of the high and low byte
; to save some cycles.
;
; The old routine took around 700 cycles for a 16bitx16bit=32bit mutiply
; This routine, at an expense of 2kB of looku tables, takes around 250
; If you reuse a term the next time this drops closer to 200
; This routine was described by Stephen Judd and found
; in The Fridge and in the C=Hacking magazine
; http://codebase64.org/doku.php?id=base:seriously_fast_multiplication
; The key thing to note is that
; (a+b)^2 (a-b)^2
; a*b = ------- - --------
; 4 4
; So if you have tables of the squares of 0..511 you can lookup and subtract
; instead of multiplying.
; Table generation: I:0..511
; square1_lo = <((I*I)/4)
; square1_hi = >((I*I)/4)
; square2_lo = <(((I-255)*(I-255))/4)
; square2_hi = >(((I-255)*(I-255))/4)
; Note: DOS3.3 starts at $9600
.ifndef square1_lo
square1_lo EQU $8E00
square1_hi EQU $9000
square2_lo EQU $9200
square2_hi EQU $9400
.endif
; for(i=0;i<512;i++) {
; square1_lo[i]=((i*i)/4)&0xff;
; square1_hi[i]=(((i*i)/4)>>8)&0xff;
; square2_lo[i]=( ((i-255)*(i-255))/4)&0xff;
; square2_hi[i]=(( ((i-255)*(i-255))/4)>>8)&0xff;
; }
init_multiply_tables:
; Build the add tables
ldx #$00
txa
.byte $c9 ; CMP #immediate - skip TYA and clear carry flag
lb1: tya
adc #$00 ; 0
ml1: sta square1_hi,x ; square1_hi[0]=0
tay ; y=0
cmp #$40 ; subtract 64 and update flags (c=0)
txa ; a=0
ror ; rotate
ml9: adc #$00 ; add 0
sta ml9+1 ; update add value
inx ; x=1
ml0: sta square1_lo,x ; square1_lo[0]=1
bne lb1 ; if not zero, loop
inc ml0+2 ; increment values
inc ml1+2 ; increment values
clc ; c=0
iny ; y=1
bne lb1 ; loop
; Build the subtract tables based on the existing one
ldx #$00
ldy #$ff
second_table:
lda square1_hi+1,x
sta square2_hi+$100,x
lda square1_hi,x
sta square2_hi,y
lda square1_lo+1,x
sta square2_lo+$100,x
lda square1_lo,x
sta square2_lo,y
dey
inx
bne second_table
rts
; Fast 16x16 bit unsigned multiplication, 32-bit result
; Input: NUM1H:NUM1L * NUM2H:NUM2L
; Result: RESULT3:RESULT2:RESULT1:RESULT0
;
; Does self-modifying code to hard-code NUM1H:NUM1L into the code
; carry=0: re-use previous NUM1H:NUM1L
; carry=1: reload NUM1H:NUM1L (58 cycles slower)
;
; clobbered: RESULT, X, A, C
; Allocation setup: T1,T2 and RESULT preferably on Zero-page.
;
; NUM1H (x_i), NUM1L (x_f)
; NUM2H (y_i), NUM2L (y_f)
; NUM1L * NUM2L = AAaa
; NUM1L * NUM2H = BBbb
; NUM1H * NUM2L = CCcc
; NUM1H * NUM2H = DDdd
;
; AAaa
; BBbb
; CCcc
; + DDdd
; ----------
; RESULT
;fixed_16x16_mul_unsigned:
multiply:
bcc num1_same_as_last_time ; 2nt/3
;============================
; Set up self-modifying code
; this changes the code to be hard-coded to multiply by NUM1H:NUM1L
;============================
lda NUM1L ; load the low byte ; 3
sta sm1a+1 ; 3
sta sm3a+1 ; 3
sta sm5a+1 ; 3
sta sm7a+1 ; 3
eor #$ff ; invert the bits for subtracting ; 2
sta sm2a+1 ; 3
sta sm4a+1 ; 3
sta sm6a+1 ; 3
sta sm8a+1 ; 3
lda NUM1H ; load the high byte ; 3
sta sm1b+1 ; 3
sta sm3b+1 ; 3
sta sm5b+1 ; 3
; sta sm7b+1 ;
eor #$ff ; invert the bits for subtractin ; 2
sta sm2b+1 ; 3
sta sm4b+1 ; 3
sta sm6b+1 ; 3
; sta sm8b+1 ;
;===========
; 52
num1_same_as_last_time:
;==========================
; Perform NUM1L * NUM2L = AAaa
;==========================
ldx NUM2L ; (low le) ; 3
sec ; 2
sm1a:
lda square1_lo,x ; 4
sm2a:
sbc square2_lo,x ; 4
; a is _aa
; sta RESULT+0 ;
sm3a:
lda square1_hi,x ; 4
sm4a:
sbc square2_hi,x ; 4
; a is _AA
sta _AA+1 ; 3
;===========
; 24
; Perform NUM1H * NUM2L = CCcc
sec ; 2
sm1b:
lda square1_lo,x ; 4
sm2b:
sbc square2_lo,x ; 4
; a is _cc
sta _cc+1 ; 3
sm3b:
lda square1_hi,x ; 4
sm4b:
sbc square2_hi,x ; 4
; a is _CC
sta _CC+1 ; 3
;===========
; 24
;==========================
; Perform NUM1L * NUM2H = BBbb
;==========================
ldx NUM2H ; 3
sec ; 2
sm5a:
lda square1_lo,x ; 4
sm6a:
sbc square2_lo,x ; 4
; a is _bb
sta _bb+1 ; 3
sm7a:
lda square1_hi,x ; 4
sm8a:
sbc square2_hi,x ; 4
; a is _BB
sta _BB+1 ; 3
;===========
; 27
;==========================
; Perform NUM1H * NUM2H = DDdd
;==========================
sec ; 2
sm5b:
lda square1_lo,x ; 4
sm6b:
sbc square2_lo,x ; 4
; a is _dd
sta _dd+1 ; 3
;sm7b:
; lda square1_hi,x ;
;sm8b:
; sbc square2_hi,x ;
; a = _DD
; sta RESULT+3 ;
;===========
; 13
;===========================================
; Add the separate multiplications together
;===========================================
clc ; 2
_AA:
lda #0 ; loading _AA ; 2
_bb:
adc #0 ; adding in _bb ; 2
sta RESULT+1 ; 3
;==========
; 9
; product[2]=_BB+_CC+c
_BB:
lda #0 ; loading _BB ; 2
_CC:
adc #0 ; adding in _CC ; 2
sta RESULT+2 ; 3
;===========
; 7
; product[3]=_DD+c
; bcc dd_no_carry1 ;
; inc RESULT+3 ;
clc ; 2
;=============
; 2
dd_no_carry1:
; product[1]=_AA+_bb+_cc
_cc:
lda #0 ; load _cc ; 2
adc RESULT+1 ; 3
sta RESULT+1 ; 3
; product[2]=_BB+_CC+_dd+c
_dd:
lda #0 ; load _dd ; 2
adc RESULT+2 ; 3
sta RESULT+2 ; 3
;===========
; 16
; product[3]=_DD+c
; bcc dd_no_carry2 ;
; inc RESULT+3 ;
;=============
; 0
dd_no_carry2:
; *z_i=product[1];
; *z_f=product[0];
; rts ; 6
;=================
; Signed multiply
;=================
;multiply:
; jsr fixed_16x16_mul_unsigned ; 6
lda NUM1H ; x_i ; 3
;===========
; 12
bpl x_positive ;^3/2nt
sec ; 2
lda RESULT+2 ; 3
sbc NUM2L ; 3
sta RESULT+2 ; 3
; lda RESULT+3 ;
; sbc NUM2H ;
; sta RESULT+3 ;
;============
; 10
x_positive:
lda NUM2H ; y_i ; 3
;============
; ; 6
bpl y_positive ;^3/2nt
sec ; 2
lda RESULT+2 ; 3
sbc NUM1L ; 3
sta RESULT+2 ; 3
; lda RESULT+3 ;
; sbc NUM1H ;
; sta RESULT+3 ;
;===========
; 10
y_positive:
ldx RESULT+2 ; *z_i=product[2]; ; 3
lda RESULT+1 ; *z_f=product[1]; ; 3
rts ; 6
;==========
; 12