Looking for an arcsin algorithm

Question

Does anyone have a simple algorithm for computing a reasonably accurate arcsine?

By "simple" I mean some sort of polynomial that requires <= 5 multiplications per output sample. And by "reasonably accurate" I mean an algorithm whose error is no more than 10% when the input argument is close to plus or minus one.

I searched the web for a while but found nothing immediately useful.

Answer 1

Here's just a polynomial version:

$$ \arcsin(x) = x + \frac{1}{2} \frac{x^3}{3} + \frac{1 \cdot 3}{2 \cdot 4} \frac{x^5}{5} + \frac{1\cdot 3 \cdot 5}{2 \cdot 4 \cdot 6} \frac{x^7}{7} $$

function y = arcsin_test3(x)
    y = x.*(1+x.*x.*(1/6+ x.*x.*(3/(2*4*5) + x.*x.*((1*3*5)/(2*4*6*7)))))
endfunction

which seems to have five multiplies (assuming you can save the result of x.*x) and three additions.

And the scilab plot is:

Top is scilab's asin vs this one, bottom is the error between the two.

---

Original Answer

The square root here might be a hassle, but I thought I'd write it up because it looks like fun. :-)

This page suggests:

from page 81 of the Handbook of Mathematical Functions, by Milton Abramowitz and Irene Stegun:

$$ \arcsin(x) = \pi/2 - \sqrt{1 - x}(a_0 + a_1*x + a_2*x^2 + a_3*x^3), $$ where $$ a_0 = 1.5707288\\ a_1 = -0.2121144\\ a_2 = 0.0742610\\ a_3 = -0.0187293 $$

I've implemented this in scilab and it works OK, except around $x= -1$. Just reflecting the $0 \le x \le 1$ over to $-1 \le x \le 0$ makes for a much better approximation.

The top plot shows scilab's asin function against the above approximation (in dashed red) against my change in green.

The bottom plot shows the error for my change (plotting that and the original on the same axes means the green looks zero everywhere).

// 25770
function y = arcsin_test(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293
    
    xx = abs(x)
    
    y = %pi/2 - sqrt(1-x).*(a0 + a1*x + a2.*x.*x + a3.*x.*x.*x)
    
endfunction

function y = arcsin_test2(x)
    a0 = 1.5707288
    a1 = -0.2121144
    a2 = 0.0742610
    a3 = -0.0187293
    
    xx = abs(x)
    
    y = %pi/2 - sqrt(1-xx).*(a0 + a1*xx + a2.*xx.*xx + a3.*xx.*xx.*xx)
    
    y = y.*sign(x); 
endfunction

x = [-1: .0100001 : 1];

clf
subplot(211)
plot(x,arcsin_test2(x),'g.');
plot(x,arcsin_test(x),'r:');
plot(x,asin(x))
subplot(212)
//plot(x,(arcsin_test(x) - asin(x)),'r:')
plot(x,(arcsin_test2(x) - asin(x)),'g.')

Answer 2

I have a pretty good implementation of $\arctan()$ here.

I think you can use the identity:

$$ \arcsin(x) = \arctan\left( \frac{x}{\sqrt{1-x^2}} \right) $$

to get what you want.

Answer 3

The central part of the curve isn't a real problem as it is fairly linear and the Taylor approximation to two or three terms is a good starting point (least squares polynomial fit slightly better).

The sides are more problematic because of the infinite slope. A way to cope is via the transform

$$\arcsin(x)=\frac\pi2-\arcsin(\sqrt{1-x^2}),$$

which involves a square root.

---

If your argument $z$ is represented with floating-point, a fast approximation of the square root is obtained by halving the exponent and applying a linear transform to the mantissa.

Let $z=m2^e$, with $1\le m<2$, then $\sqrt z=\sqrt m\,2^{e/2}$. You can approximate $\sqrt{m}$ by $(\sqrt2-1)(m+\sqrt2)$.

• take the exponent $e$ apart (clearing it yields the representation of $m$);
• if $e$ is even, compute $(\sqrt2-1)(m+\sqrt2)$;
• if $e$ is odd, compute $\sqrt2(\sqrt2-1)(m+\sqrt2)$;
• set the exponent $\lfloor e/2\rfloor$.

Answer 4

So we have two requirements:

1. Use at most 5 multiplications.

2. Errror is no more than 10%.

Requirement 2 means that we want to optimize (minimze) the relative error of the approximation.

Requiring 5 multiplications + 3 additions

I'll interpret requirement 1) such that divisions are not wanted because they are usually (much) more expensive than multiplictatios.

In order to minimize the relative error, we'll have to "get rid" of the zero of $\arcsin(x)$ at $x=0$. We achieve this by approximating $\arcsin(x) / x$. Also we want to take advantage of the symmetry of $\arcsin$. Thus, we approximate $$ f(x) = \frac{\arcsin \sqrt x}{\sqrt x} \quad\text{ over }\quad [0,1]\tag 1 $$ which has no zeros, and then use $$ \arcsin x = x\cdot f(x^2) \tag2 $$

Equation (2) costs 2 multiplications, thus we have 3 multiplications left to approximate $f$ by a polynomial of degree 3. Remez algorithm yields

$$ f(x) = 0.9678828 + 0.8698691 x - 2.166373 x^2 + 1.848968 x^3 $$ with a relative error bounded by 0.033 which is well below the 10% margin. As the relative error of $f(x)$ is the same like the relative error of $xf(x^2)$, we are done. A C-ish implementation would read:

// Approximate arcsin(x) with a relative error better than 0.033.
// Costs = 5 multiplications, 3 additions.
double approx_asin (double x)
{
   double a0 = +0.9678828;
   double a1 = +0.8698691;
   double a2 = -2.166373;
   double a3 = +1.848968;
   double q = x * x;
   return x * (a0 + q * (a1 + q * (a2 + q * a3)));
}

which costs 5 multiplications and 3 additions. Here is a plot of the relative error of $f(x)$:

Requiring 3 multiplications + 1 addition

Now $\arcsin$ has singularities at $\pm1$ which means that polynomial approximations won't converge well as the polynomial degree increases. The other way round, we don't lose much of the accuracy when we decrease the polynomial degree, and indeed a polynomial of degree 1 has a relative error of below 0.084 which is better that the required 10%:

// Approximate arcsin(x) with a relative error better than 0.084.
// Costs = 3 multiplications, 1 addition.
double approx_asin (double x)
{
   double q = x * x;
   return x * (0.916735 + 0.523269 * q);
}

The relative error of $f(x)$ looks like:

Finally, here is a Desmos plot of the two approximations. It shows $\arcsin$ in black, the approximations (dotted) and the relative errors (multiplied by 10 for display): In red the relative error for the degree-3-approximaton, and in green for the degree-1-approximaton.

Approximating arcsin and arccos to IEEE single precision

In order to achieve IEEE single precision, using polynomials or rational functions alone won't work well due to the singularities (branch points) at $\pm1$. As $\cos x\approx 1-x^2/2$ near $x=0$, we have $\arccos x\approx \sqrt{2-2x}$ near $x=1$. Doing the math, we find that $$ \arccos x = \sqrt{2-2x}\,a(1-x) $$ where $a(x)$ is analytic at $x=0$ with radius of convergence of $2$. Actually, it's $$ a(x) = \frac{\arcsin \sqrt{x/2}}{\sqrt{x/2}} = 1 + \frac1{12}x+\frac3{160}x^2+{\mathcal O}(x^3) $$ It's not advisable to approximate $a(x)$ over all of $[0,2]$ (or over $[0,1]$ for that matter). Instead, we approximate $a(x)$ over $[0,1/2]$ and do a piece-wise approximation:

• If $|x| \leqslant 1/2$, compute $A_x = x\cdot a(2x^2)$.
• If $|x| > 1/2$, compute $C_x = \sqrt{2z}\cdot a(z)$ where $z=1-|x|$.

Then use:

$$\begin{align} \arcsin x &= \begin{cases} A_x &;\text{ if } |x|\leqslant 1/2 \\ (\pi/2-C_x)\operatorname{sign}(x) &;\text{ if } |x|> 1/2 \\ \end{cases} \\ \\ \arccos x &= \begin{cases} C_x &;\text{ if } x> 1/2 \\ \pi/2-A_x &;\text{ if } |x|\leqslant 1/2 \\ \pi-C_x &;\text{ if } x<-1/2 \\ \end{cases} \end{align}$$ where we approximate $a(x)$ by a MiniMax rational function of degree [2/2]: $$ a(x)\approx\frac {45.210185257899 - 18.617417552712 x + x^2} {45.210185141956 - 22.384922725383 x + 2.0175735681637 x^2} {} $$ This approximates $\arccos$ with a relative error of less than 2.6·10⁻⁹ (which is better than 28.5 bits). And it approximates $\arcsin$ with a relative error of less than 5·10⁻⁹ (which is better than 27.5 bits).

• Desmos plot of the relative errors

If one division costs more than 2 multiplications and 1 addition, then it's more efficient to approximate $a(x)$ by a polynomial of degree 5 with coefficients

 a0 = 0.99999999552357
 a1 = 0.08333395776062
 a2 = 0.018735953473147
 a3 = 0.0056946210646921
 a4 = 0.0014889678470141
 a5 = 0.0013342813628708

evaluated using Horner's method. The resulting precision is around one bit less than with using the [2/2] rational approximation for $a$.

Answer 5

You can use the Infinite Series to compute arctan, then you can use arctan & sqrt to compute arcsin.

Note: edge-cases not intercepted - might explode at the points of discontinuity, aka divZero.

public MyRational Arctan()
{
    // https://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Infinite_series
    MyRational sum = new MyRational(0);
    MyRational epsilon = new MyRational(1, 1000000000000);

    for (int n = 0; true; ++n)
    {
        MyRational dividend = MyRational.MinusOne.Pow(n) * this.Pow(2 * n + 1);
        MyRational divisor = new MyRational(2 * n + 1);

        MyRational quotient = dividend.Divide(divisor);

        MyRational newSum = sum + quotient;

        if (newSum.Subtract(sum).Abs() < epsilon)
        {
            return newSum;
        }

        sum = newSum;
    } // Next i 

}

public MyRational Arcsin()
{
    // https://en.wikipedia.org/wiki/Inverse_trigonometric_functions#Extension_to_complex_plane
    // https://dsp.stackexchange.com/questions/25770/looking-for-an-arcsin-algorithm
    MyRational divisor = Sqrt(One.Subtract(this.Pow(2)), new MyRational(1, 10000));
    MyRational x = this.Divide(divisor);

    return x.Arctan();
}

// https://www.geeksforgeeks.org/find-root-of-a-number-using-newtons-method/
public static MyRational Sqrt(MyRational radicand, MyRational epsilon)
{
    // Assuming the sqrt of radicand as radicand only 
    MyRational x = radicand;

    // The closed guess will be stored in the root 
    MyRational root;

    // To count the number of iterations 
    int count = 0;

    while (true)
    {
        count++;

        // Calculate more closed x 
        root = MyRational.OneHalf * (x + (radicand / x));

        // Check for closeness 
        if ((root - x).Abs() < epsilon)
            break;

        // Update root 
        x = root;
    } // Whend 

    return root;
} // End Function SquareRoot