[LegionMammal978]

The coefficients given are indeed a near-optimal cubic minimax approximation for (π/2 - arcsin(x))/sqrt(1-x) on [0,1]. But those coefficients aren't actually optimal for approximating arcsin(x) itself.

For reference, the coefficients given are [1.5707288, -0.2121144, 0.0742610, -0.0187293]: if we optimize P(x) = (π/2 - arcsin(x))/sqrt(1-x) ourselves, we can extend them to double precision as [1.5707288189560218, -0.21211524058527342, 0.0742623449400704, -0.018729868776598532]. Increasing the precision reduces the max error, at x = 0, by 0.028%.

Adjusting our error function to optimize the absolute error of arcsin(x) = π/2 - P(x)*sqrt(1-x) on [0,1], we get the coefficients [1.5707583404833712, -0.2128751841625164, 0.07689738736091772, -0.02089203710669022]. The max error is reduced by 44%, from 6.75e-5 to 3.80e-5. If we plot the error function [0], we see that the new max error is achieved at five points, x = 0, 0.105, 0.386, 0.730, 0.967.

(Alternatively, adjusting our error function to optimize the relative error of arcsin(x), we get the coefficients [1.5707963267948966, -0.21441792645252514, 0.08365774237116316, -0.02732304481232744]. The max absolute error is 2.24e-4, but the max relative error is now 0.0181%, even in the vicinity of the root at x = 0. Though we'd almost certainly want to use a different formula to avoid catastrophic cancellation.)

So it goes to show, we can nearly double our accuracy, without modifying the code, just by optimizing for the right error metric.

[0] https://www.desmos.com/calculator/nj3b8rpvbs

[LegionMammal978]

Actually, we can improve this a bit further, by also adjusting the "π/2" constant in arcsin(x) = π/2 - P(x)*sqrt(1-x). We take coefficients [1.5707256467180715, -0.21298179775496026, 0.07727939759417458, -0.02132102849918157] for P(x), then take arcsin(x) = 1.570760986756484 - P(x)*sqrt(1-x). This reduces the max error by 6.97%, from 3.80e-5 to 3.53e-5.

Adjusting the "1" upward in sqrt(1-x) does not seem to help at all.

[kzrdude]

How did you find out that his optimization was done for a different equation, just by trial?

[LegionMammal978]

Just looking at the formula in the code (and the book it came from), we see that the approximation is of form arcsin(x) = π/2 - P(x)*sqrt(1-x). It is called a minimax solution in both, and the simplest form of minimax optimization is for polynomials. So we look at P(x) = (π/2 - arcsin(x))/sqrt(1-x): plotting out its error function with the original coefficients, it has the clear equioscillations that you'd expect from an optimized polynomial, i.e., each local peak has the exact same height, which is the max error. But if we look at the error curve in terms of arcsin(x), then its oscillations no longer have the same height, which indicates that the approximation can be improved.

[kzrdude]

Thank you for elaborating!