3 * Complex number arithmetic
11 * float i; imaginary part
16 * caddf( a, b, c ); c = b + a
17 * csubf( a, b, c ); c = b - a
18 * cmulf( a, b, c ); c = b * a
19 * cdivf( a, b, c ); c = b / a
21 * cmovf( b, c ); c = b
36 * c.r = b.r * a.r - b.i * a.i
37 * c.i = b.r * a.i + b.i * a.r
40 * d = a.r * a.r + a.i * a.i
41 * c.r = (b.r * a.r + b.i * a.i)/d
42 * c.i = (b.i * a.r - b.r * a.i)/d
45 * In DEC arithmetic, the test (1/z) * z = 1 had peak relative
46 * error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had
47 * peak relative error 8.3e-17, rms 2.1e-17.
49 * Tests in the rectangle {-10,+10}:
51 * arithmetic function # trials peak rms
52 * IEEE cadd 30000 5.9e-8 2.6e-8
53 * IEEE csub 30000 6.0e-8 2.6e-8
54 * IEEE cmul 30000 1.1e-7 3.7e-8
55 * IEEE cdiv 30000 2.1e-7 5.7e-8
58 * complex number arithmetic
63 Cephes Math Library Release 2.1: December, 1988
64 Copyright 1984, 1987, 1988 by Stephen L. Moshier
65 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
70 extern float MAXNUMF, MACHEPF, PIF, PIO2F;
71 #define fabsf(x) ( (x) < 0 ? -(x) : (x) )
73 float sqrtf(float), frexpf(float, int *);
74 float ldexpf(float, int);
75 float cabsf(cmplxf *), atan2f(float, float), cosf(float), sinf(float);
77 float sqrtf(), frexpf(), ldexpf();
78 float cabsf(), atan2f(), cosf(), sinf();
87 cmplxf czerof = {0.0, 0.0};
89 cmplxf conef = {1.0, 0.0};
95 register cmplxf *a, *b;
106 void csubf( a, b, c )
107 register cmplxf *a, *b;
117 void cmulf( a, b, c )
118 register cmplxf *a, *b;
123 y = b->r * a->r - b->i * a->i;
124 c->i = b->r * a->i + b->i * a->r;
132 void cdivf( a, b, c )
133 register cmplxf *a, *b;
139 y = a->r * a->r + a->i * a->i;
140 p = b->r * a->r + b->i * a->i;
141 q = b->i * a->r - b->r * a->i;
146 if( (fabsf(p) > w) || (fabsf(q) > w) || (y == 0.0f) )
150 mtherr( "cdivf", OVERFLOW );
162 register short *a, *b;
184 * Complex absolute value
205 * a = sqrt( x**2 + y**2 ).
207 * Overflow and underflow are avoided by testing the magnitudes
208 * of x and y before squaring. If either is outside half of
209 * the floating point full scale range, both are rescaled.
215 * arithmetic domain # trials peak rms
216 * IEEE -10,+10 30000 1.2e-7 3.4e-8
221 Cephes Math Library Release 2.1: January, 1989
222 Copyright 1984, 1987, 1989 by Stephen L. Moshier
223 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
234 /* square root of max and min numbers */
235 #define SMAX 1.3043817825332782216E+19
236 #define SMIN 7.6664670834168704053E-20
241 #define SMAXT (2.0f * SMAX)
242 #define SMINT (0.5f * SMIN)
247 float x, y, b, re, im;
262 /* Get the exponents of the numbers */
263 x = frexpf( re, &ex );
264 y = frexpf( im, &ey );
266 /* Check if one number is tiny compared to the other */
273 /* Find approximate exponent e of the geometric mean. */
276 /* Rescale so mean is about 1 */
277 x = ldexpf( re, -e );
278 y = ldexpf( im, -e );
280 /* Hypotenuse of the right triangle */
281 b = sqrtf( x * x + y * y );
283 /* Compute the exponent of the answer. */
284 y = frexpf( b, &ey );
287 /* Check it for overflow and underflow. */
290 mtherr( "cabsf", OVERFLOW );
296 /* Undo the scaling */
302 * Complex square root
318 * If z = x + iy, r = |z|, then
321 * Im w = [ (r - x)/2 ] ,
326 * Note that -w is also a square root of z. The solution
327 * reported is always in the upper half plane.
329 * Because of the potential for cancellation error in r - x,
330 * the result is sharpened by doing a Heron iteration
331 * (see sqrt.c) in complex arithmetic.
338 * arithmetic domain # trials peak rms
339 * IEEE -10,+10 100000 1.8e-7 4.2e-8
381 /* Approximate sqrt(x^2+y^2) - x = y^2/2x - y^4/24x^3 + ... .
382 * The relative error in the first term is approximately y^2/12x^2 .
384 if( (fabsf(y) < fabsf(0.015f*x))
399 /* Heron iteration in complex arithmetic: