1 28994509 2004-04-21 devnull #include <u.h>
 2 28994509 2004-04-21 devnull #include <libc.h>
 3 28994509 2004-04-21 devnull #include "map.h"
 4 28994509 2004-04-21 devnull 
 5 28994509 2004-04-21 devnull /*complex divide, defensive against overflow from
 6 28994509 2004-04-21 devnull  *	* and /, but not from + and -
 7 28994509 2004-04-21 devnull  *	assumes underflow yields 0.0
 8 28994509 2004-04-21 devnull  *	uses identities:
 9 28994509 2004-04-21 devnull  *	(a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c)
10 28994509 2004-04-21 devnull  *	(a + bi)/(c + di) = (b - ai)/(d - ci)
11 28994509 2004-04-21 devnull */
12 28994509 2004-04-21 devnull void
13 28994509 2004-04-21 devnull cdiv(double a, double b, double c, double d, double *u, double *v)
14 28994509 2004-04-21 devnull {
15 28994509 2004-04-21 devnull 	double r,t;
16 28994509 2004-04-21 devnull 	if(fabs(c)<fabs(d)) {
17 28994509 2004-04-21 devnull 		t = -c; c = d; d = t;
18 28994509 2004-04-21 devnull 		t = -a; a = b; b = t;
19 28994509 2004-04-21 devnull 	}
20 28994509 2004-04-21 devnull 	r = d/c;
21 28994509 2004-04-21 devnull 	t = c + r*d;
22 28994509 2004-04-21 devnull 	*u = (a + r*b)/t;
23 28994509 2004-04-21 devnull 	*v = (b - r*a)/t;
24 28994509 2004-04-21 devnull }
25 28994509 2004-04-21 devnull 
26 28994509 2004-04-21 devnull void
27 28994509 2004-04-21 devnull cmul(double c1, double c2, double d1, double d2, double *e1, double *e2)
28 28994509 2004-04-21 devnull {
29 28994509 2004-04-21 devnull 	*e1 = c1*d1 - c2*d2;
30 28994509 2004-04-21 devnull 	*e2 = c1*d2 + c2*d1;
31 28994509 2004-04-21 devnull }
32 28994509 2004-04-21 devnull 
33 28994509 2004-04-21 devnull void
34 28994509 2004-04-21 devnull csq(double c1, double c2, double *e1, double *e2)
35 28994509 2004-04-21 devnull {
36 28994509 2004-04-21 devnull 	*e1 = c1*c1 - c2*c2;
37 28994509 2004-04-21 devnull 	*e2 = c1*c2*2;
38 28994509 2004-04-21 devnull }
39 28994509 2004-04-21 devnull 
40 28994509 2004-04-21 devnull /* complex square root
41 28994509 2004-04-21 devnull  *	assumes underflow yields 0.0
42 28994509 2004-04-21 devnull  *	uses these identities:
43 28994509 2004-04-21 devnull  *	sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t))
44 28994509 2004-04-21 devnull  *	           = sqrt(r)(cos(t/2)+_isin(t/2))
45 28994509 2004-04-21 devnull  *	cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2)
46 28994509 2004-04-21 devnull  *	sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2)
47 28994509 2004-04-21 devnull */
48 28994509 2004-04-21 devnull void
49 28994509 2004-04-21 devnull csqrt(double c1, double c2, double *e1, double *e2)
50 28994509 2004-04-21 devnull {
51 28994509 2004-04-21 devnull 	double r,s;
52 28994509 2004-04-21 devnull 	double x,y;
53 28994509 2004-04-21 devnull 	x = fabs(c1);
54 28994509 2004-04-21 devnull 	y = fabs(c2);
55 28994509 2004-04-21 devnull 	if(x>=y) {
56 28994509 2004-04-21 devnull 		if(x==0) {
57 28994509 2004-04-21 devnull 			*e1 = *e2 = 0;
58 28994509 2004-04-21 devnull 			return;
59 28994509 2004-04-21 devnull 		}
60 28994509 2004-04-21 devnull 		r = x;
61 28994509 2004-04-21 devnull 		s = y/x;
62 28994509 2004-04-21 devnull 	} else {
63 28994509 2004-04-21 devnull 		r = y;
64 28994509 2004-04-21 devnull 		s = x/y;
65 28994509 2004-04-21 devnull 	}
66 28994509 2004-04-21 devnull 	r *= sqrt(1+ s*s);
67 28994509 2004-04-21 devnull 	if(c1>0) {
68 28994509 2004-04-21 devnull 		*e1 = sqrt((r+c1)/2);
69 28994509 2004-04-21 devnull 		*e2 = c2/(2* *e1);
70 28994509 2004-04-21 devnull 	} else {
71 28994509 2004-04-21 devnull 		*e2 = sqrt((r-c1)/2);
72 28994509 2004-04-21 devnull 		if(c2<0)
73 28994509 2004-04-21 devnull 			*e2 = -*e2;
74 28994509 2004-04-21 devnull 		*e1 = c2/(2* *e2);
75 28994509 2004-04-21 devnull 	}
76 28994509 2004-04-21 devnull }
77 28994509 2004-04-21 devnull 
78 28994509 2004-04-21 devnull 
79 28994509 2004-04-21 devnull void cpow(double c1, double c2, double *d1, double *d2, double pwr)
80 28994509 2004-04-21 devnull {
81 28994509 2004-04-21 devnull 	double theta = pwr*atan2(c2,c1);
82 28994509 2004-04-21 devnull 	double r = pow(hypot(c1,c2), pwr);
83 28994509 2004-04-21 devnull 	*d1 = r*cos(theta);
84 28994509 2004-04-21 devnull 	*d2 = r*sin(theta);
85 28994509 2004-04-21 devnull }