## Introduction

The fundamental algorithms used in matrix and polynominal algebra are to be presented in this article. All algorithms' implementations in polynominal are shown both in the fields of complex numbers. Algorithms' implementation in matrix (and vector) algebra are shown in the field of real numbers.

## Background

The intermediate algebra definitions and operations knowledge will be required to follow implemenatation.

These are:

- complex numbers operations;
- basic opertaions at matrices (addition, subtraction, multiplication, scaling, transpose, trace, determinant);
- matrix characteristic polynominal calculation using the
*Faddeev-Leverrier* algorithm; - inverse of matrix calculation using the
*Faddeev-Leverrier* algorithm; - polynominal derivative calculation algorithm;
- the
*Horner's* method of polynominal by binominal division and polynominal evaluation; - the
*Laguerre’s *method for calculating polynominal roots;

The complex number structure and operations are provided with `System.Numerics`

namespace.

## Using the code

### Convert Double to Complex method

The very simple methods that converts Double values to Complex structure.

public Complex Double2Complex(Double Value)
{
Complex ValueComplex = new Complex(Value, 0);
return ValueComplex;
}
public Complex[] Double2Complex(Double[] Values)
{
Complex[] ValuesComplex = new Complex[Values.GetLength(0)];
for (int m = 0; m < Values.GetLength(0); m++)
{
ValuesComplex[m] = new Complex(Values[m], 0);
}
return ValuesComplex;
}
public Complex[,] Double2Complex(Double[,] Values)
{
Complex[,] ValuesComplex = new Complex[Values.GetLength(0), Values.GetLength(1)];
for (int m = 0; m < Values.GetLength(0); m++)
{
for (int n = 0; n < Values.GetLength(1); n++)
{
ValuesComplex[m, n] = new Complex(Values[m, n], 0);
}
}
return ValuesComplex;
}

### Matrix operations

#### Identity matrix

The method returns square identity matrix of given dimension.

public Double[,] IdentityMatrix(int Dimension)
{
if (Dimension > 0)
{
Double[,] I = new Double[Dimension, Dimension];
for (int m = 0; m < Dimension; m++)
{
for (int n = 0; n < Dimension; n++)
{
if (m == n)
{
I[m, n] = 1;
}
else
{
I[m, n] = 0;
}
}
}
return I;
}
else
{
throw new Exception("Output identity matrix dimension must be positive.");
}
}

Matrix addition

The method returns a result of metrices addition. Both input matrices must have the same dimensions.

public Double[,] MatrixAddition(Double[,] MatrixA, Double[,] MatrixB)
{
int M_A = MatrixA.GetLength(0);
int N_A = MatrixA.GetLength(1);
int M_B = MatrixB.GetLength(0);
int N_B = MatrixB.GetLength(1);
if (M_A == M_B && N_A == N_B)
{
Double[,] Sum = new Double[M_A, N_A];
for (int m = 0; m < M_A; m++)
{
for (int n = 0; n < N_A; n++)
{
Sum[m, n] = MatrixA[m, n] + MatrixB[m, n];
}
}
return Sum;
}
else
{
throw new Exception("Matrices dimensions must be equal.");
}
}

#### Matrix subtraction

The method returns a result of matrices subtraction. Both input matrices must have the same dimensions.

public Double[,] MatrixSubtraction(Double[,] MatrixA, Double[,] MatrixB)
{
int M_A = MatrixA.GetLength(0);
int N_A = MatrixA.GetLength(1);
int M_B = MatrixB.GetLength(0);
int N_B = MatrixB.GetLength(1);
if (M_A == M_B && N_A == N_B)
{
Double[,] Sub = new Double[M_A, N_A];
for (int m = 0; m < M_A; m++)
{
for (int n = 0; n < N_A; n++)
{
Sub[m, n] = MatrixA[m, n] - MatrixB[m, n];
}
}
return Sub;
}
else
{
throw new Exception("Matrices dimensions must be equal.");
}
}

#### Matrix multiplication

The method returns a result of metrices multiplication. The number of columns in first multiplication component must be equal to the number of rows in second one.

public Double[,] MatrixMultiplication(Double[,] MatrixA, Double[,] MatrixB)
{
int M_A = MatrixA.GetLength(0);
int N_A = MatrixA.GetLength(1);
int M_B = MatrixB.GetLength(0);
int N_B = MatrixB.GetLength(1);
if (N_A != M_B)
{
throw new Exception("Number of columns in A matrix must be equal to number of rows in B matrix.");
}
else
{
Double[,] _MatrixAB = new Double[M_A, N_B];
for (int m_a = 0; m_a < M_A; m_a++)
{
for (int n_b = 0; n_b < N_B; n_b++)
{
_MatrixAB[m_a, n_b] = 0;
for (int n_a = 0; n_a < N_A; n_a++)
{
_MatrixAB[m_a, n_b] = _MatrixAB[m_a, n_b] + MatrixA[m_a, n_a] * MatrixB[n_a, n_b];
}
}
}
return _MatrixAB;
}
}

#### Matrix scaling

The method returns a result of matrix muliplication by constans coefficient.

public Double[,] MatrixScaling(Double Coefficient, Double[,] MatrixA)
{
int M_A = MatrixA.GetLength(0);
int N_A = MatrixA.GetLength(1);
for (int m = 0; m < M_A; m++)
{
for (int n = 0; n < N_A; n++)
{
MatrixA[m, n] = Coefficient * MatrixA[m, n];
}
}
return MatrixA;
}

#### Transpose of the matrix

The method returns a result of matrix transpose operation.

public Double[,] MatrixTranspose(Double[,] Matrix)
{
int M = Matrix.GetLength(0);
int N = Matrix.GetLength(1);
Double[,] _MatrixT = new Double[N, M];
for (int m = 0; m < M; m++)
{
for (int n = 0; n < N; n++)
{
_MatrixT[n, m] = Matrix[m, n];
}
}
return _MatrixT;
}

#### Matrix trace

The method returns a matrix trace value - the sum of all diagonal elements. To calculate trace, input matrix must be square.

public Double MatrixTrace(Double[,] Matrix)
{
int M_A = Matrix.GetLength(0);
int N_A = Matrix.GetLength(1);
if (M_A == N_A)
{
Double Trace = 0;
for (int m = 0; m < M_A; m++)
{
Trace += Matrix[m, m];
}
return Trace;
}
else
{
throw new Exception("Matrix must be square.");
}
}

#### Matrix determinant

A recursive method using *Laplace* expansion and algebraic complement calculation.

public Double MatrixDeterminant(Double[,] Matrix)
{
int M = Matrix.GetLength(0);
int N = Matrix.GetLength(1);
if (M == N)
{
Double Det = 0;
if (M == 1)
{
Det = Matrix[0, 0];
}
else if (M == 2)
{
Det = Matrix[0, 0] * Matrix[1, 1] - Matrix[0, 1] * Matrix[1, 0];
}
else
{
Double[,] AlgebraicComplement = new Double[M - 1, M - 1];
for (int m = 0; m < M; m++)
{
int a = 0;
for (int k = 1; k < M; k++)
{
int b = 0;
for (int l = 0; l < M; l++)
{
if (l != m)
{
AlgebraicComplement[a, b] = Matrix[k, l];
b++;
}
}
a++;
}
Det += Math.Pow(-1, m) * Matrix[0, m] * MatrixDeterminant(AlgebraicComplement);
}
}
return Det;
}
else
{
throw new Exception("Matrix must be square.");
}
}

#### Matrix characteristic polynominal

The method calculates matrix characteristic polynominal coefficients using *Faddeev - Leverrier* algorithm. the method returns vector with elements corresponding to next polynominal coefficients starting with free term.

For polynominal with given description:

the vector of coefficients will be:

** **

The same polynominal's coefficients structure will be used in all presented methods.

public Double[] CharacteristicPolynominalCoefficients(Double[,] Matrix)
{
int M = Matrix.GetLength(0);
int N = Matrix.GetLength(1);
if (M == N)
{
Double[] Coeffs = new Double[M + 1];
Double[] CoeffsSorted = new Double[M + 1];
Double[,] B;
Double LastCoeff;
Coeffs[0] = 1;
B = Matrix;
LastCoeff = MatrixTrace(B);
Coeffs[1] = -LastCoeff;
for(int m = 2; m < M + 1; m++)
{
B = MatrixMultiplication(Matrix, MatrixSubtraction(B, MatrixScaling(LastCoeff, IdentityMatrix(B.GetLength(0)))));
LastCoeff = MatrixTrace(B) / m;
Coeffs[m] = -LastCoeff;
}
for (int m = 0; m < M + 1; m++)
{
CoeffsSorted[m] = Coeffs[M - m];
}
return CoeffsSorted;
}
else
{
throw new Exception("Input data matrix must be square.");
}
}

#### The inverse of matrix

The method calculates the inverse of matrix using *Faddeev-Leverrier* algorithm.

public Double[,] MatrixInverse(Double[,] Matrix)
{
int M = Matrix.GetLength(0);
int N = Matrix.GetLength(1);
if (M == N)
{
if (MatrixDeterminant(Matrix) != 0)
{
Double[,] LastB;
Double[,] NextB;
Double[,] Inv = new Double[M, N];
Double LastCoeff;
Double NextCoeff;
LastB = Matrix;
LastCoeff = MatrixTrace(LastB);
for (int m = 2; m < M; m++)
{
LastB = MatrixMultiplication(Matrix, MatrixSubtraction(LastB, MatrixScaling(LastCoeff, IdentityMatrix(LastB.GetLength(0)))));
LastCoeff = MatrixTrace(LastB) / m;
}
NextB = MatrixMultiplication(Matrix, MatrixSubtraction(LastB, MatrixScaling(LastCoeff, IdentityMatrix(LastB.GetLength(0)))));
NextCoeff = MatrixTrace(NextB) / M;
Inv = MatrixScaling(1 / NextCoeff, MatrixSubtraction(LastB, MatrixScaling(LastCoeff, IdentityMatrix(LastB.GetLength(0)))));
return Inv;
}
else
{
throw new Exception("The matrix is not invertible over commutative ring.");
}
}
else
{
throw new Exception("Input data matrix must be square.");
}
}

### Polynominal operations

#### Polynominal evaluation

The method calculates polynomianl evaluation for given value.

public Complex PolynominalEvaluation(Complex[] Coefficients, Complex Value)
{
int M = Coefficients.GetLength(0);
if (M > 1)
{
Complex[] CoeffsSorted = new Complex[M];
Complex Last;
for (int m = 0; m < M; m++)
{
CoeffsSorted[m] = Coefficients[M - m - 1];
}
Last = CoeffsSorted[0];
for (int m = 1; m < M; m++)
{
Last = CoeffsSorted[m] + Last * Value;
}
return Last;
}
else
{
return Coefficients[0];
}
}

#### Polynominal addition

The method calculates the sum of two polynominals (the sum of coefficients).

public Complex[] PolynominalAddition(Complex[] Poly1, Complex[] Poly2)
{
int M = Poly1.GetLength(0);
int N = Poly2.GetLength(0);
int K = Math.Max(M, N);
Complex[] Poly1Ext = new Complex[K];
Complex[] Poly2Ext = new Complex[K];
Complex[] Add = new Complex[K];
for (int m = 0; m < M; m++)
{
Poly1Ext[m] = Poly1[m];
}
for (int n = 0; n < N; n++)
{
Poly2Ext[n] = Poly2[n];
}
for (int k = 0; k < K; k++)
{
Add[k] = Poly1Ext[k] + Poly2Ext[k];
}
return Add;
}

#### Polynominal scaling

The method calculates the polynominal multiplication by value.

public Complex[] PolynominalScaling(Complex[] Polynominal, Complex Value)
{
int M = Polynominal.GetLength(0);
for (int m = 0; m < M; m++)
{
Polynominal[m] = Value * Polynominal[m];
}
return Polynominal;
}

#### Polynominal subtraction

Polynominal subtraction can be calculated using polynominal addition and polynominal scaling methods.

Complex[] Result = PolynominalAddition(Polynominal1, PolynominalScaling(Polynominal2, -1));

Polynominal by binominal division

The method calculates a coefficients of new polynominal which is a result of polynominal and binominal division. The *Horner's* method is used for calculation. This method is not returning the remainder.

public Complex[] PolynominalByBinominalDivision(Complex[] Coefficients, Complex BinominalRoot)
{
int M = Coefficients.GetLength(0);
if (M > 1)
{
int N = M - 1;
Complex[] Quotient = new Complex[N];
Complex[] QuotientSorted = new Complex[N];
Complex[] CoeffsSorted = new Complex[M];
Complex Last;
for (int m = 0; m < M; m++)
{
CoeffsSorted[m] = Coefficients[M - m - 1];
}
Last = CoeffsSorted[0];
Quotient[0] = Last;
for (int m = 1; m < N; m++)
{
Last = CoeffsSorted[m] + BinominalRoot * Last;
Quotient[m] = Last;
}
Complex Remainder = CoeffsSorted[M - 1] + BinominalRoot * Last;
for (int n = 0; n < N; n++)
{
QuotientSorted[n] = Quotient[N - n - 1];
}
return QuotientSorted;
}
else
{
throw new Exception("Given coefficients number is not corresponding to correct polynominal structure.");
}
}

####

#### Polynominal derivative

The method calculates a derivative of given polynominal.

public Complex[] PolynominalDerivative(Complex[] Coefficients)
{
int Degree = Coefficients.GetLength(0) - 1;
int NumCoefs = Coefficients.GetLength(0);
Complex[] Derivative = new Complex[NumCoefs - 1];
if (Degree > 0)
{
for (int i = 1; i < NumCoefs; i++)
{
Derivative[i - 1] = i * Coefficients[i];
}
}
else
{
return new Complex[1] { 0 };
}
return Derivative;
}

#### Polynominal roots

The method calculates polynominal roots with given accuracy and seed root (if unknown any real number can be given) using *Laguerre’s* algorithm. The method calculation assume that every polynominal root is a complex numer.

public Complex[] PolynominalRoots(Complex[] Coefficients, Complex Seed, Double Accuracy)
{
int N = Coefficients.GetLength(0) - 1;
Complex[] Roots = new Complex[N];
Complex[] OrginalCoeffs = Coefficients;
int degree = N;
for (int n = 0; n < N; n++)
{
while (true)
{
Complex _tmp0 = PolynominalEvaluation(PolynominalDerivative(Coefficients), Seed);
Complex _tmp1 = (degree - 1) * Complex.Pow(PolynominalEvaluation(PolynominalDerivative(Coefficients), Seed), 2);
Complex _tmp2 = degree * PolynominalEvaluation(Coefficients, Seed) * PolynominalEvaluation(PolynominalDerivative(PolynominalDerivative(Coefficients)), Seed);
Complex _tmp3 = Complex.Pow((degree - 1) * (_tmp1 - _tmp2), 0.5);
Complex _tmp4 = MaximizeResultAbsolute(_tmp0, _tmp3);
Seed = Seed - (degree * PolynominalEvaluation(Coefficients, Seed)) / _tmp4;
Complex result = PolynominalEvaluation(Coefficients, Seed);
if (Complex.Abs(result) < Math.Abs(Accuracy))
{
Roots[n] = Seed;
Coefficients = PolynominalByBinominalDivision(Coefficients, Seed);
degree--;
Random _Random = new Random();
Seed = -10 + (_Random.NextDouble() * 15);
break;
}
}
}
return Roots;
}

#### Maximize absolute value of operation result

The method calculates the result of operation (addition or subtraction) which gives greater absolute value. Used in *Laguerre’s *method.

public Complex MaximizeResultAbsolute(Complex Value1, Complex Value2)
{
if (Complex.Abs(Value1 + Value2) > Complex.Abs(Value1 - Value2))
{
return Value1 + Value2;
}
else
{
return Value1 - Value2;
}
}