1 files changed, 172 insertions, 0 deletions
diff --git a/libm/ldouble/sqrtl.c b/libm/ldouble/sqrtl.c
new file mode 100644
index 000000000..a3b17175f
--- /dev/null
+++ b/libm/ldouble/sqrtl.c
@@ -0,0 +1,172 @@
+/*							sqrtl.c
+ *
+ *	Square root, long double precision
+ *
+ *
+ *
+ * SYNOPSIS:
+ *
+ * long double x, y, sqrtl();
+ *
+ * y = sqrtl( 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.
+ *
+ * Note, some arithmetic coprocessors such as the 8087 and
+ * 68881 produce correctly rounded square roots, which this
+ * routine will not.
+ *
+ * ACCURACY:
+ *
+ *
+ *                      Relative error:
+ * arithmetic   domain     # trials      peak         rms
+ *    IEEE      0,10        30000       8.1e-20     3.1e-20
+ *
+ *
+ * ERROR MESSAGES:
+ *
+ *   message         condition      value returned
+ * sqrt domain        x < 0            0.0
+ *
+ */
+
+/*
+Cephes Math Library Release 2.2:  December, 1990
+Copyright 1984, 1990 by Stephen L. Moshier
+Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+*/
+
+
+#include <math.h>
+
+#define SQRT2 1.4142135623730950488017E0L
+#ifdef ANSIPROT
+extern long double frexpl ( long double, int * );
+extern long double ldexpl ( long double, int );
+#else
+long double frexpl(), ldexpl();
+#endif
+
+long double sqrtl(x)
+long double x;
+{
+int e;
+long double z, w;
+#ifndef UNK
+short *q;
+#endif
+
+if( x <= 0.0 )
+	{
+	if( x < 0.0 )
+		mtherr( "sqrtl", DOMAIN );
+	return( 0.0 );
+	}
+w = x;
+/* separate exponent and significand */
+#ifdef UNK
+z = frexpl( x, &e );
+#endif
+
+/* Note, frexp and ldexp are used in order to
+ * handle denormal numbers properly.
+ */
+#ifdef IBMPC
+z = frexpl( x, &e );
+q = (short *)&x; /* point to the exponent word */
+q += 4;
+/*
+e = ((*q >> 4) & 0x0fff) - 0x3fe;
+*q &= 0x000f;
+*q |= 0x3fe0;
+z = x;
+*/
+#endif
+#ifdef MIEEE
+z = frexpl( x, &e );
+q = (short *)&x;
+/*
+e = ((*q >> 4) & 0x0fff) - 0x3fe;
+*q &= 0x000f;
+*q |= 0x3fe0;
+z = x;
+*/
+#endif
+
+/* approximate square root of number between 0.5 and 1
+ * relative error of linear approximation = 7.47e-3
+ */
+/*
+x = 0.4173075996388649989089L + 0.59016206709064458299663L * z;
+*/
+
+/* quadratic approximation, relative error 6.45e-4 */
+x = ( -0.20440583154734771959904L  * z
+     + 0.89019407351052789754347L) * z
+     + 0.31356706742295303132394L;
+
+/* adjust for odd powers of 2 */
+if( (e & 1) != 0 )
+	x *= SQRT2;
+
+/* re-insert exponent */
+#ifdef UNK
+x = ldexpl( x, (e >> 1) );
+#endif
+#ifdef IBMPC
+x = ldexpl( x, (e >> 1) );
+/*
+*q += ((e >>1) & 0x7ff) << 4;
+*q &= 077777;
+*/
+#endif
+#ifdef MIEEE
+x = ldexpl( x, (e >> 1) );
+/*
+*q += ((e >>1) & 0x7ff) << 4;
+*q &= 077777;
+*/
+#endif
+
+/* Newton iterations: */
+#ifdef UNK
+x += w/x;
+x = ldexpl( x, -1 );	/* divide by 2 */
+x += w/x;
+x = ldexpl( x, -1 );
+x += w/x;
+x = ldexpl( x, -1 );
+#endif
+
+/* Note, assume the square root cannot be denormal,
+ * so it is safe to use integer exponent operations here.
+ */
+#ifdef IBMPC
+x += w/x;
+*q -= 1;
+x += w/x;
+*q -= 1;
+x += w/x;
+*q -= 1;
+#endif
+#ifdef MIEEE
+x += w/x;
+*q -= 1;
+x += w/x;
+*q -= 1;
+x += w/x;
+*q -= 1;
+#endif
+
+return(x);
+}