math: move i386 sqrt to C with inline asm

diff --git a/src/math/i386/sqrt.c b/src/math/i386/sqrt.c
new file mode 100644 (file)
index 0000000..934fbcc
--- /dev/null
+++ b/src/math/i386/sqrt.c
@@ -0,0 +1,15 @@
+#include "libm.h"
+
+double sqrt(double x)
+{
+       union ldshape ux;
+       unsigned fpsr;
+       __asm__ ("fsqrt; fnstsw %%ax": "=t"(ux.f), "=a"(fpsr) : "0"(x));
+       if ((ux.i.m & 0x7ff) != 0x400)
+               return (double)ux.f;
+       /* Rounding to double would have encountered an exact halfway case.
+          Adjust mantissa downwards if fsqrt rounded up, else upwards.
+          (result of fsqrt could not have been exact) */
+       ux.i.m ^= (fpsr & 0x200) + 0x300;
+       return (double)ux.f;
+}

diff --git a/src/math/i386/sqrt.s b/src/math/i386/sqrt.s
deleted file mode 100644 (file)
index 57837e2..0000000
--- a/src/math/i386/sqrt.s
+++ /dev/null
@@ -1,21 +0,0 @@
-.global sqrt
-.type sqrt,@function
-sqrt:  fldl 4(%esp)
-       fsqrt
-       fnstsw %ax
-       sub $12,%esp
-       fld %st(0)
-       fstpt (%esp)
-       mov (%esp),%ecx
-       and $0x7ff,%ecx
-       cmp $0x400,%ecx
-       jnz 1f
-       and $0x200,%eax
-       sub $0x100,%eax
-       sub %eax,(%esp)
-       fstp %st(0)
-       fldt (%esp)
-1:     add $12,%esp
-       fstpl 4(%esp)
-       fldl 4(%esp)
-       ret