# Looking for an arcsin algorithm

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.

• This might give some ideas stackoverflow.com/questions/5920467/… Sep 12, 2015 at 0:05
• But why not just a lookup table? Sep 12, 2015 at 0:05
• 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. Sep 12, 2015 at 1:25
• 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.
– user7657
Sep 12, 2015 at 8:52
• What about CORDIC, which only takes a few additions and subtractions and no multiplications. Sep 12, 2015 at 20:10

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

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

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.')
• "Handbook of Mathematical Functions" love that book Sep 12, 2015 at 2:46
• 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.) Sep 12, 2015 at 9:05
• 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. Sep 12, 2015 at 22:42
• 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! Sep 13, 2015 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.

• 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? Sep 13, 2015 at 7:42
• 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$. Sep 13, 2015 at 19:27
• 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. Sep 14, 2015 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

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

• I look forward to experimenting with that interesting square root algorithm! Sep 13, 2015 at 7:02
• With these coefficients, the error is always positive. By a slight adjustment, we can make it symmetric and halve the maximum error.
– user7657
Sep 13, 2015 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)
{

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