`cosc`: inaccuracy on the transition between two different implementations

[nsajko]

cosc(x::Float64) is implemented using a Taylor polynomial approximation around x = 0, and a generic formula relying on sinpi(x::Float64) and cospi(x::Float64) for the general case.

The polynomial approximation around the zero at x = 0 is necessary because the generic formula for cosc has bad accuracy around the zeros of cosc. However, the transition between the two implementation pieces is awkward, see the error plot (after downsampling by taking the max error) of the approximation error in ULPs:

The transition between the two implementation happens at 0.14, so the error is a bit high on both sides of the transition. It might make sense to lower it to like 5 ULPs or less, instead of the current 30 ULPs.

Two possible improvements:

• Using a Taylor polynomial is suboptimal, a minimax polynomial (odd, like the Taylor series) should be more efficient.

• If necessary, it might make sense to introduce an additional minimax polynomial to transition between the first approximation polynomial and the generic formula-based implementation.

[nsajko]

a bit high on both sides of the transition

Here's a closeup:

The transition between the two implementations, in 0.14, is the first of the two big jumps on the plot. So, on the side closer to x = 0, the error is above 15 ULP, on the other side it's above 30 ULP, both of which seem too high, compared to the usual accuracy of the generic formula, which is almost always less than 3.5 ULP when not close to a zero.

[mikmoore]

I got nerd-sniped into seeing if I could find some "improved" arrangement of the calculation. I suspect I haven't, but I'll put what I did here anyway.

Define the following two functions:

$$ c_2(x) = \frac{\cos(x) - 1}{x^2} \qquad s_2(x) = \frac{\sin(x) - x}{x^2} $$

Note that these have the same Taylor series as $\cos$ and $\sin$ but with the first term of each removed (and the remainder shifted down two orders into their places). Computing these functions more accurately may be one path to improving cosc:

$$ \text{cosc}(x) = \frac{\pi x \cos(\pi x) - \sin(\pi x)}{\pi x^2} = \frac{\pi x (\cos(\pi x) - 1) - (\sin(\pi x) - \pi x)}{\pi x^2} = \pi \left( \pi x c_2(\pi x) - s_2(\pi x) \right) $$

Although I suspect the main advantage is that these functions (carefully implemented) can avoid branching to separately handle small values. That said, computing the difference of these functions is probably marginally less accurate than the Horner-method single polynomial of cosc at small values. And on further thought, for non-small values the subtraction is probably real killer of accuracy anyway and this just replaces one subtraction with another and I'm not certain it's any better.

---

Now a probably-not-useful tangent:

The above observation motivates a function $e_2(x) = \frac{\exp(x) - 1 - x}{x^2}$ (the Taylor series of $\exp$ with the first two terms removed and the remainder shifted down) which satisfies

$$ -e_2(ix) = c_2(x) + i s_2(x) $$

$e_2$ can be calculated with a Taylor-series-plus-doubling algorithm similar to exp (although the doubling part is uglier). Alternatively, an easy (but inefficient) calculation is via a matrix exponential:

$$ e_2(x) = \left[ \exp\left( \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & x \end{bmatrix} \right) \right]_{1,3} $$

(i.e., the (1,3) entry of that matrix exponential).

But the range reduction one can normally do with sin/cos doesn't seem applicable here, so I expect this behaves poorly with large complex arguments. Which is to say I don't think $e_2$ is a likely path to a better cosc implementation. Although it could be used to extend cosc (and a similar function for sinc) to matrices (including singular matrices). Not that I've ever seen anyone ask for that...