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, 461 deletions

diff --git a/libm/double/cmplx.c b/libm/double/cmplx.c
deleted file mode 100644
index dcd972bea..000000000
--- a/libm/double/cmplx.c
+++ /dev/null
@@ -1,461 +0,0 @@
-/*							cmplx.c
- *
- *	Complex number arithmetic
- *
- *
- *
- * SYNOPSIS:
- *
- * typedef struct {
- *      double r;     real part
- *      double i;     imaginary part
- *     }cmplx;
- *
- * cmplx *a, *b, *c;
- *
- * cadd( a, b, c );     c = b + a
- * csub( a, b, c );     c = b - a
- * cmul( a, b, c );     c = b * a
- * cdiv( a, b, c );     c = b / a
- * cneg( c );           c = -c
- * cmov( b, c );        c = b
- *
- *
- *
- * DESCRIPTION:
- *
- * Addition:
- *    c.r  =  b.r + a.r
- *    c.i  =  b.i + a.i
- *
- * Subtraction:
- *    c.r  =  b.r - a.r
- *    c.i  =  b.i - a.i
- *
- * Multiplication:
- *    c.r  =  b.r * a.r  -  b.i * a.i
- *    c.i  =  b.r * a.i  +  b.i * a.r
- *
- * Division:
- *    d    =  a.r * a.r  +  a.i * a.i
- *    c.r  = (b.r * a.r  + b.i * a.i)/d
- *    c.i  = (b.i * a.r  -  b.r * a.i)/d
- * ACCURACY:
- *
- * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
- * error 3.1e-17, rms 1.2e-17.  The test (y/z) * (z/y) = 1 had
- * peak relative error 8.3e-17, rms 2.1e-17.
- *
- * Tests in the rectangle {-10,+10}:
- *                      Relative error:
- * arithmetic   function  # trials      peak         rms
- *    DEC        cadd       10000       1.4e-17     3.4e-18
- *    IEEE       cadd      100000       1.1e-16     2.7e-17
- *    DEC        csub       10000       1.4e-17     4.5e-18
- *    IEEE       csub      100000       1.1e-16     3.4e-17
- *    DEC        cmul        3000       2.3e-17     8.7e-18
- *    IEEE       cmul      100000       2.1e-16     6.9e-17
- *    DEC        cdiv       18000       4.9e-17     1.3e-17
- *    IEEE       cdiv      100000       3.7e-16     1.1e-16
- */
-/*				cmplx.c
- * complex number arithmetic
- */
-
-
-/*
-Cephes Math Library Release 2.8:  June, 2000
-Copyright 1984, 1995, 2000 by Stephen L. Moshier
-*/
-
-
-#include <math.h>
-
-#ifdef ANSIPROT
-extern double fabs ( double );
-extern double cabs ( cmplx * );
-extern double sqrt ( double );
-extern double atan2 ( double, double );
-extern double cos ( double );
-extern double sin ( double );
-extern double sqrt ( double );
-extern double frexp ( double, int * );
-extern double ldexp ( double, int );
-int isnan ( double );
-void cdiv ( cmplx *, cmplx *, cmplx * );
-void cadd ( cmplx *, cmplx *, cmplx * );
-#else
-double fabs(), cabs(), sqrt(), atan2(), cos(), sin();
-double sqrt(), frexp(), ldexp();
-int isnan();
-void cdiv(), cadd();
-#endif
-
-extern double MAXNUM, MACHEP, PI, PIO2, INFINITY, NAN;
-/*
-typedef struct
-	{
-	double r;
-	double i;
-	}cmplx;
-*/
-cmplx czero = {0.0, 0.0};
-extern cmplx czero;
-cmplx cone = {1.0, 0.0};
-extern cmplx cone;
-
-/*	c = b + a	*/
-
-void cadd( a, b, c )
-register cmplx *a, *b;
-cmplx *c;
-{
-
-c->r = b->r + a->r;
-c->i = b->i + a->i;
-}
-
-
-/*	c = b - a	*/
-
-void csub( a, b, c )
-register cmplx *a, *b;
-cmplx *c;
-{
-
-c->r = b->r - a->r;
-c->i = b->i - a->i;
-}
-
-/*	c = b * a */
-
-void cmul( a, b, c )
-register cmplx *a, *b;
-cmplx *c;
-{
-double y;
-
-y    = b->r * a->r  -  b->i * a->i;
-c->i = b->r * a->i  +  b->i * a->r;
-c->r = y;
-}
-
-
-
-/*	c = b / a */
-
-void cdiv( a, b, c )
-register cmplx *a, *b;
-cmplx *c;
-{
-double y, p, q, w;
-
-
-y = a->r * a->r  +  a->i * a->i;
-p = b->r * a->r  +  b->i * a->i;
-q = b->i * a->r  -  b->r * a->i;
-
-if( y < 1.0 )
-	{
-	w = MAXNUM * y;
-	if( (fabs(p) > w) || (fabs(q) > w) || (y == 0.0) )
-		{
-		c->r = MAXNUM;
-		c->i = MAXNUM;
-		mtherr( "cdiv", OVERFLOW );
-		return;
-		}
-	}
-c->r = p/y;
-c->i = q/y;
-}
-
-
-/*	b = a
-   Caution, a `short' is assumed to be 16 bits wide.  */
-
-void cmov( a, b )
-void *a, *b;
-{
-register short *pa, *pb;
-int i;
-
-pa = (short *) a;
-pb = (short *) b;
-i = 8;
-do
-	*pb++ = *pa++;
-while( --i );
-}
-
-
-void cneg( a )
-register cmplx *a;
-{
-
-a->r = -a->r;
-a->i = -a->i;
-}
-
-/*							cabs()
- *
- *	Complex absolute value
- *
- *
- *
- * SYNOPSIS:
- *
- * double cabs();
- * cmplx z;
- * double a;
- *
- * a = cabs( &z );
- *
- *
- *
- * DESCRIPTION:
- *
- *
- * If z = x + iy
- *
- * then
- *
- *       a = sqrt( x**2 + y**2 ).
- * 
- * Overflow and underflow are avoided by testing the magnitudes
- * of x and y before squaring.  If either is outside half of
- * the floating point full scale range, both are rescaled.
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    DEC       -30,+30     30000       3.2e-17     9.2e-18
- *    IEEE      -10,+10    100000       2.7e-16     6.9e-17
- */
-
-
-/*
-Cephes Math Library Release 2.1:  January, 1989
-Copyright 1984, 1987, 1989 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-
-/*
-typedef struct
-	{
-	double r;
-	double i;
-	}cmplx;
-*/
-
-#ifdef UNK
-#define PREC 27
-#define MAXEXP 1024
-#define MINEXP -1077
-#endif
-#ifdef DEC
-#define PREC 29
-#define MAXEXP 128
-#define MINEXP -128
-#endif
-#ifdef IBMPC
-#define PREC 27
-#define MAXEXP 1024
-#define MINEXP -1077
-#endif
-#ifdef MIEEE
-#define PREC 27
-#define MAXEXP 1024
-#define MINEXP -1077
-#endif
-
-
-double cabs( z )
-register cmplx *z;
-{
-double x, y, b, re, im;
-int ex, ey, e;
-
-#ifdef INFINITIES
-/* Note, cabs(INFINITY,NAN) = INFINITY. */
-if( z->r == INFINITY || z->i == INFINITY
-   || z->r == -INFINITY || z->i == -INFINITY )
-  return( INFINITY );
-#endif
-
-#ifdef NANS
-if( isnan(z->r) )
-  return(z->r);
-if( isnan(z->i) )
-  return(z->i);
-#endif
-
-re = fabs( z->r );
-im = fabs( z->i );
-
-if( re == 0.0 )
-	return( im );
-if( im == 0.0 )
-	return( re );
-
-/* Get the exponents of the numbers */
-x = frexp( re, &ex );
-y = frexp( im, &ey );
-
-/* Check if one number is tiny compared to the other */
-e = ex - ey;
-if( e > PREC )
-	return( re );
-if( e < -PREC )
-	return( im );
-
-/* Find approximate exponent e of the geometric mean. */
-e = (ex + ey) >> 1;
-
-/* Rescale so mean is about 1 */
-x = ldexp( re, -e );
-y = ldexp( im, -e );
-		
-/* Hypotenuse of the right triangle */
-b = sqrt( x * x  +  y * y );
-
-/* Compute the exponent of the answer. */
-y = frexp( b, &ey );
-ey = e + ey;
-
-/* Check it for overflow and underflow. */
-if( ey > MAXEXP )
-	{
-	mtherr( "cabs", OVERFLOW );
-	return( INFINITY );
-	}
-if( ey < MINEXP )
-	return(0.0);
-
-/* Undo the scaling */
-b = ldexp( b, e );
-return( b );
-}
-/*							csqrt()
- *
- *	Complex square root
- *
- *
- *
- * SYNOPSIS:
- *
- * void csqrt();
- * cmplx z, w;
- *
- * csqrt( &z, &w );
- *
- *
- *
- * DESCRIPTION:
- *
- *
- * If z = x + iy,  r = |z|, then
- *
- *                       1/2
- * Im w  =  [ (r - x)/2 ]   ,
- *
- * Re w  =  y / 2 Im w.
- *
- *
- * Note that -w is also a square root of z.  The root chosen
- * is always in the upper half plane.
- *
- * Because of the potential for cancellation error in r - x,
- * the result is sharpened by doing a Heron iteration
- * (see sqrt.c) in complex arithmetic.
- *
- *
- *
- * ACCURACY:
- *
- *                      Relative error:
- * arithmetic   domain     # trials      peak         rms
- *    DEC       -10,+10     25000       3.2e-17     9.6e-18
- *    IEEE      -10,+10    100000       3.2e-16     7.7e-17
- *
- *                        2
- * Also tested by csqrt( z ) = z, and tested by arguments
- * close to the real axis.
- */
-
-
-void csqrt( z, w )
-cmplx *z, *w;
-{
-cmplx q, s;
-double x, y, r, t;
-
-x = z->r;
-y = z->i;
-
-if( y == 0.0 )
-	{
-	if( x < 0.0 )
-		{
-		w->r = 0.0;
-		w->i = sqrt(-x);
-		return;
-		}
-	else
-		{
-		w->r = sqrt(x);
-		w->i = 0.0;
-		return;
-		}
-	}
-
-
-if( x == 0.0 )
-	{
-	r = fabs(y);
-	r = sqrt(0.5*r);
-	if( y > 0 )
-		w->r = r;
-	else
-		w->r = -r;
-	w->i = r;
-	return;
-	}
-
-/* Approximate  sqrt(x^2+y^2) - x  =  y^2/2x - y^4/24x^3 + ... .
- * The relative error in the first term is approximately y^2/12x^2 .
- */
-if( (fabs(y) < 2.e-4 * fabs(x))
-   && (x > 0) )
-	{
-	t = 0.25*y*(y/x);
-	}
-else
-	{
-	r = cabs(z);
-	t = 0.5*(r - x);
-	}
-
-r = sqrt(t);
-q.i = r;
-q.r = y/(2.0*r);
-/* Heron iteration in complex arithmetic */
-cdiv( &q, z, &s );
-cadd( &q, &s, w );
-w->r *= 0.5;
-w->i *= 0.5;
-}
-
-
-double hypot( x, y )
-double x, y;
-{
-cmplx z;
-
-z.r = x;
-z.i = y;
-return( cabs(&z) );
-}