summaryrefslogtreecommitdiff
path: root/src/math/ef_sqrt.c
diff options
context:
space:
mode:
authorBdale Garbee <bdale@gag.com>2013-12-19 01:38:40 -0700
committerBdale Garbee <bdale@gag.com>2013-12-19 01:38:40 -0700
commit575bbaf976c5840fd0e308549c45a466fdec1352 (patch)
tree11bfb498348bf7687bffc24699c4b1a998988ee4 /src/math/ef_sqrt.c
parentb825116df173b77e2cab217a7b76112c742f9279 (diff)
parentbc3610d8cecbfed40c62d4dcb93fc9a4d2a7c9e3 (diff)
Merge branch 'branch-1.3' into debian
Conflicts: ChangeLog altoslib/AltosRecordMM.java altosui/Makefile.am altosui/altos-windows.nsi.in configure.ac debian/changelog debian/control doc/Makefile doc/altusmetrum.xsl doc/release-notes-1.2.1.xsl doc/release-notes-1.2.xsl
Diffstat (limited to 'src/math/ef_sqrt.c')
-rw-r--r--src/math/ef_sqrt.c89
1 files changed, 89 insertions, 0 deletions
diff --git a/src/math/ef_sqrt.c b/src/math/ef_sqrt.c
new file mode 100644
index 00000000..80e7f360
--- /dev/null
+++ b/src/math/ef_sqrt.c
@@ -0,0 +1,89 @@
+/* ef_sqrtf.c -- float version of e_sqrt.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 "fdlibm.h"
+
+#ifdef __STDC__
+static const float one = 1.0, tiny=1.0e-30;
+#else
+static float one = 1.0, tiny=1.0e-30;
+#endif
+
+#ifdef __STDC__
+ float __ieee754_sqrtf(float x)
+#else
+ float __ieee754_sqrtf(x)
+ float x;
+#endif
+{
+ float z;
+ __uint32_t r,hx;
+ __int32_t ix,s,q,m,t,i;
+
+ GET_FLOAT_WORD(ix,x);
+ hx = ix&0x7fffffff;
+
+ /* take care of Inf and NaN */
+ if(!FLT_UWORD_IS_FINITE(hx))
+ return x*x+x; /* sqrt(NaN)=NaN, sqrt(+inf)=+inf
+ sqrt(-inf)=sNaN */
+ /* take care of zero and -ves */
+ if(FLT_UWORD_IS_ZERO(hx)) return x;/* sqrt(+-0) = +-0 */
+ if(ix<0) return (x-x)/(x-x); /* sqrt(-ve) = sNaN */
+
+ /* normalize x */
+ m = (ix>>23);
+ if(FLT_UWORD_IS_SUBNORMAL(hx)) { /* subnormal x */
+ for(i=0;(ix&0x00800000L)==0;i++) ix<<=1;
+ m -= i-1;
+ }
+ m -= 127; /* unbias exponent */
+ ix = (ix&0x007fffffL)|0x00800000L;
+ 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 = 0x01000000L; /* 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 = one-tiny; /* trigger inexact flag */
+ if (z>=one) {
+ z = one+tiny;
+ if (z>one)
+ q += 2;
+ else
+ q += (q&1);
+ }
+ }
+ ix = (q>>1)+0x3f000000L;
+ ix += (m <<23);
+ SET_FLOAT_WORD(z,ix);
+ return z;
+}