d4r11.cpp

数值计算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;
}

void gser(double& gamser, double& a, double& x, double& gln)
{
	int itmax,n;
	double ap,sum,del,eps;
    itmax = 100;
    eps = 0.0000003;
    gln = gammln(a);
    if (x <= 0.0)
	{
        if (x < 0.0)
		{
            cout<<" pause";
            _c_exit();
		}
        gamser = 0.0;
        _c_exit();
    }
    ap = a;
    sum = 1.0 / a;
    del = sum;
    for( n = 1;n<=itmax;n++)
	{
        ap = ap + 1.0;
        del = del * x / ap;
        sum = sum + del;
        if (fabs(del) < (fabs(sum) * eps))
		{//	goto loop;
		gamser = sum * exp(-x + a * log(x) - gln);
		break;
		}
    }
}

void gcf(double& gammcf, double& a, double& x, double& gln)
{
	int itmax,n;
	double eps,a0,a1,b0,b1,fac,an,ana,anf,gold,g;
    itmax = 100;
    eps = 0.0000003;
    gln = gammln(a);
    gold = 0.0;
    a0 = 1.0;
    a1 = x;
    b0 = 0.0;
    b1 = 1.0;
    fac = 1.0;
    for (n = 1;n<=itmax;n++)
	{
        an = n;
        ana = an - a;
        a0 = (a1 + a0 * ana) * fac;
        b0 = (b1 + b0 * ana) * fac;
        anf = an * fac;
        a1 = x * a0 + anf * a1;
        b1 = x * b0 + anf * b1;
        if (a1 != 0.0)
		{
            fac = 1.0 / a1;
            g = b1 * fac;
            if (fabs((g - gold) / g) < eps )
			{
				gammcf = exp(-x + a * log(x) - gln) * g;
				break;
			}
			//goto yi;
            gold = g;
        }
    }
}

double gammq(double a, double x)
{
	double temp,gamser,gammcf,gln;
    if (x < 0.0 || a <= 0.0)
	{
        cout<< "pause";
        exit(1);
    }
    if (x < a + 1.0)
	{
        gser(gamser, a, x, gln);
        temp= 1.0 - gamser;
	}
    else
	{
        gcf(gammcf, a, x, gln);
        temp = gammcf;
	}
	return temp;   
}

double gammp(double a, double x)
{
	double temp,gamser,gln,gammcf;
    if (x < 0.0 || a <= 0.0)
	{
        cout<<"pause";
        exit(1); 
    }
    if (x < a + 1.0)
	{
        gser(gamser, a, x, gln);
        temp= gamser;
	}
    else
	{
        gcf(gammcf, a, x, gln);
        temp = 1.0 - gammcf;
	}
	return temp;
}

double erfc(double x)
{
	double temp;
    if (x < 0.0)
	{
		temp= 1.0 + gammp(0.5, x*x);
	}
    else
	{
        temp = gammq(0.5, x *x);
	}
	return temp;
}

double erf(double x)
{
	double temp;
    if (x < 0.0)
	{
		temp= -gammp(0.5, x*x);
	}
    else
	{
        temp = gammp(0.5, x *x);
	}
	return temp;
}

void main()
{
    //program d4r11
    //driver for routine erfc
    int i;
	char text[20];
    double nval,actual,x; 
    const double pi = 3.1415926;
    fstream fin;
    fin.open("d:\\vc常用数值算法集\\data\\fncval.dat",ios::in);
    while (strcmp(text,"Error")!=0)
	{
      fin>>text;
	}
	fin>>text;
    fin>>nval; 
    cout<<"Complementary Error Function "<<endl;
    fin>>text;
    cout<<endl;
    cout<<"      x          actual          erfc(x)"<<endl;

    for( i = 1;i<=nval;i++)
	{
        fin>>x;
		fin>>actual;  
		actual = 1 - actual;
		cout<<setw(8)<<x;
		cout<<setw(16)<<actual;
		cout<<setw(16)<<erfc(x)<<endl;			          
    }
}