Problem with least squares rational approximation to `asin(x)+ sqrt(1-x^2)` in [3,1] form

I’m trying to generate a decent [3,1] rational least squares polynomial approximation to asin(x)+sqrt(1-x^2) on [0,1] and failing dismally 🙁
The problem is that it has a pole for this particular combination and also a handful of others (approx. 1%) of those I tested rational approximations to this function go completely haywire. This is the first step in a Remez equal ripple code.

I have ported the original Julia code which works in arbitrary precision (so numerical instability can be ruled out) to C and get about the same results as the Julia code for modest N,D. I still find debugging easier in C/C++. The algebra all looks fine to me and it works in ~99% of cases for N, D <= 20. It uses a trick described here asin when fast sqrt available. This is of interest again because hardware sqrt is now so fast.

Unfortunately the equation is rather fickle and rational Remez solutions to it are a bit thin on the ground. When the original author tried it he listed all the solutions that were soluble in Maple/Mathematica with digits = 30 for N, D <= 20 at the end. My improvements have got past some of the obstacles. I expect Maple & Mathematica have also improved in the meantime.

The remaining failing cases in the full Julia implementation are (my C code can duplicate them up to N+D ~9):

Total failures  42/400

There is a sort of pattern visible in the failures. Notably that the first failures occur in pairs at N, D=5N and D=5N+2 [1,5][1,7], [2,10][2,12], [3,15][3,17] and the other heuristic is that failing ones tend to be at an offset at N+1, D+3 or N+1, D+5 from a previous failing one at N,D. I have no idea why this is.

My concrete question is can someone else please try and construct a rational polynomial [3,1] that will fit this function and does not have a pole in the denominator or if that isn’t possible (and I begin to think that it might not be suggestions on how else to tame asin(x)as x->1). I’m using x_n = sin(pi*n/N)^2 Chebyshev for the x coordinates which places more emphasis on fitting the ends of the range (equal spaced is worse). I’m wondering about increasing the density of points sampled at larger x values. Increasing the weights there didn’t help.

These are the (truncated coefficients for the various polynomials along the slice N+D = 4 to which [3,1] belongs.

These are the respective error graphs on the same scale for each of the approximations above (the adjusted [3,1] isn’t shown but just looks like a W). A correct solution with 6 extrema should look like W\ or /W ie HLHLHL or the opposite. Note the pole going off to infinity on the red [3,1] curve at x = 0.8897.

AOBE I would probably use [2,2] or [3,2] as the best trade off on simplicity vs accuracy. But [3,1] would be very handy if it exists.

I have omitted the customary MRE for the moment since it isn’t really minimal. I hope to get it down to under 400 lines without removing any of the crucial explanatory comments. It is a one to one code match for the ARM Remez Julia as far as is possible which makes the C rather ugly (I’m not especially proud of it but it does work and gets the same answers as Julia in the range N+D< 9. Obviously numerical stability is an issue with only double precision and the matrix solver from Julia is crude but it works well enough for [3,1] not to be a problem.

I have tried weighting the fit and various analytic transformations to move excess coefficient value from the denominator into the numerator and obtained a sort of solution at the bottom of the table but which doesn’t oscillate enough to be acceptable input for Remez. Remez also tends to restore the pole and/or create new ones (but that is another story).

My goal is to have N+D+2 extrema in the rational least squares approximation including the two end points. An interesting curiosity is that [0,4] actually has N+D+3 extrema due to serendipity. The other low order ones giving problems with poles are [1,5], [1,7], [4,4] and [6,1]. I have focused on [3,1] specifically because I had hoped that analytic methods for cubic equations might allow me to find a well behaved solution – if one exists. That approach got me past the previous gotcha which triggered on a quadratic form [0,2] and I had hoped that was the last bug in this section.

If an ideal equal ripple solution to target f(x) with error e and 6 extrema does exist we have a couple of hard algebraic constraints that can be used to decrease the search space or to sanity check putative answers. Based on the values taken at the ends of the range x = 0 and x=1.

 r(x) = (a + bx + cx^2 + dx^3)/(1-gx)

 f(0) - r(0) = e = 1 - a    IOW a = 1 - e
 f(1) - r(0) = -e = pi/2 - (a+b+c+d)/(1-g)
 IOW
 a+b+c+d = (pi/2+e)(1-g)
 so d = (pi/2+e)(1-g)-1+e+b+c
 And also to avoid having a pole in the range 0 - 1
 g < 1 

I’m now looking at tricks from solving the cubic and knowing that a final solution with the right number of extrema in the target range 0-1 must have a derivative with a numerator with all real roots in the same range. This approach is proving tough. I suspect possibly because no such solution exists (brute force on the four remaining free parameters hasn’t found one yet either). I can take a guess at a value for e ~ 0.0014 and then try all legal combos of b,c,d,g but it still isn’t yeilding…

You need to sign in to view this answers