d4r20.cpp

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

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

double bessi0(double x)
{
	double p1,p2,p3,p4,p5,p6,p7,q1,q2,q3,q4,q5,q6,q7,q8,q9;
	double y,aaa,bbb,ax,temp;
    p1 = 1.0;              p2 = 3.5156229;
    p3 = 3.0899424;        p4 = 1.2067492;
    p5 = 0.2659732;        p6 = 0.0360768;
    p7 = 0.0045813;
    q1 = 0.39894228;       q2 = 0.01328592;
    q3 = 0.00225319;       q4 = -0.00157565;
    q5 = 0.00916281;       q6 = -0.02057706;
    q7 = 0.02635537;       q8 = -0.01647633;
    q9 = 0.00392377;
    if( fabs(x) < 3.75)
	{
        y = (x / 3.75)*(x / 3.75);
        aaa = y * (p5 + y * (p6 + y * p7));
        temp = p1 + y * (p2 + y * (p3 + y * (p4 + aaa)));
	}
    else
	{
        ax = fabs(x);
        y = 3.75 / ax;
        aaa = exp(ax) / sqrt(ax);
        bbb = q4 + y * (q5 + y * (q6 + y * (q7 + y * (q8 + y * q9))));
        temp = aaa * (q1 + y * (q2 + y * (q3 + y * bbb)));
	}
	return temp;
}

void main()
{
    //program d4r20
    //driver for routine bessi0
    int i;
	char text[20];
    double nval,value,x; 
    const double pi = 3.1415926;
    fstream fin;
    fin.open("d:\\vc常用数值算法集\\data\\fncval.dat",ios::in);
    while ( strcmp(text,"I0")!=0 )
	{
      fin>>text;
	}
    fin>>nval; 
    fin>>text;
    cout<<"Modified Bessel Function I0"<<endl;
    cout<<endl;
    cout<<"    x          actual        bessi0(x)"<<endl;
    for( i = 1;i<=nval;i++)
	{
        fin>>x;
		fin>>value;           
		cout<<setprecision(8)<<setw(6)<<x;
        cout<<setprecision(8)<<setw(16)<<value;
		cout<<setprecision(8)<<setw(16)<<bessi0(x)<<endl;			          
    }
}