[goffice] Complex: fix power function.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [goffice] Complex: fix power function.
- Date: Thu, 5 Dec 2013 21:17:01 +0000 (UTC)
commit 4bb8df511efca735c233b8af92f050c93640c55e
Author: Morten Welinder <terra gnome org>
Date: Thu Dec 5 16:16:27 2013 -0500
Complex: fix power function.
This should take care of accuracy in the complex power function.
ChangeLog | 3 +-
NEWS | 4 +-
goffice/math/go-complex.c | 105 +++++++++++++++++++++++-----------
goffice/math/go-quad.c | 141 +++++++++++++++++++++++++++++++++++++++++++++
goffice/math/go-quad.h | 4 +
tests/test-quad.c | 76 ++++++++++++++++++++----
6 files changed, 285 insertions(+), 48 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index c5abb55..dec7cdb 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,9 +1,10 @@
2013-12-05 Morten Welinder <terra gnome org>
* goffice/math/go-quad.c (go_quad_pi): New constant.
- (go_quad_log): New function.
+ (go_quad_log, go_quad_atan2, go_quad_atan2pi): New functions.
* tests/test-quad.c (main): Check value of go_quad_2pi.
+ Test go_quad_atan2.
2013-12-04 Morten Welinder <terra gnome org>
diff --git a/NEWS b/NEWS
index be39466..c21d8d4 100644
--- a/NEWS
+++ b/NEWS
@@ -4,8 +4,8 @@ Jean
* Fix graph guru appearance when used with gtk+-3.10. [#719681]
Morten:
- * Improve complex math.
- * Add quad precision log function.
+ * Improve complex math, notably complex power.
+ * Add quad precision log, atan2 functions.
--------------------------------------------------------------------------
goffice 0.10.9:
diff --git a/goffice/math/go-complex.c b/goffice/math/go-complex.c
index afe7d40..bb227b5 100644
--- a/goffice/math/go-complex.c
+++ b/goffice/math/go-complex.c
@@ -44,6 +44,7 @@
#undef M_PIgo
#undef go_strto
#undef GO_const
+#define LONG_DOUBLE_VERSION
#ifdef HAVE_SUNMATH_H
#include <sunmath.h>
@@ -193,6 +194,35 @@ SUFFIX(go_complex_from_polar) (SUFFIX(GOComplex) *dst, DOUBLE mod, DOUBLE angle)
SUFFIX(go_complex_init) (dst, mod * SUFFIX(cos) (angle), mod * SUFFIX(sin) (angle));
}
+static void
+SUFFIX(go_complex_from_polar_pi) (SUFFIX(GOComplex) *dst, DOUBLE mod, DOUBLE angle)
+{
+ DOUBLE s, c;
+
+ if (SUFFIX(fabs) (angle) >= 1) {
+ angle = SUFFIX(fmod) (angle, 2);
+ if (angle > 1)
+ angle -= 2;
+ else if (angle <= -1)
+ angle += 2;
+ }
+
+ if (angle == 0)
+ s = 0, c = 1;
+ else if (angle == 0.5)
+ s = 1, c = 0;
+ else if (angle == 1)
+ s = 0, c = -1;
+ else if (angle == -0.5)
+ s = -1, c = 0;
+ else {
+ s = SUFFIX(sin) (angle * M_PIgo);
+ c = SUFFIX(cos) (angle * M_PIgo);
+ }
+
+ SUFFIX(go_complex_init) (dst, mod * c, mod * s);
+}
+
/* ------------------------------------------------------------------------- */
void
@@ -238,9 +268,9 @@ SUFFIX(go_complex_sqrt) (SUFFIX(GOComplex) *dst, SUFFIX(GOComplex) const *src)
else
SUFFIX(go_complex_init) (dst, 0, SUFFIX(sqrt) (-src->re));
} else
- SUFFIX(go_complex_from_polar) (dst,
+ SUFFIX(go_complex_from_polar_pi) (dst,
SUFFIX(sqrt) (SUFFIX(go_complex_mod) (src)),
- SUFFIX(go_complex_angle) (src) / 2);
+ SUFFIX(go_complex_angle_pi) (src) / 2);
}
/* ------------------------------------------------------------------------- */
@@ -254,42 +284,49 @@ SUFFIX(go_complex_pow) (SUFFIX(GOComplex) *dst, SUFFIX(GOComplex) const *a, SUFF
else
SUFFIX(go_complex_real) (dst, 0);
} else {
- DOUBLE res_r, res_a1, res_a2, res_a2_pi, r, rre, arg;
- SUFFIX(GOComplex) F;
+ DOUBLE res_r, res_a;
+ SUFFIX(GOQuad) qr, qa, qb, qarg;
+ void *state = SUFFIX(go_quad_start) ();
+
+ SUFFIX(go_quad_init) (&qa, a->re);
+ SUFFIX(go_quad_init) (&qb, a->im);
+ SUFFIX(go_quad_atan2pi) (&qarg, &qb, &qa);
+ SUFFIX(go_quad_mul) (&qa, &qa, &qa);
+ SUFFIX(go_quad_mul) (&qb, &qb, &qb);
+ SUFFIX(go_quad_add) (&qa, &qa, &qb);
+ SUFFIX(go_quad_sqrt) (&qr, &qa);
- SUFFIX(go_complex_to_polar) (&r, &arg, a);
/*
* This is the square root of the power we really want,
* but it is much less likely to cause overflow.
*/
- rre = SUFFIX(pow) (r, b->re / 2);
- res_r = rre * SUFFIX(exp) (-b->im * arg) * rre;
- 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);
+ SUFFIX(go_quad_init) (&qa, b->re / 2);
+ SUFFIX(go_quad_pow) (&qa, NULL, &qr, &qa);
+ SUFFIX(go_quad_init) (&qb, -b->im);
+ SUFFIX(go_quad_mul) (&qb, &qb, &qarg);
+ SUFFIX(go_quad_mul) (&qb, &qb, &SUFFIX(go_quad_pi));
+ SUFFIX(go_quad_exp) (&qb, NULL, &qb);
+ SUFFIX(go_quad_mul) (&qb, &qa, &qb);
+ SUFFIX(go_quad_mul) (&qb, &qa, &qb);
+ res_r = SUFFIX(go_quad_value) (&qb);
+
+ SUFFIX(go_quad_log) (&qa, &qr);
+ SUFFIX(go_quad_init) (&qb, b->im);
+ SUFFIX(go_quad_mul) (&qa, &qa, &qb);
+ SUFFIX(go_quad_div) (&qa, &qa, &SUFFIX(go_quad_2pi));
+ SUFFIX(go_quad_init) (&qb, b->re / 2);
+ SUFFIX(go_quad_mul) (&qb, &qb, &qarg);
+ SUFFIX(go_quad_add) (&qa, &qa, &qb);
+ SUFFIX(go_quad_init) (&qb, 0.5);
+ SUFFIX(go_quad_add) (&qb, &qb, &qa);
+ SUFFIX(go_quad_floor) (&qb, &qb);
+ SUFFIX(go_quad_sub) (&qa, &qa, &qb);
+ SUFFIX(go_quad_add) (&qa, &qa, &qa);
+ res_a = SUFFIX(go_quad_value) (&qa);
+
+ SUFFIX(go_complex_from_polar_pi) (dst, res_r, res_a);
+
+ SUFFIX(go_quad_end) (state);
}
}
@@ -445,3 +482,5 @@ void SUFFIX(go_complex_tan) (SUFFIX(GOComplex) *dst, SUFFIX(GOComplex) const *sr
SUFFIX(go_complex_cos) (&c, src);
SUFFIX(go_complex_div) (dst, &s, &c);
}
+
+/* ------------------------------------------------------------------------- */
diff --git a/goffice/math/go-quad.c b/goffice/math/go-quad.c
index 87fda74..0f27620 100644
--- a/goffice/math/go-quad.c
+++ b/goffice/math/go-quad.c
@@ -772,3 +772,144 @@ SUFFIX(go_quad_log) (QUAD *res, const QUAD *a)
*res = xi;
}
}
+
+static void
+SUFFIX(go_quad_atan_internal) (QUAD *res, const QUAD *x)
+{
+ QUAD g, gp, dk[20], dpk[20], qr, qrp;
+ int n, k;
+ gboolean converged = FALSE;
+
+ g_return_if_fail (SUFFIX(fabs) (x->h) <= 1);
+
+ qrp = SUFFIX(go_quad_zero);
+
+ dpk[0] = SUFFIX(go_quad_one);
+ SUFFIX(go_quad_mul) (&gp, x, x);
+ SUFFIX(go_quad_add) (&gp, &gp, &SUFFIX(go_quad_one));
+ SUFFIX(go_quad_sqrt) (&gp, &gp);
+
+ for (n = 1; n < (int)G_N_ELEMENTS(dk); n++) {
+ SUFFIX(go_quad_add) (&dk[0], &dpk[0], &gp);
+ dk[0].h *= 0.5; dk[0].l *= 0.5;
+
+ SUFFIX(go_quad_mul) (&g, &dk[0], &gp);
+ SUFFIX(go_quad_sqrt) (&g, &g);
+
+ for (k = 1; k <= n; k++) {
+ QUAD f;
+
+ SUFFIX(go_quad_init) (&f, SUFFIX(ldexp) (1, -(2 * k)));
+ SUFFIX(go_quad_mul) (&dk[k], &f, &dpk[k-1]);
+ SUFFIX(go_quad_sub) (&dk[k], &dk[k-1], &dk[k]);
+ SUFFIX(go_quad_sub) (&f, &SUFFIX(go_quad_one), &f);
+ SUFFIX(go_quad_div) (&dk[k], &dk[k], &f);
+ }
+
+ SUFFIX(go_quad_div) (&qr, x, &dk[n]);
+
+ SUFFIX(go_quad_sub) (&qrp, &qrp, &qr);
+ if (SUFFIX(fabs)(qrp.h) <= SUFFIX(ldexp) (SUFFIX(fabs)(qr.h), -100)) {
+ converged = TRUE;
+ break;
+ }
+
+ qrp = qr;
+ gp = g;
+ for (k = 0; k <= n; k++) dpk[k] = dk[k];
+ }
+
+ if (!converged)
+ g_warning ("go_quad_atan_internal(%.20g) failed to converge\n",
+ (double)SUFFIX(go_quad_value) (x));
+
+ *res = qr;
+}
+
+static gboolean
+SUFFIX(go_quad_atan2_special) (const QUAD *y, const QUAD *x, DOUBLE *f)
+{
+ DOUBLE dy = SUFFIX(go_quad_value) (y);
+ DOUBLE dx = SUFFIX(go_quad_value) (x);
+
+ if (dy == 0) {
+ /* This assumes all zeros are +0. */
+ *f = (dx >= 0 ? 0 : +1);
+ return TRUE;
+ }
+
+ if (dx == 0) {
+ *f = (dy >= 0 ? 0.5 : -0.5);
+ return TRUE;
+ }
+
+ /* We could do quarters too */
+
+ return FALSE;
+}
+
+/**
+ * go_quad_atan2: (skip)
+ **/
+/**
+ * go_quad_atan2l: (skip)
+ **/
+void
+SUFFIX(go_quad_atan2) (QUAD *res, const QUAD *y, const QUAD *x)
+{
+ DOUBLE f;
+ DOUBLE dy = SUFFIX(go_quad_value) (y);
+ DOUBLE dx = SUFFIX(go_quad_value) (x);
+ QUAD qr;
+
+ if (SUFFIX(go_quad_atan2_special) (y, x, &f)) {
+ res->h = f * SUFFIX(go_quad_pi).h;
+ res->l = f * SUFFIX(go_quad_pi).l;
+ return;
+ }
+
+ if (SUFFIX(fabs) (dx) >= SUFFIX(fabs) (dy)) {
+ SUFFIX(go_quad_div) (&qr, y, x);
+ SUFFIX(go_quad_atan_internal) (res, &qr);
+ } else {
+ DOUBLE f;
+ QUAD qa;
+
+ SUFFIX(go_quad_div) (&qr, x, y);
+ SUFFIX(go_quad_atan_internal) (res, &qr);
+
+ f = (qr.h >= 0) ? 0.5 : -0.5;
+ qa.h = f * SUFFIX(go_quad_pi).h;
+ qa.l = f * SUFFIX(go_quad_pi).l;
+ SUFFIX(go_quad_sub) (res, &qa, res);
+ }
+
+ if (dx < 0) {
+ /* Correct for quadrants 2 and 3. */
+ if (dy > 0)
+ SUFFIX(go_quad_add) (res, res, &SUFFIX(go_quad_pi));
+ else
+ SUFFIX(go_quad_sub) (res, res, &SUFFIX(go_quad_pi));
+ }
+}
+
+/**
+ * go_quad_atan2pi: (skip)
+ **/
+/**
+ * go_quad_atan2pil: (skip)
+ **/
+void
+SUFFIX(go_quad_atan2pi) (QUAD *res, const QUAD *y, const QUAD *x)
+{
+ DOUBLE f;
+
+ if (SUFFIX(go_quad_atan2_special) (y, x, &f)) {
+ SUFFIX(go_quad_init) (res, f);
+ return;
+ }
+
+ /* Fallback. */
+ SUFFIX(go_quad_atan2) (res, y, x);
+ SUFFIX(go_quad_div) (res, res, &SUFFIX(go_quad_pi));
+}
diff --git a/goffice/math/go-quad.h b/goffice/math/go-quad.h
index 76f40a0..912b76e 100644
--- a/goffice/math/go-quad.h
+++ b/goffice/math/go-quad.h
@@ -27,6 +27,8 @@ void go_quad_pow (GOQuad *res, double *exp2, const GOQuad *x, const GOQuad *y);
void go_quad_exp (GOQuad *res, double *exp2, const GOQuad *a);
void go_quad_expm1 (GOQuad *res, const GOQuad *a);
void go_quad_log (GOQuad *res, const GOQuad *a);
+void go_quad_atan2 (GOQuad *res, const GOQuad *y, const GOQuad *x);
+void go_quad_atan2pi (GOQuad *res, const GOQuad *y, const GOQuad *x);
void go_quad_mul12 (GOQuad *res, double x, double y);
@@ -66,6 +68,8 @@ void go_quad_powl (GOQuadl *res, long double *exp2, const GOQuadl *x, const GOQu
void go_quad_expl (GOQuadl *res, long double *exp2, const GOQuadl *a);
void go_quad_expm1l (GOQuadl *res, const GOQuadl *a);
void go_quad_logl (GOQuadl *res, const GOQuadl *a);
+void go_quad_atan2l (GOQuadl *res, const GOQuadl *y, const GOQuadl *x);
+void go_quad_atan2pil (GOQuadl *res, const GOQuadl *y, const GOQuadl *x);
void go_quad_mul12l (GOQuadl *res, long double x, long double y);
diff --git a/tests/test-quad.c b/tests/test-quad.c
index 7ae5991..943770d 100644
--- a/tests/test-quad.c
+++ b/tests/test-quad.c
@@ -85,6 +85,47 @@ pow_tests (void)
/* ------------------------------------------------------------------------- */
+#define ATAN2PI(y_,x_) (atan2((y_),(x_))/M_PI)
+
+static void
+atan2_tests (void)
+{
+#define TEST1(a_,b_) do { \
+ BINTEST1(a_,b_,go_quad_atan2,atan2,"atan2"); \
+ BINTEST1(a_,b_,go_quad_atan2pi,ATAN2PI,"atan2pi"); \
+} while (0)
+ TEST1 (0, +2);
+ TEST1 (0, -2);
+ TEST1 (3, 0);
+ TEST1 (-3, 0);
+ TEST1 (0, 0);
+ TEST1 (+1, +1);
+ TEST1 (+2.3, +1.2);
+ TEST1 (+2.3, -1.2);
+ TEST1 (+2.3, +0.2);
+ TEST1 (+2.3, -0.2);
+ TEST1 (+0.2, +2.3);
+ TEST1 (+0.2, -2.3);
+ TEST1 (+2.3, +3);
+ TEST1 (+2.3, -3);
+ TEST1 (-2.3, +3);
+ TEST1 (-2.3, -3);
+ TEST1 (+2.3, +4);
+ TEST1 (+2.3, -4);
+ TEST1 (-2.3, +4);
+ TEST1 (-2.3, -4);
+ TEST1 (2, 100);
+ TEST1 (2, -100);
+ TEST1 (1.0/256, 1.0/1024);
+ TEST1 (exp(1), 1);
+ TEST1 (exp(1), -1);
+ TEST1 (exp(1), log(2));
+ TEST1 (exp(1), -log(2));
+}
+#undef TEST1
+
+/* ------------------------------------------------------------------------- */
+
#define TEST1(a_,b_) \
do { \
BINTEST1(a_,b_,go_quad_add,ADD,"add"); \
@@ -202,6 +243,27 @@ floor_tests (void)
/* ------------------------------------------------------------------------- */
+static void
+const_tests (void)
+{
+ GOQuad a;
+
+ g_assert (fabs (go_quad_value (&go_quad_pi) - M_PI) < 1e-14);
+ g_assert (fabs (go_quad_value (&go_quad_2pi) - 2 * M_PI) < 1e-14);
+ g_assert (fabs (go_quad_value (&go_quad_e) - exp(1)) < 1e-14);
+ g_assert (fabs (go_quad_value (&go_quad_ln2) - log(2)) < 1e-14);
+ g_assert (go_quad_value (&go_quad_zero) == 0);
+ g_assert (go_quad_value (&go_quad_one) == 1);
+
+ go_quad_mul (&a, &go_quad_sqrt2, &go_quad_sqrt2);
+ go_quad_sub (&a, &a, &go_quad_one);
+ go_quad_sub (&a, &a, &go_quad_one);
+ g_assert (go_quad_value (&a) < ldexp (1.0, -100));
+}
+
+
+/* ------------------------------------------------------------------------- */
+
int
main (int argc, char **argv)
{
@@ -223,21 +285,11 @@ main (int argc, char **argv)
go_quad_sub (&c, &c, &a);
g_assert (go_quad_value (&c) == ldexp (1.0, -80));
- g_assert (fabs (go_quad_value (&go_quad_pi) - M_PI) < 1e-14);
- g_assert (fabs (go_quad_value (&go_quad_2pi) - 2 * M_PI) < 1e-14);
- g_assert (fabs (go_quad_value (&go_quad_e) - exp(1)) < 1e-14);
- g_assert (fabs (go_quad_value (&go_quad_ln2) - log(2)) < 1e-14);
- g_assert (go_quad_value (&go_quad_zero) == 0);
- g_assert (go_quad_value (&go_quad_one) == 1);
-
- go_quad_mul (&a, &go_quad_sqrt2, &go_quad_sqrt2);
- go_quad_sub (&a, &a, &go_quad_one);
- go_quad_sub (&a, &a, &go_quad_one);
- g_assert (go_quad_value (&a) < ldexp (1.0, -100));
-
+ const_tests ();
basic4_tests ();
pow_tests ();
floor_tests ();
+ atan2_tests ();
go_quad_end (state);
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
