unit2.pas

用于开发税务票据管理的软件

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls,unit1, Forms, Dialogs;
Function GASDEV:real;
Procedure SVDVAR(V:matrx2; MA, NP:integer; W:array of real;
                  var CVM:matrx2; NCVM:integer);
Procedure SVDFIT(X, Y, SIG:array of real; NDATA:integer;
    var A:array of real; MA:integer;var U, V:matrx2;var W:array of real;
     MP, NP:integer;var CHISQ:real; FUNCS:string);

implementation
Function GASDEV:real;
var
   V1,V2,FAC,R:real;
begin
    If ISET^= 0 Then
    begin
      repeat
        V1:=2 * Random - 1;
        V2:=2 * Random - 1;
        R:=Sqr(V1) + Sqr(V2);
      until (R < 1);
      FAC:=Sqrt(-2 * Ln(R) / R);
      GSET^:=V1 * FAC;
      GASDEV:=V2 * FAC;
      ISET^:=1;
    end
    Else
    begin
      GASDEV:=GSET^;
      ISET^:=0;
    end;
end;

procedure SVDCMP(A:matrx2; M, N:integer;var W:array of real; V:matrx2);
label 1,2,3;
var
    RV1:array [0..100] of real;
    I,J,L,ITS,JJ,NM,K:integer;
    G,S,SCALE1,C,F,H,X,Y,Z,Sgn,ANORM,AAAAA:real;
begin
    If M < N Then
         ShowMessage('You must augment A with extra zero rows.');
    G:=0;
    SCALE1:=0;
    ANORM:=0;
    For I:=1 To N do
    begin
        L:=I + 1;
        RV1[I]:=SCALE1 * G;
        G:=0;
        S:=0;
        SCALE1:=0;
        If I <= M Then
        begin
            For K:=I To M do
                SCALE1:=SCALE1 + Abs(A[K, I]);
            If SCALE1 <> 0 Then
            begin
                For K:=I To M do
                begin
                    A[K, I]:=A[K, I] / SCALE1;
                    S:=S + A[K, I] * A[K, I];
                end;
                F:=A[I, I];
                If F > 0 Then
                    SGN:=1
                else
                    SGN:=-1;
                G:=-Sqrt(S) * Sgn;
                H:=F * G - S;
                A[I, I]:=F - G;
                If I <> N Then
                begin
                    For J:=L To N DO
                    begin
                        S:=0;
                        For K:=I To M do
                            S:=S + A[K, I] * A[K, J];
                        F:=S / H;
                        For K:=I To M do
                            A[K, J]:=A[K, J] + F * A[K, I];
                    end;
                end;
                For K:=I To M do
                    A[K, I]:=SCALE1 * A[K, I];
            end;
        end;
        W[I]:=SCALE1 * G;
        G:=0;
        S:=0;
        SCALE1:=0;
        If (I <= M) And (I <> N) Then
        begin
            For K:=L To N do
                SCALE1:=SCALE1 + Abs(A[I, K]);
            If SCALE1 <> 0 Then
            begin
                For K:=L To N do
                begin
                    A[I, K]:=A[I, K] / SCALE1;
                    S:=S + A[I, K] * A[I, K];
                end;
                F:=A[I, L];
                If F > 0 Then
                    SGN:=1
                else
                    SGN:=-1;
                G:=-Sqrt(S) * Sgn;
                H:=F * G - S;
                A[I, L]:=F - G;
                For K:=L To N do
                    RV1[K]:=A[I, K] / H;
                If I <> M Then
                begin
                    For J:=L To M do
                    begin
                        S:=0;
                        For K:=L To N do
                            S:=S + A[J, K] * A[I, K];
                        For K:=L To N do
                            A[J, K]:=A[J, K] + S * RV1[K];
                    end;
                end;
                For K:=L To N do
                    A[I, K]:=SCALE1 * A[I, K];
            end;
        end;
        If ANORM > Abs(W[I]) + Abs(RV1[I]) Then
            ANORM:=ANORM
        Else
            ANORM:=Abs(W[I]) + Abs(RV1[I]);
    end;
    For I:=N DownTo 1 do
    begin
        If I < N Then
        begin
            If G <> 0 Then
            begin
                For J:=L To N do
                    V[J, I]:=(A[I, J] / A[I, L]) / G;
                For J:=L To N do
                begin
                    S:=0;
                    For K:=L To N do
                        S:=S + A[I, K] * V[K, J];
                    For K:=L To N do
                        V[K, J]:=V[K, J] + S * V[K, I];
                end;
            end;
            For J:=L To N do
            begin
                V[I, J]:=0;
                V[J, I]:=0;
            end;
        end;
        V[I, I]:=1;
        G:=RV1[I];
        L:=I;
    end;
    For I:=N DownTo 1 do
    begin
        L:=I + 1;
        G:=W[I];
        If I < N Then
        begin
            For J:=L To N do
                A[I, J]:=0;
        end;
        If G <> 0 Then
        begin
            G:=1 / G;
            If I <> N Then
            begin
                For J:=L To N do
                begin
                    S:=0;
                    For K:=L To M do
                        S:=S + A[K, I] * A[K, J];
                    F:=(S / A[I, I]) * G;
                    For K:=I To M do
                        A[K, J]:=A[K, J] + F * A[K, I];
                end;
            end;
            For J:=I To M do
            begin
                A[J, I]:=A[J, I] * G;
            end;
        end
        Else
            For J:=I To M do
                A[J, I]:=0;
        A[I, I]:=A[I, I] + 1;
    end;
    For K:=N DownTo 1 do
    begin
        For ITS:=1 To 30 do
        begin
            For L:=K DownTo 1 do
            begin
                NM:=L - 1;
                If Abs(RV1[L]) + ANORM = ANORM Then GoTo 2;
                If Abs(W[NM]) + ANORM = ANORM Then GoTo 1;
            end;
1:          C:=0;
            S:=1;
            For I:=L To K do
            begin
                F:=S * RV1[I];
                If Abs(F) + ANORM <> ANORM Then
                begin
                    G:=W[I];
                    H:=Sqrt(F * F + G * G);
                    W[I]:=H;
                    H:=1 / H;
                    C:=(G * H);
                    S:=-(F * H);
                    For J:=1 To M do
                    begin
                        Y:=A[J, NM];
                        Z:=A[J, I];
                        A[J, NM]:=(Y * C) + (Z * S);
                        A[J, I]:=-(Y * S) + (Z * C);
                    end;
                end;
            end;
2:          Z:=W[K];
            If L = K Then
            begin
                If Z < 0 Then
                begin
                    W[K]:=-Z;
                    For J:=1 To N do
                        V[J, K]:=-V[J, K];
                end;
                GoTo 3;
            end;
            If ITS = 30 Then
                ShowMessage('No convergence in 30 iterations');
            X:=W[L];
            NM:=K - 1;
            Y:=W[NM];
            G:=RV1[NM];
            H:=RV1[K];
            F:=((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2 * H * Y);
            G:=Sqrt(F * F + 1);
            If F > 0 Then
                    Sgn:=1
                else
                    Sgn:=-1;
            F:=((X - Z) * (X + Z) + H * ((Y / (F + Abs(G) * Sgn)) - H)) / X;
            C:=1;
            S:=1;
            For J:=L To NM do
            begin
                I:=J + 1;
                G:=RV1[I];
                Y:=W[I];
                H:=S * G;
                G:=G * C;
                Z:=Sqrt(F * F + H * H);
                RV1[J]:=Z;
                C:=F / Z;
                S:=H / Z;
                F:=(X * C) + (G * S);
                G:=-(X * S) + (G * C);
                H:=Y * S;
                Y:=Y * C;
                For NM:=1 To N do
                begin
                    X:=V[NM, J];
                    Z:=V[NM, I];
                    V[NM, J]:=(X * C) + (Z * S);
                    V[NM, I]:=-(X * S) + (Z * C);
                end;
                Z:=Sqrt(F * F + H * H);
                W[J]:=Z;
                If Z <> 0 Then
                begin
                    Z:=1 / Z;
                    C:=F * Z;
                    S:=H * Z;
                end;
                F:=(C * G) + (S * Y);
                X:=-(S * G) + (C * Y);
                For NM:=1 To M do
                begin
                    Y:=A[NM, J];
                    Z:=A[NM, I];
                    A[NM, J]:=(Y * C) + (Z * S);
                    A[NM, I]:=-(Y * S) + (Z * C);
                end;
            end;
            RV1[L]:=0;
            RV1[K]:=F;
            W[K]:=X;
        end;
3:  AAAAA:=1;
    end;
end;

procedure SVBKSB(U:matrx2; W:array of real; V:matrx2;
                   M, N:integer; B:array of real;var X:array of reaL);
var
    TMP:array[0..100] of real;
    I,J,JJ:integer;
    S:real;
begin
    For J:=1 To N do
    begin
        S:=0;
        If W[J] <> 0 Then
        begin
            For I:=1 To M do
                S:=S + U[I, J] * B[I];
            S:=S / W[J];
        end;
        TMP[J]:=S;
    end;
    For J:=1 To N do
    begin
        S:=0;
        For JJ:=1 To N do
            S:=S + V[J, JJ] * TMP[JJ];
        X[J]:=S;
    end;
end;

Procedure FPOLY(X:real; var P:array of real; NP:integer);
var
    J:integer;
begin
    P[1]:=1;
    For J:=2 To NP do
      P[J]:=P[J - 1] * X;
end;

Procedure FLEG(X:real; var PL:array of real; NL:integer);
var
    J:integer;  TWOX,F1,F2,D:real;
begin
    PL[1]:=1;
    PL[2]:=X;
    If NL > 2 Then
    begin
        TWOX:=2 * X;
        F2:=X;
        D:=1; 
        For J:=3 To NL do
        begin
            F1:=D;
            F2:=F2 + TWOX;
            D:=D + 1;
            PL[J]:=(F2 * PL[J - 1] - F1 * PL[J - 2]) / D;
        end;
    end;
end;

Procedure SVDVAR(V:matrx2; MA, NP:integer; W:array of real;
                  var CVM:matrx2; NCVM:integer);
var
    WTI:array[0..20] of real;
    I,K,J:integer;  SUM1:real;
begin
    For I:=1 To MA do
    begin
        WTI[I]:=0 ;
        If W[I] <> 0  Then WTI[I]:=1 / (W[I] * W[I]);
    end; 
    For I:=1 To MA do
    begin
        For J:=1 To I do
        begin
            SUM1:=0;
            For K:=1 To MA do
                SUM1:=SUM1 + V[I, K] * V[J, K] * WTI[K];
            CVM[I, J]:=SUM1;
            CVM[J, I]:=SUM1;
        end; 
    end; 
end;


Procedure SVDFIT(X, Y, SIG:array of real; NDATA:integer;
    var A:array of real; MA:integer;var U, V:matrx2;var W:array of real;
     MP, NP:integer;var CHISQ:real; FUNCS:string);
var
    B:array[0..1000] of real; AFUNC:array[0..50] of real;
    I,J:integer;   TMP,WMAX,THRESH,SUM1:real;
const
    TOL=0.00001;
begin
    For I:=1 To NDATA do
    begin
        If FUNCS = 'FPOLY' Then FPOLY(X[I], AFUNC, MA);
        If FUNCS = 'FLEG'  Then FLEG(X[I], AFUNC, MA);
        TMP:=1 / SIG[I];
        For J:=1 To MA do
            U[I, J]:=AFUNC[J] * TMP;
        B[I]:=Y[I] * TMP;
    end; 
    SVDCMP(U, NDATA, MA, W, V);
    WMAX:=0; 
    For J:=1 To MA do
        If W[J] > WMAX Then WMAX:=W[J];
    THRESH:=TOL * WMAX;
    For J:=1 To MA do
        If W[J] < THRESH Then W[J]:=0; 
    SVBKSB(U, W, V, NDATA, MA, B, A);
    CHISQ:=0; 
    For I:=1 To NDATA do
    begin
        If FUNCS = 'FPOLY' Then FPOLY(X[I], AFUNC, MA);
        If FUNCS = 'FLEG' Then FLEG(X[I], AFUNC, MA);
        SUM1:=0; 
        For J:= 1 To MA do
            SUM1:=SUM1 + A[J] * AFUNC[J];
        CHISQ:=CHISQ + Sqr((Y[I] - SUM1) / SIG[I]);
    end; 
end;
end.