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