[gnumeric] Math: New REDUCEPI function.
- From: Morten Welinder <mortenw src gnome org>
- To: commits-list gnome org
- Cc:
- Subject: [gnumeric] Math: New REDUCEPI function.
- Date: Wed, 5 Apr 2017 12:29:16 +0000 (UTC)
commit 4f963e1a214619de38026aa99c35737bc335d10b
Author: Morten Welinder <terra gnome org>
Date: Wed Apr 5 08:28:43 2017 -0400
Math: New REDUCEPI function.
ChangeLog | 5 ++
NEWS | 1 +
plugins/fn-math/ChangeLog | 4 ++
plugins/fn-math/functions.c | 33 ++++++++++++
plugins/fn-math/plugin.xml.in | 1 +
src/sf-trig.c | 116 +++++++++++++++++++++--------------------
src/sf-trig.h | 2 +
7 files changed, 105 insertions(+), 57 deletions(-)
---
diff --git a/ChangeLog b/ChangeLog
index 5b2a367..6d82e92 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,8 @@
+2017-04-05 Morten Welinder <terra gnome org>
+
+ * src/sf-trig.c (gnm_reduce_pi): Rename from reduce_pi_half.
+ Generalize to Pi/2^e and make public.
+
2017-03-28 Morten Welinder <terra gnome org>
* src/sf-bessel.c (y_via_j_series): Use gnm_yn if we can.
diff --git a/NEWS b/NEWS
index 9ca860f..b97231c 100644
--- a/NEWS
+++ b/NEWS
@@ -3,6 +3,7 @@ Gnumeric 1.12.35
Morten:
* Test suite improvements.
* Minor BESSELY improvements.
+ * New function REDUCEPI.
--------------------------------------------------------------------------
Gnumeric 1.12.34
diff --git a/plugins/fn-math/ChangeLog b/plugins/fn-math/ChangeLog
index 9d6c62c..389fe06 100644
--- a/plugins/fn-math/ChangeLog
+++ b/plugins/fn-math/ChangeLog
@@ -1,3 +1,7 @@
+2017-04-05 Morten Welinder <terra gnome org>
+
+ * functions.c (gnumeric_reducepi): New function.
+
2017-03-20 Morten Welinder <terra gnome org>
* Release 1.12.34
diff --git a/plugins/fn-math/functions.c b/plugins/fn-math/functions.c
index 7800812..56581af 100644
--- a/plugins/fn-math/functions.c
+++ b/plugins/fn-math/functions.c
@@ -1522,6 +1522,36 @@ gnumeric_radians (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
/***************************************************************************/
+static GnmFuncHelp const help_reducepi[] = {
+ { GNM_FUNC_HELP_NAME, F_("REDUCEPI:reduce modulo Pi divided by a power of 2")},
+ { GNM_FUNC_HELP_ARG, F_("x:number")},
+ { GNM_FUNC_HELP_ARG, F_("e:scale")},
+ { GNM_FUNC_HELP_ARG, F_("q:get lower bits of quotient, defaults to FALSE")},
+ { GNM_FUNC_HELP_EXAMPLES, "=REDUCEPI(10,1)" },
+ { GNM_FUNC_HELP_NOTE, F_("This function returns a value, xr, such that @{x}=xr+j*Pi/2^@{e} where j is
an integer and the absolute value of xr does not exceed Pi/2^(@{e}+1). If optional argument @{q} is TRUE,
returns instead the @e+1 lower bits of j. The reduction is performed as-if using an exact value of Pi.")},
+ { GNM_FUNC_HELP_NOTE, F_("The lowest valid @{e} is -1 representing reduction modulo 2*Pi; the highest
is 7 representing reduction modulo Pi/256.")},
+ { GNM_FUNC_HELP_SEEALSO, "PI"},
+ { GNM_FUNC_HELP_END}
+};
+
+static GnmValue *
+gnumeric_reducepi (GnmFuncEvalInfo *ei, GnmValue const * const *argv)
+{
+ gnm_float x = value_get_as_float (argv[0]);
+ int e = value_get_as_int (argv[1]);
+ gboolean q = argv[2] ? value_get_as_checked_bool (argv[2]) : FALSE;
+ int j;
+ gnm_float xr;
+
+ if (e < -1 || e > 7)
+ return value_new_error_VALUE (ei->pos);
+
+ xr = gnm_reduce_pi (x, (int)e, &j);
+ return q ? value_new_int (j) : value_new_float (xr);
+}
+
+/***************************************************************************/
+
static GnmFuncHelp const help_sin[] = {
{ GNM_FUNC_HELP_NAME, F_("SIN:the sine of @{x}")},
{ GNM_FUNC_HELP_ARG, F_("x:angle in radians")},
@@ -3573,6 +3603,9 @@ GnmFuncDescriptor const math_functions[] = {
{ "radians", "f", help_radians,
gnumeric_radians, NULL, NULL, NULL,
GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_EXHAUSTIVE },
+ { "reducepi", "ff|f", help_reducepi,
+ gnumeric_reducepi, NULL, NULL, NULL,
+ GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_UNIQUE_TO_GNUMERIC, GNM_FUNC_TEST_STATUS_NO_TESTSUITE },
{ "roman", "f|f", help_roman,
gnumeric_roman, NULL, NULL, NULL,
GNM_FUNC_SIMPLE, GNM_FUNC_IMPL_STATUS_COMPLETE, GNM_FUNC_TEST_STATUS_BASIC },
diff --git a/plugins/fn-math/plugin.xml.in b/plugins/fn-math/plugin.xml.in
index 9dccaef..3778e63 100644
--- a/plugins/fn-math/plugin.xml.in
+++ b/plugins/fn-math/plugin.xml.in
@@ -83,6 +83,7 @@
<function name="power"/>
<function name="quotient"/>
<function name="radians"/>
+ <function name="reducepi"/>
<function name="roman"/>
<function name="round"/>
<function name="rounddown"/>
diff --git a/src/sf-trig.c b/src/sf-trig.c
index 8660506..e1935da 100644
--- a/src/sf-trig.c
+++ b/src/sf-trig.c
@@ -79,12 +79,10 @@ gnm_acoth (gnm_float x)
/* ------------------------------------------------------------------------- */
-#ifdef GNM_REDUCES_TRIG_RANGE
-
static gnm_float
-reduce_pi_half_simple (gnm_float x, int *km4)
+reduce_pi_simple (gnm_float x, int *pk, int kbits)
{
- static const gnm_float two_over_pi = 0.63661977236758134307553505349;
+ static const gnm_float two_over_pi = GNM_const(0.63661977236758134307553505349);
static const gnm_float pi_over_two[] = {
+0x1.921fb5p+0,
+0x1.110b46p-26,
@@ -93,18 +91,18 @@ reduce_pi_half_simple (gnm_float x, int *km4)
+0x1.c1cd12p-107
};
int i;
- gnm_float k = gnm_floor (x * two_over_pi + 0.5);
+ gnm_float k = gnm_floor (x * gnm_ldexp (two_over_pi, kbits - 2) + 0.5);
gnm_float xx = 0;
g_assert (k < (1 << 26));
- *km4 = (int)k & 3;
+ *pk = (int)k;
if (k == 0)
return x;
- x -= k * pi_over_two[0];
+ x -= k * gnm_ldexp (pi_over_two[0], 2 - kbits);
for (i = 1; i < 5; i++) {
- gnm_float dx = k * pi_over_two[i];
+ gnm_float dx = k * gnm_ldexp (pi_over_two[i], 2 - kbits);
gnm_float s = x - dx;
xx += x - s - dx;
x = s;
@@ -140,10 +138,10 @@ add_at (guint32 *dst, guint64 p)
}
static gnm_float
-reduce_pi_half_full (gnm_float x, int *km4)
+reduce_pi_full (gnm_float x, int *pk, int kbits)
{
/* FIXME? This table isn't big enough for long double's range */
- static guint32 one_over_two_pi[] = {
+ static const guint32 one_over_two_pi[] = {
0x28be60dbu,
0x9391054au,
0x7f09d5f4u,
@@ -181,7 +179,7 @@ reduce_pi_half_full (gnm_float x, int *km4)
0x9afed7ecu,
0x47e35742u
};
- static guint32 pi_over_four[] = {
+ static const guint32 pi_over_four[] = {
0xc90fdaa2u,
0x2168c234u,
0xc4c6628bu,
@@ -229,21 +227,21 @@ reduce_pi_half_full (gnm_float x, int *km4)
add_at (&r[i + 1], (guint64)w2 * one_over_two_pi[i + 2 + di]);
/*
- * We're done at i==3 unless the first 29 bits, not counting
- * those ending up in sign and *km4, are all zeros or ones.
+ * We're done at i==3 unless the first 31-kbits bits, not counting
+ * those ending up in sign and *pk, are all zeros or ones.
*/
- if (i == 3 && ((r[1] + 1u) & 0x1fffffffu) > 1)
+ if (i == 3 && ((r[1] + 1u) & (0x7fffffffu >> kbits)) > 1)
break;
}
- *km4 = r[1] >> 30;
- if ((neg = ((r[1] >> 29) & 1))) {
- *km4 = (*km4 + 1) & 3;
+ *pk = kbits ? (r[1] >> (32 - kbits)) : 0;
+ if ((neg = ((r[1] >> (31 - kbits)) & 1))) {
+ (*pk)++;
/* Two-complement negation */
for (j = 1; j <= i; j++) r[j] ^= 0xffffffffu;
add_at (&r[i], 1);
}
- r[1] &= 0x3fffffffu;
+ r[1] &= (0xffffffffu >> kbits);
j = 1;
if (r[j] == 0) j++;
@@ -256,9 +254,9 @@ reduce_pi_half_full (gnm_float x, int *km4)
add_at (&r4[3], (guint64)r[j + 2] * pi_over_four[0]);
h48 = gnm_ldexp (((guint64)r4[0] << 16) | (r4[1] >> 16),
- -32 * j + 3 - 16);
+ -32 * j + (kbits + 1) - 16);
l48 = gnm_ldexp (((guint64)(r4[1] & 0xffff) << 32) | r4[2],
- -32 * j + 3 - 64);
+ -32 * j + (kbits + 1) - 64);
rh = h48 + l48;
rl = h48 - rh + l48;
@@ -268,59 +266,67 @@ reduce_pi_half_full (gnm_float x, int *km4)
rl = -rl;
}
- return rh + rl;
+ return gnm_ldexp (rh + rl, 2 - kbits);
}
-/*
- * Subtract k times Pi from @x where k is picked such that the absolute
- * value of x-k*Pi is no larger than Pi/4. k%4 is returned in @km4.
- * @x must not be negative.
- *
- * This function is used to reduce arguments to trigonometric functions
- * to [-Pi/4;+Pi/4], a range in which hardware and libm generally are
- * accurate. This avoids problems for example with Intel chips for
- * values near multiples of Pi/2 or values above 10^15.
- *
- * Note, that the Pi used has sufficient accuracy to not suffer cancellation
- * errors throughout the entire range of "double". That means 1000+ bits.
- *
- * The "no larger than Pi/4" condition may be violated by a hair. That's
- * not important since none of the trigonometric functions are critical
- * there.
- */
-static gnm_float
-reduce_pi_half (gnm_float x, int *km4)
+/**
+ * gnm_reduce_pi:
+ * @x: number of reduce
+ * @e: scale between -1 and 8, inclusive.
+ * @k: (out): location to return lower @e+1 bits of reduction count
+ *
+ * This function performs range reduction for trigonometric functions.
+ *
+ * Returns: a value, xr, such that x = xr + j * Pi/2^@e for some integer
+ * number j and |xr| <= Pi/2^(@e+1). The lower @e+1 bits of j will be
+ * returned in @k.
+ */
+gnm_float
+gnm_reduce_pi (gnm_float x, int e, int *k)
{
- void *state;
gnm_float xr;
+ void *state;
- if (gnm_isnan (x))
- return x;
+ g_return_val_if_fail (e >= -1 && e <= 8, x);
+ g_return_val_if_fail (k != NULL, x);
- g_assert (x >= 0);
+ if (!gnm_finite (x)) {
+ *k = 0;
+ return x * gnm_nan;
+ }
/*
* We aren't actually using quads, but we rely somewhat on
* proper ieee double semantics.
*/
state = gnm_quad_start ();
- if (x < (1u << 26))
- xr = reduce_pi_half_simple (x, km4);
+
+ if (gnm_abs (x) < (1u << (27 - e)))
+ xr = reduce_pi_simple (gnm_abs (x), k, e + 1);
else
- xr = reduce_pi_half_full (x, km4);
+ xr = reduce_pi_full (gnm_abs (x), k, e + 1);
+
+ if (x < 0) {
+ xr = -xr;
+ *k = -*k;
+ }
+ *k &= ((1 << (e + 1)) - 1);
+
gnm_quad_end (state);
return xr;
}
+
+
+
+#ifdef GNM_REDUCES_TRIG_RANGE
+
gnm_float
gnm_sin (gnm_float x)
{
int km4;
- gnm_float ax = gnm_abs (x);
- gnm_float xr = reduce_pi_half (ax, &km4);
-
- if (x < 0) km4 ^= 2;
+ gnm_float xr = gnm_reduce_pi (x, 1, &km4);
switch (km4) {
default:
@@ -335,8 +341,7 @@ gnm_float
gnm_cos (gnm_float x)
{
int km4;
- gnm_float ax = gnm_abs (x);
- gnm_float xr = reduce_pi_half (ax, &km4);
+ gnm_float xr = gnm_reduce_pi (x, 1, &km4);
switch (km4) {
default:
@@ -351,10 +356,7 @@ gnm_float
gnm_tan (gnm_float x)
{
int km4;
- gnm_float ax = gnm_abs (x);
- gnm_float xr = reduce_pi_half (ax, &km4);
-
- if (x < 0) xr = -xr;
+ gnm_float xr = gnm_reduce_pi (x, 1, &km4);
switch (km4) {
default:
diff --git a/src/sf-trig.h b/src/sf-trig.h
index fb1c72c..f2d4cb6 100644
--- a/src/sf-trig.h
+++ b/src/sf-trig.h
@@ -8,6 +8,8 @@ gnm_float gnm_acot (gnm_float x);
gnm_float gnm_coth (gnm_float x);
gnm_float gnm_acoth (gnm_float x);
+gnm_float gnm_reduce_pi (gnm_float x, int e, int *k);
+
#ifdef GNM_REDUCES_TRIG_RANGE
/* gnm_sin, gnm_cos, gnm_tan prototyped in numbers.h */
#endif
[
Date Prev][
Date Next] [
Thread Prev][
Thread Next]
[
Thread Index]
[
Date Index]
[
Author Index]
