/*===============================================================================================
erf.c
確率積分
積分範囲：０〜５、計算精度：小数点以下５桁まで保証
１９９９．９．１２
ｂｙ Ｔ．Ｋｕｂｏｔａ
================================================================================================*/
#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("誤差区間 ａ〜ｂ（０＜ａ＜ｂ）を入力して下さい．--->);
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 で計算しますが，精度が高々１３！までしか保証されませんので
* 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;
}