1 /* cmplx.c
2 *
3 * Complex number arithmetic
4 *
5 *
6 *
7 * SYNOPSIS:
8 *
9 * typedef struct {
10 * double r; real part
11 * double i; imaginary part
12 * }cmplx;
13 *
14 * cmplx *a, *b, *c;
15 *
16 * cadd( a, b, c ); c = b + a
17 * csub( a, b, c ); c = b - a
18 * cmul( a, b, c ); c = b * a
19 * cdiv( a, b, c ); c = b / a
20 * cneg( c ); c = -c
21 * cmov( b, c ); c = b
22 *
23 *
24 *
25 * DESCRIPTION:
26 *
27 * Addition:
28 * c.r = b.r + a.r
29 * c.i = b.i + a.i
30 *
31 * Subtraction:
32 * c.r = b.r - a.r
33 * c.i = b.i - a.i
34 *
35 * Multiplication:
36 * c.r = b.r * a.r - b.i * a.i
37 * c.i = b.r * a.i + b.i * a.r
38 *
39 * Division:
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
43 * ACCURACY:
44 *
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.
48 *
49 * Tests in the rectangle {-10,+10}:
50 * Relative error:
51 * arithmetic function # trials peak rms
52 * DEC cadd 10000 1.4e-17 3.4e-18
53 * IEEE cadd 100000 1.1e-16 2.7e-17
54 * DEC csub 10000 1.4e-17 4.5e-18
55 * IEEE csub 100000 1.1e-16 3.4e-17
56 * DEC cmul 3000 2.3e-17 8.7e-18
57 * IEEE cmul 100000 2.1e-16 6.9e-17
58 * DEC cdiv 18000 4.9e-17 1.3e-17
59 * IEEE cdiv 100000 3.7e-16 1.1e-16
60 */
61 /* cmplx.c
62 * complex number arithmetic
63 */
64
65
66 /*
67 Cephes Math Library Release 2.8: June, 2000
68 Copyright 1984, 1995, 2000 by Stephen L. Moshier
69 */
70
71
72 #include "mconf.h"
73
74 #ifdef ANSIPROT
75 extern double md_fabs ( double );
76 extern double md_cabs ( cmplx * );
77 extern double sqrt ( double );
78 extern double md_atan2 ( double, double );
79 extern double md_cos ( double );
80 extern double md_sin ( double );
81 extern double sqrt ( double );
82 extern double md_frexp ( double, int * );
83 extern double md_ldexp ( double, int );
84 int isnan ( double );
85 void cdiv ( cmplx *, cmplx *, cmplx * );
86 void cadd ( cmplx *, cmplx *, cmplx * );
87 #else
88 double md_fabs(), md_cabs(), sqrt(), md_atan2(), md_cos(), md_sin();
89 double sqrt(), md_frexp(), md_ldexp();
90 int isnan();
91 void cdiv(), cadd();
92 #endif
93
94 extern double MAXNUM, MACHEP, PI, PIO2, INFINITY, NAN;
95 /*
96 typedef struct
97 {
98 double r;
99 double i;
100 }cmplx;
101 */
102 cmplx czero = {0.0, 0.0};
103 extern cmplx czero;
104 cmplx cone = {1.0, 0.0};
105 extern cmplx cone;
106
107 /* c = b + a */
108
cadd(a,b,c)109 void cadd( a, b, c )
110 register cmplx *a, *b;
111 cmplx *c;
112 {
113
114 c->r = b->r + a->r;
115 c->i = b->i + a->i;
116 }
117
118
119 /* c = b - a */
120
csub(a,b,c)121 void csub( a, b, c )
122 register cmplx *a, *b;
123 cmplx *c;
124 {
125
126 c->r = b->r - a->r;
127 c->i = b->i - a->i;
128 }
129
130 /* c = b * a */
131
cmul(a,b,c)132 void cmul( a, b, c )
133 register cmplx *a, *b;
134 cmplx *c;
135 {
136 double y;
137
138 y = b->r * a->r - b->i * a->i;
139 c->i = b->r * a->i + b->i * a->r;
140 c->r = y;
141 }
142
143
144
145 /* c = b / a */
146
cdiv(a,b,c)147 void cdiv( a, b, c )
148 register cmplx *a, *b;
149 cmplx *c;
150 {
151 double y, p, q, w;
152
153
154 y = a->r * a->r + a->i * a->i;
155 p = b->r * a->r + b->i * a->i;
156 q = b->i * a->r - b->r * a->i;
157
158 if( y < 1.0 )
159 {
160 w = MAXNUM * y;
161 if( (md_fabs(p) > w) || (md_fabs(q) > w) || (y == 0.0) )
162 {
163 c->r = MAXNUM;
164 c->i = MAXNUM;
165 mtherr( "cdiv", OVERFLOW );
166 return;
167 }
168 }
169 c->r = p/y;
170 c->i = q/y;
171 }
172
173
174 /* b = a
175 Caution, a `short' is assumed to be 16 bits wide. */
176
cmov(a,b)177 void cmov( a, b )
178 void *a, *b;
179 {
180 register short *pa, *pb;
181 int i;
182
183 pa = (short *) a;
184 pb = (short *) b;
185 i = 8;
186 do
187 *pb++ = *pa++;
188 while( --i );
189 }
190
191
cneg(a)192 void cneg( a )
193 register cmplx *a;
194 {
195
196 a->r = -a->r;
197 a->i = -a->i;
198 }
199
200 /* md_cabs()
201 *
202 * Complex absolute value
203 *
204 *
205 *
206 * SYNOPSIS:
207 *
208 * double md_cabs();
209 * cmplx z;
210 * double a;
211 *
212 * a = md_cabs( &z );
213 *
214 *
215 *
216 * DESCRIPTION:
217 *
218 *
219 * If z = x + iy
220 *
221 * then
222 *
223 * a = sqrt( x**2 + y**2 ).
224 *
225 * Overflow and underflow are avoided by testing the magnitudes
226 * of x and y before squaring. If either is outside half of
227 * the floating point full scale range, both are rescaled.
228 *
229 *
230 * ACCURACY:
231 *
232 * Relative error:
233 * arithmetic domain # trials peak rms
234 * DEC -30,+30 30000 3.2e-17 9.2e-18
235 * IEEE -10,+10 100000 2.7e-16 6.9e-17
236 */
237
238
239 /*
240 Cephes Math Library Release 2.1: January, 1989
241 Copyright 1984, 1987, 1989 by Stephen L. Moshier
242 Direct inquiries to 30 Frost Street, Cambridge, MA 02140
243 */
244
245
246 /*
247 typedef struct
248 {
249 double r;
250 double i;
251 }cmplx;
252 */
253
254 #ifdef UNK
255 #define PREC 27
256 #define MAXEXP 1024
257 #define MINEXP -1077
258 #endif
259 #ifdef DEC
260 #define PREC 29
261 #define MAXEXP 128
262 #define MINEXP -128
263 #endif
264 #ifdef IBMPC
265 #define PREC 27
266 #define MAXEXP 1024
267 #define MINEXP -1077
268 #endif
269 #ifdef MIEEE
270 #define PREC 27
271 #define MAXEXP 1024
272 #define MINEXP -1077
273 #endif
274
275
md_cabs(z)276 double md_cabs( z )
277 register cmplx *z;
278 {
279 double x, y, b, re, im;
280 int ex, ey, e;
281
282 #ifdef INFINITIES
283 /* Note, md_cabs(INFINITY,NAN) = INFINITY. */
284 if( z->r == INFINITY || z->i == INFINITY
285 || z->r == -INFINITY || z->i == -INFINITY )
286 return( INFINITY );
287 #endif
288
289 #ifdef NANS
290 if( isnan(z->r) )
291 return(z->r);
292 if( isnan(z->i) )
293 return(z->i);
294 #endif
295
296 re = md_fabs( z->r );
297 im = md_fabs( z->i );
298
299 if( re == 0.0 )
300 return( im );
301 if( im == 0.0 )
302 return( re );
303
304 /* Get the exponents of the numbers */
305 x = md_frexp( re, &ex );
306 y = md_frexp( im, &ey );
307
308 /* Check if one number is tiny compared to the other */
309 e = ex - ey;
310 if( e > PREC )
311 return( re );
312 if( e < -PREC )
313 return( im );
314
315 /* Find approximate exponent e of the geometric mean. */
316 e = (ex + ey) >> 1;
317
318 /* Rescale so mean is about 1 */
319 x = md_ldexp( re, -e );
320 y = md_ldexp( im, -e );
321
322 /* Hypotenuse of the right triangle */
323 b = sqrt( x * x + y * y );
324
325 /* Compute the exponent of the answer. */
326 y = md_frexp( b, &ey );
327 ey = e + ey;
328
329 /* Check it for overflow and underflow. */
330 if( ey > MAXEXP )
331 {
332 mtherr( "md_cabs", OVERFLOW );
333 return( INFINITY );
334 }
335 if( ey < MINEXP )
336 return(0.0);
337
338 /* Undo the scaling */
339 b = md_ldexp( b, e );
340 return( b );
341 }
342 /* md_csqrt()
343 *
344 * Complex square root
345 *
346 *
347 *
348 * SYNOPSIS:
349 *
350 * void md_csqrt();
351 * cmplx z, w;
352 *
353 * md_csqrt( &z, &w );
354 *
355 *
356 *
357 * DESCRIPTION:
358 *
359 *
360 * If z = x + iy, r = |z|, then
361 *
362 * 1/2
363 * Im w = [ (r - x)/2 ] ,
364 *
365 * Re w = y / 2 Im w.
366 *
367 *
368 * Note that -w is also a square root of z. The root chosen
369 * is always in the upper half plane.
370 *
371 * Because of the potential for cancellation error in r - x,
372 * the result is sharpened by doing a Heron iteration
373 * (see sqrt.c) in complex arithmetic.
374 *
375 *
376 *
377 * ACCURACY:
378 *
379 * Relative error:
380 * arithmetic domain # trials peak rms
381 * DEC -10,+10 25000 3.2e-17 9.6e-18
382 * IEEE -10,+10 100000 3.2e-16 7.7e-17
383 *
384 * 2
385 * Also tested by md_csqrt( z ) = z, and tested by arguments
386 * close to the real axis.
387 */
388
389
md_csqrt(z,w)390 void md_csqrt( z, w )
391 cmplx *z, *w;
392 {
393 cmplx q, s;
394 double x, y, r, t;
395
396 x = z->r;
397 y = z->i;
398
399 if( y == 0.0 )
400 {
401 if( x < 0.0 )
402 {
403 w->r = 0.0;
404 w->i = sqrt(-x);
405 return;
406 }
407 else
408 {
409 w->r = sqrt(x);
410 w->i = 0.0;
411 return;
412 }
413 }
414
415
416 if( x == 0.0 )
417 {
418 r = md_fabs(y);
419 r = sqrt(0.5*r);
420 if( y > 0 )
421 w->r = r;
422 else
423 w->r = -r;
424 w->i = r;
425 return;
426 }
427
428 /* Approximate sqrt(x^2+y^2) - x = y^2/2x - y^4/24x^3 + ... .
429 * The relative error in the first term is approximately y^2/12x^2 .
430 */
431 if( (md_fabs(y) < 2.e-4 * md_fabs(x))
432 && (x > 0) )
433 {
434 t = 0.25*y*(y/x);
435 }
436 else
437 {
438 r = md_cabs(z);
439 t = 0.5*(r - x);
440 }
441
442 r = sqrt(t);
443 q.i = r;
444 q.r = y/(2.0*r);
445 /* Heron iteration in complex arithmetic */
446 cdiv( &q, z, &s );
447 cadd( &q, &s, w );
448 w->r *= 0.5;
449 w->i *= 0.5;
450 }
451
452
md_hypot(x,y)453 double md_hypot( x, y )
454 double x, y;
455 {
456 cmplx z;
457
458 z.r = x;
459 z.i = y;
460 return( md_cabs(&z) );
461 }
462