improve cosc(::Float32) and cosc(::Float64) accuracy

[nsajko]

Accomplished by:

• Using minimax instead of Taylor polynomials. They're more efficient. The minimax polynomials are found using Sollya.

• Using multiple polynomials for different subintervals of the domain. Known as "domain splitting". Enables good accuracy for a wider region around the origin than possible with only a single polynomial.

The polynomial degrees are kept as-is.

Fixes #59065

[nsajko]

cosc(::Float64)

Before:

After:

cosc(::Float32)

Before:

After:

[nsajko]

To clarify: the above plots show the error in ULPs. In each case, the plot happens after downsampling by taking the maximum. Also see the linked issue.

[nsajko]

The failing test is just flaky CI, issue #59073

[nsajko]

Test failure seems surely unrelated.

[oscardssmith]

I've read the pr 3 time by now, and I'm still unsure whether the struct it introduces is worth it. I feel like I n better read the more function based approach in exp.jl, but some portion of that is that I'm the author of those.

[nsajko]

still unsure whether the struct it introduces is worth it

The latest commit ("improve style") improves the style a bit. I wonder if it addresses your concerns somewhat?

[stevengj]

To be clear, is this the minimax polynomial that minimizes the maximum absolute error over the interval, or which minimizes the maximum relative error? Presumably the latter, I hope?

I'm a bit concerned that your plot seems to show 1ulp error at x == 0. What does that mean? A nonzero error at the root would correspond to an infinite relative error. Since the x == 0 root is known exactly, and is exactly representable, we should strive to have small relative error as we approach it. (Ensuring this is a nice property of using a Taylor series around a zero.)

Or am I deceived and your plot doesn't actually have the x == 0 data point, so it is just showing a small relative error as you approach the root? Also, you might want to use a semilogx scale near x == 0 to make sure nothing bad is happening close to the zero.

(I'm also not entirely clear on what precisely you are plotting. I'm imagining a plot of relative error |approx–exact|/|exact|, divided by eps(T) to get something ulp-like?)

[stevengj]

Also, you might want to check that cosc(+0.0) === -0.0 and cosc(-0.0) === +0.0; I don't remember if we have a test for this.

[nsajko]

is this the minimax polynomial that minimizes the maximum absolute error over the interval, or which minimizes the maximum relative error?

It minimizes relative error. The fpminimax commands that appear in the included Sollya script allow choosing either absolute or relative error to optimize for, but defaults to the relative option.

Anyway if I had chosen to minimize the absolute error with the polynomial, the ULP error plots would have looked quite different.

I'm a bit concerned that your plot seems to show 1ulp error at x == 0. What does that mean? A nonzero error at the root would correspond to an infinite relative error.

At x == 0, the result is exact (easy to check exhaustively, see a relevant separate PR inspired by your comments: #59265). Otherwise the relative error would not be finite.

Each point on the plot is the result of many evaluations in the same neighborhood, which are aggregated basically by taking the maximum value. Trying to plot every evaluation would result in an unreadably noisy plot, so the idea is basically to get rid of the downward spikes of the noise. Thus the plot visualizes the maximum error in the neighborhood of each point, but not the minimum error. The justification is that the worst case is exactly what is interesting.

As an additional explanation, I'll copy over a paragraph from my message on this PR:

• fix accuracy of logcosh(::Union{Float16, Float32, Float64}) JuliaStats/LogExpFunctions.jl#101:

The plotting happens after downsampling (which just takes the maximum value), which itself happens after smoothing (which also just takes the maximum, i.e., a sliding window maximum). Thus the downsampling and smoothing basically just remove the downward spikes of the noise, preserving the maximum error at the region around each point.

If you want I can also share the code I use for plotting. The ULP error itself is calculated using the same code included in this PR, see module ULPError. The downsampling is implemented using some iterators, at least some of which I want to eventually include into IterTools.jl. The initial PR (implementing the sliding window maximum, used for smoothing, which happens before downsampling proper) is:

• add sliding_window_maxima JuliaCollections/IterTools.jl#126

Since the x == 0 root is known exactly, and is exactly representable, we should strive to have small relative error as we approach it.

IMO it is fine to minimize the maximum relative error over the entire interval, there's no need to treat the immediate neighborhood of the root specially. (In fact, Sollya is not able to do that in the first place. My own package FindMinimaxPolynomial is able to do that, however it lacks Sollya's functionality for optimizing the coefficients of the polynomial in machine precision.)

Ensuring this is a nice property of using a Taylor series around a zero.

An exhaustive comparison of the maximum ULP errors between zero and nextfloat(type(0), 100_000_000) before and after this PR shows the maximum error happens to stays the same in this immediate neighborhood of zero. They are 0.58824617f0 for Float32 and 0.49999997f0 for Float64.

I'm also not entirely clear on what precisely you are plotting. I'm imagining a plot of relative error |approx–exact|/|exact|, divided by eps(T) to get something ulp-like?

See the module ULPError in test/testhelpers/ULPError.jl. In particular, the generic formula is:

abs((approximate - acc) / eps(approximate))

I believe this corresponds exactly to one of the definitions of ULP error, so it's not merely "ULP-like"?

[nsajko]

bump

[StefanKarpinski]

@stevengj, if that addresses your concerns and @oscardssmith is happy, can we merge this?

[oscardssmith]

I don't think sinc/cosc should exist in Base, but while they do, they might as well be good. Approved.