WB01432_.gif (3228 バイト)

確率積分ソースプログラム

       

/*===============================================================================================

        erf.c

                                                 確率積分

             積分範囲:0〜5、計算精度:小数点以下5桁まで保証

                                                                                                             1999.9.12

                                                                                                                             by T.Kubota

 ================================================================================================*/

#include <stdio.h>

#include <math.h>

#define    PI         3.141592653589793238462

#define    EPS    (1e-6)

double prob_inte(double);                             /* 確率積分関 */

double fact(long);                                          /* 階乗関数 */

int main(void)

{

    double    a,b;                                             /* 積分区間 */

    double    pa,pb,pab;                                 /* 確率 */

    printf("確率積分\n");

    printf("誤差区間 a〜b(0<a<b)を入力して下さい.--->);

    scanf("%lf %lf",&a,&b);

    if(a==0){

        pa=0;

    }

    else{

        pa=prob_inte(a);

    }

    pb=prob_inte(b);

    pab=pb-pa;

    printf("誤差区間 %g 〜 %g の確率は %7.5lf です.\n",a,b,pab);

    return 0;

}

/*

* 確率積分

*/

double prob_inte(double x)

{

    double    p,pp,pn,q;

    double    a,b;

    int        sign;

    int        i,j,k;

    p=a=x;

    pp=pn=0.0;

    i=1;    j=3;    k=1;

    sign=1;

    do{

        a*=(-x*x);                             /* 分子 */

        b=j*factg(k);                          /* 分母 */

        sign*=(-1);

        if(sign>0)

            pp+=(a/b);                         /* 正項の加算 */

        if(sign<0)

            pn+=(a/b);                         /* 負項の加算 */

        j+=2;

        k++;

    }while(fabs(a/b)>EPS);                 /* EPS:打ち切り誤差 */

    p+=(pp+pn);

    q=2.0*p/sqrt(PI);

    return q;

}

/*

*    階乗計算

*   階乗計算は通常 unsiged long で計算しますが,精度が高々13!までしか保証されませんので

*   double で計算しています.

*/

double fact(long n)

{

    double    f=n;

    int         i;

    inf(n==0) return 1.0;

    for(i=n-1;i>0;i--)

        f*=i;

    return f;

}

このページの先頭へ