[goffice] Trig: improve go_sinpi, go_cospi, go_tanpi, go_cotpi.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [goffice] Trig: improve go_sinpi, go_cospi, go_tanpi, go_cotpi.
- Date: Wed, 15 Apr 2015 17:06:16 +0000 (UTC)
commit 3c54f44de131951d5d32c65ebe5148acee2e439f
Author: Morten Welinder <terra gnome org>
Date: Wed Apr 15 13:05:07 2015 -0400
Trig: improve go_sinpi, go_cospi, go_tanpi, go_cotpi.
Pay more attention to -0 and other special cases.
goffice/math/go-math.c | 108 ++++++++++++++++++++++++++++++++++++++------
tests/test-math.c | 118 ++++++++++++++++++++++++++++++-----------------
2 files changed, 169 insertions(+), 57 deletions(-)
---
diff --git a/goffice/math/go-math.c b/goffice/math/go-math.c
index 2999c1d..e2bad63 100644
--- a/goffice/math/go-math.c
+++ b/goffice/math/go-math.c
@@ -990,7 +990,16 @@ double
go_sinpi (double x)
{
int k;
+ double x0 = x;
x = reduce_half (x, &k);
+
+ /*
+ * Per IEEE 754 2008:
+ * sinpi(n) == 0 with sign of n.
+ */
+ if (x == 0 && (k & 1) == 0)
+ return copysign (0, x0);
+
return do_sinpi (x, k);
}
@@ -1006,6 +1015,14 @@ go_cospi (double x)
{
int k;
x = reduce_half (x, &k);
+
+ /*
+ * Per IEEE 754 2008:
+ * cospi(n+0.5) == +0 for any integer n.
+ */
+ if (x == 0 && (k & 1) == 1)
+ return +0.0;
+
return do_sinpi (x, k + 1);
}
@@ -1019,17 +1036,40 @@ go_cospi (double x)
double
go_tanpi (double x)
{
- int k;
- x = reduce_half (x, &k);
- return do_sinpi (x, k) / do_sinpi (x, k + 1);
+ /*
+ * IEEE 754 2008 doesn't have tanpi and thus doesn't define the
+ * behaviour for -0 argument or result. crlibm has tanpi, but
+ * doesn't seem to be fully clear on these cases.
+ */
+
+ /* inf -> nan; -n -> -0; +n -> +0 */
+ x = fmod (x, 1.0);
+
+ if (x == 0)
+ return copysign (0.0, x);
+ if (fabs (x) == 0.5)
+ return copysign (go_nan, x);
+ else
+ return go_sinpi (x) / go_cospi (x);
}
double
go_cotpi (double x)
{
- int k;
- x = reduce_half (x, &k);
- return do_sinpi (x, k + 1) / do_sinpi (x, k);
+ /*
+ * IEEE 754 2008 doesn't have tanpi. Neither does crlibm. Mirror
+ * tanpi here.
+ */
+
+ /* inf -> nan; -n -> -0; +n -> +0 */
+ x = fmod (x, 1.0);
+
+ if (x == 0)
+ return copysign (go_nan, x);
+ if (fabs (x) == 0.5)
+ return copysign (0.0, x);
+ else
+ return go_cospi (x) / go_sinpi (x);
}
double
@@ -1045,7 +1085,7 @@ go_atan2pi (double y, double x)
double
go_atanpi (double x)
{
- return x < 0 ? go_atan2pi (-x, -1) : go_atan2pi (x, 1);
+ return x < 0 ? -go_atan2pi (-x, 1) : go_atan2pi (x, 1);
}
#ifdef GOFFICE_WITH_LONG_DOUBLE
@@ -1103,7 +1143,16 @@ long double
go_sinpil (long double x)
{
int k;
+ long double x0 = x;
x = reduce_halfl (x, &k);
+
+ /*
+ * Per IEEE 754 2008:
+ * sinpi(n) == 0 with sign of n.
+ */
+ if (x == 0 && (k & 1) == 0)
+ return copysign (0, x0);
+
return do_sinpil (x, k);
}
@@ -1119,6 +1168,14 @@ go_cospil (long double x)
{
int k;
x = reduce_halfl (x, &k);
+
+ /*
+ * Per IEEE 754 2008:
+ * cospi(n+0.5) == +0 for any integer n.
+ */
+ if (x == 0 && (k & 1) == 1)
+ return +0.0l;
+
return do_sinpil (x, k + 1);
}
@@ -1132,17 +1189,40 @@ go_cospil (long double x)
long double
go_tanpil (long double x)
{
- int k;
- x = reduce_halfl (x, &k);
- return do_sinpil (x, k) / do_sinpil (x, k + 1);
+ /*
+ * IEEE 754 2008 doesn't have tanpi and thus doesn't define the
+ * behaviour for -0 argument or result. crlibm has tanpi, but
+ * doesn't seem to be fully clear on these cases.
+ */
+
+ /* inf -> nan; -n -> -0; +n -> +0 */
+ x = fmodl (x, 1.0l);
+
+ if (x == 0)
+ return copysignl (0.0l, x);
+ if (fabsl (x) == 0.5l)
+ return copysignl (go_nanl, x);
+ else
+ return go_sinpil (x) / go_cospil (x);
}
long double
go_cotpil (long double x)
{
- int k;
- x = reduce_halfl (x, &k);
- return do_sinpil (x, k + 1) / do_sinpil (x, k);
+ /*
+ * IEEE 754 2008 doesn't have tanpi. Neither does crlibm. Mirror
+ * tanpi here.
+ */
+
+ /* inf -> nan; -n -> -0; +n -> +0 */
+ x = fmodl (x, 1.0l);
+
+ if (x == 0)
+ return copysignl (go_nanl, x);
+ if (fabsl (x) == 0.5)
+ return copysignl (0.0l, x);
+ else
+ return go_cospil (x) / go_sinpil (x);
}
long double
@@ -1158,7 +1238,7 @@ go_atan2pil (long double y, long double x)
long double
go_atanpil (long double x)
{
- return x < 0 ? go_atan2pil (-x, -1) : go_atan2pil (x, 1);
+ return x < 0 ? -go_atan2pil (-x, 1) : go_atan2pil (x, 1);
}
#endif
diff --git a/tests/test-math.c b/tests/test-math.c
index fc737c4..c39da3b 100644
--- a/tests/test-math.c
+++ b/tests/test-math.c
@@ -1,54 +1,79 @@
#include <goffice/goffice.h>
-#define REFTEST(a_,f_,r_, txt_) \
-do { \
- double a = (a_); \
- double r = (r_); \
- double fa = f_(a); \
- g_printerr ("%s(%g) = %g [%g]\n", \
- txt_, a, fa, r); \
- if (r == floor (r)) \
- g_assert (r == fa); \
- else \
- g_assert (fabs (r - fa) / r < 1e-14); \
-} while (0)
-
-#define REFTEST2(a_,b_,f_,r_, txt_) \
-do { \
- double a = (a_); \
- double b = (b_); \
- double r = (r_); \
- double fab = f_(a, b); \
- g_printerr ("%s(%g,%g) = %g [%g]\n", \
- txt_, a, b, fab, r); \
- if (r * 256 == floor (r * 256)) \
- g_assert (r == fab); \
- else \
- g_assert (fabs (r - fab) / r < 1e-14); \
-} while (0)
+#define REFTEST(a_,f_,r_, txt_) \
+ do { \
+ double a = (a_); \
+ double r = (r_); \
+ double fa = f_(a); \
+ g_printerr ("%s(%g) = %g [%g]\n", \
+ txt_, a, fa, r); \
+ g_assert (copysign (1,r) == copysign (1,fa)); \
+ if (isnan (r)) \
+ g_assert (isnan (fa)); \
+ else if (r == floor (r)) \
+ g_assert (r == fa); \
+ else \
+ g_assert (fabs (r - fa) / fabs (r) < 1e-14); \
+ } while (0)
+
+#define REFTEST2(a_,b_,f_,r_, txt_) \
+ do { \
+ double a = (a_); \
+ double b = (b_); \
+ double r = (r_); \
+ double fab = f_(a, b); \
+ g_printerr ("%s(%g,%g) = %g [%g]\n", \
+ txt_, a, b, fab, r); \
+ g_assert (copysign (1,r) == copysign(1,fab)); \
+ if (r * 256 == floor (r * 256)) \
+ g_assert (r == fab); \
+ else \
+ g_assert (fabs (r - fab) / fabs (r) < 1e-14); \
+ } while (0)
#ifdef GOFFICE_WITH_LONG_DOUBLE
-#define REFTEST2l(a_,b_,f_,r_, txt_) \
-do { \
- long double a = (a_); \
- long double b = (b_); \
- long double r = (r_); \
- long double fab = f_(a, b); \
- g_printerr ("%s(%Lg,%Lg) = %Lg [%Lg]\n", \
- txt_, a, b, fab, r); \
- if (r * 256 == floorl (r * 256)) \
- g_assert (r == fab); \
- else \
- g_assert (fabsl (r - fab) / r < 1e-14); \
-} while (0)
+
+#define REFTESTl(a_,f_,r_, txt_) \
+ do { \
+ long double a = (a_); \
+ long double r = (r_); \
+ long double fa = f_(a); \
+ g_printerr ("%s(%Lg) = %Lg [%Lg]\n", \
+ txt_, a, fa, r); \
+ g_assert (copysignl (1,r) == copysignl (1,fa)); \
+ if (isnanl (r)) \
+ g_assert (isnanl (fa)); \
+ else if (r == floorl (r)) \
+ g_assert (r == fa); \
+ else \
+ g_assert (fabsl (r - fa) / fabsl (r) < 1e-14); \
+ } while (0)
+
+#define REFTEST2l(a_,b_,f_,r_, txt_) \
+ do { \
+ long double a = (a_); \
+ long double b = (b_); \
+ long double r = (r_); \
+ long double fab = f_(a, b); \
+ g_printerr ("%s(%Lg,%Lg) = %Lg [%Lg]\n", \
+ txt_, a, b, fab, r); \
+ g_assert (copysignl (1,r) == copysignl(1,fab)); \
+ if (r * 256 == floorl (r * 256)) \
+ g_assert (r == fab); \
+ else \
+ g_assert (fabsl (r - fab) / fabsl (r) < 1e-14); \
+ } while (0)
#else
+#define REFTESTl(a_,f_,r_, txt_)
#define REFTEST2l(a_,b_,f_,r_, txt_)
#endif
/* ------------------------------------------------------------------------- */
-static double fake_sinpi (double x) { return x == floor (x) ? 0 : sin (x * M_PI); }
-static double fake_cospi (double x) { return fake_sinpi (x + 0.5); }
+static double fake_sinpi (double x) { return x == floor (x) ? copysign(0,x) : sin (x * M_PI); }
+static double fake_cospi (double x) { return x + 0.5 == floor (x + 0.5) ? +0 : cos (x * M_PI); }
+static double fake_tanpi (double x) { x = fmod(x,1); return x + 0.5 == floor (x + 0.5) ? copysign (go_nan,x)
: tan (x * M_PI); }
+static double fake_cotpi (double x) { x = fmod(x,1); return x == 0 ? copysign (go_nan,x) : fake_cospi (x) /
fake_sinpi (x); }
static double fake_atanpi (double x) { return atan (x) / M_PI; }
static double fake_atan2pi (double y, double x) { return atan2 (y,x) / M_PI; }
@@ -56,13 +81,20 @@ static double fake_atan2pi (double y, double x) { return atan2 (y,x) / M_PI; }
#define TEST1(a_) do { \
REFTEST(a_,go_sinpi,fake_sinpi(a_),"sinpi"); \
REFTEST(a_,go_cospi,fake_cospi(a_),"cospi"); \
- REFTEST(a_,go_tanpi,(fake_sinpi(a_)/fake_cospi(a_)),"tanpi"); \
- REFTEST(a_,go_cotpi,(fake_cospi(a_)/fake_sinpi(a_)),"cotpi"); \
+ REFTEST(a_,go_tanpi,fake_tanpi(a_),"tanpi"); \
+ REFTEST(a_,go_cotpi,fake_cotpi(a_),"cotpi"); \
REFTEST(a_,go_atanpi,fake_atanpi(a_),"atanpi"); \
+ \
+ REFTEST(a_,go_sinpil,fake_sinpi(a_),"sinpil"); \
+ REFTEST(a_,go_cospil,fake_cospi(a_),"cospil"); \
+ REFTEST(a_,go_tanpil,fake_tanpi(a_),"tanpil"); \
+ REFTEST(a_,go_cotpil,fake_cotpi(a_),"cotpil"); \
+ REFTEST(a_,go_atanpil,fake_atanpi(a_),"atanpil"); \
} while (0)
#define TEST2(a_,b_) do { \
REFTEST2(a_,b_,go_atan2pi,fake_atan2pi(a_,b_),"atan2pi"); \
+ REFTEST2(a_,b_,go_atan2pil,fake_atan2pi(a_,b_),"atan2pil"); \
} while (0)
static void
@@ -70,7 +102,7 @@ trig_tests (void)
{
double d;
- for (d = 0; d < 10; d += 0.125) {
+ for (d = 0; d < 5; d += 0.125) {
TEST1(d);
TEST1(-d);
}
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
