ref: 88ccea37f65b87312292b2b1a107c39927c6a2c3
dir: /sys/src/cmd/map/libmap/complex.c/
#include <u.h> #include <libc.h> #include "map.h" /*complex divide, defensive against overflow from * * and /, but not from + and - * assumes underflow yields 0.0 * uses identities: * (a + bi)/(c + di) = ((a + bd/c) + (b - ad/c)i)/(c + dd/c) * (a + bi)/(c + di) = (b - ai)/(d - ci) */ void cdiv(double a, double b, double c, double d, double *u, double *v) { double r,t; if(fabs(c)<fabs(d)) { t = -c; c = d; d = t; t = -a; a = b; b = t; } r = d/c; t = c + r*d; *u = (a + r*b)/t; *v = (b - r*a)/t; } void cmul(double c1, double c2, double d1, double d2, double *e1, double *e2) { *e1 = c1*d1 - c2*d2; *e2 = c1*d2 + c2*d1; } void csq(double c1, double c2, double *e1, double *e2) { *e1 = c1*c1 - c2*c2; *e2 = c1*c2*2; } /* complex square root * assumes underflow yields 0.0 * uses these identities: * sqrt(x+_iy) = sqrt(r(cos(t)+_isin(t)) * = sqrt(r)(cos(t/2)+_isin(t/2)) * cos(t/2) = sin(t)/2sin(t/2) = sqrt((1+cos(t)/2) * sin(t/2) = sin(t)/2cos(t/2) = sqrt((1-cos(t)/2) */ void csqrt(double c1, double c2, double *e1, double *e2) { double r,s; double x,y; x = fabs(c1); y = fabs(c2); if(x>=y) { if(x==0) { *e1 = *e2 = 0; return; } r = x; s = y/x; } else { r = y; s = x/y; } r *= sqrt(1+ s*s); if(c1>0) { *e1 = sqrt((r+c1)/2); *e2 = c2/(2* *e1); } else { *e2 = sqrt((r-c1)/2); if(c2<0) *e2 = -*e2; *e1 = c2/(2* *e2); } } void cpow(double c1, double c2, double *d1, double *d2, double pwr) { double theta = pwr*atan2(c2,c1); double r = pow(hypot(c1,c2), pwr); *d1 = r*cos(theta); *d2 = r*sin(theta); }