Index: fusl/src/math/sqrtf.c |
diff --git a/fusl/src/math/sqrtf.c b/fusl/src/math/sqrtf.c |
new file mode 100644 |
index 0000000000000000000000000000000000000000..28cb4ad3713e6447843631a8d4e686f0f1dfdbb9 |
--- /dev/null |
+++ b/fusl/src/math/sqrtf.c |
@@ -0,0 +1,84 @@ |
+/* origin: FreeBSD /usr/src/lib/msun/src/e_sqrtf.c */ |
+/* |
+ * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. |
+ */ |
+/* |
+ * ==================================================== |
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
+ * |
+ * Developed at SunPro, a Sun Microsystems, Inc. business. |
+ * Permission to use, copy, modify, and distribute this |
+ * software is freely granted, provided that this notice |
+ * is preserved. |
+ * ==================================================== |
+ */ |
+ |
+#include "libm.h" |
+ |
+static const float tiny = 1.0e-30; |
+ |
+float sqrtf(float x) |
+{ |
+ float z; |
+ int32_t sign = (int)0x80000000; |
+ int32_t ix,s,q,m,t,i; |
+ uint32_t r; |
+ |
+ GET_FLOAT_WORD(ix, x); |
+ |
+ /* take care of Inf and NaN */ |
+ if ((ix&0x7f800000) == 0x7f800000) |
+ return x*x + x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf, sqrt(-inf)=sNaN */ |
+ |
+ /* take care of zero */ |
+ if (ix <= 0) { |
+ if ((ix&~sign) == 0) |
+ return x; /* sqrt(+-0) = +-0 */ |
+ if (ix < 0) |
+ return (x-x)/(x-x); /* sqrt(-ve) = sNaN */ |
+ } |
+ /* normalize x */ |
+ m = ix>>23; |
+ if (m == 0) { /* subnormal x */ |
+ for (i = 0; (ix&0x00800000) == 0; i++) |
+ ix<<=1; |
+ m -= i - 1; |
+ } |
+ m -= 127; /* unbias exponent */ |
+ ix = (ix&0x007fffff)|0x00800000; |
+ if (m&1) /* odd m, double x to make it even */ |
+ ix += ix; |
+ m >>= 1; /* m = [m/2] */ |
+ |
+ /* generate sqrt(x) bit by bit */ |
+ ix += ix; |
+ q = s = 0; /* q = sqrt(x) */ |
+ r = 0x01000000; /* r = moving bit from right to left */ |
+ |
+ while (r != 0) { |
+ t = s + r; |
+ if (t <= ix) { |
+ s = t+r; |
+ ix -= t; |
+ q += r; |
+ } |
+ ix += ix; |
+ r >>= 1; |
+ } |
+ |
+ /* use floating add to find out rounding direction */ |
+ if (ix != 0) { |
+ z = 1.0f - tiny; /* raise inexact flag */ |
+ if (z >= 1.0f) { |
+ z = 1.0f + tiny; |
+ if (z > 1.0f) |
+ q += 2; |
+ else |
+ q += q & 1; |
+ } |
+ } |
+ ix = (q>>1) + 0x3f000000; |
+ ix += m << 23; |
+ SET_FLOAT_WORD(z, ix); |
+ return z; |
+} |