fix accuracy of logcosh(::Union{Float32, Float64})
#101
Conversation
Plots of the error in ULPs after this change
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. |
oscardssmith has been reviewing the `ULPError` code on JuliaLang/julia#59087 So sync the code with the changes there.
|
|
||
| The kernel of `logcosh`. | ||
|
|
||
| The polynomial coefficients were found using Sollya: |
There was a problem hiding this comment.
Why Sollya and not Remez.jl as the kernel for log1pmx? Would Remez.jl give a different result?
There was a problem hiding this comment.
Last time I checked out Remez.jl, it only implemented the Remez algorithm. That is, Remez.jl returns the coefficients in multiple precision, while Sollya's fpminimax further takes care of rounding the coefficients to machine precision without unnecessary loss of accuracy. See https://hal.science/inria-00119513
Therefore, we see that the general situation for L∞ approximation by real polynomials can be considered quite satisfying. The problem for the scientist that implements in software or hardware such approximations is that he uses finite-precision arithmetic and unfortunately, most of the time, the minimax approximation given by Chebyshev’s theorem and computed by Remez’ algorithm has coefficients which are transcendental (or at least irrational) numbers, hence not exactly representable with a finite number of bits.
Thus, the coefficients of the approximation usually need to be rounded according to the requirements of the application targeted (for example, in current software implementations, one often uses FP numbers in IEEE single or double precision for storing the coefficients of the polynomial approximation). But this rounding, if carelessly done, may lead to an important loss of accuracy. For instance, if we choose to round to the nearest each coefficient of the minimax approximation to the required format (this yields a polynomial that we will call rounded minimax in the sequel of the paper), the quality of the approximation we get can be very poor.
As suggested by devmotion on PR JuliaStats#99.
As suggested by devmotion on PR JuliaStats#99.
Specialization on the callable argument types should happen anyway.
Fixes #100