/* 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 "complex.h" #include "mconf.h" 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(float complex); float atan2f(float, float); float cosf(float); float sinf(float); #else float sqrtf(), frexpf(), ldexpf(); float cabsf(), atan2f(), cosf(), sinf(); #endif /* typedef struct { float r; float i; }cmplxf; */ float complex czerof = 0.0f + 0.0f * I; float complex conef = 1.0f + 0.0f * I; /* c = b + a */ float complex caddf( float complex a, float complex b ) { return (creal (b) + creal (a) + (cimag (b) + cimag (a)) * I); } /* c = b - a */ float complex csubf( float complex a, float complex b ) { return (creal (b) - creal (a) + (cimag (b) - cimag (a)) * I); } /* c = b * a */ float complex cmulf( float complex a, float complex b ) { return ((creal (b) * creal (a) - cimag (b) * cimag (a)) + (creal (b) * cimag (a) + cimag (b) * creal (a)) * I); } /* c = b / a */ float complex cdivf( float complex a, float complex b ) { float y, p, q, w; y = creal (a) * creal (a) + cimag (a) * cimag (a); p = creal (b) * creal (a) + cimag (b) * cimag (a); q = cimag (b) * creal (a) - creal (b) * cimag (a); if( y < 1.0f ) { w = MAXNUMF * y; if( (fabsf(p) > w) || (fabsf(q) > w) || (y == 0.0f) ) { mtherr( "cdivf", OVERFLOW ); return (MAXNUMF + MAXNUMF * I); } } return (p/y + (q/y) * 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.0 * SMAX) #define SMINT (0.5 * SMIN) float cabsf( float complex z ) { float x, y, b, re, im; int ex, ey, e; re = fabsf( (float) creal(z) ); im = fabsf( (float) cimag(z) ); if( re == 0.0f ) return( im ); if( im == 0.0 ) 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.0); /* Undo the scaling */ b = ldexpf( b, e ); return( b ); } /* csqrtf() * * Complex square root * * * * SYNOPSIS: * * float complex csqrtf(); * float complex z, w; * * w = csqrtf( z ); * * * * DESCRIPTION: * * * If z = x + iy, r = |z|, then * * 1/2 * Re w = [ (r + x)/2 ] , * * 1/2 * Im w = [ (r - x)/2 ] . * * Cancellation error in r-x or r+x is avoided by using the * identity 2 Re w Im w = y. * * Note that -w is also a square root of z. The root chosen * is always in the right half plane and Im w has the same sign as y. * * * * ACCURACY: * * * Relative error: * arithmetic domain # trials peak rms * IEEE -10,+10 1,000,000 1.8e-7 3.5e-8 * */ float complex csqrtf( float complex z ) { float complex w; float x, y, r, t, scale; x = creal(z); y = cimag(z); if( y == 0.0f ) { if( x < 0.0f ) { w = 0.0f + sqrtf(-x) * I; return (w); } else if (x == 0.0f) { return (0.0f + y * I); } else { w = sqrtf(x) + y * I; return (w); } } if( x == 0.0f ) { r = fabsf(y); r = sqrtf(0.5f*r); if( y > 0 ) w = r + r * I; else w = r - r * I; return (w); } /* Rescale to avoid internal overflow or underflow. */ if ((fabsf(x) > 4.0f) || (fabsf(y) > 4.0f)) { x *= 0.25f; y *= 0.25f; scale = 2.0f; } else { #if 1 x *= 6.7108864e7f; /* 2^26 */ y *= 6.7108864e7f; scale = 1.220703125e-4f; /* 2^-13 */ #else x *= 4.0f; y *= 4.0f; scale = 0.5f; #endif } w = x + y * I; r = cabsf(w); if( x > 0 ) { t = sqrtf( 0.5f * r + 0.5f * x ); r = scale * fabsf( (0.5f * y) / t ); t *= scale; } else { r = sqrtf( 0.5f * r - 0.5f * x ); t = scale * fabsf( (0.5f * y) / r ); r *= scale; } if (y < 0) w = t - r * I; else w = t + r * I; return (w); }