summaryrefslogtreecommitdiff
path: root/libm/double/sqrt.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/double/sqrt.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/double/sqrt.c')
-rw-r--r--libm/double/sqrt.c178
1 files changed, 0 insertions, 178 deletions
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 <math.h>
-#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);
-}