Does anyone have a simple algorithm for computing a reasonably accurate arcsine? By "simple" I mean some sort of polynomial that requires <= 5 multiplies per output sample. And by "reasonably accurate" I mean an algo 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.

  • $\begingroup$ This might give some ideas stackoverflow.com/questions/5920467/… $\endgroup$ – geometrikal Sep 12 '15 at 0:05
  • $\begingroup$ But why not just a lookup table? $\endgroup$ – geometrikal Sep 12 '15 at 0:05
  • $\begingroup$ I'm thinking of an implementation where available memory is painfully limited. So I didn't consider any 'look-up table' solution. Thanks for your thoughts. $\endgroup$ – Richard Lyons Sep 12 '15 at 1:25
  • $\begingroup$ Do you allow square roots ? Due to the behavior of the function close to $\pm1$ (infinite slope), a polynomial approximation doesn't work well there. $\endgroup$ – Yves Daoust Sep 12 '15 at 8:52
  • $\begingroup$ What about CORDIC, which only takes a few additions and subtractions and no multiplications. $\endgroup$ – mattgately Sep 12 '15 at 20:10

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)))))

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

And the scilab plot is:

enter image description here

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).

enter image description here

// 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)


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); 

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

//plot(x,(arcsin_test(x) - asin(x)),'r:')
plot(x,(arcsin_test2(x) - asin(x)),'g.')
  • 2
    $\begingroup$ "Handbook of Mathematical Functions" love that book $\endgroup$ – geometrikal Sep 12 '15 at 2:46
  • 1
    $\begingroup$ Oh shoot. I saw that wiki 'inverse trig functions' web page during my web search but I didn't scroll down far enough to see the 'infinite series' material! Shame on me. Peter K., that's another one I owe ya'. (My original problem was to improve the performance of a central-difference digital differentiator which, I believe, can be done by performing an arcsine operation.) $\endgroup$ – Richard Lyons Sep 12 '15 at 9:05
  • $\begingroup$ yeah, but Rick, you can't do an infinite series. if you're gonna make it finite, then the optimal coefficients won't be exactly what you get from truncating the infinite series. if you have MATLAB going (they have a relatively cheap "home use" license now), i can send you MATLAB code for doing Remez exchange on the function of your heart's desire. $\endgroup$ – robert bristow-johnson Sep 12 '15 at 22:42
  • 1
    $\begingroup$ UPDATE: I was able to solve my original problem (creating a very simple digital differentiator that has improved performance over that of a central-difference differentiator) without using an arcsin() function. I may post a blog at dsprelated.com describing my results. I thank everyone here for their help! $\endgroup$ – Richard Lyons Sep 13 '15 at 7:34

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.

  • $\begingroup$ Your link to various functions is interesting Robert. Just for giggles I tried to implement the sqrt(1+x) function. Instead of using the correct limits of 0 -to- 4, I screwed up and used 1 -to- 5. Of course I computed an incorrect result. However, when I multiplied my "incorrect" result by two I ended up with the correct result. Interesting, huh? $\endgroup$ – Richard Lyons Sep 13 '15 at 7:42
  • $\begingroup$ Rick, i'm pretty sure the functions are "correct" (or reasonably accurate) stated as they are with the limits of $x$ as stated. for $\sqrt{x}$, it's only good for $1 \le x \le 2$, so if you're between $2$ and $4$, then you'll have to have a little constant (the $\sqrt{2}$) stored in there and you'll need to know the difference between an even exponent of $2$ and an odd exponent of $2$. $\endgroup$ – robert bristow-johnson Sep 13 '15 at 19:27
  • $\begingroup$ I too believe your functions are correct. I was merely commenting on a silly 'summation limits' mistake I made, and in making my mistake I computed an incorrect result. But I noticed my incorrect result was exactly one half the correct result. I was just saying that my incorrect result had an "interesting" relationship to the correct result, that's all. Sorry for the confusion Robert. $\endgroup$ – Richard Lyons Sep 14 '15 at 11:26

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


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$.

enter image description here

  • $\begingroup$ I look forward to experimenting with that interesting square root algorithm! $\endgroup$ – Richard Lyons Sep 13 '15 at 7:02
  • $\begingroup$ With these coefficients, the error is always positive. By a slight adjustment, we can make it symmetric and halve the maximum error. $\endgroup$ – Yves Daoust Sep 13 '15 at 8:52

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)

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

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

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

    return root;
} // End Function SquareRoot 

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.