unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics,UNIT1, Controls, Forms, Dialogs;

Function DBRENT(AX, BX, CX, TOL:real;var XMIN:real):real;
Function BESSJ0(X:real):real;
Function BESSJ1(X:real):real;
Procedure MNBRAK(var AX, BX, CX, FA, FB, FC:real);

implementation
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;


Function BESSJ0(X:real):real;
var
   AAA,BBB,CCC,Y,AX,Z,DDD,EEE,XX:real;
const
  P1=1;                    P2=-0.001098628627;
  P3=0.2734510407e-4;      P4=-0.2073370639e-5;
  P5=2.093887211E-07;
  Q1=-0.1562499995e-1;     Q2=0.1430488765e-3;
  Q3=-0.6911147651e-5;     Q4=7.621095161E-07;
  Q5=-9.34945152E-08;
  R1=57568490574;          R2=-13362590354;
  R3=651619640.7;          R4=-11214424.18;
  R5=77392.33017;          R6=-184.9052456;
  S1=57568490411;          S2=1029532985;
  S3=9494680.718;          S4=59272.64853;
  S5=267.8532712;          S6=1;
begin
  If Abs(X) < 8 Then
    begin
      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)));
      BESSJ0:= AAA / (S1+ Y* (S2+ CCC));
    end
  Else
    begin
      AX:=Abs(X);
      Z:=8/ 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));
      BESSJ0:= Sqrt(0.636619772 / AX) * (Cos(XX) * AAA- EEE);
    End;
End;

Function BESSJ1(X:real):real;
VAR
    AAA,BBB,CCC,AX,Z,Y,XX,SGN:real;
const
    R1 = 72362614232;         R2 = -7895059235;
    R3 = 242396853.1;         R4 = -2972611.439;
    R5 = 15704.4826;          R6 = -30.16036606;
    S1 = 144725228442;        S2 = 2300535178;
    S3 = 18583304.74;         S4 = 99447.43394;
    S5 = 376.9991397;         S6 = 1;
    P1 = 1;                   P2 = 0.00183105;
    P3 = -0.3516396496e-4;    P4 = 0.000002457520174;
    P5 = -0.240337019e-6;
    Q1 = 0.04687499995;       Q2 = -0.2002690873e-3;
    Q3 = 0.8449199096e-5;     Q4 = -0.88228987e-6;
    Q5 = 0.105787412e-6;
begin
    If Abs(X) < 8 Then
      begin
        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))));
        BESSJ1:=X * AAA / BBB;
      end
    Else
      begin
        AX:=Abs(X);
        Z:=8 / 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);
        if X > 0 THEN
          SGN:= 1
        ELSE
          SGN:= -1;
        BESSJ1:=CCC * (Cos(XX) * AAA - Z * Sin(XX) * BBB * Sgn);
      End;
End;

Procedure MNBRAK(var AX, BX, CX, FA, FB, FC:real);
const
    GOLD = 1.618034;  GLIMIT = 100;  TINY = 1E-20;
var 
    DUM,R,Q,U,ULIM,FU:real;   DONE:boolean;
begin
    FA:=FUNC(AX);
    FB:=FUNC(BX);
    If FB > FA Then
    begin
        DUM:=AX;
        AX:=BX;
        BX:=DUM;
        DUM:=FB;
        FB:=FA;
        FA:=DUM;
    end;
    CX:=BX + GOLD * (BX - AX);
    FC:=FUNC(CX);
    repeat
        If FB < FC Then break;
        DONE:=true;
        R:=(BX - AX) * (FB - FC);
        Q:=(BX - CX) * (FB - FA);
        DUM:=Q - R;
        If Abs(DUM) < TINY Then DUM:=TINY;
        U:=BX - ((BX - CX) * Q - (BX - AX) * R) / (2 * DUM);
        ULIM:=BX + GLIMIT * (CX - BX);
        If (BX - U) * (U - CX) > 0  Then
        begin
            FU:=FUNC(U);
            If FU < FC Then
            begin
                AX:=BX;
                FA:=FB;
                BX:=U;
                FB:=FU;
                Exit;
            end
            Else If FU > FB Then
            begin
                CX:=U;
                FC:=FU;
                Exit;
            end;
            U:=CX + GOLD * (CX - BX);
            FU:=FUNC(U);
        end
        Else If (CX - U) * (U - ULIM) > 0  Then
        begin
            FU:=FUNC(U);
            If FU < FC Then
            begin
                BX:=CX;
                CX:=U;
                U:=CX + GOLD * (CX - BX);
                FB:=FC;
                FC:=FU;
                FU:=FUNC(U);
            end;
        end
        Else If (U - ULIM) * (ULIM - CX) >= 0  Then
        begin
            U:=ULIM;
            FU:=FUNC(U);
        end
        Else
        begin
            U:=CX + GOLD * (CX - BX);
            FU:=FUNC(U);
        end;
        If DONE Then
        begin
            AX:=BX;
            BX:=CX;
            CX:=U;
            FA:=FB;
            FB:=FC;
            FC:=FU;
        end
        Else
            DONE:=false;
    until not DONE
end;

end.