0

次のコードを使用して、単純な二分アルゴリズムを使用して関数の根を見つけます

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

typedef float (*continous_function)(float);

static const float epsilon = 0.0000001;

#if __STDC__
static __inline float fabsf(float x)
{
        return x >=0 ? x : -1*x;
}
#else
#include <math.h>
#endif

float poly_1(float x)
{
        return x*x*x+2*x*x+3;
}

float poly_2(float x)
{
        return 2*x + 3;
}

float poly_3(float x)
{
        return 5*x*x*x*x*x+3*x*x*x+2*x+7;
}

static __inline void swap_in_place(float *a, float *b)
{
        float tmp = *a;
        *a = *b;
        *b = tmp;
}

float bisection_root(continous_function f, float a, float b)
{

        float neg = f(a);
        float pos = f(b);
        float c = 0.5 * (a + b);
        float mid;

        if (neg * pos > 0) {
                /* neg and pos should have different sizes */
                abort();
        }

        if (neg > 0) {
                /* Ensure f(a)=neg is negative */
                swap_in_place(&neg, &pos);
                swap_in_place(&a, &b);
        }

        mid = f(c);

        if (fabsf(mid) < epsilon) {
                return c;
        } else {
                if (mid > 0)
                        return bisection_root(f, a, c);
                else
                        return bisection_root(f, c, b);
        }

}

int main()
{
        {
                float a = -5;
                float b = 2;
                printf("Root of x^3+2*x^2+3 is %f\n", bisection_root(poly_1,a,b));
        }

        {
                float a = -2;
                float b = 0;
                printf("Root of 2*x+3 is %f\n", bisection_root(poly_2,a,b));
        }

        {
                float a = -1;
                float b = 0;
                printf("Root of 5*x^5+3*x^3+2*x+7 is %f\n", bisection_root(poly_3,a,b));
        }

        return 0;
}

このプログラムは、mingw gcc を使用して Windows XP (32 ビット) でコンパイルすると、セグメンテーション違反を引き起こします。

小数点以下の桁数がepsilon減ると、セグメンテーション違反を回避できます。したがって、セグメンテーション違反はオーバーフローと関係があると結論付けています。

イプシロンを設定したり、問題を引き起こす可能性のある他のエラーを修正したりするための信頼できる方法を見つけることができるように、このエラーが発生する理由と正確な方法を知りたいです。

4

1 に答える 1

1 
static const float epsilon = 0.0000001;

多項式の値が 0 から離れている浮動小数点数が存在する必要はありませんepsilon。特に、ルートが の大きな値の場合、多項式は1 つの浮動小数点数とその後継の浮動小数点数の間で、より小さい値からより大きい値にxジャンプする可能性があります。 .-epsilon+epsilon

この場合、および他の多くの場合、アルゴリズムはループします。再帰を使用して実装し、C コンパイラは通常末尾呼び出しの最適化を保証しないため、この無限ループはセグメンテーション違反を引き起こす可能性があります (スタックがいっぱいになると)。

再帰を使用するか単純なwhileループを使用するかに関係なく適用できる解決策は、反復回数を制限することです。

注:多項式を計算するためのホーナーのスキームについて読む必要があります。

于 2015-03-24T12:26:59.147 に答える