Consider the following function:
This is an entire function with a Taylor series at 0 of: 1/1!, -1/3!, 1/5!, -1/7!, ... that I want to implement in floating point (think IEEE-754 single, IEEE-754 double or arbitrary precision like with MPFR.
My current C-ish implementation with a floating-point type F
looks something like:
F func (F x)
{
F q = sqrt (fabs (x));
if (x * x == 0)
return 1 - x / 6;
else
return x > 0
? sin (q) / q
: sinh (q) / q;
}
From what I can tell the implementation works reasonably well in an environment of 0, but I am wondering if there are loophole values where it does not work (well).
Loophole values where it does not work
Infinities: When x
is +infinity, the expected math result is 0.0, yet OP's code may return Not-a-number due to sin(infinity)
.
When x
is -infinity, the expected math result is infinity, yet OP's code may return Not-a-number due to sinh(infinity)/infinity
.
Code could test for these non-finite x
values.
Overflow: When x
is modestly negative, sinh(q)
may overflow into infinity, yet the extend math quotient sinh(q)/q
is still finite.
// float example
float yn = sinhf(q) / q;
// If a non-finite value resulted, try a slower, less precise, yet wider range approach.
if (!isfinite(yn)) {
yn = expf(q - logf(q) - logf(2));
}
@njuffa offers a good alternative.
Another improvement:
float yn = sinhf(q) / q;
if (!isfinite(yn)) {
// I expect this to work for all `x` where the math `sinhf(q) / q <= FLT_MAX`.
yn = expf(q/2);
yn *= yn/(2*q);
}
Modest magnitude x
: I found that, using float
, q = sqrtf(fabsf(x)); return sinhf(q) / q;
loses about 6 bits of correctness before sinhf(q)
hits overflow. This is because the error is forming the square root q
and then q
applied to sinhf(g)
incurs an error proportional to q
.
To really reduce this function's error, a better, TBD, approach is needed as x
grows negative (about x < -4
).
Trig errors
q = sqrtf(x)
has up to 0.5 ULP error in the result. When q < 99.99%*π
, sinf(q)
is well behaved. Yet at/near every following multiple of π, the ULP error of sinf(q)
can be terrible, even if the absolution value is very reasonable. As q
grows and its least significant significant digit nears 1.0, sign and magnitude of sinf(q)
becomes meaningless.
Additional code is needed to compensate for square root errors. Considerations start at x > 1.57
for ULP errors.
There are no for absolute errors concerns for finite x
.
OP likely area of concern.
By using if (x * x == 0) return 1 - x / 6;
code is well behaved near zero. The end-points of the region is about x = [-sqrt(DBL_TRUE_MIN) ... sqrt(DBL_TRUE_MIN)]
. Now let us look at the end points as x
is just within this region and just outside (e.g. near +/-sqrt(DBL_TRUE_MIN)
)
Taylor's series for sin(x)
is x - x^3/3! + x^5/5! +/- ...
Taylor's series for sinh(x)
is x + x^3/3! + x^5/5! + ...
Mathematically, all looks fine as the function is continuous, the 1st,2nd derivative also.
We could look close and compute with double x
values near this area, chart that and report loophole values where the y
steps are a tad jagged. Yet doing this examination is dubious as there is scant benefit over the simpler:
F func2 (F x)
{
if (x == 0)
return 1;
else
F q = sqrt (fabs (x));
return x > 0
? sin (q) / q
: sinh (q) / q;
}
Here, the nature of x
near 0.0 for sin(x), sinh(x)
is well researched and tested.
Still, to test, consider using float
as it is feasible to test all float
.
Find first difference
#include <float.h>
#include <math.h>
#include <stdio.h>
int main() {
float x = sqrtf(FLT_TRUE_MIN);
printf("sqrtf(FLT_TRUE_MIN):%24a\n", sqrtf(FLT_TRUE_MIN));
x = FLT_TRUE_MIN;
for (int i = 0; 1; i++) {
float q = sqrtf(x);
float yp = sinf(q) / q;
float yn = sinhf(q) / q;
float y0 = 1.0f - x / 6.0f;
if (y0 != yp || y0 != yn) {
printf("x:%a, q:%a, sin(q)/q:%a, sinh(q)/q:%a, 1-x/6:%a\n", //
x, q, yp, yn, y0);
break;
}
x = nextafterf(x, 2 * x);
}
return 0;
}
Sample output:
sqrtf(FLT_TRUE_MIN): 0x1.6a09e6p-75
x: 0x1p-46, q:0x1p-23, sin(q)/q:0x1p+0, sinh(q)/q:0x1.fffffep-1, 1-x/6:0x1p+0
It is clear that OP's original code has no discernable issues in the +/-sqrtf(FLT_TRUE_MIN)
range as all 3 sub-function result in the same y
values for x
until much larger than sqrtf(FLT_TRUE_MIN)
- as expected.
Still, I do not see a benefit (unless most usage values are near 0.0) in doing the if (x*x == 0)
test versus the simple if (x==0)
.
Certainly OP was most interesting is dealing with x
at or near 0.0 and the "loopholes" there about x*x == 0
. So was I when I first saw this post. Later I realized that many x
values in the range [-8828.3 ... -7995.2] failed.
It highlights the importance of functional testing across the spectrum of x
before attempting optimization in speed or low error performance.