У вас есть разные неправильные представления о вашем коде.Сначала давайте посмотрим на формулу, которую вы предоставили:
sin(x) = sum((−1)^k * x^(2*k + 1) / (2*k + 1)! for x ∈ R; k = 0, ..., infinity
Функция синуса принимает действительное и возвращает действительное.Следовательно, вы должны использовать тип с плавающей точкой для x и sin (x).Используйте double.Давайте также напишем функцию, которая эмулирует sin из <math.h>:
double my_sin(double x);
Вышеприведенная серия точна, когда существует бесконечно много терминов.Конечно, мы не можем рассчитать это количество, и это также будет пустой тратой времени, потому что термины становятся все меньше, пока они больше не могут быть представлены double.Итак, давайте выберем максимальное количество терминов, скажем
enum {
nTerms = 8
};
Факториалы быстро растут.Обычный 32-битный int может содержать 12!= 479 001 60064-битный int может содержать 20!= 2 432 902 008 176 640 000.Поскольку мы собираемся использовать эти факториалы в расчете double, мы также можем использовать здесь double.Это даже позволит нам представлять 22!= 1 124 000 727 777 607 680 000 точно.
Ваша силовая функция также должна иметь double базу.Показатель является целым числом.(Но, пожалуйста, используйте более естественный порядок power(base, exp).
Наконец, (−1)^k - это просто знакопеременный знак. Он положителен, когда k четный, а в противном случае нечетный.
Помещение всехэто вместе:
double fact(int n)
{
double result = 1.0;
while (n > 0) {
result *= n;
n--;
}
return result;
}
double power(double a, int n)
{
double result = 1.0;
while (n > 0) {
result *= a;
n--;
}
return result;
}
enum {
nTerms = 8
};
double my_sin(double x)
{
double result = 0.0;
double sign = 1.0;
for(int k = 0; k < nTerms; k++)
{
double num = power(x, 2*k + 1);
double denom = fact(2*k + 1);
double term = sign * num / denom;
result = result + term;
sign = -sign;
}
return result;
}
Если мы напишем программу-драйвер для распечатки некоторых тестовых значений и сравним их с реализацией стандартной математической библиотеки sin:
int main(void)
{
for (int i = 0; i < 15; i++) {
double x = 0.1 * i;
double m = my_sin(x); // series approximation
double s = sin(x); // <math.h> implementation
printf("%16g%16g%16g%16g\n", x, m, s, m - s);
}
return 0;
}
, мы можем увидетьчто мы не так уж плохо:
x my_sin(x) sin(x) difference
-------- ------------ ------------ ------------
0 0 0 0
0.1 0.0998334 0.0998334 1.38778e-17
0.2 0.198669 0.198669 2.77556e-17
0.3 0.29552 0.29552 0
0.4 0.389418 0.389418 -5.55112e-17
0.5 0.479426 0.479426 0
0.6 0.564642 0.564642 0
0.7 0.644218 0.644218 0
0.8 0.717356 0.717356 0
0.9 0.783327 0.783327 -4.44089e-16
1 0.841471 0.841471 -2.77556e-15
1.1 0.891207 0.891207 -1.43219e-14
1.2 0.932039 0.932039 -6.20615e-14
1.3 0.963558 0.963558 -2.42029e-13
1.4 0.98545 0.98545 -8.52318e-13
(но чем дальше мы идем от нуля, тем хуже. Попробуйте другие значения для nTerms.)
I 'Вы сказали в комментарии выше, что вам не нужно вычислять факториалы и степени, и это правда. Если вы посмотрите на условия ряда, вы увидите, что:
s[n] = -1 * s[n - 1] * x^2 / (2*n * (2*n +1))
s[0] = x
s[1] = x^3 / (1 * 2 * 3) = x * x^2 / (2 * 3)
s[2] = x^5 / (1 * 2 * 3 * 4 * 5) = x^3 / (1 * 2 * 3) * x^2 / (4 * 5)
s[3] = ...
Вот функция, которая реализуетЭто вычисляет термины, пока их добавление к сумме не изменит их, потому что они слишком малы:
double sin_r(double x)
{
double sum = x;
double a = x;
int n;
for (n = 1; ; n++) {
double was = sum;
a = -a * x*x / (2*n) / (2*n + 1);
sum += a;
if (was == sum) break;
}
return sum;
}
Сложение по-прежнему теряет некоторую точность, суммируя сначала первые слагаемые, но оно имеет преимуществочто он не должен рассчитывать факторииLS и полномочия.Вам даже не нужно <math.h>.