jyapayne/musl
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

math: change the formula used for acos.s

Browse source 
old: 2*atan2(sqrt(1-x),sqrt(1+x))
new: atan2(fabs(sqrt((1-x)*(1+x))),x)
improvements:
* all edge cases are fixed (sign of zero in downward rounding)
* a bit faster (here a single call is about 131ns vs 162ns)
* a bit more precise (at most 1ulp error on 1M uniform random
samples in [0,1), the old formula gave some 2ulp errors as well)
This commit is contained in:
nsz 2012-05-05 01:11:56 +02:00
commit f697d66b81
2 changed files with 18 additions and 22 deletions
Download patch file Download diff file Expand all files Collapse all files

20
src/math/i386/acos.s
View file

@ -1,3 +1,5 @@
# use acos(x) = atan2(fabs(sqrt((1-x)*(1+x))), x)
.global acosf .global acosf
.type acosf,@function .type acosf,@function
acosf: acosf:
@ -14,17 +16,13 @@ acosl:
.type acos,@function .type acos,@function
acos: acos:
fldl 4(%esp) fldl 4(%esp)
1: fld1 1: fld %st(0)
fld %st(1)
fld1 fld1
fsubp fsub %st(0),%st(1)
fsqrt fadd %st(2)
fxch %st(2)
faddp
fsqrt
fpatan
fld1
fld1
faddp
fmulp fmulp
fsqrt
fabs # fix sign of zero (matters in downward rounding mode)
fxch %st(1)
fpatan
ret ret

20
src/math/x86_64/acosl.s
View file

@ -1,18 +1,16 @@
# see ../i386/acos.s
.global acosl .global acosl
.type acosl,@function .type acosl,@function
acosl: acosl:
fldt 8(%rsp) fldt 8(%rsp)
1: fld %st(0)
fld1 fld1
fld %st(1) fsub %st(0),%st(1)
fld1 fadd %st(2)
fsubp
fsqrt
fxch %st(2)
faddp
fsqrt
fpatan
fld1
fld1
faddp
fmulp fmulp
fsqrt
fabs
fxch %st(1)
fpatan
ret ret