d4r18.cpp

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

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

int sgn(double x)
{
	int temp;
	if(x>0)temp=1;
	if(x=0)temp=0;
	if(x<0)temp=-1;
	return temp;
}

double bessj1(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double bbb,ccc,aaa,temp,ax,xx,z,y;
    r1 = 72362614232.0;           r2 = -7895059235.0;
    r3 = 242396853.1;             r4 = -2972611.439;
    r5 = 15704.4826;              r6 = -30.16036606;
    s1 = 144725228442.0;          s2 = 2300535178.0;
    s3 = 18583304.74;             s4 = 99447.43394;
    s5 = 376.9991397;             s6 = 1.0;
    p1 = 1.0;                     p2 = 0.00183105;
    p3 = -0.00003516396496;       p4 = 0.000002457520174;
    p5 = -0.000000240337019;
    q1 = 0.04687499995;           q2 = -0.0002002690873;
    q3 = 0.000008449199096;       q4 = -0.00000088228987;
    q5 = 0.000000105787412;
    if (fabs(x) < 8.0 )
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
        temp = x * aaa / bbb;
	}
    else
	{
        ax = fabs(x);
        z = 8.0 / ax;
        y = z*z;;
        xx = ax - 2.356194491;
        aaa = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        bbb = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        ccc = sqrt(0.636619772 / ax);
        temp = ccc * (cos(xx) * aaa - z * sin(xx) * bbb * sgn(x));
	}  
	return temp;
}

double bessy1(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6,s7;
	double bbb,ccc,aaa,temp,xx,z,y;
    p1 = 1.0;                     p2 = 0.00183105;
    p3 = -0.00003516396496;       p4 = 0.000002457520174;
    p5 = -0.000000240337019;
    q1 = 0.04687499995;           q2 = -0.0002002690873;
    q3 = 0.000008449199096;       q4 = -0.00000088228987;
    q5 = 0.000000105787412;
    r1 = -4900604943000.0;        r2 = 1275274390000.0;
    r3 = -51534381390.0;          r4 = 734926455.1;
    r5 = -4237922.726;            r6 = 8511.937935;
    s1 = 24995805700000.0;        s2 = 424441966400.0;
    s3 = 3733650367.0;            s4 = 22459040.02;
    s5 = 102042.605;              s6 = 354.9632885;
    s7 = 1.0;
    if (x < 8.0)
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s4 + y * (s5 + y * (s6 + y * s7));
        bbb = s1 + y * (s2 + y * (s3 + y * bbb));
        ccc = bessj1(x) * log(x) - 1.0 / x;
        temp = x * aaa / bbb + 0.636619772 * ccc;
	}
    else
	{
        z = 8.0 / x;
        y = z*z;
        xx = x - 2.356194491;
        aaa = sqrt(0.636619772 / x);
        bbb = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        ccc = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        temp = aaa * (sin(xx) * bbb + z * cos(xx) * ccc);
	}
	return temp;
}

double bessj0(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double y,bbb,ccc,aaa,temp,eee,ddd,ax,xx,z;
    p1 = 1.0;                  p2 = -0.001098628627;
    p3 = 0.00002734510407;     p4 = -0.000002073370639;
    p5 = 2.093887211e-07;
    q1 = -0.01562499995;       q2 = 0.0001430488765;
    q3 = -0.000006911147651;   q4 = 7.621095161e-07;
    q5 = -9.34945152e-08;
    r1 = 57568490574.0;        r2 = -13362590354.0;
    r3 = 651619640.7;          r4 = -11214424.18;
    r5 = 77392.33017;          r6 = -184.9052456;
    s1 = 57568490411.0;        s2 = 1029532985.0;
    s3 = 9494680.718;          s4 = 59272.64853;
    s5 = 267.8532712;          s6 = 1.0;
    if (fabs(x) < 8.0)
	{
       y = x * x;
       bbb = y * (r4 + y * (r5 + y * r6));
       aaa = r1 + y * (r2 + y * (r3 + bbb));
       ccc = y * (s3 + y * (s4 + y * (s5 + y * s6)));
       temp = aaa / (s1 + y * (s2 + ccc));
	}
    else
	{
       ax = fabs(x);
       z = 8.0 / ax;
       y = z * z;
       xx = ax - 0.785398164;
       ccc = y * (p3 + y * (p4 + y * p5));
       aaa = p1 + y * (p2 + ccc);
       ddd = y * (q3 + y * (q4 + y * q5));
       eee = z * sin(xx) * (q1 + y * (q2 + ddd));
       temp = sqrt(0.636619772 / ax) * (cos(xx) * aaa - eee);
	}
	return temp;
}

double bessy0(double x)
{
	double p1,p2,p3,p4,p5,q1,q2,q3,q4,q5;
	double r1,r2,r3,r4,r5,r6,s1,s2,s3,s4,s5,s6;
	double y,bbb,ccc,aaa,temp,xx,z;
    p1 = 1.0;                  p2 = -0.001098628627;
    p3 = 0.00002734510407;     p4 = -0.000002073370639;
    p5 = 2.093887211e-07;
    q1 = -0.01562499995;       q2 = 0.0001430488765;
    q3 = -0.000006911147651;   q4 = 7.621095161e-07;
    q5 = -9.34945152e-08;
    r1 = -2957821389.0;        r2 = 7062834065.0;
    r3 = -512359803.6;         r4 = 10879881.29;
    r5 = -86327.92757;         r6 = 228.4622733;
    s1 = 40076544269.0;        s2 = 745249964.8;
    s3 = 7189466.438;          s4 = 47447.2647;
    s5 = 226.1030244;          s6 = 1.0;
    if (x < 8.0 )
	{
        y = x*x;
        aaa = r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6))));
        bbb = s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))));
        temp = aaa / bbb + 0.636619772 * bessj0(x) * log(x);
	}
    else
	{
        z = 8.0 / x;
        y = z*z;
        xx = x - 0.785398164;
        aaa = sqrt(0.636619772 / x);
        bbb = p1 + y * (p2 + y * (p3 + y * (p4 + y * p5)));
        ccc = q1 + y * (q2 + y * (q3 + y * (q4 + y * q5)));
        temp = aaa * (sin(xx) * bbb + z * cos(xx) * ccc);
	}
	return temp;
}

double bessy(int n, double x)
{
	double tox,by,bym,byp,temp;
	int j;
    if (n < 2)
	{
        cout<<"bad argument n in bessy";
        _c_exit();
	}
    tox = 2.0 / x;
    by = bessy1(x);
    bym = bessy0(x);
    for( j = 1;j<=n-1;j++)
	{
        byp = j * tox * by - bym;
        bym = by;
        by = byp;
    }
    temp = by;
	return temp;
}

void main()
{
    //program d4r18
    //driver for routine bessy
    int i,n;
	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,"Yn,")!=0) 
	{
      fin>>text;
	}
	fin>>text;
    fin>>nval; 
    cout<<"Bessel Function Yn "<<endl;
    fin>>text;
    cout<<endl;
    cout<<"     n       x          actual          bessy(n,x)"<<endl;

    for( i = 1;i<=nval;i++)
	{
		fin>>n;
        fin>>x;
		fin>>value;           
		cout<<setw(6)<<n;
		cout<<setw(8)<<x;
        cout<<setw(18)<<value;
		cout<<setw(18)<<bessy(n,x)<<endl;			          
    }
}