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

• 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

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−9 (which is better than 28.5 bits). And it approximates $$\arcsin$$ with a relative error of less than 5·10−9 (which is better than 27.5 bits).

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

• Nice stuff! How did you compute the rational minimax approximation? Aug 29 at 16:10
• To compute MiniMax polynomials / rational functions I am using my own go on the Remez algorithm. I didn't have a Computer Algebra System handy back then and implemented it all by myself (Python, then C++); including Linear Algebra based on GNU MPFR. I think every major CAS or Numerics Package like Octave or Mathematica would come with a decent implementation of Remez. Aug 29 at 17:54
• Yes, I know the Remez exchange algorithm, but I mostly know it for polynomial Chebyshev approximation. I haven't tried it for rational functions, and I wasn't sure if it's any different in that case and if convergence is guaranteed. Aug 29 at 18:34
• From my experience Remez works reasonable for rational functions provided the target function is smooth. And when optimizing for relative error like with the current question, one has to get rid of zeros, usually by dividing them out. For the [2/2] $a(x)$ from my answer, it takes 15 iterations until the change in height ($\sum |a_i|$) is so small that any further improvements will no more lead to better approximations. Aug 29 at 19:25

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