summaryrefslogtreecommitdiff
path: root/libm/float/sqrtf.c
diff options
context:
space:
mode:
Diffstat (limited to 'libm/float/sqrtf.c')
-rw-r--r--libm/float/sqrtf.c140
1 files changed, 140 insertions, 0 deletions
diff --git a/libm/float/sqrtf.c b/libm/float/sqrtf.c
new file mode 100644
index 000000000..bc75a907b
--- /dev/null
+++ b/libm/float/sqrtf.c
@@ -0,0 +1,140 @@
+/* sqrtf.c
+ *
+ * Square root
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * float x, y, sqrtf();
+ *
+ * y = sqrtf( x );
+ *
+ *
+ *
+ * DESCRIPTION:
+ *
+ * Returns the square root of x.
+ *
+ * Range reduction involves isolating the power of two of the
+ * argument and using a polynomial approximation to obtain
+ * a rough value for the square root. Then Heron's iteration
+ * is used three times to converge to an accurate value.
+ *
+ *
+ *
+ * ACCURACY:
+ *
+ *
+ * Relative error:
+ * arithmetic domain # trials peak rms
+ * IEEE 0,1.e38 100000 8.7e-8 2.9e-8
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ * message condition value returned
+ * sqrtf domain x < 0 0.0
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2: June, 1992
+Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+/* Single precision square root
+ * test interval: [sqrt(2)/2, sqrt(2)]
+ * trials: 30000
+ * peak relative error: 8.8e-8
+ * rms relative error: 3.3e-8
+ *
+ * test interval: [0.01, 100.0]
+ * trials: 50000
+ * peak relative error: 8.7e-8
+ * rms relative error: 3.3e-8
+ *
+ * Copyright (C) 1989 by Stephen L. Moshier. All rights reserved.
+ */
+#include <math.h>
+
+#ifdef ANSIC
+float frexpf( float, int * );
+float ldexpf( float, int );
+
+float sqrtf( float xx )
+#else
+float frexpf(), ldexpf();
+
+float sqrtf(xx)
+float xx;
+#endif
+{
+float f, x, y;
+int e;
+
+f = xx;
+if( f <= 0.0 )
+ {
+ if( f < 0.0 )
+ mtherr( "sqrtf", DOMAIN );
+ return( 0.0 );
+ }
+
+x = frexpf( f, &e ); /* f = x * 2**e, 0.5 <= x < 1.0 */
+/* If power of 2 is odd, double x and decrement the power of 2. */
+if( e & 1 )
+ {
+ x = x + x;
+ e -= 1;
+ }
+
+e >>= 1; /* The power of 2 of the square root. */
+
+if( x > 1.41421356237 )
+ {
+/* x is between sqrt(2) and 2. */
+ x = x - 2.0;
+ y =
+ ((((( -9.8843065718E-4 * x
+ + 7.9479950957E-4) * x
+ - 3.5890535377E-3) * x
+ + 1.1028809744E-2) * x
+ - 4.4195203560E-2) * x
+ + 3.5355338194E-1) * x
+ + 1.41421356237E0;
+ goto sqdon;
+ }
+
+if( x > 0.707106781187 )
+ {
+/* x is between sqrt(2)/2 and sqrt(2). */
+ x = x - 1.0;
+ y =
+ ((((( 1.35199291026E-2 * x
+ - 2.26657767832E-2) * x
+ + 2.78720776889E-2) * x
+ - 3.89582788321E-2) * x
+ + 6.24811144548E-2) * x
+ - 1.25001503933E-1) * x * x
+ + 0.5 * x
+ + 1.0;
+ goto sqdon;
+ }
+
+/* x is between 0.5 and sqrt(2)/2. */
+x = x - 0.5;
+y =
+((((( -3.9495006054E-1 * x
+ + 5.1743034569E-1) * x
+ - 4.3214437330E-1) * x
+ + 3.5310730460E-1) * x
+ - 3.5354581892E-1) * x
+ + 7.0710676017E-1) * x
+ + 7.07106781187E-1;
+
+sqdon:
+y = ldexpf( y, e ); /* y = y * 2**e */
+return( y);
+}