diff options
Diffstat (limited to 'src/cmd/map/libmap/complex.c')
-rw-r--r-- | src/cmd/map/libmap/complex.c | 85 |
1 files changed, 85 insertions, 0 deletions
diff --git a/src/cmd/map/libmap/complex.c b/src/cmd/map/libmap/complex.c new file mode 100644 index 00000000..b3099fd4 --- /dev/null +++ b/src/cmd/map/libmap/complex.c @@ -0,0 +1,85 @@ +#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); +} |