Собственная функция asin () (с рядом Тейлора) неточна

Мне нужно написать свою собственную функцию asin () без библиотеки math.h с использованием серии Taylor. Он отлично работает для чисел между , но когда я близок к ограничениям, он останавливается с 1604 итерациями и поэтому является неточным.

Я не знаю, как сделать его более аккуратным. Любые предложения очень ценятся!

Код выглядит следующим образом:

#include  #include  #include  #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) printf("nan\n"); else { printf("%.10e\n", my_arcsin(x)); //printf("%.10e\n", asin(x)); } return 0; } } 

А также короткий список моих ценностей и ожидаемых:

 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 

    ПОЖАЛУЙСТА, УВЕДОМЛЕНИЕ: В этом случае я настоятельно рекомендую метод @ Bence, так как вы не можете ожидать медленного конвергентного метода с низкой точностью данных для получения произвольной точности.

    Однако я готов показать вам, как улучшить результат, используя ваш текущий алгоритм.

    Основная проблема заключается в том, что a и b растут слишком быстро и вскоре становятся inf (после примерно 150 итераций). Другая аналогичная проблема – my_pow(xp,n) быстро растет, когда n растет, однако в этом случае это не имеет большого значения, поскольку мы можем предположить, что входные данные попадают в диапазон [-1, 1] .

    Поэтому я только что изменил метод, который вы имеете дело с a/b , введя ab_ratio , см. Мой отредактированный код:

     #include  #include  #include  #define EPS 0.000000000001 #include  #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; } } 

    Для 0.99 (и даже 0.9999999 ) он скоро дает правильные результаты с более чем 10 значащими цифрами. Однако он приближается к 1 .
    На самом деле процесс работает почти 12 минут на моем ноутбуке, вычисляя --asin 1 , и текущий результат равен 1.570786871 после 3560000000 итераций.

    ОБНОВЛЕНО: прошло 1h51min, а результат 1.570792915 и счетчик итераций – 27340000000 .

    Несмотря на то, что радиус сходимости расширения серии, который вы используете, равен 1, поэтому серия в конечном счете сходится при -1

    Я предлагаю вам

    • используйте свой оригинальный алгоритм для | x | <= 1 / sqrt (2),
    • используйте arcsin (x) = pi / 2 – arcsin (sqrt (1-x ^ 2)) для 1 / sqrt (2)
    • используйте arcsin (x) = -pi / 2 + arcsin (sqrt (1-x ^ 2)) для -1.0 <= x <-1 / sqrt (2).

    Таким образом, вы можете преобразовать свой вход x в [-1 / sqrt (2), 1 / sqrt (2)], где конвергенция относительно быстрая.