Approximation of arcsin in C

Question

I've got a program that calculates the approximation of an arcsin value based on Taylor's series.

My friend and I have come up with an algorithm which has been able to return the almost "right" values, but I don't think we've done it very crisply. Take a look:

double my_asin(double x)
{
    double a = 0;
    int i = 0;
    double sum = 0;
    a = x;
    for(i = 1; i < 23500; i++)
    {
        sum += a;
        a = next(a, x, i);
    }
}

double next(double a, double x, int i)
{
    return a*((my_pow(2*i-1, 2)) / ((2*i)*(2*i+1)*my_pow(x, 2)));
}

I checked if my_pow works correctly so there's no need for me to post it here as well. Basically I want the loop to end once the difference between the current and next term is more or equal to my EPSILON (0.00001), which is the precision I'm using when calculating a square root.

---

This is how I would like it to work:

while(my_abs(prev_term - next_term) >= EPSILON)

But the function double next is dependent on i, so I guess I'd have to increment it in the while statement too. Any ideas how I should go about doing this?

---

Example output for -1:

$ -1.5675516116e+00

Instead of:

$ -1.5707963268e+00

---

Thanks so much guys.

Answer 1

using Taylor series for arcsin is extremly imprecise as the stuff converge very badly and there will be relatively big differencies to the real stuff for finite number of therms. Also using pow with integer exponents is not very precise and efficient.

However using arctan for this is OK

arcsin(x) = arctan(x/sqrt(1-(x*x)));

as its Taylor series converges OK on the <0.0,0.8> range all the other parts of the range can be computed through it (using trigonometric identities). So here my C++ implementation (from my arithmetics template):

T atan    (const T &x)                                              // = atan(x)
    {
    bool _shift=false;
    bool _invert=false;
    bool _negative=false;
    T z,dz,x1,x2,a,b; int i;
    x1=x; if (x1<0.0) { _negative=true; x1=-x1; }
    if (x1>1.0) { _invert=true; x1=1.0/x1; }
    if (x1>0.7) { _shift=true; b=::sqrt(3.0)/3.0; x1=(x1-b)/(1.0+(x1*b)); }
    x2=x1*x1;
    for (z=x1,a=x1,b=1,i=1;i<1000;i++)  // if x1>0.8 convergence is slow
        {
        a*=x2; b+=2; dz=a/b; z-=dz;
        a*=x2; b+=2; dz=a/b; z+=dz;
        if (::abs(dz)<zero) break;
        }
    if (_shift) z+=pi/6.0;
    if (_invert) z=0.5*pi-z;
    if (_negative) z=-z;
    return z;
    }
T asin    (const T &x)                                              // = asin(x)
    {
    if (x<=-1.0) return -0.5*pi;
    if (x>=+1.0) return +0.5*pi;
    return ::atan(x/::sqrt(1.0-(x*x)));
    }

Where T is any floating point type (float,double,...). As you can see you need sqrt(x), pi=3.141592653589793238462643383279502884197169399375105, zero=1e-20 and +,-,*,/ operations implemented. The zero constant is the target precision.

So just replace T with float/double and ignore the :: ...

Answer 2

Issues with your code and question include:

1. Your image file showing the Taylor series for arcsin has two errors: There is a minus sign on the x⁵ term instead of a plus sign, and the power of x is shown as xⁿ but should be x²ⁿ⁺¹.
2. The x factor in the terms of the Taylor series for arcsin increases by x² in each term, but your formula a*((my_pow(2*i-1, 2)) / ((2*i)*(2*i+1)*my_pow(x, 2))) divides by x² in each term. This does not matter for the particular value -1 you ask about, but it will produce wrong results for other values, except 1.
3. You ask how to end the loop once the difference in terms is “more or equal to” your epsilon, but, for most values of x, you actually want less than (or, conversely, you want to continue, not end, while the difference is greater than or equal to, as you show in code).
4. The Taylor series is a poor way to evaluate functions because its error increases as you get farther from the point around which the series is centered. Most math library implementations of functions like this use a minimax series or something related to it.
5. Evaluating the series from low-order terms to high-order terms causes you to add larger values first, then smaller values later. Due to the nature of floating-point arithmetic, this means that accuracy from the smaller terms is lost, because it is “pushed out” of the width of the floating-point format by the larger values. This effect will limit how accurate any result can be.
6. Finally, to get directly to your question, the way you have structured the code, you directly update a, so you never have both the previous term and the next term at the same time. Instead, create another double b so that you have an object b for a previous term and an object a for the current term, as shown below.

Example:

double a = x, b, sum = a;
int i = 0;
do
{
    b = a;
    a = next(a, x, ++i);
    sum += a;
} while (abs(b-a) > threshold);

Answer 3

so I guess I'd have to increment it in the while statement too

Yes, this might be a way. And what stops you?

int i=0;
while(condition){
   //do something
   i++;
}

Another way would be using the for condition:

for(i = 1; i < 23500 && my_abs(prev_term - next_term) >= EPSILON; i++)

Answer 4

Your formula is wrong. Here is the correct formula: http://scipp.ucsc.edu/~haber/ph116A/taylor11.pdf.

P.S. also note that your formula and your series are not correspond to each other.

---

You can use while like this:

while( std::abs(sum_prev - sum) < 1e-15 )
    {
        sum_prev = sum;
        sum += a;
        a = next(a, x, i);
    }