diff options
Diffstat (limited to 'libm/ldouble/sqrtl.c')
-rw-r--r-- | libm/ldouble/sqrtl.c | 172 |
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); -} |