From 28994509cc11ac6a5443054dfae1fedfb69039bc Mon Sep 17 00:00:00 2001
From: rsc <devnull@localhost>
Date: Wed, 21 Apr 2004 22:19:33 +0000
Subject: Why not?

---
 src/cmd/map/libmap/complex.c | 85 ++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 85 insertions(+)
 create mode 100644 src/cmd/map/libmap/complex.c

(limited to 'src/cmd/map/libmap/complex.c')

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);
+}
-- 
cgit v1.2.3