d4r13.cpp

Visual C++ 常用数值算法集 源代码

#include <iostream.h>
#include <math.h>
#include <iomanip.h>
#include <stdlib.h>
#include <fstream.h>
#include <string>
#include <process.h>

double gammln(double xx)
{
	int j;
	float temp;
    double cof[6],stp,half,one,fpf,x,tmp,ser;
    cof[1] = 76.18009173;
    cof[2] = -86.50532033;
    cof[3] = 24.01409822;
    cof[4] = -1.231739516;
    cof[5] = 0.00120858003;
    cof[6] = -0.00000536382;
    stp = 2.50662827465;
    half = 0.5;
    one = 1.0;
    fpf = 5.5;
    x = xx - one;
    tmp = x + fpf;
    tmp = (x + half) * log(tmp) - tmp;
    ser = one;
    for (j = 1;j<=6;j++)
	{
        x = x + one;
        ser = ser + cof[j] / x;
    }
    temp = tmp + log(stp * ser);
	return temp;
}

double betacf(double a,double b,double x)
{
	double itmax,eps,am,bm,az,qab,qap,qam,bz,em,tem,ap,bp,aap,bpp,aold,temp,d;
	int m;
    itmax = 100;
    eps = 0.0000003;
    am = 1.0;
    bm = 1.0;
    az = 1.0;
    qab = a + b;
    qap = a + 1.0;
    qam = a - 1.0;
    bz = 1.0 - qab * x / qap;
    for( m = 1;m<=itmax;m++)
	{
        em = m;
        tem = em + em;
        d = em * (b - m) * x / ((qam + tem) * (a + tem));
        ap = az + d * am;
        bp = bz + d * bm;
        d = -(a + em) * (qab + em) * x / ((a + tem) * (qap + tem));
        aap = ap + d * az;
        bpp = bp + d * bz;
        aold = az;
        am = ap / bpp;
        bm = bp / bpp;
        az = aap / bpp;
        bz = 1.0;
        if( fabs(az - aold) < eps * fabs(az)) 
			goto yi;
    }
    cout<<"a or b too big, or itmax too small";
yi: temp = az;
	return temp;
}

double betai(double a,double b,double x)
{
	double bt,aaa,temp;
    if (x < 0.0 || x > 1.0) cout<<"bad argument x in betai";
    if (x == 0.0 || x == 1.0)
	{
        bt = 0.0;
	}
    else
	{
        aaa = gammln(a + b) - gammln(a) - gammln(b);
        bt = exp(aaa + a * log(x) + b * log(1.0 - x));
	}
    if (x < ((a + 1.0) / (a + b + 2.0)))
	{
        temp = bt * betacf(a, b, x) / a;
        _c_exit();
	}
    else
	{
        temp = 1.0 - bt * betacf(b, a, 1.0 - x) / b;
        _c_exit();
	}
	return temp;
}

void main()
{
    //program d4r13
    //driver for routine betai betacf
    int i;
	char text[20],text1[20];
    double nval,a,b,x,value; 
    const double pi = 3.1415926;
    fstream fin;
    fin.open("d:\\vc常用数值算法集\\data\\fncval.dat",ios::in);
    while ( (strcmp(text,"Incomplete")!=0)||(strcmp(text1,"Beta")!=0) )
	{
      fin>>text;
	  if(strcmp(text,"Incomplete")==0)
		  fin>>text1;
	}
	fin>>text;
    fin>>nval; 
    cout<<"Incomplete Beta Function "<<endl;
    fin>>text;
    cout<<endl;
    cout<<"      a       b       x        actual       betai(x)"<<endl;
    for( i = 1;i<=nval;i++)
	{
        fin>>a;
		fin>>b;
		fin>>x;
		fin>>value;
		cout<<setprecision(3)<<setw(8)<<a;
		cout<<setprecision(3)<<setw(8)<<b;
		cout<<setprecision(3)<<setw(8)<<x;
		cout<<setprecision(6)<<setw(14)<<value;
		cout<<setprecision(6)<<setw(14)<<betai(a, b, x)<<endl;			          
    }
}