dbrent.txt

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

Function DBRENT(AX, BX, CX, TOL:real;var XMIN:real):real;
const
    ITMAX = 100;  ZEPS = 0.1e-9;
var
    ITER:integer;
    A,B,D,V,W,X,E,FX,FV1,DV,DW,DX,FW,XM:real;
    TOL1,TOL2,D1,D2,U1,U2,U,FU,DU,OLDE,ZZ:real;
    DONE,OK1,OK2:boolean;
begin
    A:=AX;
    If CX < AX Then A:=CX;
    B:=AX;
    If CX > AX Then B:=CX;
    V:=BX;
    W:=V;
    X:=V;
    E:=0; 
    FX:=FUNC(X);
    FV1:=FX;
    FW:=FX;
    DX:=DERIV(X);
    DV:=DX;
    DW:=DX;
    For ITER:=1 To ITMAX do
    begin
        XM:=0.5 * (A + B);
        TOL1:=TOL * Abs(X) + ZEPS;
        TOL2:=2  * TOL1;
        If Abs(X - XM) <= (TOL2 - 0.5 * (B - A)) Then
        begin
          DONE:=true;
          break;
        end
        else
          DONE:=false;
        If Abs(E) > TOL1 Then
        begin
            D1:=2 * (B - A);
            D2:=D1;
            If DW <> DX Then D1:=(W - X) * DX / (DX - DW);
            If DV <> DX Then D2:=(V - X) * DX / (DX - DV);
            U1:=X + D1;
            U2:=X + D2;
            OK1:=((A - U1) * (U1 - B) > 0 ) And (DX * D1 <= 0 );
            OK2:=((A - U2) * (U2 - B) > 0 ) And (DX * D2 <= 0 );
            OLDE:=E;
            E:=D;
            If OK1 Or OK2 Then
               If OK1 And OK2 Then
                   D:=D1
               Else
                   D:=D2
            Else
               If OK1 Then
                   D:=D1
               Else
                   D:=D2;
            If Abs(D) > Abs(0.5 * OLDE) Then
            begin
              U:=X + D;
              If (U - A < TOL2) Or (B - U < TOL2) Then
              Begin
                  If XM-X>=0 then
                      ZZ:=1
                  else
                      ZZ:=-1; 
                 D:=Abs(TOL1) * ZZ;
              end;
            end;
        end;
        If DX >= 0  Then
            E:=A - X
        Else
            E:=B - X;
        D:=0.5 * E;
        If Abs(D) >= TOL1 Then
        begin
            U:=X + D;
            FU:=FUNC(U);
        end
        Else
        begin
            if d >=0 then
                ZZ:=1
            else
                ZZ:=-1;
            U:=X + Abs(TOL1) * ZZ;
            FU:=FUNC(U);
            If FU > FX Then
            begin
              DONE:=true;
              break;
            end
            Else
              DONE:=false;
        end;
        DU:=DERIV(U);
        If FU <= FX Then
        begin
            If U >= X Then
                A:=X
            Else
                B:=X;
            V:=W;
            FV1:=FW;
            DV:=DW;
            W:=X;
            FW:=FX;
            DW:=DX;
            X:=U;
            FX:=FU;
            DX:=DU;
        end
        Else
        begin
            If U < X Then
                A:=U
            Else
                B:=U;
            If (FU <= FW) Or (W = X) Then
            begin
                V:=W;
                FV1:=FW;
                DV:=DW;
                W:=U;
                FW:=FU;
                DW:=DU;
            end
            Else If (FU <= FV1) Or (V = X) Or (V = W) Then
            begin
                V:=U;
                FV1:=FU;
                DV:=DU;
            end;
        end;
    end; 
    If Not DONE Then
        ShowMessage('DBRENT exceeded maximum iterations.')
    Else
    begin
        XMIN:=X;
        DBRENT:=FX;
    end;
end;