text1.txt

数值计算c++源代码,包括各种算法。很有用的。

double CABS(double A1, double A2)
{
	double X,Y,t1;
    X = fabs(A1);
    Y = fabs(A2);
    if (X == 0.0)
        t1 = Y;
    else if (Y == 0)
        t1 = X;
    else if (X > Y)
        t1 = X * sqrt(1 + sqrt(Y / X));
    else
        t1 = Y * sqrt(1 + sqrt(X / Y));
	return t1;
    
}

 double CDIV1(double A1,double A2,double B1, double B2)
 {
	 double R,DEN,t2;
    if (fabs(B1) >= fabs(B2)) 
	{
        R = B2 / B1;
        DEN = B1 + R * B2;
        t2 = (A1 + A2 * R) / DEN;
	}
    else
	{
        R = B1 / B2;
        DEN = B2 + R * B1;
        t2 =(A1 * R + A2) / DEN;
    }
	return t2;
 }

double CDIV2(double A1, double A2, double B1,double B2)
{
	double R,DEN,t3;
    if (fabs(B1) >= fabs(B2)) 
	{
        R = B2 / B1;
        DEN = B1 + R * B2;
        t3 = (A2 - A1 * R) / DEN;
	}
    else
	{
        R = B1 / B2;
        DEN = B2 + R * B1;
       t3 = (A2 * R - A1) / DEN;
    }
	return t3;
}

double  CSQR1(double X,double Y)

{
	double tt,U,R,W,V;

    if ((X == 0) &&( Y == 0))
        U = 0.0;
    else
	{
        if (fabs(X) >= fabs(Y))
            W = sqrt(fabs(X)) * sqrt(0.5 * (1 + sqrt(1 + sqrt(fabs(Y / X)))));
        else
		{
            R = fabs(X / Y);
            W = sqrt(fabs(Y)) * sqrt(0.5 * (R + sqrt(1 + sqrt(R))));
		}
        
        if (X >= 0)
		{
            U = W;
            V = Y / (2 * U);
		}
        else
		{
            if( Y >= 0)
                V = W;
            else
                V = -W;
        U = Y / (2 * V);
        }
    }
    tt = U;
	return tt;
}

double  CSQR2(double X,double Y)
{
	double ttt,U,R,W,V;
    if ((X = 0) && (Y = 0))
        V = 0;
    else
	{
        if (fabs(X) >= fabs(Y))
            W = sqrt(fabs(X)) * sqrt(0.5 * (1 + sqrt(1 + sqrt(fabs(Y / X)))));
        else
		{
            R = fabs(X / Y);
            W = sqrt(fabs(Y)) * sqrt(0.5 * (R + sqrt(1 + sqrt(R))));
        }
        if (X >= 0)
		{
            U = W;
            V = Y / (2 * U);
		}
        else
		{
            if (Y >= 0 )
				V = W;
			else
				V = -W;
        
        U = Y / (2 * V);
        }
    }
    ttt = V;
	return ttt;
}

void LAGUER(double A[3][102],int  M, double X[3], double EPS,int & POLISH)
{
    double ZERO[3], B[3], D[3], F[3], G[3], H[3];
    double G2[3], SQ[3], GP[3], GM[3], DX[3], X1[3];
	int ITER,J,MAXIT;
	double EPSS,DXOLD,ERQ,ABX,DUM,DUM1,DUM2,CDX;
    ZERO[1] = 0.0;
    ZERO[2] = 0.0;
    EPSS = 0.00000006;
    MAXIT = 100;
    DXOLD = CABS(X[1], X[2]);
    for (ITER = 1;ITER<=MAXIT;ITER++)
	{
        B[1] = A[1][ M + 1];
        B[2] = A[2][ M + 1];
        ERQ = CABS(B[1], X[2]);
        D[1] = ZERO[1];
        D[2] = ZERO[2];
        F[1] = ZERO[1];
        F[2] = ZERO[2];
        ABX = CABS(X[1], X[2]);
        for (J = M;J>=1;J--)
		{
            DUM = X[1] * F[1] - X[2] * F[2] + D[1];
            F[2] = X[2] * F[1] + X[1] * F[2] + D[2];
            F[1] = DUM;
            DUM = X[1] * D[1] - X[2] * D[2] + B[1];
            D[2] = X[2] * D[1] + X[1] * D[2] + B[2];
            D[1] = DUM;
            DUM = X[1] * B[1] - X[2] * B[2] + A[1][ J];
            B[2] = X[2] * B[1] + X[1] * B[2] + A[2][ J];
            B[1] = DUM;
            ERQ = CABS(B[1], B[2]) + ABX * ERQ;
        }
        ERQ = EPSS * ERQ;
        if (CABS(B[1], B[2]) <= ERQ )
		{
            value( X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO);
            return;
		}
        else
		{
            G[1] = CDIV1(D[1], D[2], B[1], B[2]);
            G[2] = CDIV2(D[1], D[2], B[1], B[2]);
            G2[1] = G[1] * G[1] - G[2] * G[2];
            G2[2] = 2 * G[1] * G[2];
            H[1] = G2[1] - 2 * CDIV1(F[1], F[2], B[1], B[2]);
            H[2] = G2[2] - 2 * CDIV2(F[1], F[2], B[1], B[2]);
            DUM1 =(M - 1) *(M * H[1] - G2[1]);
            DUM2 = (M - 1) * (M * H[2] - G2[2]);
            SQ[1] = CSQR1(DUM1, DUM2);
            SQ[2] = CSQR2(DUM1, DUM2);
            GP[1] = G[1] + SQ[1];
            GP[2] = G[2] + SQ[2];
            GM[1] = G[1] - SQ[1];
            GM[2] = G[2] - SQ[2];
            if (CABS(GP[1], GP[2]) < CABS(GM[1], GM[2]))
			{
                GP[1] = GM[1];
                GP[2] = GM[2];
            }
            DX[1] = CDIV1(M, 0, GP[1], GP[2]);
            DX[2] = CDIV2(M, 0, GP[1], GP[2]);
        }
        X1[1] = X[1] - DX[1];
        X1[2] = X[2] - DX[2];
        if ((X[1] = X1[1]) && (X[2] = X1[2]))
		{
            value(X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO);
            return;
        }
        X[1] = X1[1];
        X[2] = X1[2];
        CDX = CABS(DX[1], DX[2]);
        DXOLD = CDX;
        if (! POLISH)
		{
            if (CDX <= EPS * CABS(X[1], X[2]))
			{
               value( X1, DX, GM, GP, SQ, G2, H, F, D, B, ZERO);
               return;
            }
        }
    }
    cout<< "too many iterations"<<endl;
}