From 818707ed8c88fe42573f3af0376bc1a15c80bbe9 Mon Sep 17 00:00:00 2001 From: Stephen Heumann Date: Tue, 21 Dec 2021 19:11:18 -0600 Subject: [PATCH] Use a more accurate implementation of cbrt(). The previous simple one could be wrong in several low-order digits due to the inaccuracy in the representation of the exponent (1/3). This version effectively breaks the number up into the form a*8^b, computes the cube root of 8^b exactly (i.e. 2^b), and uses the slightly inaccurate exponentiation only for a. --- math2.asm | 99 ++++++++++++++++++++++++++++++++++++++-------------- math2.macros | 6 +++- 2 files changed, 78 insertions(+), 27 deletions(-) diff --git a/math2.asm b/math2.asm index 68fcf80..7f87d48 100644 --- a/math2.asm +++ b/math2.asm @@ -340,42 +340,89 @@ cbrt start cbrtf entry cbrtl entry using MathCommon2 - - phb place the number in a work area - plx (except exponent/sign word) - ply +scale equ 1 + + csubroutine (10:x),2 + + phb phk plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla get exponent/sign word - phy - phx - + + stz scale scale by 0 by default (for inf/nan) + + lda x+8 pha save original sign and #$7FFF - sta t1+8 force sign to + + sta x+8 force sign to + + cmp #32767 skip scaling for inf/nan + beq do_calc + + ldx x+6 if number is denormalized + bmi div_exp + bne normaliz + ldx x+4 + bne normaliz + ldx x+2 + bne normaliz + ldx x + beq div_exp - ph4 #onethird compute abs(x)^(1/3) +normaliz dec a normalize it and adjust exponent + asl x + rol x+2 + rol x+4 + rol x+6 + bpl normaliz + +div_exp pha calculate exponent/3 + pha + pha + pea 3 + _SDivide + pla a = quotient + plx x = remainder + cpx #2 adjust remainder of 2 to -1 + bne setscale + ldx #-1 + inc a + +setscale sec calculate amount to scale result by + sbc #16383/3 + sta scale + txa use remainder as exponent for calc. + clc + adc #16383 +do_calc sta t1+8 + + lda x place mantissa in work area + sta t1 + lda x+2 + sta t1+2 + lda x+4 + sta t1+4 + lda x+6 + sta t1+6 + + ph4 #onethird compute val^(1/3) ph4 #t1 FXPWRY - pla if sign of x was - - bpl ret + clc apply scaling lda t1+8 - ora #$8000 - sta t1+8 set sign of result to - + adc scale + sta t1+8 -ret plb - ldx #^t1 return a pointer to the result - lda #t1 - rtl + asl t1+8 set sign of result to orig. sign of x + pla + asl a + ror t1+8 + + plb + lda #t1 return t1 + sta x + lda #^t1 + sta x+2 + creturn 4:x onethird dc e'0.33333333333333333333' end diff --git a/math2.macros b/math2.macros index c666b3f..375f0de 100644 --- a/math2.macros +++ b/math2.macros @@ -616,4 +616,8 @@ LDX #$0B0A JSL $E10000 MEND - + MACRO +&lab _SDivide +&lab ldx #$0A0B + jsl $E10000 + MEND