Какой более точный алгоритм я могу использовать для вычисления синуса числа? - PullRequest
Новая
       24

Какой более точный алгоритм я могу использовать для вычисления синуса числа?

2 голосов
/ 21 сентября 2019

У меня есть этот код, который вычисляет предположение для синуса и сравнивает его с результатом стандартной библиотеки C (glibc в моем случае): 

#include <stdio.h>
#include <math.h>

double double_sin(double a)
{
    a -= (a*a*a)/6;

    return a;
}

int main(void)
{
    double clib_sin = sin(.13),
             my_sin = double_sin(.13);
    printf("%.16f\n%.16f\n%.16f\n", clib_sin, my_sin, clib_sin-my_sin);
    return 0;
}

Точность для double_sin низкая (около 5-6 цифр).Вот мой вывод: 

0.1296341426196949
0.1296338333333333
0.0000003092863615

Как видите, после .12963 результаты отличаются.

Некоторые примечания:

  • Я понимаюНе думаю, что ряд Тейлора будет работать для этой конкретной ситуации, факториалы, необходимые для большей точности, не могут быть сохранены внутри unsigned long long.

  • Таблицы поиска не являются опциейони занимают слишком много места и, как правило, не предоставляют никакой информации о том, как рассчитать результат.

  • Если вы используете магические числа, объясните их (хотя я бы предпочел, чтобы онине использовались).

  • Я бы предпочел, чтобы алгоритм был легко понятен и мог использоваться в качестве справочного по сравнению с тем, который не является.

  • Результат не должен быть абсолютно точным.Минимумом будут требования IEEE 754, C и / или POSIX.

  • Я использую формат IEEE-754 double, на который можно положиться.

  • Поддерживаемый диапазон должен быть как минимум от -2*M_PI до 2*M_PI.Было бы неплохо, если бы было включено уменьшение диапазона.

Какой более точный алгоритм я могу использовать для вычисления синуса числа?

У меня былоидея о чем-то похожем на Ньютона-Рафсона, но вместо этого для вычисления синуса. Однако я ничего не смог найти на нем и исключаю эту возможность.

1 Ответ

3 голосов
/ 21 сентября 2019

С серией Тейлора вы можете довольно близко подойти.Хитрость заключается не в том, чтобы рассчитать полный факториал на каждой итерации.

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

sin(x) = x^1/1! - x^3/3! + x^5/5! - x^7/7!

Глядя на термины, вы вычисляете следующий член, умножая числитель наx ^ 2, умножая знаменатель на следующие два числа в факториале и меняя знак.Затем вы останавливаетесь, когда добавление следующего термина не меняет результат.

Таким образом, вы можете закодировать его так: 

double double_sin(double x)
{
    double result = 0;
    double factor = x;
    int i;

    for (i=2; result+factor!=result; i+=2) {
        result += factor;
        factor *= -(x*x)/(i*(i+1));
    }
    return result;
}

Мой вывод: 

0.1296341426196949
0.1296341426196949
-0.0000000000000000

РЕДАКТИРОВАТЬ:

Точность может быть дополнительно увеличена, если термины добавляются в обратном направлении, однако это означает вычисление фиксированного количества терминов: 

#define FACTORS 30

double double_sin(double x)
{
    double result = 0;
    double factor = x;
    int i, j;
    double factors[FACTORS];

    for (i=2, j=0; j<FACTORS; i+=2, j++) {
        factors[j] = factor;
        factor *= -(x*x)/(i*(i+1));
    }
    for (j=FACTORS-1;j>=0;j--) {
        result += factors[j];
    }
    return result;
}

Эта реализация теряет точность, если x выходит за пределы диапазона от 0 до 2 * PI.Это можно исправить, вызвав x = fmod(x, 2*M_PI); в начале функции для нормализации значения.

...