floating-pointapproximation

floating-point implementation of a function with removable singularity

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).

Solution

  • 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).

    Soap box

    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.