From 7ce331c01ce6eb7b3f5c715a38a24359da9c6ee2 Mon Sep 17 00:00:00 2001 From: Eric Andersen Date: Thu, 22 Nov 2001 14:04:29 +0000 Subject: 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 --- libm/double/sqrt.c | 178 ----------------------------------------------------- 1 file changed, 178 deletions(-) delete mode 100644 libm/double/sqrt.c (limited to 'libm/double/sqrt.c') diff --git a/libm/double/sqrt.c b/libm/double/sqrt.c deleted file mode 100644 index 92bbce53b..000000000 --- a/libm/double/sqrt.c +++ /dev/null @@ -1,178 +0,0 @@ -/* sqrt.c - * - * Square root - * - * - * - * SYNOPSIS: - * - * double x, y, sqrt(); - * - * y = sqrt( 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 - * DEC 0, 10 60000 2.1e-17 7.9e-18 - * IEEE 0,1.7e308 30000 1.7e-16 6.3e-17 - * - * - * ERROR MESSAGES: - * - * message condition value returned - * sqrt domain x < 0 0.0 - * - */ - -/* -Cephes Math Library Release 2.8: June, 2000 -Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier -*/ - - -#include -#ifdef ANSIPROT -extern double frexp ( double, int * ); -extern double ldexp ( double, int ); -#else -double frexp(), ldexp(); -#endif -extern double SQRT2; /* SQRT2 = 1.41421356237309504880 */ - -double sqrt(x) -double x; -{ -int e; -#ifndef UNK -short *q; -#endif -double z, w; - -if( x <= 0.0 ) - { - if( x < 0.0 ) - mtherr( "sqrt", DOMAIN ); - return( 0.0 ); - } -w = x; -/* separate exponent and significand */ -#ifdef UNK -z = frexp( x, &e ); -#endif -#ifdef DEC -q = (short *)&x; -e = ((*q >> 7) & 0377) - 0200; -*q &= 0177; -*q |= 040000; -z = x; -#endif - -/* Note, frexp and ldexp are used in order to - * handle denormal numbers properly. - */ -#ifdef IBMPC -z = frexp( x, &e ); -q = (short *)&x; -q += 3; -/* -e = ((*q >> 4) & 0x0fff) - 0x3fe; -*q &= 0x000f; -*q |= 0x3fe0; -z = x; -*/ -#endif -#ifdef MIEEE -z = frexp( 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 approximation = 7.47e-3 - */ -x = 4.173075996388649989089E-1 + 5.9016206709064458299663E-1 * z; - -/* adjust for odd powers of 2 */ -if( (e & 1) != 0 ) - x *= SQRT2; - -/* re-insert exponent */ -#ifdef UNK -x = ldexp( x, (e >> 1) ); -#endif -#ifdef DEC -*q += ((e >> 1) & 0377) << 7; -*q &= 077777; -#endif -#ifdef IBMPC -x = ldexp( x, (e >> 1) ); -/* -*q += ((e >>1) & 0x7ff) << 4; -*q &= 077777; -*/ -#endif -#ifdef MIEEE -x = ldexp( x, (e >> 1) ); -/* -*q += ((e >>1) & 0x7ff) << 4; -*q &= 077777; -*/ -#endif - -/* Newton iterations: */ -#ifdef UNK -x = 0.5*(x + w/x); -x = 0.5*(x + w/x); -x = 0.5*(x + w/x); -#endif - -/* Note, assume the square root cannot be denormal, - * so it is safe to use integer exponent operations here. - */ -#ifdef DEC -x += w/x; -*q -= 0200; -x += w/x; -*q -= 0200; -x += w/x; -*q -= 0200; -#endif -#ifdef IBMPC -x += w/x; -*q -= 0x10; -x += w/x; -*q -= 0x10; -x += w/x; -*q -= 0x10; -#endif -#ifdef MIEEE -x += w/x; -*q -= 0x10; -x += w/x; -*q -= 0x10; -x += w/x; -*q -= 0x10; -#endif - -return(x); -} -- cgit v1.2.3