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/float/cmplxf.c | 407 ---------------------------------------------------- 1 file changed, 407 deletions(-) delete mode 100644 libm/float/cmplxf.c (limited to 'libm/float/cmplxf.c') diff --git a/libm/float/cmplxf.c b/libm/float/cmplxf.c deleted file mode 100644 index 949b94e3d..000000000 --- a/libm/float/cmplxf.c +++ /dev/null @@ -1,407 +0,0 @@ -/* cmplxf.c - * - * Complex number arithmetic - * - * - * - * SYNOPSIS: - * - * typedef struct { - * float r; real part - * float i; imaginary part - * }cmplxf; - * - * cmplxf *a, *b, *c; - * - * caddf( a, b, c ); c = b + a - * csubf( a, b, c ); c = b - a - * cmulf( a, b, c ); c = b * a - * cdivf( a, b, c ); c = b / a - * cnegf( c ); c = -c - * cmovf( 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 - * IEEE cadd 30000 5.9e-8 2.6e-8 - * IEEE csub 30000 6.0e-8 2.6e-8 - * IEEE cmul 30000 1.1e-7 3.7e-8 - * IEEE cdiv 30000 2.1e-7 5.7e-8 - */ - /* cmplx.c - * complex number arithmetic - */ - - -/* -Cephes Math Library Release 2.1: December, 1988 -Copyright 1984, 1987, 1988 by Stephen L. Moshier -Direct inquiries to 30 Frost Street, Cambridge, MA 02140 -*/ - - -#include -extern float MAXNUMF, MACHEPF, PIF, PIO2F; -#define fabsf(x) ( (x) < 0 ? -(x) : (x) ) -#ifdef ANSIC -float sqrtf(float), frexpf(float, int *); -float ldexpf(float, int); -float cabsf(cmplxf *), atan2f(float, float), cosf(float), sinf(float); -#else -float sqrtf(), frexpf(), ldexpf(); -float cabsf(), atan2f(), cosf(), sinf(); -#endif -/* -typedef struct - { - float r; - float i; - }cmplxf; -*/ -cmplxf czerof = {0.0, 0.0}; -extern cmplxf czerof; -cmplxf conef = {1.0, 0.0}; -extern cmplxf conef; - -/* c = b + a */ - -void caddf( a, b, c ) -register cmplxf *a, *b; -cmplxf *c; -{ - -c->r = b->r + a->r; -c->i = b->i + a->i; -} - - -/* c = b - a */ - -void csubf( a, b, c ) -register cmplxf *a, *b; -cmplxf *c; -{ - -c->r = b->r - a->r; -c->i = b->i - a->i; -} - -/* c = b * a */ - -void cmulf( a, b, c ) -register cmplxf *a, *b; -cmplxf *c; -{ -register float 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 cdivf( a, b, c ) -register cmplxf *a, *b; -cmplxf *c; -{ -float 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.0f ) - { - w = MAXNUMF * y; - if( (fabsf(p) > w) || (fabsf(q) > w) || (y == 0.0f) ) - { - c->r = MAXNUMF; - c->i = MAXNUMF; - mtherr( "cdivf", OVERFLOW ); - return; - } - } -c->r = p/y; -c->i = q/y; -} - - -/* b = a */ - -void cmovf( a, b ) -register short *a, *b; -{ -int i; - - -i = 8; -do - *b++ = *a++; -while( --i ); -} - - -void cnegf( a ) -register cmplxf *a; -{ - -a->r = -a->r; -a->i = -a->i; -} - -/* cabsf() - * - * Complex absolute value - * - * - * - * SYNOPSIS: - * - * float cabsf(); - * cmplxf z; - * float a; - * - * a = cabsf( &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 - * IEEE -10,+10 30000 1.2e-7 3.4e-8 - */ - - -/* -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 - { - float r; - float i; - }cmplxf; -*/ -/* square root of max and min numbers */ -#define SMAX 1.3043817825332782216E+19 -#define SMIN 7.6664670834168704053E-20 -#define PREC 12 -#define MAXEXPF 128 - - -#define SMAXT (2.0f * SMAX) -#define SMINT (0.5f * SMIN) - -float cabsf( z ) -register cmplxf *z; -{ -float x, y, b, re, im; -int ex, ey, e; - -re = fabsf( z->r ); -im = fabsf( z->i ); - -if( re == 0.0f ) - { - return( im ); - } -if( im == 0.0f ) - { - return( re ); - } - -/* Get the exponents of the numbers */ -x = frexpf( re, &ex ); -y = frexpf( 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 = ldexpf( re, -e ); -y = ldexpf( im, -e ); - -/* Hypotenuse of the right triangle */ -b = sqrtf( x * x + y * y ); - -/* Compute the exponent of the answer. */ -y = frexpf( b, &ey ); -ey = e + ey; - -/* Check it for overflow and underflow. */ -if( ey > MAXEXPF ) - { - mtherr( "cabsf", OVERFLOW ); - return( MAXNUMF ); - } -if( ey < -MAXEXPF ) - return(0.0f); - -/* Undo the scaling */ -b = ldexpf( b, e ); -return( b ); -} - /* csqrtf() - * - * Complex square root - * - * - * - * SYNOPSIS: - * - * void csqrtf(); - * cmplxf z, w; - * - * csqrtf( &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 solution - * reported 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 - * IEEE -10,+10 100000 1.8e-7 4.2e-8 - * - */ - - -void csqrtf( z, w ) -cmplxf *z, *w; -{ -cmplxf q, s; -float x, y, r, t; - -x = z->r; -y = z->i; - -if( y == 0.0f ) - { - if( x < 0.0f ) - { - w->r = 0.0f; - w->i = sqrtf(-x); - return; - } - else - { - w->r = sqrtf(x); - w->i = 0.0f; - return; - } - } - -if( x == 0.0f ) - { - r = fabsf(y); - r = sqrtf(0.5f*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( (fabsf(y) < fabsf(0.015f*x)) - && (x > 0) ) - { - t = 0.25f*y*(y/x); - } -else - { - r = cabsf(z); - t = 0.5f*(r - x); - } - -r = sqrtf(t); -q.i = r; -q.r = 0.5f*y/r; - -/* Heron iteration in complex arithmetic: - * q = (q + z/q)/2 - */ -cdivf( &q, z, &s ); -caddf( &q, &s, w ); -w->r *= 0.5f; -w->i *= 0.5f; -} - -- cgit v1.2.3