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

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