创建自定义正弦函数

Creating a custom sine function

我尝试使用 cTaylor Series 创建一个自定义 sine 函数来计算系列中有 10 个项的 sin,但我得到了当我尝试在 x > 6.

中查找 sine(x) 时出现错误结果

它适用于 -5 < x < 5,但超出该范围的任何内容都不会产生正确的结果。

我希望 sin(10) 到 return 接近 -0.5440 的值,但得到 1418.0269775391

我已将所有内容放在一个文件中,这样更容易。

#include <stdio.h>
#include <stdlib.h>

double factorial(double n);
double power(double n, double pow);
double sine(double n);

// This is supposed to all go in a .c file and reference the .h stuff above
// This is the actual implementation of the functions declared above
double factorial(double n) {
    // 0! = 1 so just return it
    if(n == 0) {
        return 1;
    }
    // Recursively call factorial with n-1 until n == 0
    return n * (factorial(n - 1)); 
}


double power(double n, double power) {
    double result = n;
    // Loop as many times as the power and just multiply itself power amount of times
    for(int i = 1; i < power; i++) {
        result = n * result;
    }
    return result;
}

double sine(double n) {
    double result = n;
    double coefficent = 3; // Increment this by 2 each loop
    for(int i = 0; i < 10; i++) { // Change 10 to go out to more/less terms
        double pow = power(n, coefficent);
        double frac = factorial(coefficent);
        printf("Loop %d:\n%2.3f ^ %2.3f = %2.3f\n", i, n, coefficent, pow);
        printf("%2.3f! = %2.3f\n", coefficent, frac);

        // Switch between adding/subtracting
        if(i % 2 == 0) { // If the index of the loop is divided by 2, the index is even, so subtract
            result = result - (pow/frac); // x - ((x^3)/(3!)) - ((x^5)/(5!))...
        } else {
            result = result + (pow/frac); // x - ((x^3)/(3!)) + ((x^5)/(5!))...
        }
        coefficent = coefficent + 2;
        printf("Result = %2.3f\n\n", result);
    }
    return result;
}

// main starting point. This is suppossed to #include "functions.c" which contain the above functions in it
int main(int argc, char** argv) {

    double number = atof(argv[1]); // argv[1] = "6"
    double sineResult = sine(number);
    printf("%1.10f", sineResult);
    return (0);
}

进行更正后,如我对问题的评论中所列,建议的代码如下所示:

#include <stdio.h>
#include <stdlib.h>

double factorial(double n);
double power(double n, double pow);
double sine(double n);

// This is supposed to all go in a .c file and reference the .h stuff above
// This is the actual implementation of the functions declared above
double factorial(double n) {
    // 0! = 1 so just return it
    if(n == 0) {
        return 1;
    }
    // Recursively call factorial with n-1 until n == 0
    return n * (factorial(n - 1));
}


double power(double n, double power) {
    double result = n;
    // Loop as many times as the power and just multiply itself power amount of times
    for(int i = 1; i < power; i++) {
        result = n * result;
    }
    return result;
}

double sine(double n) {
    double result = n;
    double coefficent = 3.0; // Increment this by 2 each loop

    for(int i = 0; i < 10; i++) { // Change 10 to go out to more/less terms
        double pow = power(n, coefficent);
        double frac = factorial(coefficent);
        printf("Loop %d:\n%2.3f ^ %2.3f = %2.3f\n", i, n, coefficent, pow);
        printf("%2.3f! = %2.3f\n", coefficent, frac);

        // Switch between adding/subtracting
        if(i % 2 == 0) { // If the index of the loop is divided by 2, the index is even, so subtract
            result = result - (pow/frac); // x - ((x^3)/(3!)) - ((x^5)/(5!))...
        } else {
            result = result + (pow/frac); // x - ((x^3)/(3!)) + ((x^5)/(5!))...
        }
        coefficent = coefficent + 2;
        printf("Result = %2.3f\n\n", result);
    }
    return result;
}


// main starting point. This is suppossed to #include "functions.c" which contain the above functions in it
int main( void )
{
    double number = atof("6");
    double sineResult = sine(number);
    printf("%1.10f", sineResult);
    return (0);
}

结果输出如下:

Loop 0:
6.000 ^ 3.000 = 216.000
3.000! = 6.000
Result = -30.000

Loop 1:
6.000 ^ 5.000 = 7776.000
5.000! = 120.000
Result = 34.800

Loop 2:
6.000 ^ 7.000 = 279936.000
7.000! = 5040.000
Result = -20.743

Loop 3:
6.000 ^ 9.000 = 10077696.000
9.000! = 362880.000
Result = 7.029

Loop 4:
6.000 ^ 11.000 = 362797056.000
11.000! = 39916800.000
Result = -2.060

Loop 5:
6.000 ^ 13.000 = 13060694016.000
13.000! = 6227020800.000
Result = 0.037

Loop 6:
6.000 ^ 15.000 = 470184984576.000
15.000! = 1307674368000.000
Result = -0.322

Loop 7:
6.000 ^ 17.000 = 16926659444736.000
17.000! = 355687428096000.000
Result = -0.275

Loop 8:
6.000 ^ 19.000 = 609359740010496.000
19.000! = 121645100408832000.000
Result = -0.280

Loop 9:
6.000 ^ 21.000 = 21936950640377856.000
21.000! = 51090942171709440000.000
Result = -0.279

-0.2793866930

泰勒展开有一个错误,这取决于参数范围和泰勒展开的顺序。我相信你已经越界了。有关更多示例,请参见此处:www.dotancohen.com/eng/taylor-sine.php

正如我在

中所说的

The real Taylor expansion centered in x0 is:

where Rn is the Lagrange Remainder

Note that Rn grows fast as soon as x moves away from the center x0.

Since you are implementing the Maclaurin series (Taylor series centered in 0) and not the general Taylor series, your function will give really wrong results when trying to calculate sin(x) for big values of x.

因此,在 sine() 函数中的 for 循环之前,您必须将域至少减少到 [-pi, pi]... 更好如果将其减少到 [0, pi] 并利用正弦奇偶校验。

要修复您的代码,您需要 math.h 中的 fmod(),因此您可以:

#include <math.h>

// Your code

double sine (double n) {
    // Define PI
    const double my_pi = 3.14159265358979323846;
    // Sine's period is 2*PI
    n = fmod(n, 2 * my_pi);
    // Any negative angle can be brought back
    // to it's equivalent positive angle
    if (n < 0) {
        n = 2 * my_pi - n;
    }
    // Sine is an odd function...
    // let's take advantage of it.
    char sign = 1;
    if (n > my_pi) {
        n -= my_pi;
        sign = -1;
    }
    // Now n is in range [0, PI].

    // The rest of your function is fine

    return sign * result;
}

现在如果你真的讨厌 math.h 模块,你可以像这样实现你自己的 fmod()

double fmod(double a, double b)
{
    double frac = a / b;
    int floor = frac > 0 ? (int)frac : (int)(frac - 0.9999999999999999);
    return (a - b * floor);
}

Try it online!