d4r19.h

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 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 BESSJ(int N, double X)
{
	int IACC;
	double temp,BIGNO,BIGNI,AX,BESSJ,TOX,BJM,BJ,BJP,BJM,BESJ,JSUM,Sum,BJP;
    IACC = 40;
    BIGNO = 10000000000.0;
    BIGNI = 0.0000000001;
    if (N < 2 )
	{
        cout<<"bad argument N in BASSJ";
        _c_exit();
    }
    AX = fabs(X);
    if( AX = 0) 
	{
        temp = 0.0;
	}
    else if (AX > float(N))
	{
        TOX = 2.0 / AX;
        BJM = BESSJ0(AX);
        BJ = BESSJ1(AX);
        for( J = 1;J<=N-1;J++)
		{
            BJP = J * TOX * BJ - BJM;
            BJM = BJ;
            BJ = BJP;
        }
        temp = BJ;
	}
    else
	{
        TOX = 2.0 / AX;
        M = 2 * Int(((N + Int(Sqr(IACC * N)))) / 2);
        BESJ = 0.0
        JSUM = 0;
        Sum = 0.0;
        BJP = 0.0;
        BJ = 1.0;
        for (J = M;J>=1;J--)
		{
            BJM = J * TOX * BJ - BJP;
            BJP = BJ;
            BJ = BJM;
            if (fabs(BJ) > BIGNO)
			{
                BJ = BJ * BIGNI;
                BJP = BJP * BIGNI;
                BESJ = BESJ * BIGNI;
                Sum = Sum * BIGNI;
            }
            if (JSUM != 0) Sum = Sum + BJ;
            JSUM = 1 - JSUM;
            if (J = N) BESJ = BJP;
        }
        Sum = 2.0 * Sum - BJ;
        temp = BESJ / Sum;		
    }
	return temp;
}

void main()
{
    //PROGRAM D4R2
    //Driver for routine FACTRL
    int I,N;
	char Text[20],Text1[20],Text2[20];
    double NVAL,Value,X; 
    const double PI = 3.1415926;
    fstream fin;
    fin.open("e:\\VB常用数值算法集\\DATA\\FNCVAL.DAT",ios::in);
    while ( (strcmp(Text,"Bessel")!=0)||(strcmp(Text1,"Function")!=0)||(strcmp(Text2,"Jn,")!=0) )
	{
      //Text[0]=0; 
      fin>>Text;
	  if(strcmp(Text,"Bessel")==0)
	  {
		  fin>>Text1;
		  fin>>Text2;
	  }
	}
	fin>>Text;
    fin>>NVAL; 
    cout<<Text<<" "<<endl;
    fin>>Text;
    cout<<endl;
    cout<<"   N      X           Actual        BESSJ(N,X)"<<endl;

    for( I = 1;I<=NVAL;I++)
	{
		fin>>N;
        fin>>X;
		fin>>Value;           
		cout<<setw(1)<<N;
		cout<<setw(18)<<X;
        cout<<setw(18)<<Value;
		cout<<setw(18)<<BESSJ(N,X)<<endl;			          
    }
}