d4r26.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;
}

double t(double y)
{
	double aaa,bbb,z,p,temp;
    if (y <= 10.0)
	{
        aaa = 0.000057941 * y - 0.00176148 * y + 0.0208645;
        bbb = 1;
        temp = ((aaa * y - 0.129013) * y + 0.85777) * y + 1.0125;
	}
    else
	{
        z = log(y) - 0.775;
        p = (0.775 - log(z)) / (1.0 + z);
        temp = y / (p + 1.0) / z;
    }
	return temp;
}

void jap(double &x, double &a, double &nmax, double f[])
{
	int n,m;
    double fa[10], rr[10],mmax,e,eps,sum,d1,r,s,nu,al,li,am;
    mmax = 10;
    e = 4;
    eps = 0.5 * pow(10 , (-e));    
    for (n = 0;n<=nmax;n++)
	{
        fa[n] = 0.0;
    }
    sum = exp(gammln(1.0 + a));
    sum = exp(a * log(x / 2.0)) / sum;
    d1 = 2.3026 * e + 1.3863;
    if( nmax > 0 )
	{
        r = t(0.5 * d1 / nmax) * nmax;
	}
    else
	{
        r = 0.0;
    }
    s = t(0.73576 * d1 / x) * 1.3591 * x;
    if (r <= s)
	{
        nu = 1 + int(s);
	}
    else
	{
        nu = 1 + int(r);
	}
yi: m = 0;
    al = 1.0;
    li = int(nu / 2);
er: m = m + 1;
    al = al * (m + a) / (m + 1);
    if (m < li) 
		goto er;
    n = 2 * m;
    r = 0.0;
    s = 0.0;
san: r = 1.0 / (2.0 * (a + n) / x - r);
    if (int(n / 2) != n / 2.0)
	{
        am = 0.0;
	}
    else
	{
        al = al * (n + 2) / (n + 2.0 * a);
        am = al * (n + a);
	}
    s = r * (am + s);
    if (n <= nmax)
		rr[n - 1] = r;
    n = n - 1;
    if (n >= 1) 
		goto san;
    f[0] = sum / (1.0 + s);
    for (n = 0;n<=nmax - 1;n++)
        f[n + 1] = rr[n] * f[n];
    for (n = 0 ;n<=nmax;n++)
	{
        if (fabs((f[n] - fa[n]) / f[n]) > eps)
		{
            for( m = 0;m<=nmax;m++)
                fa[m] = f[m];            
            nu = nu + 5;
            goto yi;
        }
    }
}

void besjan(double x,double a,double nm, double ih,double f[])
{
	int i;
	double temp;
	temp=1;
    if (ih == 1)
	{
        jap(x, a, nm, f);
        _c_exit(); 
	}
    else
	{
        jap(x, a, temp, f);
        f[1] = 2.0 * a * f[0] / x - f[1];
        for (i = 1;i<=nm;i++)
            f[i + 1] = 2.0 * (a - i) * f[i] / x - f[i - 1];
        _c_exit();
	}
}

void main()
{
    //program d4r26
    //driver for routine besjan
    int i,n[19],ih;
	char text[20];
    double nval,a[19],value[19],x[19],f[19],temp; 
    const double pi = 3.1415926;
	temp=4;
    fstream fin;
    fin.open("d:\\vc常用数值算法集\\data\\fncval.dat",ios::in);
    while ( strcmp(text,"Ja+n")!=0) 
	{
      fin>>text;
	}
    fin>>nval; 
    fin>>text;
	cout<<"Bessel Function Ja+n"<<endl;
    cout<<endl;
    cout<<"     n    a     ih    x          actual         besjan"<<endl;
    for (i = 0 ; i<=nval - 1; i++)
	{
        fin>>n[i];
        fin>>a[i];
        fin>>x[i];
        fin>>value[i];
    }
    ih = 1;
    besjan(x[0], a[0], temp, ih, f);
    for( i = 0;i<=4;i++)
	{
		cout<<setw(6)<<n[i];
        cout<<setw(6)<<a[i];
		cout<<setw(6)<<ih;
		cout<<setw(6)<<x[i];
		cout<<setw(16)<<value[i];
		cout<<setw(16)<<f[i]<<endl;
	}
    besjan(x[5], a[5], temp, ih, f);
    for( i = 5;i<=9;i++)
	{
		cout<<setw(6)<<n[i];
        cout<<setw(6)<<a[i];
		cout<<setw(6)<<ih;
		cout<<setw(6)<<x[i];
		cout<<setw(16)<<value[i];
		cout<<setw(16)<<f[i-5]<<endl;
    }
    ih = -1;
    besjan(x[10], a[10], temp, ih, f);
    for( i = 10;i<=13;i++)
	{
		cout<<setw(6)<<n[i];
        cout<<setw(6)<<a[i];
		cout<<setw(6)<<ih;
		cout<<setw(6)<<x[i];
		cout<<setw(16)<<value[i];
		cout<<setw(16)<<f[i-9]<<endl;
    }
    besjan(x[14], a[14], temp, ih, f);
    for( i = 14;i<=17;i++)
	{
		cout<<setw(6)<<n[i];
        cout<<setw(6)<<a[i];
		cout<<setw(6)<<ih;
		cout<<setw(6)<<x[i];
		cout<<setw(16)<<value[i];
		cout<<setw(16)<<f[i-13]<<endl;
    }
}