[gnumeric] IMPOWER: fix problem on real axis and improve accuracy.

[Morten Welinder <mortenw src gnome org>]

commit db1824f94cc6b2206a1e99a536637b2e502917f2
Author: Morten Welinder <terra gnome org>
Date:   Sun May 2 22:37:44 2010 -0400

    IMPOWER: fix problem on real axis and improve accuracy.

 ChangeLog     |    2 +
 NEWS          |    1 +
 src/complex.c |   62 +++++++++++++++++++++++++++++++++++++++++++++++---------
 3 files changed, 55 insertions(+), 10 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 62c7fdc..71b6ee4 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,7 @@
 2010-05-02  Morten Welinder  <terra gnome org>
 
+	* src/complex.c (complex_pow): Fix problem on negative real axis.
+
 	* src/complex.h (complex_invalid, complex_invalid_p): New
 	functions.
 
diff --git a/NEWS b/NEWS
index b4b6db1..db72b98 100644
--- a/NEWS
+++ b/NEWS
@@ -26,6 +26,7 @@ Morten:
 	* Fix advanced-filter problem with multiple criteria.  [#164169]
 	* Fix problem parsing comples numbers.
 	* Improve error handle for complex functions.
+	* Improve IMPOWER.
 
 --------------------------------------------------------------------------
 Gnumeric 1.10.2
diff --git a/src/complex.c b/src/complex.c
index 7642dc6..7e8a4cc 100644
--- a/src/complex.c
+++ b/src/complex.c
@@ -13,7 +13,7 @@
 
 #include <stdlib.h>
 #include <errno.h>
-
+#include <mathfunc.h>
 
 /* ------------------------------------------------------------------------- */
 
@@ -206,18 +206,60 @@ complex_sqrt (complex_t *dst, complex_t const *src)
 
 /* ------------------------------------------------------------------------- */
 
+/* Like complex_angle, but divide result by pi.  */
+static gnm_float
+complex_angle_pi (complex_t 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 complex_angle (src) / M_PIgnum;
+}
+
+
 void
 complex_pow (complex_t *dst, complex_t const *a, complex_t const *b)
 {
-	complex_t lna, b_lna;
-
-	/* ln is not defined for reals less than or equal to zero.  */
-	if (complex_real_p (a) && complex_real_p (b))
-		complex_init (dst, gnm_pow (a->re, b->re), 0);
-	else {
-		complex_ln (&lna, a);
-		complex_mul (&b_lna, b, &lna);
-		complex_exp (dst, &b_lna);
+	if (complex_zero_p (a) && complex_zero_p (b)) {
+		complex_invalid (dst);
+	} else {
+		gnm_float res_r, res_a1, res_a2, res_a2_pi, r, arg;
+		complex_t F;
+
+		complex_to_polar (&r, &arg, a);
+		res_r = gnm_pow (r, b->re) * gnm_exp (-b->im * arg);
+		res_a1 = b->im * gnm_log (r);
+		res_a2 = b->re * arg;
+		res_a2_pi = b->re * complex_angle_pi (a);
+
+		res_a2_pi = gnm_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;
+			complex_init (&F, 0, 1);
+		} else if (res_a2_pi == 1) {
+			res_a2 = 0;
+			complex_real (&F, -1);
+		} else if (res_a2_pi == 1.5) {
+			res_a2 = 0;
+			complex_init (&F, 0, -1);
+		} else
+			complex_real (&F, 1);
+
+		complex_from_polar (dst, res_r, res_a1 + res_a2);
+		complex_mul (dst, dst, &F);
 	}
 }