Модуль реализации матричных вычислений для массивов больших размеров

Материал из DRKB



В этом модуле «осели» все операции с матрицами и векторами, которые я использовал для работы. Но есть алгоритмы, которые многие, наверняка, увидят впервые: Divide - алгоритм прямого деления, MSqrt - квадратный корень, MAbs - абсолютная величина. Поскольку модуль содержит все, от элементарных операций до матричных, разобраться будет несложно:

Например, решение системы ЛУ (консольное приложение)

var
  N: Integer;
  A: Matrix;
  b, x: Vector;
begin
  N := {..};
  A.Init(N, N);
  b.Init(N);
  x.Init(N); // или x.Init(B); или x.InitRow(A);
  {...}
  { формирование A и b }
  {...}
  x.Divide( b, A );
  x.Print;
  {...}
end.


Некоторые алгоритмы требуют пояснения, например:

Matrix.E(i, j: LongWord)

или

Vector.E(i: Integer): RealPtr

(RealPtr = ^Real) функция для вычисления адреса элемента матрицы/вектора. Перешла из ДОС когда в модуле использовался алгоритм управления виртуальной памятью для больших размерностей.

Matrix.Multiple(X, Y: Vector)

Результатом, которого является произведение вектора X на транспонированный вектор Y - матрица ранга 1.


Matrix.Invert(A: Matrix)

- если A[N,M], и N <> M то результат - матрица размера [M,N] - псевдообратная = A+.


Matrix.Addition(A: Matrix; B: Real)

- добавление числа в главную диагональ.


Matrix.Diag(r: Real)

- присваивание значения главной диагонали.

HMatrixW.pas
unit HMatrixW;

{
        Матричная алгебра 1996-1998,
        HNV 1996 .
}
interface

uses
  SysUtils,
  HNVString,
  Math,
  // VHeapLow,
  Windows,
  Dialogs;

const
  MaxRow = 214748364;
  { Максимальное к-во элементов в векторе }

  NumbSpline: longword = 1;
  Indicator: word = $118F;
  Condition: Boolean = True;
  InitCount: Integer = 0;
  InitSize: Integer = 0;

var
  NIteration: longword; { Число итераций процесса }

  Temp: Real;
  MathEps: Real; { Относительная машинная точность         }
  SqrtEps: Real; { и корень из нее                         }

  Eps: Real;
  MinEps: Real;

type

  VirtualPtr = record
    Ptr: Pointer;
    Size: longword;
  end;

  RealPtr = ^Real;

  ArrayType =
    array[1..1] of Real;

  ArrayPtr = ^ArrayType;

  VectorPrim = { Вектор }  object
    VAddr: VirtualPtr; { Виртуальный адрес  (tail of DOS)  }
    FInit: word; { Флаг инициализации (tail of DOS)  }
    n: longword; { Размерность вектора  }
    function Addr: ArrayPtr; { Адрес начала вектора }
    function E(I: longword): RealPtr; { Адрес элемента }
    {                                       _ }
    procedure Copy(X: VectorPrim); {   = X }
    {                                          _   _ }
    procedure Substrac(X, Y: VectorPrim); { X - Y }
    {                                          _   _ }
    procedure Addition(X, Y: VectorPrim); { X + Y }
    procedure Init(Size: longword); overload;
    procedure Init(X: VectorPrim); overload;
    procedure Free;
    function Summ: Real; { Сумма элементов вектора   }
    function SummSqr: Real; { Сумма квадратов элементов }

    procedure Sort(a1: VectorPrim);
    { Сортировка вектора по возрастанию }
    procedure SortAbs(a1: VectorPrim); { -<>- по модулю                }
    procedure SortBy(a1: VectorPrim);
    { Сортировка вектора по возрастанию }
    procedure Clear; { Очистить нулями }
    procedure Print;
    procedure PrintF(M: Byte);
    procedure ScMultiple(X: Real; y: VectorPrim);
    { Умножение на скаляр }
    procedure ScDivide(X: Real; y: VectorPrim);
    { Деление на скаляр }
    procedure SetAll(X: Real); { Все эл-ты = X }
    procedure Multiple(x, y: VectorPrim); overload;
    { поэлементное умножение }
  end;

  MatrixPrimArray =
    array[1..1] of VectorPrim;

  MatrixVPrimPtr = ^MatrixPrimArray;

  Matrix = { Матрица }  object
    VAddr: VirtualPtr; { Виртуальный адрес массива векторов
      матрица хранится по строкам. Максималь-
      ная размерность MaxRow * MaxRow (tail of DOS)}
    FInit: word; { Флаг инициализации (tail of DOS) }
    n, m: longword; { Размерности матрицы }
    function Addr: MatrixVPrimPtr;
    { Адрес массива векторов-строк }
    function E(i, j: longword): RealPtr; { Адрес элемента }
    procedure Copy(B: Matrix); {   = B }
    procedure Multiple(A, B: Matrix); overload; { A * B }
    procedure Multiple(A, B: VectorPrim); overload;
    procedure Substrac(A, B: Matrix); { A - B }
    procedure Addition(A, B: Matrix); overload; { A + B }
    procedure Addition(A: Matrix; B: VectorPrim); overload;
    procedure Addition(A: Matrix; B: Real); overload;
    procedure Divide(Bx, Ax: Matrix); { B / A }
    procedure Diag(B: VectorPrim); overload; { Diag = b    }
    procedure Diag(r: Real); overload; { Diag(*) = r    }
    procedure Col(i: longword; B: VectorPrim); overload; { Столбец = b }
    procedure Row(i: longword; B: VectorPrim); overload; { Строка = b  }
    procedure Col(i: longword; r: Real); overload; { Столбец = r }
    procedure Row(i: longword; r: Real); overload; { Строка = r  }
    procedure Clear; { Обнуление эл-тов }
    procedure Invert(A: Matrix);
    { Вычисление обратной/псевдообратной матрицы }
    procedure Print;
    procedure PrintF(M2: Byte); { M2 - к-во знач. цифр }
    procedure ScMultiple(X: Real; Y: Matrix);
    { умножение на скаляр }
    procedure Trans(A: Matrix); { Транспонирование }
    procedure MSqrt(AX: Matrix);
    { Вычисление корня для              }
    {                                 положительно-определенной матрицы }
    procedure MAbs(AX: Matrix); { Матричная функция ABS }
    function Cond: Real; { числo обусловленности (1-Ok)}
    function Norm: Real; { Норма матрицы }
    function Alpha: Real; { Параметр регуляризации }
    procedure Init(Size1, Size2: longword); overload;
    { Размещение матрицы }
    procedure Init(A: Matrix); overload;
    procedure InitTrans(A: Matrix);
    procedure InitSqr(A: Matrix);
    procedure ReadFromFile(FN: string); { Читать из файла TXT,
      Init не требуется }
    procedure MSqr(A: Matrix); { = Trans( A ) * A }
    procedure Simmetry(A: Matrix); { = ( Trans( A ) + A ) / 2 }
    function Asimm: Real; { Ассимметрия матрицы }
    function Difference(A: Matrix): Real;
    procedure Free; { Освобождение памяти }
  end;

  Vector = object(VectorPrim)
    procedure Divide(x: Vector; AB: Matrix); overload; { =x / AB }
    procedure Divide(x, y: Vector); overload;
    procedure Multiple(A: Matrix; x: Vector); overload; { =A * x  }
    procedure Diag(A: Matrix); { =Диаголаль }
    procedure Col(i: longword; A: Matrix); { =Столбец }
    procedure Row(i: longword; A: Matrix); { =Строка }
    function Max: Real; { Максимальный элемент }
    function Min: Real; { Минимальный элемент  }
    procedure Abs(x: Vector); { Абсолютная величина  }
    function MinIndex: longword; { Индекс мин.эл-та     }
    function MaxIndex: longword; { Индекс макс. эл-та   }
    function Poly(x: Real): Real; { Полином x }
    procedure IntPower(x: Vector; P: Integer);
    { Целая степенть(поэлементно) }
    procedure SubVector(x: Vector; From, Count: longword);
    procedure InitRow(A: Matrix); overload;
    procedure InitRow(i: Integer; A: Matrix); overload;
    procedure InitCol(A: Matrix); overload;
    procedure InitCol(i: Integer; A: Matrix); overload;
    procedure FillBodyTo(XMax: Real);
    procedure FillBody;
    function Difference(x: Vector): Real;
    procedure TableDiff(x, y: Vector); // Две программы из
    procedure Smoothe(y: Vector); // Фортран-пакета(старые)
  end;

  SplineCoeffType = record
    NPoint: longword;
    h, d, bb, d1, dg, d2: Vector;
  end;

  SplineCoeffPtr = ^SplineCoeffType;

  Spline = { Кубические сплайны }  object
    VAddr: VirtualPtr;
    function Addr: SplineCoeffPtr;
    procedure Init(X, Y: Vector);
    function F(X: Real): Real; { Значение сплайна в точке X }
    function Integral(A, B: Real): Real;
      overload; { Интеграл от А до В  }
    function Integral: Real; overload; { Интеграл  }
    function FMin(A, B: Real): Real;
    function FMax(A, B: Real): Real;
    function ZeroIn(A, B: Real): Real;
    function DF(X: Real): Real;
    { Значение первой производной }
    procedure Free;
  private
    function Integr(A, B: Real): Real; { Интеграл от А до В  }
  end;

  PolyApprox = { Аппроксимация полиномом }  object
    Coeff: Vector;
    dCoeff: Vector;
    PlanM: Matrix;
    procedure Init(xt, yt: Vector; NP: Integer);
    function PX(xr: Real): Real;
    function DPX(xr: Real): Real;
    procedure Free;
  private
    XMin, XMax: Real;
  end;

  Eigen = { Собственные значения ( для симметричных матриц ) }
    object
  private
    function F(Alpha: Real): Real;
  public
    Values: Vector;
    Vectors: Matrix;
    AbsValues: Vector;
    Cond: Real;
    D, R, P: Matrix; { D, R и P - факторы : D = R * R; A = D + R + P; }
    Gamma: Matrix; { Матрица ошибок извлеч. из A }
    procedure Init(A: Matrix);
    procedure Free;
  end;

{*--------------------------------------------------------------------*}

implementation

uses
  HMinZero; // одномерный поиск

procedure Error(E: word);
begin
  ShowMessage('Internal Erorr : ' + IntToStr(E));
  Halt(E);
end;

procedure GetMem(var P: Pointer; Size: Integer); { (tail of DOS) }
begin
  Inc(InitSize, Size);
  Inc(InitCount);
  System.GetMem(P, Size);
end;

procedure FreeMem(P: Pointer; Size: Integer); { (tail of DOS) }
begin
  Dec(InitSize, Size);
  Dec(InitCount);
  System.FreeMem(P, Size);
end;

procedure Vector.TableDiff(x, y: Vector);
{
C        SUBROUTINE DDGT3
C
C        PURPOSE
C           TO COMPUTE A VECTOR OF DERIVATIVE VALUES GIVEN VECTORS OF
C           ARGUMENT VALUES AND CORRESPONDING FUNCTION VALUES.
C
C        USAGE
C           CALL DDGT3(X,Y,Z,NDIM,IER)
C
C        DESCRIPTION OF PARAMETERS
C           X     -  GIVEN VECTOR OF DOUBLE PRECISION ARGUMENT VALUES
C                    (DIMENSION NDIM)
C           Y     -  GIVEN VECTOR OF DOUBLE PRECISION FUNCTION VALUES
C                    CORRESPONDING TO X (DIMENSION NDIM)
C           Z     -  RESULTING VECTOR OF DOUBLE PRECISION DERIVATIVE
C                    VALUES (DIMENSION NDIM)
C           NDIM  -  DIMENSION OF VECTORS X,Y AND Z
C           IER   -  RESULTING ERROR PARAMETER
C                    IER  = -1  - NDIM IS LESS THAN 3
C                    IER  =  0  - NO ERROR
C                    IER POSITIVE  - X(IER) = X(IER-1) OR X(IER) =
C                                    X(IER-2)
C
C        REMARKS
C           (1)   IF IER = -1,2,3, THEN THERE IS NO COMPUTATION.
C           (2)   IF IER =  4,...,N, THEN THE DERIVATIVE VALUES Z(1)
C                 ,..., Z(IER-1) HAVE BEEN COMPUTED.
C           (3)   Z CAN HAVE THE SAME STORAGE ALLOCATION AS X OR Y.  IF
C                 X OR Y IS DISTINCT FROM Z, THEN IT IS NOT DESTROYED.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           EXCEPT AT THE ENDPOINTS X(1) AND X(NDIM), Z(I) IS THE
C           DERIVATIVE AT X(I) OF THE LAGRANGIAN INTERPOLATION
C           POLYNOMIAL OF DEGREE 2 RELEVANT TO THE 3 SUCCESSIVE POINTS
C           (X(I+K),Y(I+K)) K = -1,0,1. (SEE HILDEBRAND, F.B.,
C           INTRODUCTION TO NUMERICAL ANALYSIS, MC GRAW-HILL, NEW YORK/
C           TORONTO/LONDON, 1956, PP. 64-68.)
C
C     ..................................................................
C
      SUBROUTINE DDGT3(X,Y,Z,NDIM,IER)}
var
  DY1, DY2, DY3, A, B: Real;
  i: Integer;
begin
  { TEST OF DIMENSION AND ERROR EXIT IN CASE NDIM IS LESS THAN 3}
  if x.n <> y.n then
    Error(8001);
  { PREPARE DIFFERENTIATION LOOP}
  A := X.e(1)^;
  B := Y.e(1)^;
  I := 2;
  DY2 := X.e(2)^ - A;
  if DY2 = 0 then
    Error(8002);

  DY2 := (Y.e(2)^ - B) / DY2;
  { START DIFFERENTIATION LOOP}
  for I := 3 to x.n do
  begin
    A := X.e(I)^ - A;
    if A = 0 then
      Error(8002);
    A := (Y.e(I)^ - B) / A;
    B := X.e(I)^ - X.e(I - 1)^;
    if B = 0 then
      Error(8002);
    DY1 := DY2;
    DY2 := (Y.e(I)^ - Y.e(I - 1)^) / B;
    DY3 := A;
    A := X.e(I - 1)^;
    B := Y.e(I - 1)^;

    if (I - 3) <= 0 then
      e(1)^ := DY1 + DY3 - DY2;
    e(I - 1)^ := DY1 + DY2 - DY3;
  end;

  e(n)^ := DY2 + DY3 - DY1;

end;

procedure Vector.Smoothe(y: Vector);
{
C        SUBROUTINE DSE15
C
C        PURPOSE
C           TO COMPUTE A VECTOR OF SMOOTHED FUNCTION VALUES GIVEN A
C           VECTOR OF FUNCTION VALUES WHOSE ENTRIES CORRESPOND TO
C           EQUIDISTANTLY SPACED ARGUMENT VALUES.
C
C        USAGE
C           CALL DSE15(Y,Z,NDIM,IER)
C
C        DESCRIPTION OF PARAMETERS
C           Y     -  GIVEN VECTOR OF DOUBLE PRECISION FUNCTION VALUES
C                    (DIMENSION NDIM)
C           Z     -  RESULTING VECTOR OF DOUBLE PRECISION SMOOTHED
C                    FUNCTION VALUES (DIMENSION NDIM)
C           NDIM  -  DIMENSION OF VECTORS Y AND Z
C           IER   -  RESULTING ERROR PARAMETER
C                    IER = -1  - NDIM IS LESS THAN 5
C                    IER =  0  - NO ERROR
C
C        REMARKS
C           (1)  IF IER=-1 THERE HAS BEEN NO COMPUTATION.
C           (2)   Z CAN HAVE THE SAME STORAGE ALLOCATION AS Y.  IF Y IS
C                 DISTINCT FROM Z, THEN IT IS NOT DESTROYED.
C
C        SUBROUTINE AND FUNCTION SUBPROGRAMS REQUIRED
C           NONE
C
C        METHOD
C           IF X IS THE (SUPPRESSED) VECTOR OF ARGUMENT VALUES, THEN
C           EXCEPT AT THE POINTS X(1),X(2),X(NDIM-1) AND X(NDIM), EACH
C           SMOOTHED VALUE Z(I) IS OBTAINED BY EVALUATING AT X(I) THE
C           LEAST-SQUARES POLYNOMIAL OF DEGREE 1 RELEVANT TO THE 5
C           SUCCESSIVE POINTS (X(I+K),Y(I+K)) K = -2,-1,...,2.  (SEE
C           HILDEBRAND, F.B., INTRODUCTION TO NUMERICAL ANALYSIS,
C           MC GRAW-HILL, NEW YORK/TORONTO/LONDON, 1956, PP. 295-302.)
C
C     ..................................................................
C
      SUBROUTINE DSE15(Y,Z,NDIM,IER)
}
var
  A, B, C: Real;
  I, NDIM: Integer;
begin
  {        TEST OF DIMENSION }
  NDIM := y.n;

  if (NDIM < 5) then
    Error(8004);

  { PREPARE LOOP }
  A := Y.e(1)^ + Y.e(1)^;
  C := Y.e(2)^ + Y.e(2)^;
  B := 0.2 * (A + Y.e(1)^ + C + Y.e(3)^ - Y.e(5)^);
  C := 0.1 * (A + A + C + Y.e(2)^ + Y.e(3)^ + Y.e(3)^ + Y.e(4)^);
  { START LOOP }
  for I := 5 to NDIM do
  begin
    A := B;
    B := C;
    C := 0.2E0 * (Y.e(I - 4)^ + Y.e(I - 3)^ + Y.e(I - 2)^
      + Y.e(I - 1)^ + Y.e(I)^);
    e(I - 4)^ := A;
  end;
  { UPDATE LAST FOUR COMPONENTS}
  A := Y.e(NDIM)^ + Y.e(NDIM)^;
  A := 0.10 * (A + A + Y.e(NDIM - 1)^ + Y.e(NDIM - 1)^ +
    Y.e(NDIM - 1)^ + Y.e(NDIM - 2)^ + Y.e(NDIM - 2)^ + Y.e(NDIM - 3)^);
  e(NDIM - 3)^ := B;
  e(NDIM - 2)^ := C;
  e(NDIM - 1)^ := A;
  e(NDIM)^ := A + A - C;
end;

procedure Jacobi(AT, S: Matrix; L: Vector);
var
  i, j, k, p, q: word;

  SinT, CosT, Sin2T, Cos2T, Lambda, Mu, Omega, AMax, App,
  Aqq, Apq, Temp: Real;

  T: Vector;
  A, B, R: Matrix;

begin
  A.Init(AT);
  B.Init(A);

  S.Clear;
  for i := 1 to S.n do
    S.E(i, i)^ := 1;

  R.Init(S);
  T.InitRow(R);

  NIteration := 0;

  repeat

    Inc(NIteration);

    B.Copy(A);
    R.Copy(S);

    p := 1;
    q := 2;
    AMax := A.E(p, q)^;

    T.Diag(A);
    Temp := T.Max;

    for i := 1 to A.n do
      for j := i + 1 to A.n do
        if Abs(A.E(i, j)^) > AMax then
        begin
          p := i;
          q := j;
          AMax := Abs(A.E(p, q)^);
        end;

    if (AMax + Temp) = Temp then
      Break;

    App := A.E(p, p)^;
    Aqq := A.E(q, q)^;
    Apq := A.E(p, q)^;

    Lambda := -Apq;
    Mu := (App - Aqq) / 2.0;
    Omega := Lambda / Sqrt(Sqr(Lambda) + Sqr(Mu));

    if Mu < 0 then
      Omega := -Omega;

    SinT := Omega / Sqrt(2 * (1 + Sqrt(1 - Sqr(Omega))));
    Sin2T := Sqr(SinT);
    CosT := Sqrt(1 - Sin2T);
    Cos2T := Sqr(CosT);

    for k := 1 to A.n do
    begin
      B.E(p, k)^ := A.E(p, k)^ * CosT - A.E(q, k)^ * SinT;
      B.E(q, k)^ := A.E(p, k)^ * SinT + A.E(q, k)^ * CosT;
    end;

    for i := 1 to A.n do
    begin
      B.E(i, p)^ := A.E(i, p)^ * CosT - A.E(i, q)^ * SinT;
      B.E(i, q)^ := A.E(i, p)^ * SinT + A.E(i, q)^ * CosT;

      R.E(i, p)^ := S.E(i, p)^ * CosT - S.E(i, q)^ * SinT;
      R.E(i, q)^ := S.E(i, p)^ * SinT + S.E(i, q)^ * CosT;
    end;

    B.E(p, p)^ := App * Cos2T + Aqq * Sin2T - 2 * Apq * SinT * CosT;
    B.E(q, q)^ := App * Sin2T + Aqq * Cos2T + 2 * Apq * SinT * CosT;
    B.E(p, q)^ := 0;
    B.E(q, p)^ := 0;

    A.Copy(B);
    S.Copy(R);

  until False;

  L.Diag(A);

  B.Free;
  R.Free;
  A.Free;
  T.Free;

end;

function Eigen.F(Alpha: Real): Real;
var
  Q, T, X: Vector;
  St, S: Matrix;
  Z: Real;
begin
  Inc(NIteration);

  S.Init(Gamma);
  St.InitTrans(Vectors);
  Q.InitRow(Gamma);
  T.Init(Q);
  X.Init(Q);

  S.Clear;
  S.Diag(Exp(Alpha));

  S.Multiple(Vectors, S);
  S.Multiple(S, St);

  T.Diag(Gamma);
  T.SetAll(1);

  Q.Multiple(Gamma, T);

  S.Addition(Gamma, S);

  X.Divide(Q, S);

  T.Substrac(X, T);

  Z := T.SummSqr;

  Result := Ln(Z);

  Q.Free;
  T.Free;
  S.Free;
  ST.Free;
  X.Free;

end;

procedure Eigen.Init(A: Matrix);
var
  Iter: longword;
  Q, T, X: Vector;
  St, S: Matrix;
  xZ: Real;
begin

  if A.n <> A.m then { Матрица квадратная ? }
    Error(4001);

  if A.n < 2 then
    Error(4003);

  Values.InitRow(A);
  Vectors.Init(A);

  Jacobi(A, Vectors, Values);
  AbsValues.Init(Values);
  AbsValues.Abs(Values);

  if AbsValues.Min <> 0 then
    Cond := AbsValues.Max / AbsValues.Min
  else
    Cond := 1E+38;

  D.Init(A);
  R.Init(D);
  P.Init(R);

  Gamma.Init(A);
  Q.InitRow(A);
  X.Init(Q);
  T.Init(Q);
  St.Init(A);
  S.Init(St);

  Iter := NIteration;

  D.Clear;
  D.Diag(1);
  R.ScMultiple(4, R);
  R.Addition(R, D);
  R.MSqrt(R);
  R.Substrac(R, D);
  R.ScMultiple(0.5, R);
  D.MSqr(R);
  P.Substrac(A, R);
  P.Substrac(P, D);

  NIteration := Iter + NIteration;

  xZ := Exp(Ln(MathEps) * (Ln(1 / MathEps) /
    (Ln(Cond) + Ln(1 / MathEps))));

  Q.SetAll(Exp(xZ));

  S.Clear;
  S.Diag(Q);

  St.Trans(Vectors);
  S.Multiple(Vectors, S);
  S.Multiple(S, St);

  Gamma.Copy(S);

  Q.Free;
  T.Free;
  X.Free;
  St.Free;
  S.Free;

end;

procedure Vector.Divide(x, y: Vector);
var
  i: longword;
begin
  if x.n <> y.n then
    Error(4010);
  for i := 1 to x.n do
    E(i)^ := x.E(i)^ / y.E(i)^;
end;

procedure Eigen.Free;
begin
  Values.Free;
  Vectors.Free;
  AbsValues.Free;
  D.Free;
  R.Free;
  P.Free;
  Gamma.Free;
end;

procedure Matrix.Simmetry(A: Matrix);
var
  T: Matrix;
begin
  if A.n <> A.m then
    Error(5001);

  if n <> A.n then
    Error(5002);

  T.InitTrans(A);
  T.Addition(T, A);

  ScMultiple(0.5, T);

  T.Free;
end;

procedure Matrix.MAbs(AX: Matrix);
var
  T: Matrix;
begin
  T.InitSqr(AX);
  MSqrt(T);
  T.Free;
end;

function Spline.FMin(A, B: Real): Real;
begin
  Result := HMinZero.FMin(SELF, A, B);
end;

function Spline.FMax(A, B: Real): Real;
begin
  Result := HMinZero.FMax(SELF, A, B);
end;

function Spline.ZeroIn(A, B: Real): Real;
begin
  Result := FZeroIn(SELF, A, B);
end;

function Spline.Integr(A, B: Real): Real;

var
  x1, x3, x5, F1, F3, F5: Real;

  function Int(x1, x3, x5, F1, F3, F5, Eps: Real): Real;

  var
    h1, h2, f2, f4, x2, x4, s1, s2, s3: Real;
    q1, q2, st: Real;

  begin
    h1 := X3 - X1;
    h2 := h1 / 2.;

    x2 := (X3 + X1) / 2.;
    x4 := (X5 + X3) / 2.;

    f2 := F(x2);
    f4 := F(x4);

    {  Формула Симпсона }
    s1 := h1 * (f1 + 4. * f2 + f3) / 6.;
    s2 := h1 * (f3 + 4. * f4 + f5) / 6.;

    { Формула Ньютона-Коттеса 5-го порядка }
    s3 := h2 * 2. * (7. * f1 + 32. * f2 + 12. * f3 + 32. * f4 + 7. * f5) / 45.;

    { Условие окончания }

    st := h2 / Abs(B - A);
    q1 := Abs(s1 + s2 - s3);
    q2 := 0.5 * Eps * Abs(s1 + s2 + s3) * st + Eps;

    if (q1 > q2) and (h2 > Eps + Abs(B - A) * Eps) then
      { Рекурсия - если точность плохая }
      Result := Int(x1, x2, x3, f1, f2, f3, Eps) + Int(
        x3, x4, x5, f3, f4, f5, Eps)
    else

      { За решения принять значение формулы 5-го порядка }
      Result := s3;
  end;

begin

  Result := 0;

  if A = B then
    exit;

  x1 := A;
  x3 := B - (B - A) * 0.5;
  x5 := B;

  F1 := F(x1);
  F3 := F(x3);
  F5 := F(x5);

  Result := Int(x1, x3, x5, F1, F3, F5, Eps);

end;

function Spline.Integral(A, B: Real): Real;
begin
  Result := Integr(A, B);
end;

function Spline.Integral: Real;
begin
  Result := Integr(Addr^.d1.e(1)^, Addr^.d1.e(Addr^.d1.n)^);
end;

procedure Vector.FillBodyTo(XMax: Real);
var
  x, S: Real;
  i: Integer;
begin
  x := Addr^[1];
  S := (XMax - x) / (n - 1);
  for i := 1 to n do
  begin
    Addr^[i] := x;
    x := x + S;
  end;
end;

procedure Vector.FillBody;
begin
  FillBodyTo(Addr^[n]);
end;

procedure Vector.SubVector(x: Vector; From, Count: longword);
var
  i, j: longword;
begin
  if (n < Count) or ((Count + From - 1) > x.n) then
    Error(2004);

  if (x.n < From) then
    Error(2005);

  j := 0;
  for i := From to From + Count - 1 do
  begin
    Inc(j);
    Addr^[j] := x.e(i)^;
  end;
end;

procedure PolyApprox.Init(xt, yt: Vector; NP: Integer);
{
   Аппроксимация зависимости Y( x )
   полиномом степени NP - 1

   Вход :
      X, Y;
}
var
  A: Matrix; { Матрица плана }
  AT: Matrix; { --//-- Транспонированная }
  B: Matrix;
  BY: Vector;
  t: Vector;
  i, j: longword;
  Z: Real;

begin
  Coeff.Init(NP);

  A.Init(xt.n, Coeff.n);
   {
      Формирование матрицы плана
   }
  t.Init(xt.n);

  XMin := xt.Min;
  XMax := xt.Max;

  for i := 1 to Coeff.n do
  begin
    for j := 1 to xt.n do
      t.E(j)^ := xt.e(j)^;

    t.IntPower(t, i - 1);
    A.Col(i, t);
  end;

  { Обычный метод НК }
  PlanM.Init(A); // Запомним матрицу плана
  AT.InitTrans(A);
  B.InitSqr(A); //  B = AT * A
  Coeff.Multiple(AT, yt); // Coeff = AT * y
  Coeff.Divide(Coeff, B); { Coeff = Coeff / B }
  { но, можно и так : }
  { вычисление псевдообратной матрицы }
  //   AT.Invert( A );

  //   Coeff.Multiple( AT, yt ); // Coeff = A+ * yt

  if NP > 1 then
  begin
    dCoeff.Init(Coeff.n - 1);

    for i := 2 to Coeff.n do
      dCoeff.e(i - 1)^ := Coeff.e(i)^ * (i - 1);
  end;

  t.Free;
  A.Free;
  AT.Free;

end;

function PolyApprox.PX(xr: Real): Real;
begin
  Result := Coeff.Poly(xr);
end;

function PolyApprox.DPX(xr: Real): Real;
{
   Производная полинома в "xr"
}
begin
  Result := dCoeff.Poly(xr);
end;

procedure PolyApprox.Free;
begin
  if Coeff.n > 1 then
    dCoeff.Free;

  Coeff.Free;
  PlanM.Free;
end;

function Matrix.Difference(A: Matrix): Real;
begin
  Result := SELF.Norm - A.Norm;
end;

function Vector.Difference(x: Vector): Real;
begin
  Result := Sqrt(SELF.SummSqr) - Sqrt(x.SummSqr);
end;

function Matrix.Asimm: Real;
{
   Ассимметрия матрицы
}
var
  i, j: longword;
  S, T: Real;
begin

  if N <> M then
    Error(2002);

  S := 0;
  T := S;
  for i := 1 to N do
    for j := 1 to N do
      if i <> j then
      begin
        S := S + Math.IntPower(E(i, j)^ - E(j, i)^, 2);
        {               T := T + ( E(i,j)^ + E(j,i)^ ) / 2; }
      end;

  S := Sqrt(S) {/ ( T + 1 )};
  Result := S;

end;

procedure Vector.IntPower(x: Vector; P: Integer);
{
   Поэлементное возведение в целую степень
}
var
  i: longword;
  t: Vector;
begin
  t.Init(x.N);

  for i := 1 to x.n do
    t.E(i)^ := Math.IntPower(x.E(i)^, P);

  Copy(t);
  t.Free;
end;

procedure Matrix.MSqr(A: Matrix);
{     T
   = A  * A ( квадрат матрицы )
}
var
  T: Matrix;

begin
  if (N <> M) or (M <> A.M) then
    Error(2001);

  T.InitTrans(A);
  Multiple(T, A);

  T.Free;

end;

procedure Matrix.InitSqr(A: Matrix);
{     T
   = A  * A ( квадрат матрицы )
}
begin
  Init(A.m, A.m);
  MSqr(A);
end;

function Vector.Poly(x: Real): Real; { Полином x }
{
   Poly = .e(1) + x * ( .e(2) + x * ( .e(3) + x * ( .e(4) + ...x * .e(n))...
}
var
  i: Integer;
  S: Real;
begin
  S := e(n)^;

  for i := n - 1 downto 1 do
    S := e(i)^ + S * x;

  Result := S;
end;

procedure Matrix.ReadFromFile(FN: string);
{
   Чтение матрицы из текстового файла
   ( по строкам )
}
var
  f: Text;
  S: string;
  i, j: longword;
  Y: Matrix;
const
  Delim: CharSet =
    [' ', ';', #9];
begin
  if FileExists(FN) then
  begin
    Assign(f, FN);
    Reset(f);
    n := 0;
    j := n;

    while not EOF(f) do
    begin
      ReadLn(f, S);

      if (Trim(S) = '') then
        Continue;

      if (S[1] = '''') then
        Continue;

      Inc(n);
      m := WordCount(S, Delim);

      if j = 0 then
        j := m
      else
      if j <> m then
      begin
        ShowMessage('Error in input file.(1)');
        Halt;
      end;
    end;

    Y.Init(n, m);
    Y.Clear;
    Reset(f);
    n := 0;
    while not EOF(f) do
    begin
      ReadLn(f, S);

      if (Trim(S) = '') then
        Continue;

      if (S[1] = '''') then
        Continue;

      Inc(n);

      for i := 1 to WordCount(S, Delim) do
        Y.E(n, i)^ := StrToFloat(ExtractWord(i, S, Delim));
    end;

    Init(n, m);
    Copy(Y);
    Y.Free;
    Close(f);
  end
  else
  begin
    ShowMessage('File : "' + FN + '" - not found.');
    Halt;
  end;
end;

function Vector.Max: Real;
{
   Максимальный элемент вектора
}
var
  i: longword;
  P: ArrayPtr;
  X: Real;
begin
  P := Addr;
  X := P^[1];

  for i := 1 to n do
    if P^[i] > X then
      X := P^[i];

  Result := X;
end;

function Vector.Min: Real;
{
   Минимальный элемент вектора
}
var
  i: longword;
  P: ArrayPtr;
  X: Real;
begin
  P := Addr;

  X := P^[1];

  for i := 1 to n do
    if P^[i] < X then
      X := P^[i];

  Result := X;
end;

procedure Vector.Abs(x: Vector);
{
   .e(i)^ = Abs( x.e(i)^ )
}
var
  i: longword;
  P1: ArrayPtr;
  P2: ArrayPtr;
begin
  P1 := Addr;
  P2 := x.Addr;

  for i := 1 to n do
    P1^[i] := System.Abs(P2^[i]);

end;

function Vector.MaxIndex: longword;
{
   Индекс максимального элемента
}
var
  i: longword;
  j: longword;
  P: ArrayPtr;
  X: Real;
begin
  P := Addr;
  X := P^[1];
  j := 1;

  for i := 1 to n do
    if P^[i] > X then
    begin
      X := P^[i];
      j := i;
    end;

  Result := j;
end;

function Vector.MinIndex: longword;
{
   Индекс минимального элемента
}
var
  i: longword;
  j: longword;
  P: ArrayPtr;
  X: Real;
begin
  P := Addr;
  X := P^[1];
  j := 1;

  for i := 1 to n do
    if P^[i] < X then
    begin
      X := P^[i];
      j := i;
    end;

  Result := j;
end;

function Matrix.Cond: Real;
{
   Число обусловленности матрицы( произведение нормы обратной и нормы исходной )
              -1
   Cond = ||.|| * ||.||.
}
var
  B, A: Matrix;
  n0, n1: Real;
begin
  n0 := SELF.Norm;

  if n = m then
  begin
    B.Init(SELF);
    A.Init(B);
    B.Clear;
    B.Diag(1);
    B.Divide(B, A);
  end
  else
    Error(1011);

  n1 := B.Norm;

  Result := n1 * n0;

  B.Free;
  A.Free;

end;

function Matrix.Norm: Real;
{
   Норма Фробениуса :
   = Sqrt( Summ( Sqr( .e( *, * ) ) ) )
}
var
  i: Integer;
  S: Real;
begin
  S := 0;

  for i := 1 to n do
    S := S + Addr^[i].SummSqr;

  Result := Sqrt(S);
end;

procedure Matrix.MSqrt(AX: Matrix);
{
   HNV 1997

   Квадратный корень матрицы A.
   ( предполагается, что матрица A -
     положительно определенная. )

   R такая, что   T
                 R * R = A

   Симметричное/асимметричное разложение матрицы.
   Ассиметрия характерна для плохо обусловленных.

}
var
  Y, X, R, T: Matrix;
  Diff, DiffOld: Real;
  Alphax, EAlpha: Real;
  Reject: Boolean;

begin

  Y.Init(n, n);
  T.Init(Y);
  X.Init(T);
  R.Init(X);

  X.ScMultiple(0.5, AX);
  X.Addition(X, 1); { B = AX / 2 + 1 }

  NIteration := 0;

  Y.MSqr(X);
  Y.Substrac(Y, AX);
  DiffOld := Y.Norm;

  Alphax := 0.5;
  { Обычная реализация метода Ньютона }
  repeat
    Inc(NIteration);

    Y.MSqr(X);
    Y.Substrac(Y, AX);

    T.Trans(X);
    Y.Divide(Y, T);

    repeat
      T.ScMultiple(Alphax, Y);
      R.Substrac(X, T);

      T.MSqr(R);
      T.Substrac(T, AX);

      Diff := T.Norm;

      Reject := Diff > DiffOld;

      if Reject then
        Alphax := -Alphax / 17
      else
        Break;

      EAlpha := Alphax * Alphax + 1;
    until EAlpha = 1;

    if (not Reject) and (Abs(Alphax) * (DiffOld - Diff) >
      (DiffOld + 1) * Eps) then
    begin
      X.Copy(R);
      DiffOld := Diff;
      if Abs(Alphax) < 0.5 then
        Alphax := Alphax * 2;
    end
    else
      Break;

  until NIteration > 1023;
  // Обычно сходится за 20-30 итераций...

  Trans(X); // - результат

  Y.Free;
  X.Free;
  R.Free;
  T.Free;
end;

procedure Matrix.Init(Size1, Size2: longword);
var
  i: longword;
begin
  if Size1 = 0 then
    Error(2018);

  m := Size2;
  n := Size1;

  VAddr.Size := n * SizeOf(VectorPrim);
  GetMem(VAddr.Ptr, VAddr.Size);
  //   VAddr := GetVMem( n * SizeOf( VectorPrim ) );

  { Разпределить массив векторов-строк }
  for i := 1 to n do
    Addr^[i].Init(m);

  FInit := Indicator; { Флаг инициализации }
end;

procedure Matrix.Init(A: Matrix);
begin
  Init(A.n, A.m);
  Copy(A);
end;

procedure Matrix.InitTrans(A: Matrix);
begin
  Init(A.m, A.n);
  Trans(A);
end;

procedure Matrix.Trans(A: Matrix);
{
   Транспонирование матрицы A
}
var
  i: longword;
  x: Vector;
  B: Matrix;
begin
  if (n <> A.m) or (m <> A.n) then
    Error(116);

  x.Init(A.n);
  B.Init(A.m, A.n);

  for i := 1 to A.m do
  begin
    x.Col(i, A);
    B.Row(i, x);
  end;

  Copy(B);
  B.Free;
  x.Free;

end;

procedure VectorPrim.Init(Size: longword);
begin
  if Size = 0 then
    Error(2017);

  n := Size;
  { Распределить память для вектора }
  VAddr.Size := n * SizeOf(real);
  GetMem(VAddr.Ptr, VAddr.Size);
  //   VAddr := GetVMem( n * SizeOf( Real ) );
  FInit := Indicator;
end;

procedure Vector.InitRow(A: Matrix);
begin
  Init(A.m);
end;

procedure Vector.InitCol(A: Matrix);
begin
  Init(A.n);
end;

procedure Vector.InitRow(i: Integer; A: Matrix);
begin
  InitRow(A);
  Row(i, A);
end;

procedure Vector.InitCol(i: Integer; A: Matrix);
begin
  InitCol(A);
  Col(i, A);
end;

procedure VectorPrim.Init(X: VectorPrim);
begin
  n := X.n;

  { Распределить память для вектора }
  VAddr.Size := n * SizeOf(real);
  GetMem(VAddr.Ptr, VAddr.Size);
  //   VAddr := GetVMem( n * SizeOf( Real ) );
  FInit := Indicator;

  Copy(X);
end;

procedure VectorPrim.SetAll(X: Real);
var
  i: longword;
  t: ArrayPtr;
begin
  t := Addr;
  for i := 1 to n do
    t^[i] := X;
end;

procedure Spline.Init(X, Y: Vector);

  procedure SplineXY;
      {
        Вычисление коэффициентов кубического сплайна функции
        заданной таблицей значений Y( x ).
      }
  var
    i, j: longword;
    f: Real;
  begin
    NumbSpline := 1;
    with Addr^ do
    begin

      for i := 1 to NPoint - 1 do
      begin
        h.E(i)^ := x.E(i + 1)^ - x.E(i)^;
        d.E(i)^ := (y.E(i + 1)^ - y.E(i)^) / h.E(i)^;
      end;

      dg.E(1)^ := 2.;
      dg.E(NPoint)^ := 2.;
      d2.E(1)^ := 1.;
      d1.E(NPoint - 1)^ := 1.;
      bb.E(1)^ := d.E(1)^ * 3.;
      bb.E(NPoint)^ := d.E(NPoint - 1)^ * 3.;

      for i := 2 to NPoint - 1 do
      begin
        bb.E(i)^ :=
          (d.E(i)^ * h.E(i - 1)^ + d.E(i - 1)^ * h.E(i)^) * 3.;
        d1.E(i - 1)^ := h.E(i)^;
        d2.E(i)^ := h.E(i - 1)^;
        dg.E(i)^ := 2. * (h.E(i)^ + h.E(i - 1)^);
      end;
            {
                Решение системы линейных уравнений с 3-х диагональной
                матрицей коэффициентов
            }
      for i := 1 to NPoint - 1 do
      begin
        f := d1.E(i)^ / dg.E(i)^;
        dg.E(i + 1)^ := dg.E(i + 1)^ - d2.E(i)^ * f;
        bb.E(i + 1)^ := bb.E(i + 1)^ - bb.E(i)^ * f;
      end;

      bb.E(NPoint)^ := bb.E(NPoint)^ / dg.E(NPoint)^;

      for i := 1 to NPoint - 1 do
      begin
        j := NPoint - i;
        bb.E(j)^ :=
          (bb.E(j)^ - d2.E(j)^ * bb.E(j + 1)^) / dg.E(j)^;
      end;

      for i := 1 to NPoint do
      begin
        d1.Copy(x); { d1 - исходные значения X }
        d2.Copy(y); { d2 - ---:---- Y          }
      end;
    end;
  end;

begin
  VAddr.Size := SizeOf(SplineCoeffType);
  GetMem(VAddr.Ptr, SizeOf(SplineCoeffType));
  //   VAddr := GetVMem( SizeOf( SplineCoeffType ) );

  with Addr^ do
  begin
    NPoint := X.n;
    H.Init(X.n);
    D.Init(X.n);
    bb.Init(X.n);
    d1.Init(X.n);
    dg.Init(X.n);
    d2.Init(X.n);
  end;

  SplineXY;

end;

function Spline.F(X: Real): Real;
{
   Вычисление значения сплайн-функции в точке X.
}

  function SplineIndex: longword;
      {
         Вычисляет номер сплайна для точки 'X'.
      }
  var
    NS: longword;

    procedure IndRecurs(K, L: longword);
    var
      M: longword;
    begin
      with Addr^ do
      begin
        NS := K;

        if (L - K) <= 1 then
          Exit;

        M := ((K + L) div 2);

        if ((X - D1.E(K)^) * (X - D1.E(M)^)) <= 0.0 then
          IndRecurs(K, M)
        else
          IndRecurs(M, L);
      end;
    end;

  begin
    with Addr^ do
    begin

      Result := 1;

      if X <= d1.E(1)^ then
        exit;

      Result := NPoint - 1;

      if X >= d1.E(NPoint)^ then
        exit;

      IndRecurs(1, NPoint);
      Result := NS;

    end;
  end;

var
  j: longword;
  f1, t, s: Real;
begin
  with Addr^ do
  begin

    j := SplineIndex;

    t := (X - d1.E(j)^) / h.E(j)^;
    f1 := 1. - t;
    s := t * d2.E(j + 1)^ + f1 * d2.E(j)^;
    s := s + h.E(j)^ * f1 * t * ((bb.E(j)^ - d.E(j)^) *
      f1 - (bb.E(j + 1)^ - d.E(j)^) * t);

    Result := s;
  end;
end;

function Spline.DF(X: Real): Real;
{
   Вычисление первой производной от сплайн-функции в точке 'X'.
}
var
  Delta: Real;
  Temp: Real;
  PA: PolyApprox;
  xS, yS: Vector;

  procedure PolyCalc(n: Integer);
  begin
    xS.Init(5);
    yS.Init(5);
    xS.SubVector(Addr^.D1, n, 5);
    yS.SubVector(Addr^.D2, n, 5);
    PA.Init(xS, yS, 5);
    Result := PA.DPX(X);
    PA.Free;
    yS.Free;
    xS.Free;
  end;

begin
  Delta := Abs(X) * SqrtEps + SqrtEps;

  with Addr^ do
  begin

    if X - Delta <= D1.E(1)^ then
      Temp := -3 * F(X) + 4 * F(X + Delta) - F(X + 2 * Delta)
    else
    if X + Delta >= D1.E(Npoint)^ then
      Temp := 3 * F(X) - 4 * F(X - Delta) + F(X - 2 * Delta)
    else
      Temp := (F(X + Delta) - F(X - Delta));

    Result := Temp / Delta / 2;
  end;
end;

procedure VectorPrim.Free;
begin
  FInit := 0;
  FreeMem(VAddr.Ptr, VAddr.Size);
  VAddr.Ptr := nil;
  //   FreeVMem( VAddr );
  //   VAddr := 0;
end;

procedure Matrix.Free;
var
  i: longword;
begin
  FInit := 0;

  for i := 1 to n do
    Addr^[i].Free;

  FreeMem(VAddr.Ptr, VAddr.Size);
  VAddr.Ptr := nil;
  //FreeVMem( VAddr );
  //   VAddr := 0;
end;

procedure Spline.Free;
begin
  with Addr^ do
  begin
    H.Free;
    D.Free;
    bb.Free;
    d1.Free;
    dg.Free;
    d2.Free;
  end;
  FreeMem(VAddr.Ptr, VAddr.Size);
  VAddr.Ptr := nil;
  //FreeVMem( VAddr );
  //   VAddr := 0;
end;

function VectorPrim.Addr: ArrayPtr;
begin
  Result := VAddr.Ptr;
  //Addr := VirtualToPtr( VAddr );
end;

function Matrix.Addr: MatrixVPrimPtr;
begin
  Result := VAddr.Ptr;
  //Addr := VirtualToPtr( VAddr );
end;

function Matrix.E(I, J: longword): RealPtr;
begin
  if (i >= 1) and (i <= n) and (j >= 1) and (j <= m) then
    Result := System.Addr(Addr^[i].E(j)^) { Адрес элемента I,J }
  else
    Error(4);
end;

function VectorPrim.E(I: longword): RealPtr;
begin
  if (i >= 1) and (i <= n) then
    Result := System.Addr(Addr^[i]) { Адрес элемента I }
  else
    Error(5);
end;

procedure VectorPrim.Copy(X: VectorPrim);
begin
  if FInit <> Indicator then
    Error(100);
  Move(X.Addr^, Addr^, X.n * SizeOf(real)); { Быстрая пересылка }
end;

procedure Matrix.Copy(B: Matrix);
var
  i: longword;
begin
  if FInit <> Indicator then
    Error(101);
  for i := 1 to n do
    Addr^[i].Copy(B.Addr^[i]); { Векторное копирование }
end;

procedure Matrix.Multiple(A, B: Matrix);
{
   Вычисление произведения двух матриц
      = A * B
}
var
  i, j: longword;
  TX: Matrix;
  t1, t2: Vector;
begin

  if FInit <> Indicator then
    Error(102);

  TX.Init(A.n, B.m);
  t1.Init(A.m);

  t2.Init(B.n);

  for i := 1 to A.n do
  begin
    t1.Row(i, A);
    for j := 1 to B.m do
    begin
      t2.Col(j, B);
      t2.Multiple(t2, t1);
      TX.E(i, j)^ := t2.Summ;
    end;
  end;

  Copy(TX);

  TX.Free;
  t1.Free;
  t2.Free;

end;

procedure Matrix.Multiple(A, B: VectorPrim);
{
                                              T
   Вычисление произведения двух векторов A * B
}
var
  i, j: longword;
begin

  if FInit <> Indicator then
    Error(102);

  if A.n <> B.n then
    Error(4001);

  if n <> m then
    Error(4002);

  if n <> A.n then
    Error(4003);

  for i := 1 to n do
    for j := 1 to n do
      E(i, j)^ := A.E(i)^ * B.E(j)^;

end;

procedure Vector.Multiple(A: Matrix; x: Vector);
{
   Вычисление произведения матрицы на вектор
      = A * x
}
var
  i: longword;
  T1, T2: Vector;
begin

  if x.n <> A.m then
    Error(6);

  if FInit <> Indicator then
    Error(103);

  T1.Init(A.n);
  T2.Init(A.m);

  for i := 1 to A.n do
  begin
    T2.Row(i, A); { Быстрая пересылка }
    T2.Multiple(T2, x); { Умножение         }
    T1.Addr^[i] := T2.Summ; { Суммирование      }
  end;

  Copy(T1);

  T1.Free;
  T2.Free;

end;

procedure VectorPrim.Multiple(x, y: VectorPrim);
var
  i: longword;
  yp, xp, Sp: ArrayPtr;
begin
  if x.n <> y.n then
    Error(121);

  Sp := Addr;
  xp := x.Addr;
  yp := y.Addr;

  for i := 1 to x.n do
    Sp^[i] := xp^[i] * yp^[i];

end;

procedure Vector.Divide(x: Vector; AB: Matrix);
{
   Решение системы линейных алгебраических уравнений
   вида:
        A * y = b
   методом прямого деления.
}
var
  b: Matrix;
begin
  if (x.n <> AB.m) or (AB.m <> AB.n) then
    Error(7);

  if FInit <> Indicator then { Новый вектор }
    Error(104);

  b.Init(AB.n, 1);
  b.Col(1, x);

  b.Divide(b, AB);

  Col(1, b);
  b.Free;

end;

procedure Matrix.Divide(Bx, Ax: Matrix);
{
   Решение матричной системы вида :
                A * Y = B, где
      A[n,n], Y[n,m], B[n,m] m<=n
   может быть применена для обращения матрицы
   если В = Е.
}
var
  i, j, k: longword;
  d, t: Real;
  a, b: Matrix;
  FBreak: Boolean;
begin

  if FInit <> Indicator then { Новая матрица }
    Error(105);

  a.Init(Ax.n, Ax.m);
  b.Init(Bx.n, Bx.m);

  a.Copy(Ax);
  b.Copy(Bx);

  Condition := True;
  MinEps := 1.0;

  FBreak := False;

  for i := 1 to a.n do
  begin
    if System.Abs(a.E(i, i)^) = 0.0 then
    begin
      Condition := Condition and (i = 1);

      for j := 1 to a.n do
      begin
        if System.Abs(a.E(j, i)^) <> 0.0 then
        begin
          a.Addr^[i].Addition(a.Addr^[i], a.Addr^[j]);
          b.Addr^[i].Addition(b.Addr^[i], b.Addr^[j]);
          Break;
        end;
      end;

      FBreak := True;

    end;

    if FBreak then
      Continue;

    for k := 1 to b.m do
      b.E(i, k)^ := b.E(i, k)^ / a.E(i, i)^;

    for k := a.n downto i do
      a.E(i, k)^ := a.E(i, k)^ / a.E(i, i)^;

    for j := 1 to a.n do
    begin
      d := a.E(j, i)^;

      if (j <> i) and (System.Abs(d) <> 0.0) then
      begin
        for k := i + 1 to a.n do
        begin
          t := a.E(j, k)^ / d - a.E(i, k)^;

          a.E(j, k)^ := t;
        end;

        for k := 1 to b.m do
          b.E(j, k)^ := b.E(j, k)^ / d - b.E(i, k)^;

        if j < i then
          a.E(j, j)^ := a.E(j, j)^ / d;
      end;
    end;
  end;

  for i := 1 to a.n do
    for k := 1 to b.m do
      if a.E(i, i)^ <> 0 then
        E(i, k)^ := b.E(i, k)^ / a.E(i, i)^
      else
        E(i, k)^ := 0;

  a.Free;
  b.Free;

end;

procedure Vector.Diag(A: Matrix);
{
   = Diag( A );
}
var
  i: longword;
  P: ArrayPtr;
begin
  if FInit <> Indicator then { Новый вектор }
    Error(106);

  P := Addr;

  for i := 1 to n do
    P^[i] := A.E(i, i)^;
end;

procedure Vector.Col(i: longword; A: Matrix);
{
   = A( *, i )
}
var
  j: longword;
  P: ArrayPtr;
begin
  if FInit <> Indicator then { Проверк инициализации }
    Error(107);

  P := Addr;

  for j := 1 to A.n do
    P^[j] := A.E(j, i)^;
end;

procedure Vector.Row(i: longword; A: Matrix);
{
   = A( i, * );
}
begin
  Copy(A.Addr^[i]);
end;

procedure VectorPrim.Addition(x, y: VectorPrim);
{
   = x.e(*)^ + y.e(*)^;
}
var
  j: longword;
  xP, yP, sP: ArrayPtr;
begin

  if (x.n <> y.n) then
    Error(9);

  if FInit <> Indicator then
    Error(108);

  xP := x.Addr;
  yP := y.Addr;
  sP := Addr;

  for j := 1 to n do
    sP^[j] := xP^[j] + yP^[j];
end;

procedure VectorPrim.Substrac(x, y: VectorPrim);
{
   = x.e(*)^ - y.e(*)^;
}
var
  j: longword;
  xP, yP, sP: ArrayPtr;
begin

  if (x.n <> y.n) then
    Error(10);

  if FInit <> Indicator then
    Error(109);

  xP := x.Addr;
  yP := y.Addr;
  sP := Addr;

  for j := 1 to n do
    sP^[j] := xP^[j] - yP^[j];
end;

function VectorPrim.Summ: Real;
{
   Сумма эл-тов вектора
}
var
  i: longword;
  S: Real;
  T: VectorPrim;
  AP: ArrayPtr;
begin
  S := 0.0;
  T.Init(n);

  AP := T.Addr;
  Move(Addr^, AP^, n * SizeOf(real));

  T.SortAbs(T); { Сортировка по возрастанию }

  for i := 1 to n do { Суммируем начиная с меньших эл-тов }
    S := S + AP^[i];

  Result := S;

  T.Free;
end;

function VectorPrim.SummSqr: Real;
{
   Сумма квадратов эл-тов вектора
}
var
  i: longword;
  S: Real;
  T: VectorPrim;
  TP: ArrayPtr;
begin
  S := 0.0;

  T.Init(n);
  TP := T.Addr;

  Move(Addr^, TP^, SizeOf(real) * n);

  T.SortAbs(T); { сортировка по возрастанию }

  for i := 1 to n do
    S := S + Sqr(TP^[i]);

  Result := S;

  T.Free;
end;

procedure Matrix.Addition(A, B: Matrix);
{
   = A( *, * ) + B( *, * );
}
var
  i: longword;
begin
  if (A.n <> B.n) or (A.m <> B.m) then
    Error(11);

  if FInit <> Indicator then
    Error(110);

  for i := 1 to n do
    Addr^[i].Addition(A.Addr^[i], B.Addr^[i]);
end;

procedure Matrix.Addition(A: Matrix; B: VectorPrim);
{
   = A( *, * ) + B( *, * );
}
var
  i: longword;
begin

  if (A.n <> n) or (n <> m) then
    Error(11);

  if FInit <> Indicator then
    Error(110);

  for i := 1 to n do
    E(i, i)^ := A.E(i, i)^ + B.E(i)^;
end;

procedure Matrix.Addition(A: Matrix; B: Real);
{
   = A( *, * ) + B( *, * );
}
var
  i: longword;
begin

  if (A.n <> A.m) then
    Error(11);

  if FInit <> Indicator then
    Error(110);

  Copy(A);
  for i := 1 to n do
    E(i, i)^ := A.E(i, i)^ + B;
end;

procedure Matrix.Substrac(A, B: Matrix);
{
   = A( *, * ) - B( *, * );
}
var
  i: longword;
begin
  if (A.n <> B.n) or (A.m <> B.m) then
    Error(11);

  if FInit <> Indicator then
    Error(111);

  for i := 1 to n do
    Addr^[i].Substrac(A.Addr^[i], B.Addr^[i]);
end;

procedure Matrix.Diag(b: VectorPrim);
{
   Diag(*) = b(*);
}
var
  i: longword;
begin
  if m <> b.n then
    Error(12);
  for i := 1 to n do
    E(i, i)^ := b.E(i)^;
end;

procedure Matrix.Diag(r: Real);
{
   Diag(*) = r;
}
var
  i: longword;
begin
  for i := 1 to n do
    E(i, i)^ := r;
end;

procedure Matrix.Col(i: longword; b: VectorPrim);
{
   .e( *, i ) = b( * );
}
var
  j: longword;
begin
  if n <> b.n then
    Error(13);
  for j := 1 to n do
    E(j, i)^ := b.E(j)^;
end;

procedure Matrix.Row(i: longword; b: VectorPrim);
{
   .e( i, * ) = b( * );
}
begin
  if m <> b.n then
    Error(13);
  Addr^[i].Copy(b);
end;

procedure Matrix.Col(i: longword; r: Real);
{
   .e( *, i ) = r;
}
var
  j: longword;
begin
  for j := 1 to n do
    E(j, i)^ := r;
end;

procedure Matrix.Row(i: longword; r: Real);
{
   .e( i, * ) = r;
}
var
  j: longword;
begin
  for j := 1 to m do
    E(i, j)^ := r;
end;

function Spline.Addr: SplineCoeffPtr;
begin
  Result := VAddr.Ptr;
  //Addr := VirtualToPtr( VAddr );
end;

procedure VectorPrim.Sort(a1: VectorPrim);
{
   Quick-сортировка по возрастанию.
}
var
  a: ArrayPtr;
  t: Vector;

  procedure SortPrim(l, r: longword);
  var
    i, j: longword;
    x, y: Real;
  begin
    i := l;
    j := r;
    x := a^[(l + r) div 2];
    repeat

      while a^[i] < x do
        i := i + 1;

      while x < a^[j] do
        j := j - 1;

      if i <= j then
      begin
        y := a^[i];
        a^[i] := a^[j];
        a^[j] := y;

        i := i + 1;
        j := j - 1;

      end;
    until i > j;

    if l < j then
      SortPrim(l, j);
    if i < r then
      SortPrim(i, r);

  end;

begin {quicksort}
  ;

  t.Init(a1.n);
  t.Copy(a1);
  a := t.Addr;

  SortPrim(1, N);

  Copy(t);
  t.Free;
end;

procedure VectorPrim.SortBy(a1: VectorPrim);
{
   Quick-сортировка по возрастанию.
}
var
  a, b: ArrayPtr;
  t: Vector;

  procedure SortPrim(l, r: longword);
  var
    i, j: longword;
    x, y: Real;
  begin
    i := l;
    j := r;
    x := a^[(l + r) div 2];
    repeat

      while a^[i] < x do
        i := i + 1;

      while x < a^[j] do
        j := j - 1;

      if i <= j then
      begin
        y := a^[i];
        a^[i] := a^[j];
        a^[j] := y;

        y := b^[i];
        b^[i] := b^[j];
        b^[j] := y;

        i := i + 1;
        j := j - 1;

      end;
    until i > j;

    if l < j then
      SortPrim(l, j);

    if i < r then
      SortPrim(i, r);

  end;

begin {quicksort}
  t.Init(a1.n);
  t.Copy(a1);
  a := t.Addr;
  b := Addr;

  SortPrim(1, N);

  {   Copy( t );}
  t.Free;
end;

procedure VectorPrim.SortAbs(a1: VectorPrim);
{
   Quick-сортировка по возрастанию ABS.
}
var
  a: ArrayPtr;
  t: Vector;

  procedure SortPrim(l, r: longword);
  var
    i, j: longword;
    x, y: Real;
  begin
    i := l;
    j := r;
    x := Abs(a^[(l + r) div 2]);
    repeat

      while Abs(a^[i]) < x do
        i := i + 1;

      while x < Abs(a^[j]) do
        j := j - 1;

      if i <= j then
      begin
        y := a^[i];
        a^[i] := a^[j];
        a^[j] := y;

        i := i + 1;
        j := j - 1;

      end;
    until i > j;

    if l < j then
      SortPrim(l, j);
    if i < r then
      SortPrim(i, r);

  end;

begin {quicksort}
  t.Init(a1.n);
  t.Copy(a1);
  a := t.Addr;

  SortPrim(1, N);

  Copy(t);
  t.Free;
end;

procedure VectorPrim.Clear;
begin
  SetAll(0);
end;

procedure Matrix.Clear;
var
  i: longword;
begin
  for i := 1 to n do
    Addr^[i].Clear;
end;

function Matrix.Alpha: Real;
var
  QCond, xZ, HNorm: Real;
  Temp: Vector;
begin
  Temp.Init(n);
  QCond := SELF.Cond;
  Temp.Diag(SELF);
  HNorm := Sqrt(Temp.SummSqr);
  xZ := Ln(MathEps) * (Ln(1 / MathEps) / (Ln(QCond) + Ln(1 / MathEps)));
  Result := Exp(xZ) * HNorm;
  Temp.Free;
end;

procedure Matrix.Invert(A: Matrix);
{
   Метод деления для вычисления обратной/псевдообратной матрицы 'A'.
}
var
  B: Matrix;
  Q: Matrix;
  R: Matrix;
  QCond, RCond: Real;
begin

  Q.InitSqr(A);

  //  Комбинированный метод корня и понижениея числа обусловленности

  Q.Addition(Q, Q.Alpha);

  R.Init(Q);

  B.InitTrans(A);

  QCond := Q.Cond;

  //   WriteLn('Cond(Q)=', QCond );

  R.MSqrt(Q); // квадратный корень Q
  RCond := R.Cond;

  //   WriteLn('Cond(R)=', RCond );

  //  Вычисление псевдообратной матрицы

  B.Divide(B, R);
  R.Trans(R);
  Divide(B, R);

  Q.Free;
  R.Free;
  B.Free;

end;

procedure VectorPrim.Print;
{ (tail of DOS) }
var
  i: longword;
  l: Integer;
begin
  l := Trunc((79 - n) / n);
  if l < 9 then
  begin
    WriteLn('Warning, length of output digit < 9. ');
    for i := 1 to n do
      WriteLn(i: 5, ' ', E(i)^);
    exit;
  end;

  for i := 1 to n do
    Write(E(i)^: l, ' ');

  WriteLn;

end;

procedure VectorPrim.PrintF(M: Byte);
{
  (tail of DOS)
}
var
  i: longword;
  l: Byte;
begin
  l := Trunc((79 - n) / n);

  if l - M < 2 then
  begin
    WriteLn('Warning, invalid format');
    exit;
  end;

  if l < 7 then
  begin
    WriteLn('Warning, length of fixed digit < 7. ');
    exit;
  end;

  if l > 18 then
    l := 18;

  for i := 1 to n do
    Write(E(i)^: l: M, ' ');
  WriteLn;

end;

procedure Matrix.Print;
{
   (tail of DOS)
}
var
  i: longword;
begin
  for i := 1 to n do
    Addr^[i].Print;
  WriteLn;
end;

procedure Matrix.PrintF(M2: Byte);
{
   (tail of DOS)
}
var
  i: longword;
begin
  for i := 1 to n do
    Addr^[i].PrintF(M2);
  WriteLn;
end;

procedure VectorPrim.ScMultiple(X: Real; y: VectorPrim);
var
  i: longword;
  t: ArrayPtr;
  yp: ArrayPtr;
begin
  t := Addr;
  yp := y.Addr;
  for i := 1 to n do
    t^[i] := X * yp^[i]; { Быстрое умножение }
end;

procedure VectorPrim.ScDivide(X: Real; y: VectorPrim);
var
  i: longword;
  t: ArrayPtr;
  yp: ArrayPtr;
begin
  t := Addr;
  yp := y.Addr;
  for i := 1 to n do
    t^[i] := yp^[i] / X; { Быстрое умножение }
end;

procedure Matrix.ScMultiple(X: Real; Y: Matrix);
var
  i: longword;
begin
  for i := 1 to n do
    Addr^[i].ScMultiple(X, Y.Addr^[i]);
end;

begin
   {
      Вычисление точности:
          MathEps = относительная машинная точность.
                ( -log( MathEps ) ~= к-во знаков в мантиссе. )
          SqrtEps = MathEps ^ 1/2;
                ( - точность половина знаков мантиссы. )
          Eps = MathEps ^ 3/4;
                ( --//-- 3/4 знаков мантиссы. )
   }
  MathEps := 1;
  repeat
    MathEps := MathEps / 2;
    Temp := 1 + MathEps;
  until Temp = 1;

  SqrtEps := Sqrt(MathEps);

  Eps := Sqrt(SqrtEps) * SqrtEps;

  MinEps := 1;

  Randomize;
  Indicator := Random($7FFF);
  //InitVHeap( '$$Matr.vhp', true );
end.

Этот модуль используется почти во всех реализованных мной численных алгоритмах и методах. Те части, которые писал не я - приводятся без изменений(по возможности) стиля и комментариев.

HNVString.pas
unit HNVString;

interface

{*********************************************************}
{*                  OPSTRING.PAS 1.20                    *}
{*     Copyright (c) TurboPower Software 1987, 1992.     *}
{*                 All rights reserved.                  *}
{*********************************************************}

type
  CharSet =
    set of char;

{ Given a set of word delimiters, return number of words in S}
function WordCount(S: string; WordDelims: CharSet): Byte;

{ Given a set of word delimiters, return start position of N'th word in S}
function WordPosition(N: Byte; S: string; WordDelims: CharSet): Byte;

{ Given a set of word delimiters, return the N'th word in S}
function ExtractWord(N: Byte; S: string; WordDelims: CharSet): string;

function RightJust(S: string; N: word): string;

function HexB(B: Byte): string;

function HexText(var Buff; N: Integer): string;

implementation

const
  Dig: array[0..$F] of char = '0123456789ABCDEF';

{ Convert BYTE to Hex String }
function HexB(B: Byte): string;
var
  T: string;
begin
  T := '';
  T := T + Dig[B shr 4];
  T := T + Dig[B and $0F];
  HexB := T;
end;

{ Convert BUFFER to HEX String }
function HexText(var Buff; N: Integer): string;
type
  ByteArray =
    array of Byte;
var
  i: Integer;
  S: string;
  B: ByteArray Absolute Buff;
begin
  S := '';
  for i := 1 to N do
    S := S + HexB(B[i]);
  HexText := S;
end;


function RightJust(S: string; N: word): string;
var
  T: string;
  i: Integer;
  j: Integer;
begin
  T := '';
  j := N;
  if Length(S) > j then
    T := S
  else
  begin
    j := j - Length(S);
    for i := 1 to j do
      T := T + ' ';
    T := T + S;
  end;
  RightJust := T;
end;

function WordCount(S: string; WordDelims: CharSet): Byte;
  {-Given a set of word delimiters, return number of words in S}
var
  Count: Byte;         {!!.02}
  I: word;                  {!!.02}
  SLen: Integer;
begin

  SLen := Length(S);
  Count := 0;
  I := 1;

  while I <= SLen do
  begin
    {skip over delimiters}
    while (I <= SLen) and (S[I] in WordDelims) do
      Inc(I);

    {if we're not beyond end of S, we're at the start of a word}
    if I <= SLen then
      Inc(Count);

    {find the end of the current word}
    while (I <= SLen) and not (S[I] in WordDelims) do
      Inc(I);
  end;

  WordCount := Count;
end;

function WordPosition(N: Byte; S: string; WordDelims: CharSet): Byte;
  {-Given a set of word delimiters, return start position of N'th word in S}
var
  Count: Byte;         {!!.02}
  I: word;                  {!!.02}
  SLen: word;
begin
  Count := 0;
  I := 1;

  SLen := Length(S);

  WordPosition := 0;

  while (I <= SLen) and (Count <> N) do
  begin
    {skip over delimiters}
    while (I <= SLen) and (S[I] in WordDelims) do
      Inc(I);

    {if we're not beyond end of S, we're at the start of a word}
    if I <= SLen then
      Inc(Count);

    {if not finished, find the end of the current word}
    if Count <> N then
      while (I <= SLen) and not (S[I] in WordDelims) do
        Inc(I)
    else
      WordPosition := I;
  end;
end;

function ExtractWord(N: Byte; S: string; WordDelims: CharSet): string;
  {-Given a set of word delimiters, return the N'th word in S}
var
  I: word; {!!.12}
  SLen: Integer;
  Temp: string;
begin

  SLen := Length(S);
  Temp := '';
  I := WordPosition(N, S, WordDelims);

  if I <> 0 then
    {find the end of the current word}
    while (I <= SLen) and not (S[I] in WordDelims) do
    begin
      {add the I'th character to result}
      Temp := Temp + S[i];
      Inc(I);
    end;

  ExtractWord := Temp;

end;

end.


Author: Andrey
Source: http://delphikingdom.com
ID: 03965