acos.s fix: use the formula acos(x) = atan2(sqrt(1-x),sqrt(1+x))

the old formula atan2(1,sqrt((1+x)/(1-x))) was faster but
could give nan result at x=1 when the rounding mode is
FE_DOWNWARD (so 1-1 == -0 and 2/-0 == -inf), the new formula
gives -0 at x=+-1 with downward rounding.

diff --git a/src/math/i386/acos.s b/src/math/i386/acos.s
index 4f0168f..bfff0c5 100644 (file)
--- a/src/math/i386/acos.s
+++ b/src/math/i386/acos.s
@@ -18,12 +18,10 @@ acos:
        fld %st(1)
        fld1
        fsubp
+       fsqrt
        fxch %st(2)
        faddp
-       fdivp
        fsqrt
-       fld1
-       fxch %st(1)
        fpatan
        fld1
        fld1

diff --git a/src/math/x86_64/acosl.s b/src/math/x86_64/acosl.s
index 92c2987..db68d2d 100644 (file)
--- a/src/math/x86_64/acosl.s
+++ b/src/math/x86_64/acosl.s
@@ -6,12 +6,10 @@ acosl:
        fld %st(1)
        fld1
        fsubp
+       fsqrt
        fxch %st(2)
        faddp
-       fdivp
        fsqrt
-       fld1
-       fxch %st(1)
        fpatan
        fld1
        fld1