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

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);
-}