laguer.txt

最近开发的,比较使用,曲线拟合,谢谢下载 ,没问题

Function CABS(A1, A2:real):real;
var
    X,Y:real;
begin
    X:=Abs(A1);
    Y:=Abs(A2);
    If X = 0 Then
        CABS:=Y
    Else If Y = 0 Then
        CABS:=X
    Else If X > Y Then
        CABS:=X * Sqrt(1 + Sqrt(Y / X))
    Else
        CABS:=Y * Sqrt(1 + Sqrt(X / Y));
End;

Function CDIV1(A1, A2, B1, B2:real):real;
var
    R,DEN:real;
begin
    If Abs(B1) >= Abs(B2) Then
    begin
        R:=B2 / B1;
        DEN:=B1 + R * B2;
        CDIV1:=(A1 + A2 * R) / DEN;
    end
    Else
    begin
        R:=B1 / B2;
        DEN:=B2 + R * B1;
        CDIV1:=(A1 * R + A2) / DEN;
    end;
end;
Function CDIV2(A1, A2, B1, B2:real):real;
var
    R,DEN:real;
begin
    If Abs(B1) >= Abs(B2) Then
    begin
        R:=B2 / B1;
        DEN:=B1 + R * B2;
        CDIV2:=(A2 - A1 * R) / DEN;
    end
    Else
    begin
        R:=B1 / B2;
        DEN:=B2 + R * B1;
        CDIV2:=(A2 * R - A1) / DEN;
    end;
end;
Function CSQR1(X, Y:real):real;
var
    U,W,V,R:real;  
begin
    If (X = 0) And (Y = 0) Then
        U:=0
    Else
    begin
        If Abs(X) >= Abs(Y) Then
            W:=Sqrt(Abs(X)) * Sqrt(0.5 * (1 + Sqrt(1 + Sqrt(Abs(Y / X)))))
        Else
        begin
            R:=Abs(X / Y);
            W:=Sqrt(Abs(Y)) * Sqrt(0.5 * (R + Sqrt(1 + Sqrt(R))));
        end;
        If X >= 0 Then
        begin
            U:=W;
            V:=Y / (2 * U);
        end
        Else
        begin
            If Y >= 0 Then
                V:=W
            Else
                V:=-W;
            U:=Y / (2 * V);
        end;
    end;
    CSQR1:=U;
end;
Function CSQR2(X, Y:real):real;
var
    V,W,U,R:real;
begin
    If (X = 0) And (Y = 0) Then
        V:=0
    Else
    begin
        If Abs(X) >= Abs(Y) Then
            W:=Sqrt(Abs(X)) * Sqrt(0.5 * (1 + Sqrt(1 + Sqrt(Abs(Y / X)))))
        Else
        begin
            R:=Abs(X / Y);
            W:=Sqrt(Abs(Y)) * Sqrt(0.5 * (R + Sqrt(1 + Sqrt(R))));
        end; 
        If X >= 0 Then
        begin
            U:=W;
            V:=Y / (2 * U);
        end
        Else
        begin
            If Y >= 0 Then
                V:=W
            Else
                V:=-W;
            U:=Y / (2 * V);
        end;
    end;
    CSQR2:=V;
end;

Procedure LAGUER(A:matrx2;M:integer;var X:array of real;
                                         EPS:real;POLISH:boolean);
const
    EPSS = 0.6e-7;    MAXIT = 100;
var
    ZERO,B,D,F,G,H:array[0..2] of real;
    G2,SQ,GP,GM,DX,X1:array[0..2] of real;
    ITER,J:integer;
    DXOLD,ERQ,ABX,DUM,DUM1,DUM2,CDX:real;
begin
    ZERO[1]:=0; 
    ZERO[2]:=0; 
    DXOLD:=CABS(X[1], X[2]);
    For ITER:=1 To MAXIT do
    begin
        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 DownTo 1 do
        begin
            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;
        End; 
        ERQ:=EPSS * ERQ;
        If CABS(B[1], B[2]) <= ERQ Then
            Exit
        Else
        begin
            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]) Then
            begin
                GP[1]:=GM[1];
                GP[2]:=GM[2];
            end;
            DX[1]:=CDIV1(M, 0, GP[1], GP[2]);
            DX[2]:=CDIV2(M, 0, GP[1], GP[2]);
        end;
        X1[1]:=X[1] - DX[1];
        X1[2]:=X[2] - DX[2];
        If (X[1] = X1[1]) And (X[2] = X1[2]) Then Exit;
        X[1]:=X1[1];
        X[2]:=X1[2];
        CDX:=CABS(DX[1], DX[2]);
        DXOLD:=CDX;
        If Not POLISH Then
        begin
            If CDX <= EPS * CABS(X[1], X[2]) Then
            begin
                Exit;
            end;
        end;
    end;
    ShowMessage('too many iterations');
End;