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