The asin () native function (with Taylor series) is inaccurate

Question

I need to write my own asin () function without the math.h library using the Taylor series. It works great for numbers between <-0.98; 0.98>, but when I am close to the limits, it stops with 1604 iterations and therefore is inaccurate.

I do not know how to make it more accurete. Any suggestions are greatly appreciated!

The following code:

#include <stdio.h> #include <stdlib.h> #include <string.h> #define EPS 0.000000000001 double my_arcsin(double x) { long double a, an, b, bn; a = an = 1.0; b = bn = 2.0; long double n = 3.0; double xn; double xs = x; double xp = x; int iterace = 0; xn = xs + (a/b) * (my_pow(xp,n) / n); while (my_abs(xn - xs) >= EPS) { n += 2.0; an += 2.0; bn += 2.0; a = a * an; b = b * bn; xs = xn; xn = xs + (a/b) * (my_pow(xp,n) / n); iterace++; } //printf("%d\n", iterace); return xn; } int main(int argc, char* argv[]) { double x = 0.0; if (argc > 2) x = strtod(argv[2], NULL); if (strcmp(argv[1], "--asin") == 0) { if (x < -1 || x > 1) printf("nan\n"); else { printf("%.10e\n", my_arcsin(x)); //printf("%.10e\n", asin(x)); } return 0; } } 

And also a short list of my values ​​and expected:

 My values Expected values my_asin(x) 5.2359877560e-01 5.2359877560e-01 0.5 1.5567132089e+00 1.5707963268e+00 1 //problem 1.4292568534e+00 1.4292568535e+00 0.99 //problem 1.1197695150e+00 1.1197695150e+00 0.9 1.2532358975e+00 1.2532358975e+00 0.95 

Answer 1

PLEASE NOTE: In this case, I strongly recommend the @Bence method, since you cannot expect a slow converged method with low data accuracy to get arbitrary accuracy.

However, I am ready to show you how to improve the result using your current algorithm.

The main problem is that a and b grow too fast and soon become inf (after about 150 iterations). Another similar problem: my_pow(xp,n) grows rapidly when n grows, but in this case it does not really matter, since we can assume that the input data falls into the range [-1, 1] .

So, I just changed the method you work with a/b by typing ab_ratio , see my edited code:

 #include <stdio.h> #include <stdlib.h> #include <string.h> #define EPS 0.000000000001 #include <math.h> #define my_pow powl #define my_abs fabsl double my_arcsin(double x) { #if 0 long double a, an, b, bn; a = an = 1.0; b = bn = 2.0; #endif unsigned long _n = 0; long double ab_ratio = 0.5; long double n = 3.0; long double xn; long double xs = x; long double xp = x; int iterace = 0; xn = xs + ab_ratio * (my_pow(xp,n) / n); long double step = EPS; #if 0 while (my_abs(step) >= EPS) #else while (1) /* manually stop it */ #endif { n += 2.0; #if 0 an += 2.0; bn += 2.0; a = a * an; b = b * bn; #endif _n += 1; ab_ratio *= (1.0 + 2.0 * _n) / (2.0 + 2.0 * _n); xs = xn; step = ab_ratio * (my_pow(xp,n) / n); xn = xs + step; iterace++; if (_n % 10000000 == 0) printf("%lu %.10g %g %g %g %g\n", _n, (double)xn, (double)ab_ratio, (double)step, (double)xn, (double)my_pow(xp, n)); } //printf("%d\n", iterace); return xn; } int main(int argc, char* argv[]) { double x = 0.0; if (argc > 2) x = strtod(argv[2], NULL); if (strcmp(argv[1], "--asin") == 0) { if (x < -1 || x > 1) printf("nan\n"); else { printf("%.10e\n", my_arcsin(x)); //printf("%.10e\n", asin(x)); } return 0; } } 

For 0.99 (and even 0.9999999 ), he will soon give the correct results with more than 10 significant digits. However, it is approaching 1 . Actually, the process runs for almost 12 minutes on my laptop, calculating --asin 1 , and the current result is 1.570786871 after 3560000000 iterations.

UPDATED: Now it was 1h51min, and the result was 1.570792915 and the iteration counter was 27340000000 .

Answer 2

Despite the fact that the radius of convergence of the extension of the series that you use is 1, so the series will eventually converge for -1 <x <1, the convergence is really painfully slow near the limits of this interval. The solution is to somehow avoid these parts of the interval.

I suggest you

• use the original algorithm for | x | <= 1 / sqrt (2),
• use the identifier arcsin (x) = pi / 2 - arcsin (sqrt (1-x ^ 2)) for 1 / sqrt (2) <x <= 1.0,
• use the identifier arcsin (x) = -pi / 2 + arcsin (sqrt (1-x ^ 2)) for -1.0 <= x <-1 / SQRT (2).

This way you can convert your input x to [-1 / sqrt (2), 1 / sqrt (2)], where the convergence is relatively fast.