我应该如何处理 return 错误答案的函数？

Question

我们无法使用math.h

如果我的正弦函数输入大于 54 而我的 exp 函数输入大于 37，我会得到错误的答案；我想它溢出了，我该怎么办？我想自己写正弦函数和exp函数，我用的是泰勒展开。我只需要分数中的 6 位数字；

    //Start of exp function
float exp(float x)
{
    double expreturn=0;
    long long int fctrl=1;
    for(int n=0;n<=12;n++){

        fctrl=1;
        for(int i=2;i<=n;i++)
            fctrl *= i;
        expreturn += (pow(x,n)/fctrl);

    }

    return (float)expreturn;
}//End of exp function



//Start of sin function
float sin(float x)
{
    double sinus=0;
    long long int fctrl=1;
   while(!(0<=x&&x<6.3))
    {
        if(x<0)
            x += PI*2;
        else
            x -= PI*2;

    }
    for(int n=0;n<=8;n++)
    {
        fctrl=1;
        for(int i=2;i<=(2*n+1);i++)
            fctrl *= i;
       sinus += ((pow(-1,n)*pow(x,2*n+1))/fctrl);
    }
    return (float)sinus;
}//End of sin function


//Start of pow function
float pow(float x,int y)
{
    double pow = 1;
    for(int i=0;i<y;i++)
        pow *= x;
    return (float)pow;
}//End of pow Function

这里有一些例子：

输入

sin(200)

期望输出

-0.873297

我的函数输出

-0.872985

但它适用于较小的值

我用这个：

我现在可以做什么？

Answer 1

您在 sin() 中计算 'x modulus 2 PI' 的方式不会产生准确的结果，例如x = 200 你得到 5.221206 而 fmod() returns 5.221255.

将 x 改为 double 而不是 float 为我解决了这个错误。

这是我测试过的代码。它在没有优化或在 gcc (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0

上使用 -O3 产生了相同的结果

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

#define PI  3.14159265358979323846264338

int main( int argc, char *argv[] )
{
    double v = strtod( argv[1], 0 );
    float  f;
    double d;

    f = v;
    while( f > 6.3 ) {
        f -= 2*PI;
    }

    d = v;
    while( d > 6.3 ) {
        d -= 2*PI;
    }

    printf( "%f\n", f );
    printf( "%lf\n", d );
    printf( "%lf\n", fmod( v, 2*PI ) );
    return( 0 );
}

编辑

正如@Eric Postpischil@ 所指出的那样，这仍然不能解释错误的结果，所以我已经开始使用 sin() 本身，并且似乎增加迭代次数可以提高准确性。使 fctrl 成为 unsigned long long 您可以迭代到 n=10 而不会溢出，对于 x=200 或 x=5.221255 您现在将得到结果 -0.873260 至少更近了。

这是实际代码：

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

#define PI  3.14159265358979323846264338

int iter = 8;
//Start of sin function
float mysin(double x)
{
    double sinus=0;
    unsigned long long int fctrl=1;
    while(!(0<=x&&x<6.3))
    {
        if(x<0)
            x += PI*2;
        else
            x -= PI*2;

    }
    for(int n=0;n<=iter;n++)
    {
        fctrl=1;
        for(int i=2;i<=(2*n+1);i++)
            fctrl *= i;

        double pr = (pow(-1,n)*pow(x,2*n+1));
        printf( "%lf / %llu\n", pr, fctrl );
        sinus += pr / fctrl;
    }
    return (float)sinus;
}//End of sin function

int main( int argc, char *argv[] )
{
    double v = strtod( argv[1], 0 );
    if( argc > 2 ) {
        iter = strtol( argv[2], 0, 10 );
    }
    float  f;
    double d;

    f = v;
    while( f > 6.3 ) {
        f -= 2*PI;
    }

    d = v;
    while( d > 6.3 ) {
        d -= 2*PI;
    }

/*
    printf( "%f\n", f );
    printf( "%lf\n", d );
    printf( "%lf\n", fmod( v, 2*PI ) );
*/

    printf( "%f\n", mysin( v ) );
    return( 0 );
}

编辑 2

另一个改进是避免计算x^(2n+1)和n！首先（这将产生巨大的价值）然后划分但使用这样的循环：

for(int n=0;n<=iter;n++)
{
    double v = n % 2 == 0 ? 1.0 : -1.0;
    for( int j=1; j <= 2*n+1; j++ ) {
            v *= x;
            v /= j;
    }
    sinus += v;
}

n=10，mysin(200) returns -0.873296 现在

Answer 2

有几种方法可以提高发布代码的准确性。

使用 double 作为首选类型，而不是 int 和 float，用于存储中间值、参数和结果的变量。
利用函数的对称性，将计算限制在数值逼近更准确的范围内。
例如，sin 函数可以使用两个函数实现，一个入口点转换作为参数传递的角度并校正另一个函数实现级数计算的结果：

const double PI = 3.14159265358979323846264338; // It will be rounded, anyway.
const double TAU = PI * 2;
const double HALF_PI = PI / 2;
const double THREE_HALF_PI = PI * 1.5;

// This is called only for x in [-Pi/2, Pi/2]
double my_sin_impl(double x);

double my_sin(double x)
{
    // I want to transform all the angles in the range [-pi/2, 3pi/2]
    if ( x < -HALF_PI )
    {
        long long y = (THREE_HALF_PI - x) / TAU;
        x += y * TAU;
    }
    if ( x > THREE_HALF_PI )
    {
        long long y = (x + HALF_PI) / TAU;
        x -= y * TAU;
    }

    if ( x < HALF_PI )
        return my_sin_impl(x);
    else                              //    When x is in [pi/2, 3pi/2],
        return -my_sin_impl(x - PI);  // <- the result should be transformed
}

除非赋值是强制性的，否则避免直接使用 pow 并在每次迭代时计算阶乘。它效率低下，可能导致整数类型溢出。看，例如下面实现不同的方法。
在贴出的代码中，计算sin的迭代次数固定为8，但是我们可以在达到类型的数值精度时停止。在下面的代码片段中，我还将带有交替符号的项的总和拆分为两个部分总和，以（稍微）限制 loss of significance。

double my_sin_impl(double x)
{
    double old_positive_sum, positive_sum = x;
    double old_negative_sum, negative_sum = 0.0;
    double x2 = x * x;
    double term = x;
    int n = 2; 
    do
    {
        term *= x2 / (n * (n + 1));
        old_negative_sum = negative_sum;
        negative_sum += term;
        term *= x2 / ((n + 2) * (n + 3));
        old_positive_sum = positive_sum;
        positive_sum += term;
        n += 4;
    }
    while (positive_sum != old_positive_sum  &&  negative_sum != old_negative_sum);
    //                  ^^ The two values are still numerically different 
    return positive_sum - negative_sum;
}

否则，注意到前面代码片段中的循环最多执行 5 次（给或取），我们可以决定保持滴定次数固定并预先计算系数。像这样（可测试 here）。

double my_sin_impl(double x)
{
    static const double k_neg[] = {1.0/ 6.0, 1.0/42.0, 1.0/110.0, 1.0/210.0, 1.0/342.0};
    static const double k_pos[] = {1.0/20.0, 1.0/72.0, 1.0/156.0, 1.0/272.0, 1.0/420.0};

    double positive_sum = 0.0;
    double negative_sum = 0.0;
    double x2 = x * x;
    double term = 1.0;
    for (int i = 0; i < 5; ++i) 
    {
        term *= x2 * k_neg[i];
        negative_sum += term;
        term *= x2 * k_pos[i];
        positive_sum += term;
    }
    return x * (1.0 + positive_sum - negative_sum);
}

我应该如何处理 return 错误答案的函数？

What should I do with my functions that return wrong answers?

c

trigonometry

exp