Sources
Delphi Russian Knowledge Base
DRKB - это самая большая и удобная в использовании база знаний по Дельфи в рунете, составленная Виталием Невзоровым

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

01.01.2007
Автор: Andrey

WEB-сайт: http://delphikingdom.com

В этом модуле «осели» все операции с матрицами и векторами, которые я использовал для работы. Но есть алгоритмы, которые многие, наверняка, увидят впервые: 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 )

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

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;
 
   F := 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
   FMin := HMinZero.FMin( SELF, A, B );
end;
 
function Spline.FMax( A, B : Real ) : Real;
begin
   FMax := HMinZero.FMax( SELF, A, B );
end;
 
function Spline.ZeroIn( A, B : Real ) : Real;
begin
   ZeroIn := 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
         { Рекурсия - если точность плохая }
         Int := Int( x1, x2, x3, f1, f2, f3, Eps )
            + Int( x3, x4, x5, f3, f4, f5, Eps )
      else
 
         { За решения принять значение формулы 5-го порядка }
         Int := s3;
   end;
 
begin
 
   Integr := 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 );
 
   Integr := Int( x1, x3, x5, F1, F3, F5, Eps );
 
end;
 
function Spline.Integral( A, B : Real ) : Real;
begin
   Integral := Integr( A, B );
end;
 
function Spline.Integral : Real;
begin
   Integral := 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
   PX := Coeff.Poly( xr );
end;
 
function PolyApprox.DPX( xr : Real ) : Real;
{
   Производная полинома в "xr"
}
begin
   DPX := 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
   Difference := SELF.Norm - A.Norm;
end;
 
function Vector.Difference( x : Vector ) : Real;
begin
   Difference := 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 )};
   Asimm := 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;
 
   Poly := 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 ];
 
   Max := 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 ];
 
   Min := 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;
 
   MaxIndex := 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;
 
   MinIndex := 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;
 
   Cond := 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;
 
   Norm := 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
 
            SplineIndex := 1;
 
            if X <= d1.E( 1 )^ then
               exit;
 
            SplineIndex := NPoint - 1;
 
            if X >= d1.E( NPoint )^ then
               exit;
 
            IndRecurs( 1, NPoint );
            SplineIndex := 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 );
 
         F := 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 );
      DF := 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 ) );
 
         DF := 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
   Addr := VAddr.Ptr;
   //Addr := VirtualToPtr( VAddr );
end;
 
function Matrix.Addr : MatrixVPrimPtr;
begin
   Addr := 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
      E := 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
      E := 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 ];
 
   Summ := 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 ] );
 
   SummSqr := 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
   Addr := 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 ) ) );
   Alpha := 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.
 

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

unit HNVString;
 
interface
 
{*********************************************************}
{*                  OPSTRING.PAS 1.20                    *}
{*     Copyright (c) TurboPower Software 1987, 1992.     *}
{*                 All rights reserved.                  *}
{*********************************************************}
 
Type
   CharSet =
      Set of Char;
 
function WordCount(S : string; WordDelims : CharSet) : Byte;
    {-Given a set of word delimiters, return number of words in S}
 
function WordPosition(N : Byte; S : string; WordDelims : CharSet) : Byte;
    {-Given a set of word delimiters, return start position of N'th word in S}
 
function ExtractWord(N : Byte; S : string; WordDelims : CharSet) : string;
    {-Given a set of word delimiters, return the N'th word in S}
 
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.

https://delphiworld.narod.ru/

DelphiWorld 6.0