3

When I try to calculate the function $f(x)=1-\operatorname{sinc}x$ for small values of $x$ I get large relative errors due to catastrophic cancellation. I want an accurate way to calculate $f(x)$ without using a series expansion or an iterative method. (For my purposes, accurate means always within about 5 ULPs).

I can easily calculate $\operatorname{sinc}x$. The problem is calculating $1-\operatorname{sinc}x$. I am looking for a way that does not simply compute a truncated series like this:

$$ \begin{aligned} f(x) &= 1-\operatorname{sinc}x \\ &= 1-\frac{\sin x}{x} \\ &= 1-\left(1-\frac{x^2}{3!}+\frac{x^4}{5!}-...\right) \\ &= \frac{x^2}{3!}-\frac{x^4}{5!}+... \\ \end{aligned} $$

Notice the problematic $1-1$ in the second-last line above.

I tried to reformulate $f(x)$ in terms of standard functions by applying trigonometric identities. After finding the common denominator $x$, there are no identities that bridge the $x$ and $\sin x$ terms in the numerator.

I also tried an approach similar to Kahan's $\operatorname{log1p}()$ and $\operatorname{expm1}()$ which rely on matching the loss rates between numerator and denominator, but I still get large errors. Maybe that only works for functions involving $\log$ and $\exp$.

Ken Whatmough
  • 243
  • 1
    What if you define $ f(x) = \frac{{x^2 }}{{3!}} - \frac{{x^4 }}{{5!}} + \cdots $ and use it for small $x$ without any preceeding steps involving $1-1$? You can also try a Padé approximation of this convergent series. – Gary Apr 22 '21 at 13:49
  • 2
    Presumably, the reason you want to avoid series is that you think it takes too long. Even if your language supports $\sin(x)$ it can take a long time compared to a multiply. I would just set a threshold and compute $\sinc(x)$ if $x$ greater than the threshold, then use the series when it is less. If the threshold is small, it won't take many terms of the series to get the accuracy you want. – Ross Millikan Apr 22 '21 at 13:52
  • @RossMillikan You are correct, yes, but first I want to know if there's a way that uses algebraic manipulation or something more clever like Kahan's trick. – Ken Whatmough Apr 22 '21 at 14:06
  • 1
    Instead of the direct series, you can use a Padé approximation, which often converges faster. It is the ratio of polynomials. The fact that sinc doesn't have poles makes me suspect that it won't help much. – Ross Millikan Apr 22 '21 at 14:09
  • 1
    @Gary Good points also, but I was hoping for a clever identity or perhaps an error "balancing trick" like in Kahan. – Ken Whatmough Apr 22 '21 at 14:10
  • 1
    Disappointing. It looks just like 1-cos(x) = x^2/2! - x^4/4! + x^6/6! - ... which does have a nice formula 2*(sin(x/2)^2). Is 1-cos(x) just lucky, and 1-sinc(x) is unlucky? – Don Hatch Jun 08 '22 at 07:17
  • @DonHatch Indeed. The $x$ denominator seems to prevent a closed-form solution and Kahan's error-matching tricks seem relatively ineffective here. Padé is nice for small inputs but I was hoping for something that holds for all inputs. BTW your page "The Right Way to Calculate Stuff" is the inspiration behind this. I've been trying to expand upon your original work. Thanks for commenting. – Ken Whatmough Jun 09 '22 at 13:19

1 Answers1

2

As other contributors, $[2n+2,2n]$ Padé approximants are probably the best approximation.

For example, the simplest $$f(x)= 1-\frac{\sin x}{x} \sim\frac{x^2(420 -11 x^2)}{60 \left(x^2+42\right)}$$ gives an error of $1.10\times 10^{-8}$ for $x=\frac \pi 6$ and $1.70\times 10^{-12}$ for $x=\frac \pi{24} $.

Claude Leibovici
  • 289,558