complex.c (1592B)
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 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 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 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 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 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 }