d4r30.cpp

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

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

double bs(int n, double a)
{
	int m;
    double b[20],temp,del,r,eps,s,ak,a2,omeg,da,ajp,ajm,q1,q2;
    if (a == 0.0)
	{
       if ((n / 2) == (n / 2.0))
	   {
         b[n] = 0.0;
         goto yw;
       }
       b[n] = 2.0 / n;
       goto yw;
    }
    if (fabs(a) < 8.0)
	{
       del = 0.00000001;
       if ((n / 2) != (n / 2.0))
	   {
         r = 2.0 / n;
         eps = r * del;
         s = r;
         ak = 0.0;
         a2 = a * a;
w:        ak = ak + 2.0;
         r = r * a2 * (n + ak - 2.0) / (ak * (ak - 1.0) * (n + ak));
         s = s + r;
         if (r > eps) 
			 goto w;
         b[n] = s + r;
         goto yw;
       }
       r = 2.0 * a / (n + 1.0);
       omeg = fabs(r * del);
       s = r;
       ak = 1.0;
       a2 = a * a;
l:      ak = ak + 2.0;
       r = r * a2 * (n + ak - 2.0) / (ak * (ak - 1.0) * (n + ak));
       s = s + r;
       if (fabs(r) > omeg)
		   goto l;
       b[n] = -(s + r);
       goto yw;
    }
    da = 1.0 / a;
    ajp = da * exp(a);
    ajm = (da * da) / ajp;
    b[1] = ajp - ajm;
    if (n == 1)
		goto yw;
    q1 = -1.0;
    q2 = 1.0;
    for (m = 2;m<=n;m++)
	{
        b[m] = q1 * ajp - ajm + q2 * da * b[m - 1];
        q1 = -q1;
        q2 = q2 + 1.0;
    }
yw:  temp = b[n];
	 return temp;
}

void main()
{
    //program d4r30
    //driver for routine bs
    int i;
	char text[20];
    double nval,value,b,x; 
    const double pi = 3.1415926;
    fstream fin;
    fin.open("d:\\vc常用数值算法集\\data\\fncval.dat",ios::in);
    while ( strcmp(text,"Bs")!=0 )
	{
      fin>>text;
	}
    fin>>nval; 
    cout<<"Exponential Integral Bs "<<endl;
    fin>>text;
    cout<<endl;
    cout<<"     n      b            actual          bs(n,b)"<<endl;
    for( i = 1;i<=nval;i++)
	{
        fin>>x;
		fin>>b;
		fin>>value;           
		cout<<setw(6)<<x;
        cout<<setw(8)<<b;
        cout<<setw(18)<<value;
		cout<<setw(18)<<bs(x,b)<<endl;			          
    }
}