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