[goffice] Complex power: fewer special cases.

[Morten Welinder <mortenw src gnome org>]

commit eeb428bbc622a28663f0421e3f0ec570cf13ac32
Author: Morten Welinder <terra gnome org>
Date:   Mon Dec 23 13:56:11 2013 -0500

    Complex power: fewer special cases.
    
    The generic case is more accurate now.

 ChangeLog                 |    4 ++++
 goffice/math/go-complex.c |   44 ++++++++++++++++++++------------------------
 goffice/math/go-quad.c    |   14 +++++++++++++-
 3 files changed, 37 insertions(+), 25 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index af347bb..55c4292 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,7 @@
+2013-12-23  Morten Welinder  <terra gnome org>
+
+       * goffice/math/go-quad.c (go_quad_atan2_special): Do quarters too.
+
 2013-12-18  Morten Welinder  <terra gnome org>
 
        * goffice/math/go-quad.c (go_quad_sin, go_quad_sinpi)
diff --git a/goffice/math/go-complex.c b/goffice/math/go-complex.c
index 8f8e91b..f3f4a74 100644
--- a/goffice/math/go-complex.c
+++ b/goffice/math/go-complex.c
@@ -327,47 +327,41 @@ SUFFIX(mulmod1) (SUFFIX(GOQuad) *dst, SUFFIX(GOQuad) const *qa_, DOUBLE b)
 void
 SUFFIX(go_complex_pow) (COMPLEX *dst, COMPLEX const *a, COMPLEX const *b)
 {
-       if (SUFFIX(go_complex_real_p) (a) && a->re >= 0 &&
-           SUFFIX(go_complex_real_p) (b)) {
-               SUFFIX(go_complex_init) (dst, pow (a->re, b->re), 0);
-       } else if (SUFFIX(go_complex_real_p) (b)) {
-               DOUBLE p = SUFFIX(fabs) (b->re);
-               COMPLEX t;
-
-               if (p == 0) {
+       if (b->im == 0) {
+               if (SUFFIX(go_complex_real_p) (a) && a->re >= 0) {
+                       SUFFIX(go_complex_init) (dst, pow (a->re, b->re), 0);
+                       return;
+               }
+
+               if (b->re == 0) {
                        SUFFIX(go_complex_init) (dst, 1, 0);
                        return;
                }
 
-               if (p == 1) {
-                       t = *a;
-               } else {
-                       DOUBLE arg = SUFFIX(go_complex_angle_pi) (a);
-                       DOUBLE r = (p >= 2)
-                               ? SUFFIX(pow) (a->re * a->re +
-                                              a->im * a->im, p / 2)
-                               : SUFFIX(pow) (SUFFIX(hypot) (a->re, a->im), p);
-                       SUFFIX(go_complex_from_polar_pi) (&t, r, arg * p);
+               if (b->re == 1) {
+                       *dst = *a;
+                       return;
                }
 
-               if (b->re > 0)
-                       *dst = t;
-               else {
-                       COMPLEX one;
-                       one.re = 1; one.im = 0;
-                       SUFFIX(go_complex_div) (dst, &one, &t);
+               if (b->re == 2) {
+                       SUFFIX(go_complex_mul) (dst, a, a);
+                       return;
                }
-       } else {
+       }
+
+       {
                DOUBLE e1, e2;
                int e;
                SUFFIX(GOQuad) qr, qa, qb, qarg, qrr, qra;
                void *state = SUFFIX(go_quad_start) ();
 
+               /* Convert a to polar coordinates.  */
                SUFFIX(go_quad_init) (&qa, a->re);
                SUFFIX(go_quad_init) (&qb, a->im);
                SUFFIX(go_quad_atan2pi) (&qarg, &qb, &qa);
                SUFFIX(go_quad_hypot) (&qr, &qa, &qb);
 
+               /* Compute result modulus.  */
                SUFFIX(go_quad_init) (&qa, b->re);
                SUFFIX(go_quad_pow) (&qa, &e1, &qr, &qa);
                SUFFIX(go_quad_init) (&qb, -b->im);
@@ -379,6 +373,7 @@ SUFFIX(go_complex_pow) (COMPLEX *dst, COMPLEX const *a, COMPLEX const *b)
                qrr.h = SUFFIX(ldexp) (qrr.h, e);
                qrr.l = SUFFIX(ldexp) (qrr.l, e);
 
+               /* Compute result angle.  */
                SUFFIX(go_quad_log) (&qa, &qr);
                SUFFIX(go_quad_div) (&qa, &qa, &SUFFIX(go_quad_2pi));
                SUFFIX(mulmod1) (&qa, &qa, b->im);
@@ -386,6 +381,7 @@ SUFFIX(go_complex_pow) (COMPLEX *dst, COMPLEX const *a, COMPLEX const *b)
                SUFFIX(go_quad_add) (&qa, &qa, &qb);
                SUFFIX(go_quad_add) (&qra, &qa, &qa);
 
+               /* Convert to rectangular coordinates.  */
                SUFFIX(go_quad_sinpi) (&qa, &qra);
                SUFFIX(go_quad_mul) (&qa, &qa, &qrr);
                SUFFIX(go_quad_cospi) (&qb, &qra);
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index 0b2f883..db4ecc4 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -921,7 +921,19 @@ SUFFIX(go_quad_atan2_special) (const QUAD *y, const QUAD *x, DOUBLE *f)
                return TRUE;
        }
 
-       /* We could do quarters too */
+       if (SUFFIX(fabs) (SUFFIX(fabs)(dx) - SUFFIX(fabs)(dy)) < 1e-10) {
+               QUAD d;
+               SUFFIX(go_quad_sub) (&d, x, y);
+               if (d.h == 0) {
+                       *f = (dy >= 0 ? 0.25 : -0.75);
+                       return TRUE;
+               }
+               SUFFIX(go_quad_add) (&d, x, y);
+               if (d.h == 0) {
+                       *f = (dy >= 0 ? -0.25 : +0.75);
+                       return TRUE;
+               }
+       }
 
        return FALSE;
 }