unit xVectRN;
{ Feel free to use this unit as you want.
  First developped by Bernard Gisin

  This unit helps manipulating real vectors in a N-dimensional space.
  use  xMatRN.pp  to do matrix manipulation

  Unfortunately, the assignment operator  ":="
  cannot be not well defined.
  I would like the assignment operator := to be defined as any comparaison
  operator, but without returning any results. For example :
  operator := (var rOut:RVectRNFloat; rIn:RVectRNFloat);
  rOut.blabla:=rIn.blabla
  Thus I will use the  "<="  operator as an assignment operator,
  ignoring the boolean result.

  If you add new procedure or function using vectors  RVectRNFloat,
  If you pass by value these vectors, then at the end of the procedure,
  for each vector passed by value you should add the line :
  if (vecV.fIsTemp) then vecV.pOrg^.fFreeToUse:=TRUE;
  // If  vecV was a temporary result. It is free to use now
  ----------------------------------------------------------------}

INTERFACE
//"""""""

uses wGlobals2;

const
  gglCTEVecRNIndMaxFamilly = 10;
  // Up to 10 (from 1 to 10)  different length of vectors can be defined,
  // to allocate each time the right amount of memory for each vector.
  // Then a shorter length can also be given to the vector.
  // I think that usually, only one or two familly will be defined.
  // Two if the length difference between the two groupes of vectors
  // is important.

  gglwVectRNError : word32 = 0;
  // If a division by zero has been tempted, this value is incremented by 1.
  // It should be set to zero, before trying to test it.

  gglvVNEpsilon : Float = 1.0e-12;
  // Small number. Number having absolute value under it can be considered as zero.
  // normally, the number taken as limite should be base on  vMachEps, defined
  // as the smallest number such that  1.0 + vMachEps <> 1.0

type
  PRVectRNFloat = ^RVectRNFloat;
  RVectRNFloat = record
  //===================
  bFamilly : byte;  // in [1..gglCTEVecRNIndMaxFamilly] ( =[1..10] )
  // Up to gglCTEVecRNIndMaxFamilly different kind of vectors can be defined.
  // gglanVectRNIndMax[bFamilly] determine that the memory allocation
  // has been done for index from 1 to gglanVectRNIndMax[bFamilly]
  // Do not change it.

  fIsTemp : boolean;
  // Normally is FALSE.
  // Is set to TRUE for temporary vectors.

  nIndMax : int32;
  // The vector has length nIndMax

  pavV : PTypeArrayFloat;
  // The component of the vector.
  // pavV[1..nIndMax] is the vector.

// Private
  // Only used for vectors which are used for temporary storage
  // Normally, for internal use.
  // for intermediate results.
  fFreeToUse : boolean;
  // defines if this temporary can be used, or is holding an intermediate
  // result.

  pvecNext : PRVectRNFloat;
  // Pointer on the next temporary vector of the same familly

  pOrg : PRVectRNFloat;
  // Pointer on the original vector gglapvecRNTemp  used by this temporary vector
  end; // RVectRNFloat = record

type
  PArrayRVectRNFloat = ^ArrayRVectRNFloat;
  ArrayRVectRNFloat = array[0..100000000] of RVectRNFloat;

var
  gglanVectRNIndMax : array[1..gglCTEVecRNIndMaxFamilly] of int32;
  // Memorise the maximum index allocated in memory for each familly of vectors

  gglavecRNTemp : array[1..gglCTEVecRNIndMaxFamilly] of RVectRNFloat;
  // When a familly is created, one temporary vector is created by default.
  // If more are needed, they will be allocated automatically

procedure VNCreateNewVecRNFamilly(bFamillyP:byte; nIndMaxP:int32);
//================================================================
// This creates a new familly of vectors, which has all the same
// amount of memory allocated.
// This must be called before using the vectors.
// bFamilly, in [1..10] define which groupe of vectors we want to define.
// for example,  bFamilly=1  for vectors of 10 elements, and
// bFamilly=2  for vectors of 1000 elements.
// nIndMaxP  defines the maximum index allocated for  pavV
// in practice  pavV^[0] is allocated, but not used.

procedure VNDisposeVecRNFamilly(bFamillyP:byte);
//==============================================
// If  CreateNewVecRNFamilly(bFamillyP)  has been called then memory has been
// allocated, and must be disposed by calling this procedure.

procedure VNInit(var vecV:RVectRNFloat; bFamillyP:byte);
// Any new vector should be initialised with this procedure,
// after  CreateNewVecRNFamilly(bFamillyP)  has been called once
// for the  bFamillyP

procedure VNDone(var vecV:RVectRNFloat);
// When a vector is no more needed, this should be called, to
// free the memory allocated for this vector.

procedure VNSetToZero(var vecV:RVectRNFloat; nIndMaxP:int32);
// Define that the vector has length  nIndMaxP
// and fills it with zeros.

procedure VNSetToBasis(var vecV:RVectRNFloat; nIndMaxP,nBasisIndexP:int32);
// Define that the vector has length  nIndMaxP
// and fills it with zeros, except the nBasisIndexP which is set to 1.
// This defines the  nBasisIndexP  vector of the standard basis of
// the  Reals^nIndMaxP

procedure VNAssign( var vecV:RVectRNFloat; nIndMinP,nIndMaxP:int32;
                    vVal1,vVal2,vVal3,vVal4,vVal5:Float);
// Assign the indexes nIndMinP to nIndMaxP of the vector to the values
// vVal1 to vVal5.
// If nIndMaxP >= nIndMinP+5,
// then the values from nIndMinP+5 to nIndMaxP are set to 0

procedure VNAssign( var vecV:RVectRNFloat; nIndMinP,nIndMaxP:int32;
                    vVal1,vVal2,vVal3,vVal4,vVal5,
                    vVal6,vVal7,vVal8,vVal9,vVal10:Float);
// Assign the indexes nIndMinP to nIndMaxP of the vector to the values
// vVal1 to vVal10.
// If nIndMaxP >= nIndMinP+10,
// then the values from nIndMinP+10 to nIndMaxP are set to 0

procedure VNCopyFromTo(var vecSrc,vecRes:RVectRNFloat);
// Copy the content of vecSrc into the vector vecRes

function Norme2(vecV:RVectRNFloat):Float;
// Square of the norme of the vector

function Norme(vecV:RVectRNFloat):Float;
// Norme of the vector

procedure VNNormalise(var vecV:RVectRNFloat);
// Normalise the vector
// If the vector = 0  then return  (1,0,0,...)

function VNAngleWith(vecV1, vecV2:RVectRNFloat):Float;
// Angle in radian between the two vectors. in [0..Pi]
// No sign can be given, since no orientation is given.

function VNVectToStr(vecV:RVectRNFloat; nInd1,nInd2:int32):ShortString;
// Convert the indexes from  nInd1  to  nInd2  of the vector into a string.
// usefull to display the vector in a debugging window.

//operator := (vecIn:RVectRNFloat) vecOut:RVectRNFloat;
// The assignment operator doesn't work as I would like.
// I would like the assignment operator := to be defined as any comparaison
// operator, but without returning any results. For example :
// operator := (var rOut:RVectRNFloat; rIn:RVectRNFloat);
// rOut.blabla:=rIn.blabla
// Thus I will use the  <=  operator as an assignment operator,
// ignoring the boolean result.

procedure SetTemporaryVect(var vecV1,vecV2,vecV3:RVectRNFloat);
//============================================================
// vecV3 = temporary vector.
// free the vectors vecV1 and vecV2 if they are temporary
// Used in vector additions, substraction and mult. by scalar.
// Normally not used by an external unit, but used by  xMatrixRN.pp

operator <= (var vecRes:RVectRNFloat; vecSrc:RVectRNFloat) fRes:boolean;
//======================================================================
// It is used as an assignment operator, vecRes:=vecSrc,
// ignoring the result.
// The :  operator := (vecSrc:TVectRNFloat) vecRes:TVectRNFloat;
// cannot be redefined, and should not be used to assign one vector to another.
// DO NOT USE vec1:=vec2
//---------------*****----
// vec1:=vec2;  would assign the pointer  vec2.pavV  to  pvec1.pavV
// the  vec1.pavV  pointer would be lost.
// If vecSrc is the result of an operation which has been stored in
// gglapvecRNTemp,  then this exanges the two pointers
// vecSrc.pavV and vecRes.pavV
// Else the content of vecSrc is copied in vecRes, but not its bFamilly.

operator + (vecV1, vecV2:RVectRNFloat) vecV3:RVectRNFloat;
operator - (vecV1, vecV2:RVectRNFloat) vecV3:RVectRNFloat;
operator * (vecV1, vecV2:RVectRNFloat) vScalarProduct:Float; // scalar product
operator * (vScal:Float; vecIn:RVectRNFloat) vecOut:RVectRNFloat; // Scalar time a Vector

operator = (vecV1, vecV2:RVectRNFloat) fRes:boolean;
// True if Norme(pvecV1 - pvecV2) <= gglvVNEpsilon
// Use  not (vecV1 = vecV2)  to test that  vecV1 <> vecV2

IMPLEMENTATION
//""""""""""""

Uses wUtiles0, wDebug0;

const fFirstInit : boolean = TRUE; // Should the first initialisation be done ?

var nNotUsed : integer; // only to avoid an error message while compliling in DevPascal 1.9

procedure VNCreateNewVecRNFamilly(bFamillyP:byte; nIndMaxP:int32);
//================================================================
var nn : int32;
begin
// For the first init, inidcate that the temporary vectors gglavecRNTemp
// has not been initialised
if fFirstInit then begin
   for nn:=1 to gglCTEVecRNIndMaxFamilly do begin
      gglavecRNTemp[nn].pvecNext :=NIL;
      gglavecRNTemp[nn].pOrg     :=NIL;
      gglanVectRNIndMax[nn]      :=0;
      gglavecRNTemp[nn].pavV     :=NIL; // Indicate that not init. yet
      end;
   fFirstInit:=FALSE;
   end;

// If the familly has already be defined, must do some checking
if (gglanVectRNIndMax[bFamillyP] = nIndMaxP) then exit;
                       // The familly has alredy be defined, nothing to do

if (gglanVectRNIndMax[bFamillyP] > 0) then begin
   // The familly has alredy be defined, with another size.
   // Desallocate the memory, to reallocate it with the new size
   VNDisposeVecRNFamilly(bFamillyP);
   end;

// no range checking...
// Memorise the memoy which has to be allocated for a new vector
// belonging to  bFamillyP
gglanVectRNIndMax[bFamillyP]:=nIndMaxP;

// Vector of length zero are useless. => no memory allocation
if (nIndMaxP = 0) then exit;

// Allocation for a temporary vector.
// More temporary vectors will be created, if necessary
with gglavecRNTemp[bFamillyP] do begin
   bFamilly   :=bFamillyP;
   fIsTemp    :=FALSE;
   nIndMax    :=0;
   fFreeToUse :=TRUE;
   pvecNext   :=NIL;
   pOrg       :=@gglavecRNTemp[bFamillyP];
   GetMem(pavV, (gglanVectRNIndMax[bFamilly]+1)*SizeOf(Float));
   end;
end;  // VNCreateNewVecRNFamilly

procedure VNDisposeVecRNFamilly(bFamillyP:byte);
//==============================================
// Dispose all the temporary vectors of  bFamillyP
var
  pvecV,
  pvecVLast : PRVectRNFloat; // used to search a free temp. vector
begin
// Get the first vector of the familly
pvecV:=@gglavecRNTemp[bFamillyP];

// Loop on all chained vectors of the same familly
while (pvecV <> NIL) do begin
   pvecVLast :=pvecV;
   pvecV     :=pvecV^.pvecNext;
   FreeMem(pvecVLast^.pavV, (gglanVectRNIndMax[bFamillyP]+1)*SizeOf(Float));
   end;

// Just to indicate that this vector has not been initialised
gglavecRNTemp[bFamillyP].pOrg :=NIL;
gglavecRNTemp[bFamillyP].pavV :=NIL;
gglanVectRNIndMax[bFamillyP]  :=0; // this familly has not been defined
end;  // VNDisposeVecRNFamilly

procedure VNInit(var vecV:RVectRNFloat; bFamillyP:byte);
//======================================================
begin
vecV.bFamilly:=bFamillyP;
GetMem(vecV.pavV, (gglanVectRNIndMax[bFamillyP]+1)*SizeOf(Float));

vecV.fIsTemp    :=FALSE;
vecV.nIndMax    :=0;;
vecV.fFreeToUse :=FALSE;
vecV.pvecNext   :=NIL;
vecV.pOrg       :=NIL;
end;  //  VNInit

procedure VNDone(var vecV:RVectRNFloat);
//======================================
begin
FreeMem(vecV.pavV, (gglanVectRNIndMax[vecV.bFamilly]+1)*SizeOf(Float));
vecV.pavV:=NIL;
end;  // VNDone

procedure VNSetToZero(var vecV:RVectRNFloat; nIndMaxP:int32);
//===========================================================
// Define that the vector has length  nIndMaxP
// and fills it with zeros.
var nn : int32;
begin
vecV.nIndMax:=nIndMaxP;
for nn:=1 to nIndMaxP do vecV.pavV^[nn]:=0.0;
end;  // VNSetToZero

procedure VNSetToBasis(var vecV:RVectRNFloat; nIndMaxP,nBasisIndexP:int32);
//=========================================================================
var nn : int32;
begin
vecV.nIndMax:=nIndMaxP;
for nn:=1 to nIndMaxP do vecV.pavV^[nn]:=0.0;
if (nBasisIndexP <= nIndMaxP) then vecV.pavV^[nBasisIndexP]:=1.0;
end;  // VNSetToBasis

procedure VNAssign( var vecV:RVectRNFloat; nIndMinP,nIndMaxP:int32;
                    vVal1,vVal2,vVal3,vVal4,vVal5:Float);
//=======================================================
var nn : int32;
begin
vecV.pavV^[nIndMinP]:=vVal1;
if (nIndMinP+1 <= nIndMaxP) then vecV.pavV^[nIndMinP + 1]:=vVal2;
if (nIndMinP+2 <= nIndMaxP) then vecV.pavV^[nIndMinP + 2]:=vVal3;
if (nIndMinP+3 <= nIndMaxP) then vecV.pavV^[nIndMinP + 3]:=vVal4;
if (nIndMinP+4 <= nIndMaxP) then vecV.pavV^[nIndMinP + 4]:=vVal5;
for nn:=nIndMinP+5 to nIndMaxP do vecV.pavV^[nn]:=0.0;
end;  // VNAssign

procedure VNAssign( var vecV:RVectRNFloat; nIndMinP,nIndMaxP:int32;
{=================} vVal1,vVal2,vVal3,vVal4,vVal5,
                    vVal6,vVal7,vVal8,vVal9,vVal10:Float);
var nn : int32;
begin
vecV.pavV^[nIndMinP]:=vVal1;
if (nIndMinP+1 <= nIndMaxP) then vecV.pavV^[nIndMinP + 1]:=vVal2;
if (nIndMinP+2 <= nIndMaxP) then vecV.pavV^[nIndMinP + 2]:=vVal3;
if (nIndMinP+3 <= nIndMaxP) then vecV.pavV^[nIndMinP + 3]:=vVal4;
if (nIndMinP+4 <= nIndMaxP) then vecV.pavV^[nIndMinP + 4]:=vVal5;
if (nIndMinP+5 <= nIndMaxP) then vecV.pavV^[nIndMinP + 5]:=vVal6;
if (nIndMinP+6 <= nIndMaxP) then vecV.pavV^[nIndMinP + 6]:=vVal7;
if (nIndMinP+7 <= nIndMaxP) then vecV.pavV^[nIndMinP + 7]:=vVal8;
if (nIndMinP+8 <= nIndMaxP) then vecV.pavV^[nIndMinP + 8]:=vVal9;
if (nIndMinP+9 <= nIndMaxP) then vecV.pavV^[nIndMinP + 9]:=vVal10;
for nn:=nIndMinP+10 to nIndMaxP do vecV.pavV^[nn]:=0.0;
end;  // VNAssign

procedure VNCopyFromTo(var vecSrc,vecRes:RVectRNFloat);
//=====================================================
// Copy the content of vecSrc into the vector vecRes
var nn : int32;
begin
// Copy the content of the VecIn vector in the vecOut vector
for nn:=1 to vecSrc.nIndMax do vecRes.pavV^[nn]:=vecSrc.pavV^[nn];
vecRes.nIndMax:=vecSrc.nIndMax;

// If vecSrc was passed by reference, then the following test should be done
// If  vecV was a temporary result. It is free to use now
//if (vecSrc.fIsTemp) then vecSrc.pOrg^.fFreeToUse:=TRUE;
end; // VNCopyFromTo

function Norme2(vecV:RVectRNFloat):Float;
//=======================================
// Return the square of the norme.
var
  nn          : int32;
  vNorme2     : Float;
begin
vNorme2:=0;
for nn:=1 to vecV.nIndMax do vNorme2:=vNorme2 + sqr(vecV.pavV^[nn]);
Norme2:=vNorme2;

// If  vecV was a temporary result. It is free to use now
if (vecV.fIsTemp) then vecV.pOrg^.fFreeToUse:=TRUE;
end;  // Norme2

function Norme(vecV:RVectRNFloat):Float;
//======================================
// Return the norme.
var
  nn          : int32;
  vNorme2     : Float;
begin
vNorme2:=0;
for nn:=1 to vecV.nIndMax do vNorme2:=vNorme2 + sqr(vecV.pavV^[nn]);
Norme:=sqrt(vNorme2);

// If  vecV was a temporary result. It is free to use now
if (vecV.fIsTemp) then vecV.pOrg^.fFreeToUse:=TRUE;
end;  // Norme

procedure VNNormalise(var vecV:RVectRNFloat);
//===========================================
// Normalise the vector  vecV
// If  vecV = 0  then return  vecV = (1,0,0,...)
var
  nn          : int32;
  vNorme      : Float;
begin
vNorme:=Norme(vecV); // Return the norme of the vector

if (vNorme <= 1e-15) then begin
   inc(gglwVectRNError); // signal that a division by zero has been tempted
   for nn:=2 to vecV.nIndMax do vecV.pavV^[nn]:=0.0;
   vecV.pavV^[1]:=1.0;
   exit;
   end;

for nn:=1 to vecV.nIndMax do vecV.pavV^[nn]:=vecV.pavV^[nn] / vNorme;

// If  vecV was a temporary result. It is free to use now
if (vecV.fIsTemp) then vecV.pOrg^.fFreeToUse:=TRUE;
end;  // VNNormalise

function VNAngleWith(vecV1, vecV2:RVectRNFloat):Float;
//========================================================
// Angle in radian between the two vectors. in [0..Pi]
// No sign can be given, since no orientation is given.
var
  vNormeProd : Float;
begin
vNormeProd:=Norme(vecV1) * Norme(vecV2);

if (vNormeProd <= 1e-30) then begin
   // Angle not defines
   inc(gglwVectRNError); // signal that a division by zero has been tempted
   VNAngleWith:=0;
   end;

VNAngleWith:=ArcCos( (vecV1 * vecV2) / vNormeProd);

// If  vecV1 was a temporary result. It is free to use now
// If  vecV2 was a temporary result. It is free to use now
if (vecV1.fIsTemp) then vecV1.pOrg^.fFreeToUse:=TRUE;
if (vecV2.fIsTemp) then vecV2.pOrg^.fFreeToUse:=TRUE;
end;  // VNAngleWith

function VNVectToStr(vecV:RVectRNFloat; nInd1,nInd2:int32):ShortString;
//=====================================================================
// Convert the indexes from  nInd1  to  nInd2  of the vector into a string.
// usefull to display the vector in a debugging window.
var nn : int32;
begin
VNVectToStr:='(' + FloatToStr(vecV.pavV^[nInd1],8,3);

for nn:=nInd1+1 to nInd2 do
   VNVectToStr:=VNVectToStr + ',' + FloatToStr(vecV.pavV^[nn],8,3);

VNVectToStr:=VNVectToStr + ')';

// If  vecV was a temporary result. It is free to use now
if (vecV.fIsTemp) then vecV.pOrg^.fFreeToUse:=TRUE;
end;  // TVectRNFloat.VectToStr

//****************************************************************
// Operators definitions
//****************************************************************

procedure SetTemporaryVect(var vecV1,vecV2,vecV3:RVectRNFloat);
//=============================================================
// vecV3 = temporary vector.
// free the vectors vecV1 and vecV2 if they are temporary
var
  pvecNew,  // If a new temporary vector should be created
  pvecVLast : PRVectRNFloat; // used to search a free temp. vector
begin
// Initialise the new vector vecV3
vecV3.bFamilly  :=vecV1.bFamilly;  // = vecV2.bFamilly  normally
vecV3.fIsTemp   :=TRUE;

// if vecV1 or vecV2 are temporary vectors then they are no more used
if (vecV1.fIsTemp) then begin
   // Use vecV1 for the temporary result
   vecV3.pavV   :=vecV1.pavV;
   vecV3.pOrg   :=vecV1.pOrg; // Pointer on  gglapvecRNTemp[..] or one of the list

   if (vecV2.fIsTemp) then begin
      // vecV2 was also a temporary result. Its free to use now
      vecV2.pOrg^.fFreeToUse:=TRUE;

      // it could be that vecV2 = vecV1 = vecV3, thus this is necessary
      vecV3.pOrg^.fFreeToUse:=FALSE;
      end;
   end
else if (vecV2.fIsTemp) then begin
   // Use vecV2 for the temporary result
   vecV3.pavV  :=vecV2.pavV;
   vecV3.pOrg  :=vecV2.pOrg; // Pointer on  gglapvecRNTemp[..]  or one of the list
   end
else if gglavecRNTemp[vecV1.bFamilly].fFreeToUse then begin
   vecV3.pavV  :=gglavecRNTemp[vecV1.bFamilly].pavV;
   vecV3.pOrg  :=@gglavecRNTemp[vecV1.bFamilly];
   gglavecRNTemp[vecV1.bFamilly].fFreeToUse:=FALSE;
   end
else begin
   // Need a temporary vector to store the result
   // look for the first temporary free vector of the familly
   pvecVLast :=@gglavecRNTemp[vecV1.bFamilly];
   pvecNew   :=pvecVLast^.pvecNext;
   while (pvecNew <> NIL) and (not pvecNew^.fFreeToUse) do begin
         pvecVLast :=pvecNew;
         pvecNew   :=pvecNew^.pvecNext;
         end;

   if (pvecNew = NIL) then begin
      // No free temporary vector, a new one must be created
      // This is seldom necessary
      new(pvecNew);
      VNInit(pvecNew^, vecV1.bFamilly);
//      pvecNew^.fIsTemp    :=FALSE; used by temporary vectors, but is not temporary
      pvecNew^.nIndMax    :=0;
      pvecNew^.fFreeToUse :=FALSE;
      pvecNew^.pvecNext   :=NIL;
      pvecNew^.pOrg       :=pvecNew;

      // Add the new vector to the chained list of temporary vectors of this familly
      pvecVLast^.pvecNext:=pvecNew;
      end; // if (pvecNew = NIL) then begin

   // Use for vecV3, the array of the temporary vector pavecNew
   vecV3.pavV  :=pvecNew^.pavV;
   vecV3.pOrg  :=pvecNew;

   // Mark that  pvecNew  is not free to use
   pvecNew^.fFreeToUse:=FALSE;
   end; // if .. Need a temporary vector to store the result
end;  // SetTemporaryVect

operator <= (var vecRes:RVectRNFloat; vecSrc:RVectRNFloat) fRes:boolean;
//======================================================================
// This defines the assignment operator, since the normal assignment
// operator cannot be defined.
begin
// For assignment of the kind :
// glf:=vec1 <= vec2 + vec3, the test will be TRUE.
// glf is a side effect product to ignore.
if (vecSrc.fIsTemp) then begin
   // we know that the  vecSrc  vector use the data of
   // vecSrc.pOrg^.pavV
   // Exchange simply  vecSrc.pOrg^.pavV  with  vecRes.pavV
   // Exchange also the familly, in case they are different
   vecSrc.pOrg^.pavV     :=vecRes.pavV;
   vecRes.pavV           :=vecSrc.pavV;
   vecSrc.pOrg^.bFamilly :=vecRes.bFamilly;
   vecRes.bFamilly       :=vecSrc.bFamilly;
   vecRes.nIndMax        :=vecSrc.nIndMax;

   // It is free for another use
   vecSrc.pOrg^.fFreeToUse:=TRUE;

   // vecSrc itself is a temporary variable, which will disapear.

   fRes:=FALSE; // not important
   exit;
   end;

// Here we copy simply one vector into another

// Copy the content of the vecSrc vector in the vecRes vector
VNCopyFromTo(vecSrc, vecRes);
fRes:=TRUE; // not important
end; // operator <=

operator + (vecV1, vecV2:RVectRNFloat) vecV3:RVectRNFloat;
//========================================================
var
   nn,nMin,nMax : int32;
begin
// vecV3 = temporary vector.
// free the vectors vecV1 and vecV2 if they are temporary
SetTemporaryVect(vecV1,vecV2,vecV3);

// Selecting the shortest vector length
if (vecV1.nIndMax < vecV2.nIndMax) then begin
   nMin:=vecV1.nIndMax;
   nMax:=vecV2.nIndMax;

   // The shortest vector is taken as if it was completed with zeros.
   for nn:=nMin+1 to nMax do vecV3.pavV^[nn]:=vecV2.pavV^[nn];
   end
else begin
   nMin:=vecV2.nIndMax;
   nMax:=vecV1.nIndMax;

   // The shortest vector is taken as if it was completed with zeros.
   for nn:=nMin+1 to nMax do vecV3.pavV^[nn]:=vecV1.pavV^[nn];
   end;

// Addition of the vector
for nn:=1 to nMin do vecV3.pavV^[nn]:=vecV1.pavV^[nn] + vecV2.pavV^[nn];

vecV3.nIndMax  :=nMax; // Length of the resulting vector
end; // operator +

operator - (vecV1, vecV2:RVectRNFloat) vecV3:RVectRNFloat;
//========================================================
var
   nn,nMin,nMax : int32;
begin
// vecV3 = temporary vector.
// free the vectors vecV1 and vecV2 if they are temporary
SetTemporaryVect(vecV1,vecV2,vecV3);

// Selecting the shortest vector length
if (vecV1.nIndMax < vecV2.nIndMax) then begin
   nMin:=vecV1.nIndMax;
   nMax:=vecV2.nIndMax;

   // The shortest vector is taken as if it was completed with zeros.
   for nn:=nMin+1 to nMax do vecV3.pavV^[nn]:=-vecV2.pavV^[nn];
   end
else begin
   nMin:=vecV2.nIndMax;
   nMax:=vecV1.nIndMax;

   // The shortest vector is taken as if it was completed with zeros.
   for nn:=nMin+1 to nMax do vecV3.pavV^[nn]:=vecV1.pavV^[nn];
   end;

// Addition of the vector
for nn:=1 to nMin do vecV3.pavV^[nn]:=vecV1.pavV^[nn] - vecV2.pavV^[nn];

vecV3.nIndMax:=nMax; // Length of the resulting vector
end; // operator -

operator * (vecV1, vecV2:RVectRNFloat) vScalarProduct:Float; // scalar product
//==========================================================
var
   nn,nMin : int32;
begin
// Selecting the shortest vector length
if (vecV1.nIndMax < vecV2.nIndMax) then nMin:=vecV1.nIndMax
else                                    nMin:=vecV2.nIndMax;

// Computing the scalar product
vScalarProduct:=0;
for nn:=1 to nMin do
   vScalarProduct:=vScalarProduct + vecV1.pavV^[nn] * vecV2.pavV^[nn];

// If  vecV1 was a temporary result. It is free to use now
// If  vecV2 was a temporary result. It is free to use now
if (vecV1.fIsTemp) then vecV1.pOrg^.fFreeToUse:=TRUE;
if (vecV2.fIsTemp) then vecV2.pOrg^.fFreeToUse:=TRUE;
end; // operator *

operator * (vScal:Float; vecIn:RVectRNFloat) vecOut:RVectRNFloat;
//===============================================================
// Scalar time a Vector
var nn : int32;
begin
// vecV3 = temporary vector.
// free the vectors vecV1 and vecV2 if they are temporary
SetTemporaryVect(vecIn,vecIn,vecOut);

for nn:=1 to vecIn.nIndMax do vecOut.pavV^[nn]:=vScal * vecIn.pavV^[nn];

vecOut.nIndMax:=vecIn.nIndMax; // Length of the resulting vector
end; // operator *

operator = (vecV1, vecV2:RVectRNFloat) fRes:boolean;
//==================================================
// True if Norme(vecV1 - vecV2) <= gglvVNEpsilon
begin
fRes:=Norme(vecV1 - vecV2) <= gglvVNEpsilon;

// If  vecV1 was a temporary result. It is free to use now
// If  vecV2 was a temporary result. It is free to use now
if (vecV1.fIsTemp) then vecV1.pOrg^.fFreeToUse:=TRUE;
if (vecV2.fIsTemp) then vecV2.pOrg^.fFreeToUse:=TRUE;
end; // operator =

end. // unit xVectRN;
