summaryrefslogtreecommitdiff
path: root/libm/ldouble/sqrtl.c
diff options
context:
space:
mode:
authorEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
committerEric Andersen <andersen@codepoet.org>2001-11-22 14:04:29 +0000
commit7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 (patch)
tree3a7e8476e868ae15f4da1b7ce26b2db6f434468c /libm/ldouble/sqrtl.c
parentc117dd5fb183afb1a4790a6f6110d88704be6bf8 (diff)
Totally rework the math library, this time based on the MacOs X
math library (which is itself based on the math lib from FreeBSD). -Erik
Diffstat (limited to 'libm/ldouble/sqrtl.c')
-rw-r--r--libm/ldouble/sqrtl.c172
1 files changed, 0 insertions, 172 deletions
diff --git a/libm/ldouble/sqrtl.c b/libm/ldouble/sqrtl.c
deleted file mode 100644
index a3b17175f..000000000
--- a/libm/ldouble/sqrtl.c
+++ /dev/null
@@ -1,172 +0,0 @@
-/* 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);
-}