[goffice] Complex: improve accuracy of power.

[Morten Welinder <mortenw src gnome org>]

commit 6f49aadc32dae2e77fb135421be07d9b526e8b26
Author: Morten Welinder <terra gnome org>
Date:   Mon May 23 10:56:11 2011 -0400

    Complex: improve accuracy of power.

 ChangeLog                 |    6 +++
 NEWS                      |    1 +
 goffice/math/go-complex.c |   75 +++++++++++++++++++++++++++++++++++++++-----
 goffice/math/go-complex.h |    4 ++
 4 files changed, 77 insertions(+), 9 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index a3a2e70..ad44258 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,9 @@
+2011-05-23  Morten Welinder  <terra gnome org>
+
+	* goffice/math/go-complex.c (go_complex_invalid)
+	(go_complex_angle_pi): New functions.
+	(go_complex_pow): Improve accuracy.
+
 2011-05-22  Morten Welinder  <terra gnome org>
 
 	* goffice/math/go-complex.h (go_complexl): Only defined when
diff --git a/NEWS b/NEWS
index 88c1554..8108fa9 100644
--- a/NEWS
+++ b/NEWS
@@ -2,6 +2,7 @@ goffice 0.8.16:
 
 Morten:
 	* Fix long-double compilation issue.
+	* Improve accuracy of complex power.
 
 --------------------------------------------------------------------------
 goffice 0.8.15:
diff --git a/goffice/math/go-complex.c b/goffice/math/go-complex.c
index c96c543..52e10bf 100644
--- a/goffice/math/go-complex.c
+++ b/goffice/math/go-complex.c
@@ -18,6 +18,7 @@
 
 #define DOUBLE double
 #define SUFFIX(_n) _n
+#define M_PIgo    3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117
 #define go_strto go_strtod
 #define GO_const(_n) _n
 
@@ -25,6 +26,7 @@
 #include "go-complex.c"
 #undef DOUBLE
 #undef SUFFIX
+#undef M_PIgo
 #undef go_strto
 #undef GO_const
 
@@ -33,6 +35,7 @@
 #endif
 #define DOUBLE long double
 #define SUFFIX(_n) _n ## l
+#define M_PIgo    3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117L
 #define go_strto go_strtold
 #define GO_const(_n) _n ## L
 #endif
@@ -230,15 +233,43 @@ SUFFIX(go_complex_sqrt) (SUFFIX(go_complex) *dst, SUFFIX(go_complex) const *src)
 void
 SUFFIX(go_complex_pow) (SUFFIX(go_complex) *dst, SUFFIX(go_complex) const *a, SUFFIX(go_complex) const *b)
 {
-	SUFFIX(go_complex) lna, b_lna;
-
-	/* ln is not defined for reals less than or equal to zero.  */
-	if (SUFFIX(go_complex_real_p) (a) && SUFFIX(go_complex_real_p) (b))
-		SUFFIX(go_complex_init) (dst, SUFFIX(pow) (a->re, b->re), 0);
-	else {
-		SUFFIX(go_complex_ln) (&lna, a);
-		SUFFIX(go_complex_mul) (&b_lna, b, &lna);
-		SUFFIX(go_complex_exp) (dst, &b_lna);
+	if (SUFFIX(go_complex_zero_p) (a) && SUFFIX(go_complex_real_p) (b)) {
+		if (b->re <= 0)
+			SUFFIX(go_complex_invalid) (dst);
+		else
+			SUFFIX(go_complex_real) (dst, 0);
+	} else {
+		DOUBLE res_r, res_a1, res_a2, res_a2_pi, r, arg;
+		SUFFIX(go_complex) F;
+
+		SUFFIX(go_complex_to_polar) (&r, &arg, a);
+		res_r = SUFFIX(pow) (r, b->re) * SUFFIX(exp) (-b->im * arg);
+		res_a1 = b->im * SUFFIX(log) (r);
+		res_a2 = b->re * arg;
+		res_a2_pi = b->re * SUFFIX(go_complex_angle_pi) (a);
+
+		res_a2_pi = SUFFIX(fmod) (res_a2_pi, 2);
+		if (res_a2_pi < 0) res_a2_pi += 2;
+
+		/*
+		 * Problem: sometimes res_a2 is a nice fraction of pi.
+		 * Actually adding it will introduce pointless rounding
+		 * errors.
+		 */
+		if (res_a2_pi == 0.5) {
+			res_a2 = 0;
+			SUFFIX(go_complex_init) (&F, 0, 1);
+		} else if (res_a2_pi == 1) {
+			res_a2 = 0;
+			SUFFIX(go_complex_real) (&F, -1);
+		} else if (res_a2_pi == 1.5) {
+			res_a2 = 0;
+			SUFFIX(go_complex_init) (&F, 0, -1);
+		} else
+			SUFFIX(go_complex_real) (&F, 1);
+
+		SUFFIX(go_complex_from_polar) (dst, res_r, res_a1 + res_a2);
+		SUFFIX(go_complex_mul) (dst, dst, &F);
 	}
 }
 
@@ -252,6 +283,12 @@ void SUFFIX(go_complex_init) (SUFFIX(go_complex) *dst, DOUBLE re, DOUBLE im)
 
 /* ------------------------------------------------------------------------- */
 
+void SUFFIX(go_complex_invalid) (SUFFIX(go_complex) *dst)
+{
+	dst->re = SUFFIX(go_nan);
+	dst->im = SUFFIX(go_nan);
+}
+
 void SUFFIX(go_complex_real) (SUFFIX(go_complex) *dst, DOUBLE re)
 {
 	dst->re = re;
@@ -287,6 +324,26 @@ DOUBLE SUFFIX(go_complex_angle) (SUFFIX(go_complex) const *src)
 }
 
 /* ------------------------------------------------------------------------- */
+/*
+ * Same as go_complex_angle, but divided by pi (which occasionally produces
+ * nice round numbers not suffering from rounding errors).
+ */
+
+DOUBLE SUFFIX(go_complex_angle_pi) (SUFFIX(go_complex) const *src)
+{
+	if (src->im == 0)
+		return (src->re >= 0 ? 0 : -1);
+
+	if (src->re == 0)
+		return (src->im >= 0 ? 0.5 : -0.5);
+
+	/* We could do quarters too */
+
+	/* Fallback.  */
+	return SUFFIX(go_complex_angle) (src) / M_PIgo;
+}
+
+/* ------------------------------------------------------------------------- */
 
 void SUFFIX(go_complex_conj) (SUFFIX(go_complex) *dst, SUFFIX(go_complex) const *src)
 {
diff --git a/goffice/math/go-complex.h b/goffice/math/go-complex.h
index 1dc372f..2409913 100644
--- a/goffice/math/go-complex.h
+++ b/goffice/math/go-complex.h
@@ -31,11 +31,13 @@ void go_complex_div  (go_complex *dst, go_complex const *a, go_complex const *b)
 void go_complex_pow  (go_complex *dst, go_complex const *a, go_complex const *b);
 void go_complex_sqrt (go_complex *dst, go_complex const *src);
 void go_complex_init (go_complex *dst, double re, double im);
+void go_complex_invalid (go_complex *dst);
 void go_complex_real (go_complex *dst, double re);
 int go_complex_real_p (go_complex const *src);
 int go_complex_zero_p (go_complex const *src);
 double go_complex_mod (go_complex const *src);
 double go_complex_angle (go_complex const *src);
+double go_complex_angle_pi (go_complex const *src);
 void go_complex_conj (go_complex *dst, go_complex const *src);
 void go_complex_scale_real (go_complex *dst, double f);
 void go_complex_add (go_complex *dst, go_complex const *a, go_complex const *b);
@@ -63,11 +65,13 @@ void go_complex_divl  (go_complexl *dst, go_complexl const *a, go_complexl const
 void go_complex_powl  (go_complexl *dst, go_complexl const *a, go_complexl const *b);
 void go_complex_sqrtl (go_complexl *dst, go_complexl const *src);
 void go_complex_initl (go_complexl *dst, long double re, long double im);
+void go_complex_invalidl (go_complexl *dst);
 void go_complex_reall (go_complexl *dst, long double re);
 int go_complex_real_pl (go_complexl const *src);
 int go_complex_zero_pl (go_complexl const *src);
 long double go_complex_modl (go_complexl const *src);
 long double go_complex_anglel (go_complexl const *src);
+long double go_complex_angle_pil (go_complexl const *src);
 void go_complex_conjl (go_complexl *dst, go_complexl const *src);
 void go_complex_scale_reall (go_complexl *dst, long double f);
 void go_complex_addl (go_complexl *dst, go_complexl const *a, go_complexl const *b);