Classes of matrices and vectors in Delphi

This article discusses type design for working with objects of linear algebra: vectors, matrices, quaternions. Shows the classic use of the mechanism for overloading standard operations, using the “Copy On Write” method and annotations.

Working in the field of mathematical modeling, I often have to deal with computational algorithms that use operations on matrices, vectors, and quaternions. Surprisingly, he found that, despite the possibilities of a modern development environment, fellow programmers often use a procedural approach in solving such problems. So, to calculate the product of a matrix and a vector, types and functions like these are described:

TVec3 = array[1..3] of Extended; TMatrix3x3 = array[1..3, 1..3] of Extended; function MVMult(M: TMatrix3x3; V: TVec3): TVec3; 

It is proposed to use an object-based approach, which, in turn, involves the joint distribution of data and methods for their processing. Let's see what features Delphi provides for solving this class of tasks in an object form. Designing the structure of objects we will proceed from the following requirements:

• Version of the development environment Delphi XE7
• For numeric data, use the Extended type as the most accurate;
• Use dynamic arrays to store data, because sizes of vectors and matrices can be any and convenient to look at in the debugger;
• The elements of vectors and matrices are numbered with 1, the elements of quaternions with 0;
• For calculations, use the operations +, -, *, /;
• Provide the ability to transfer by value and copy vectors and matrices with the ": =" operation;
• To provide the possibility of automated initialization of vectors and matrices by predetermined sizes.

Design

For convenience, we define auxiliary types:
')

 TAbstractVector = array of Extended; TAbstractMatrix = array of array of Extended; 

Now we define the quaternion structure, vector and matrix:

 TQuaternion = record private FData: array[0..3] of Extended; procedure SetElement(Index: Byte; Value: Extended); function GetElement(Index: Byte): Extended; public property Element[Index: Byte]: Extended read GetElement write SetElement; default; end; TVector = record private FData: TAbstractVector; FCount: Word; procedure SetElement(Index: Word; Value: Extended); function GetElement(Index: Word): Extended; public constructor Create(ElementsCount: Word); property Count: Word read FCount; property Elements[Index: Word]: Extended read GetElement write SetElement; default; end; TMatrix = record private FData: TAbstractMatrix; FRowsCount: Word; FColsCount: Word; procedure SetElement(Row, Col: Word; Value: Extended); function GetElement(Row, Col: Word): Extended; public constructor Create(RowsCount, ColsCount: Word); property RowCount: Word read FRowsCount; property ColCount: Word read FColsCount; property Elements[Row, Col: Word]: Extended read GetElement write SetElement; default; end; 

We use exactly record , since overloading operations for class construction in Delphi is not allowed. In addition, record objects have a useful property — their data is deployed in memory at the place of declaration, in other words, the record object is not a reference to an instance in dynamic memory.
However, in our case, the elements of vectors and matrices will be stored in a dynamic array, the object of which is a reference. Therefore, it will be convenient to use explicit constructors. They perform the initialization of internal fields, allocating memory for the required number of elements:

 constructor TVector.Create(ElementsCount: Word); begin FCount := ElementsCount; FData := nil; SetLength(FData, FCount); end; constructor TMatrix.Create(RowsCount, ColsCount: Word); begin FRowsCount := RowsCount; FColsCount := ColsCount; FData := nil; SetLength(FData, FRowsCount, FColsCount); end; 

A quaternion at this stage does not require a constructor, since it stores data in a static array and is deployed in memory at the place of its declaration.
To access the elements, here are the properties-indexers, they are convenient to make default to omit the name. Access to the requested item occurs after checking its index for valid values. The implementation for TVector is shown:

 function TVector.GetElement(Index: Word): Extended; begin {$R+} Result := FData[Pred(Index)]; end; procedure TVector.SetElement(Index: Word; Value: Extended); begin {$R+} FData[Pred(Index)] := Value; end; 

At this stage, to create our objects, you will have to use the following code:

 var V: TVector; . . . V := TVector.Create(3); V[1] := 1; V[2] := 2; V[3] := 3; 

Practice has shown that it is useful to have the means to use a more concise syntax for creating a vector or matrix. To do this, add additional constructors, as well as implement the implicit casting operation, which will allow to overload the ": =".

 TQuaternion = record public . . . constructor Create(Q: TAbstractVector); class operator Implicit(V: TAbstractVector): TQuaternion; end; TVector = record public . . . constructor Create(V: TAbstractVector); overload; class operator Implicit(V: TAbstractVector): TVector; end; TMatrix = record public . . . constructor Create(M: TAbstractMatrix); overload; class operator Implicit(M: TAbstractMatrix): TMatrix; end; 

And implementation:

 constructor TQuaternion.Create(Q: TAbstractVector); begin if Length(Q) <> 4 then raise EMathError.Create(WRONG_SIZE); Move(Q[0], FData[0], SizeOf(FData)); end; class operator TQuaternion.Implicit(V: TAbstractVector): TQuaternion; begin Result.Create(V); end; constructor TVector.Create(V: TAbstractVector); begin FCount := Length(V); FData := Copy(V); end; class operator TVector.Implicit(V: TAbstractVector): TVector; begin Result.Create(V); end; constructor TMatrix.Create(M: TAbstractMatrix); var I: Integer; begin FRowsCount := Length(M); FColsCount := Length(M[0]); FData := nil; SetLength(FData, FRowsCount, FColsCount); for I := 0 to Pred(FRowsCount) do FData[I] := Copy(M[I]); end; class operator TMatrix.Implicit(M: TAbstractMatrix): TMatrix; begin Result.Create(M); end; 

Now, to create and initialize a vector or matrix, it suffices to write:

 var V: TVector; M: TMatrix; . . . V := [4, 5, 6]; // M := [[1, 2, 3], [4, 5, 6], [7, 8, 9]]; 

Operation overload

Here, for example, only the operation overload * for multiplying the matrix by the vector will be implemented. The remaining operations can be viewed in the file attached to the article. A full list of overload capabilities is here .

 TMatrix = record public . . . class operator Multiply(M: TMatrix; V: TVector): TVector; end; class operator TMatrix.Multiply(M: TMatrix; V: TVector): TVector; var I, J: Integer; begin if (M.FColsCount <> V.FCount) then raise EMathError.Create(WRONG_SIZE); Result.Create(M.FRowsCount); for I := 0 to M.FRowsCount - 1 do for J := 0 to M.FColsCount - 1 do Result.FData[I] := Result.FData[I] + M.FData[I, J] * V.FData[J]; end; 

The first argument of the Multiply () method is the matrix to the left of the * sign, the second argument is the column vector to the right of the * sign. The result of the product is a new vector whose object is created in the process of calculation. In case of mismatch of the number of columns of the matrix and the number of elements of the vector, an exception is raised. Here is the use of this operation in the program:

 var V, VResult: TVector; M: TMatrix; . . . VResult := M * V; 

It is convenient to use wrapper functions to construct anonymous vectors and matrices from array literals on the fly:

 function TVec(V: TAbstractVector): TVector; begin Result.Create(V); end; function TMat(M: TAbstractMatrix): TMatrix; begin Result.Create(M); end; function TQuat(Q: TAbstractVector): TQuaternion; begin Result.Create(Q); end; 

The use of wrappers is as follows. The equivalent of the expression from the previous example is shown:

  V := TMat([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) * TVec([4, 5, 6]); 

In addition to standard operations, it is useful to add specific methods to our object types, such as transposing or referencing. The following is an example of a matrix inversion (inversion) method. Despite its size, it is the fastest I've ever seen (in high-level languages).

 TMatrix = record public . . . function Inv: TMatrix; end; function TMatrix.Inv: TMatrix; var Ipiv, Indxr, Indxc: array of Integer; DimMat, I, J, K, L, N, ICol, IRow: Integer; Big, Dum, Pivinv: Extended; begin //  . if (FRowsCount <> FColsCount) then raise EMathError.Create(NOT_QUAD); Result := Self; DimMat := FRowsCount; SetLength(Ipiv, DimMat); SetLength(Indxr, DimMat); SetLength(Indxc, DimMat); IRow := 1; ICol := 1; for I := 1 to DimMat do begin Big := 0; for J := 1 to DimMat do if (Ipiv[J - 1] <> 1) then for K := 1 to DimMat do if (Ipiv[K - 1] = 0) then if (Abs(Result[J, K]) >= Big) then begin Big := Abs(Result[J, K]); IRow := J; ICol := K; end; Ipiv[ICol - 1] := Ipiv[ICol - 1] + 1; if (IRow <> ICol) then for L := 1 to DimMat do begin Dum := Result[IRow, L]; Result[IRow, L] := Result[ICol, L]; Result[ICol, L] := Dum; end; Indxr[I - 1] := IRow; Indxc[I - 1] := ICol; if Result[ICol, ICol] = 0 then raise EMathError.Create(SINGULAR); Pivinv := 1.0 / Result[ICol, ICol]; Result[ICol, ICol] := 1.0; for L := 1 to DimMat do Result[ICol, L] := Result[ICol, L] * Pivinv; for N := 1 to DimMat do if (N <> ICol) then begin Dum := Result[N, ICol]; Result[N, ICol] := 0.0; for L := 1 to DimMat do Result[N, L] := Result[N, L] - Result[ICol, L] * Dum; end; end; for L := DimMat downto 1 do if (Indxr[L - 1] <> Indxc[L - 1]) then for K := 1 to DimMat do begin Dum := Result[K, Indxr[L - 1]]; Result[K, Indxr[L - 1]] := Result[K, Indxc[L - 1]]; Result[K, Indxc[L - 1]] := Dum; end; end; 

Copy by value

Using dynamic arrays to store elements of vectors and matrices leads to the fact that when you try to copy them entirely in the receiving object (the one to the left of ": ="), a copy of the link to this dynamic array is created.
For example, an attempt to preserve the value of the matrix M after evaluating the expression will also result in the inversion of the MStore matrix.

 var M, MStore: TMatrix; . . . MStore := M; M := M.Inv; 

In order to correctly implement copying by value, we use the fact that, along with the negative offset from the address of the first element of the dynamic array, along with the length value, the counter of references to this array is stored. If the value of the counter is 0, then the memory manager frees this array. If the value of the counter is 1, this means that there is only one reference to the array instance in memory.
Therefore, when copying, we must analyze the value of the counter and, if it is greater than 1, then create a full-fledged copy of the array by copying it into the receiving object elementwise . Below is the function code, which returns True only in the case when the value of the reference counter passed in the input parameter of the dynamic array exceeds 1.

 {$POINTERMATH ON} function NotUnique(var Arr): Boolean; begin Result := (PCardinal(Arr) - 2)^ > 1; end; 

At what point should a full copy be performed? This is a rather expensive operation, so there is no point in performing it when accessing an element of a vector / matrix for reading. If we have at least a thousand references to the original, if he himself does not undergo any changes, then they all remain the same. Therefore, you need to copy only when accessing the item for writing. To do this, we modify the SetElement () methods for vectors and matrices, adding at the beginning a check for the uniqueness of an instance of the FData array:

 procedure TVector.SetElement(Index: Word; Value: Extended); begin {$R+} CheckUnique; FData[Pred(Index)] := Value; end; procedure TVector.CheckUnique; begin if NotUnique(FData) then FData := Copy(FData); end; procedure TMatrix.SetElement(Row, Col: Word; Value: Extended); begin {$R+} CheckUnique; FData[Pred(Row), Pred(Col)] := Value; end; procedure TMatrix.CheckUnique; var I: Integer; begin if NotUnique(FData) then begin FData := Copy(FData); for I := 0 to Pred(FRowsCount) do FData[i] := Copy(FData[i]); end; end; 

Thus, when you try to change the value of an element, a check will be made on the uniqueness of the link, and, if it is not confirmed, an element-by-element copy will be created, to which the change will be made.

Annotations and automatic initialization

Elements of vectors and matrices should be accessed only after allocating memory for them. Element values ​​are stored in a dynamic array, the dimensions of which are set in the object's constructor. An implicit constructor call can occur during object initialization, or in the process of evaluating an expression.

 var V: TVector; M: TMatrix; begin // V[1] := 1; // :    V := TVector.Create(4); //    M := TMatrix.Create(4, 4); //    // V := [1, 0, 0, 0]; //    // V := M * TVec([1, 0, 0, 0]); //    V[1] := 1; //    :   

The use of implicit constructors can lead to errors when, sooner or later, it will be allowed to refer to an element of an uncreated object. According to the rules of good tone, the constructor should be called explicitly.
But what if there are hundreds and thousands of vectors and matrices in our program? Consider the description of a class that uses a set of vectors and matrices as its fields.

 TMovement = record R: TVector; V: TVector; W: TVector; Color: TVector; end; TMovementScheme = class private FMovement: array[1..100] of TMovement; FOrientation: TMatrix; end; 

It is required to develop a method for automated initialization of all fields of the TVector and TMatrix types: to allocate memory for vectors and matrices in accordance with the necessary element numbers and sizes. The annotation mechanism (or attributes, in Delphi terms) will help us in this - a tool that allows you to add types with arbitrary metadata. So, for each vector, the number of its elements must be known in advance, for the matrix - the number of rows and columns.
Create a class that encapsulates data about dimensions, according to the rules for creating attribute classes.

 TDim = class(TCustomAttribute) private FRowCount: Integer; FColCount: Integer; public constructor Create(ARowCount: Integer; AColCount: Integer = 0); overload; property RowCount: Integer read FRowCount; property ColCount: Integer read FColCount; end; constructor TDim.Create(ARowCount: Integer; AColCount: Integer = 0); begin FRowCount := ARowCount; FColCount := AColCount; end; 

The constructor gets the number of rows and columns, and in the case of a vector, we can do with only the number of rows. Now we add the definition of the types from the previous listing with new annotations:

 TMovement = record [TDim(3)] R: TVector; [TDim(3)] V: TVector; [TDim(3)] W: TVector; [TDim(4)] Color: TVector; end; TMovementScheme = class private FMovement: array[1..100] of TMovement; [TDim(3, 3)] FOrientation: TMatrix; end; 

Below is the code that initializes objects of type TVector and TMatrix based on information taken from attributes.

 procedure Init(Obj, TypeInfoOfObj: Pointer; Offset: Integer = 0); const DefaultRowCount = 3; DefaultColCount = 3; VectorTypeName = 'TVector'; MatrixTypeName = 'TMatrix'; var RTTIContext: TRttiContext; Field : TRttiField; ArrFld: TRttiArrayType; I: Integer; Dim: TCustomAttribute; RowCount, ColCount: Integer; OffsetFromArray: Integer; begin for Field in RTTIContext.GetType(TypeInfoOfObj).GetFields do begin if Field.FieldType <> nil then begin RowCount := DefaultRowCount; ColCount := DefaultColCount; for Dim in Field.GetAttributes do begin RowCount := (Dim as TDim).RowCount; ColCount := (Dim as TDim).ColCount; end; if Field.FieldType.TypeKind = tkArray then begin ArrFld := TRttiArrayType(Field.FieldType); if ArrFld.ElementType.TypeKind = tkRecord then begin for I := 0 to ArrFld.TotalElementCount - 1 do begin OffsetFromArray := I * ArrFld.ElementType.TypeSize; if ArrFld.ElementType.Name = VectorTypeName then PVector(Integer(Obj) + Field.Offset + OffsetFromArray + Offset)^ := TVector.Create(RowCount) else if ArrFld.ElementType.Name = MatrixTypeName then PMatrix(Integer(Obj) + Field.Offset + OffsetFromArray + Offset)^ := TMatrix.Create(RowCount, ColCount) else Init(Obj, ArrFld.ElementType.Handle, Field.Offset + OffsetFromArray); end; end; end else if Field.FieldType.TypeKind = tkRecord then begin if Field.FieldType.Name = VectorTypeName then PVector(Integer(Obj) + Field.Offset + Offset)^ := TVector.Create(RowCount) else if Field.FieldType.Name = MatrixTypeName then PMatrix(Integer(Obj) + Field.Offset + Offset)^ := TMatrix.Create(RowCount, ColCount) else Init(Obj, Field.FieldType.Handle, Field.Offset) end; end; end; end; 

The Init () procedure receives as input the address of the container object and its RTTI data. Further, a recursive traversal of all container fields occurs, and for all counter fields with the names of the types “TVector” and “TMatrix” their constructors will be explicitly called.
Let's modify the TMovementScheme class using the Init () procedure:

 TMovementScheme = class . . . public constructor Create; end; constructor TMovementScheme.Create; begin Init(Self, Self.ClassInfo); end; 

Init () call option for arbitrary writing:

 var Movement: TMovement; . . . Init(@Movement, TypeInfo(TMovement)); 

By default, Init () creates vectors with three elements, and 3x3 matrices, so in the declaration of the types TMovement and TMovementScheme, the attributes [TDim (3)] and [TDim (3, 3)] can be omitted, leaving only [TDim (4) ].

Attached to the article is a file in which the implementation of the described ideas is given in full.