Собственная функция asin() (с рядом Тейлора) не точна
Мне нужно написать свою собственную функцию asin() без библиотеки math.h с использованием ряда Тейлора. Он отлично работает для чисел между <-0,98; 0,98>, но когда я близок к пределам, он останавливается с 1604 итерациями и поэтому является неточным.
Я не знаю, как сделать его более аккуратным. Любые предложения очень ценятся!
Код следующий:
#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;
}
}
А также краткий список моих ценностей и ожидаемых:
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
2 ответа
ПОЖАЛУЙСТА, ОБРАТИТЕ ВНИМАНИЕ: в этом случае я настоятельно рекомендую метод @Bence, поскольку вы не можете ожидать, что медленно сходящийся метод с низкой точностью данных получит произвольную точность.
Однако я хочу показать вам, как улучшить результат, используя ваш текущий алгоритм.
Основная проблема заключается в том, что a а также b растет слишком быстро и скоро станет inf (после примерно 150 итераций). Еще одна похожая проблема my_pow(xp,n) быстро растет, когда n растет, однако это не имеет большого значения в этом самом случае, так как мы можем предположить, что входные данные попадают в диапазон [-1, 1],
Так что я просто изменил метод, с которым вы имеете дело a/b введя ab_ratioсм. мой отредактированный код:
#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;
}
}
За 0.99 (и даже 0.9999999) это скоро дает правильные результаты с более чем 10 значащими цифрами. Однако это становится медленным, когда приближается к 1,
На самом деле этот процесс был запущен в течение почти 12 минут на моем ноутбуке расчета --asin 1и текущий результат 1.570786871 после 3560000000 итераций.
ОБНОВЛЕНО: это было 1h51min сейчас и результат 1.570792915 и счетчик итераций 27340000000,
Даже несмотря на то, что радиус сходимости раскрываемого ряда, который вы используете, равен 1, следовательно, ряд в конечном итоге будет сходиться при -1 Я предлагаю вам Таким образом, вы можете преобразовать свои входные данные x в [-1/sqrt(2),1/sqrt(2)], где конвергенция относительно быстрая.