# NEW FEATURES IN MATHPARTNER 2021

###### Abstract

We introduce new features in the MathPartner service that have recently become
available to users. We highlight the functions for calculating both
arithmetic-geometric mean and geometric-harmonic mean. They allow calculating
complete elliptic integrals of the first kind. They are useful for solving
many physics problems, for example, one can calculate the period of a simple
pendulum. Next, one can calculate the modified arithmetic-geometric mean
proposed by Semjon Adlaj. Consequently, one can calculate the complete
elliptic integrals of the second kind as well as the circumference of an
ellipse. Furthermore, one can also calculate the Sylvester matrices of the
first and the second kind. Thus, by means of a few strings, one can calculate
the resultant of two polynomials as well as the discriminant of a binary form.
Some new matrix functions are also added. So, today the list of matrix
functions includes the transpose, adjugate, conjugate, inverse, generalized
inverse, and pseudo inverse of a matrix, the matrix determinant, the kernel,
the echelon form, the characteristic polynomial, the Bruhat decomposition, the
triangular LSU decomposition, which is an exact block recursive LU
decomposition, the QR block recursive decomposition, and the singular value
decomposition. In addition, two block-recursive functions have been
implemented for calculating the Cholesky decomposition of symmetric positive-
definite matrices: one function for sparse matrices with the standard
multiplication algorithm and another function for dense matrices with
multiplication according to the Winograd–Strassen algorithm. The linear
programming problems can be solved too. So, the MathPartner service has become
better and handy. It is freely available at http://mathpar.ukma.edu.ua/ as
well as at http://mathpar.com/ .

Keywords: computer algebra, arithmetic-geometric mean, geometric-harmonic
mean, complete elliptic integral, pendulum, Sylvester matrix, Bruhat
decomposition, LSU decomposition, QR decomposition, Cholesky decomposition,
modern teaching technologies

Citation: Malaschonok G.I., Seliverstov A.V.. . Computer assisted mathematics,
2021. № 3. P. LABEL:first@page–NEW FEATURES IN MATHPARTNER 2021. .

##  1 INTRODUCTION

The MathPartner service is useful at school, university, and work [15, 16]. It
can help you to solve problems in mathematical analysis, algebra, geometry,
physics, and more. You can operate with functions and functional matrices, to
obtain the exact numerical and analytical solutions and solutions in which the
numerical coefficients have a required accuracy. Today it is available at
http://mathpar.ukma.edu.ua/ as well as at http://mathpar.com/.

We present some new features and improvements. In particular, you can
calculate the arithmetic-geometric mean and its modification to calculate the
complete elliptic integrals of the first and the second kind [1, 2]. Thus, you
can calculate the circumference of an ellipse as well as the period of a
pendulum. Another application is also proposed to compute packing properties
[13, 20, 23]. One can also calculate the Sylvester matrix [3, 4] as well as
the resultant of two polynomials. Some new matrix functions have been
implemented too [8, 14, 17, 18].

##  2 SIX MEANS AND THE COMPLETE ELLIPTIC INTEGRALS

Given two non-negative numbers $x$ and $y$, one can define their arithmetic,
geometric and harmonic means as $\frac{x+y}{2}$, $\sqrt{xy}$, and
$\frac{2xy}{x+y}$, respectively. Moreover, $\mathbf{AGM}{(x,y)}$ denotes the
arithmetic-geometric mean of $x$ and $y$. It was defined by Johann Carl
Friedrich Gauss at the end of the 18th century. $\mathbf{GHM}{(x,y)}$ denotes
the geometric-harmonic mean of $x$ and $y$. At last, $\mathbf{MAGM}{(x,y)}$
denotes the modified arithmetic-geometric mean of $x$ and $y$. It is defined
by Semjon Adlaj [1, 2]. Every mean is a symmetric homogeneous function in
their two variables $x$ and $y$. In contrast to well-known means,
$\mathbf{AGM}{(x,y)}$, $\mathbf{GHM}{(x,y)}$, and $\mathbf{MAGM}{(x,y)}$ are
calculated iteratively.

The arithmetic-geometric mean $\mathbf{AGM}{(x,y)}$ is equal to the limit of
both sequences $x_{n}$ and $y_{n}$, where $x_{0}=x$, $y_{0}=y$,
$x_{n+1}=\frac{1}{2}(x_{n}+y_{n})$, and $y_{n+1}=\sqrt{x_{n}y_{n}}$.

In the same way, the geometric-harmonic mean $\mathbf{GHM}{(x,y)}$ is equal to
the limit of both sequences $x_{n}$ and $y_{n}$, where $x_{0}=x$, $y_{0}=y$,
$x_{n+1}=\sqrt{x_{n}y_{n}}$, and $y_{n+1}=\frac{2x_{n}y_{n}}{x_{n}+y_{n}}$.
Note that $\mathbf{AGM}{(x,y)}\mathbf{GHM}{(x,y)}=xy$.

The modified arithmetic-geometric mean $\mathbf{MAGM}{(x,y)}$ is equal to the
limit of the sequence $x_{n}$, where $x_{0}=x$, $y_{0}=y$, $z_{0}=0$,
$x_{n+1}=\frac{x_{n}+y_{n}}{2}$,
$y_{n+1}=z_{n}+\sqrt{(x_{n}-z_{n})(y_{n}-z_{n})}$, and
$z_{n+1}=z_{n}-\sqrt{(x_{n}-z_{n})(y_{n}-z_{n})}$.

For example, let us run the commands, where functions begin with the symbol
$\backslash$, SPACE denotes the ring of coefficients, and FLOATPOS denotes the
number of decimal places

SPACE=R64[]; FLOATPOS=3; a=$\backslash$AGM(1,5); g=$\backslash$GHM(1,5);
m=$\backslash$MAGM(1,5); [a,g,m];

The output is equal to $[2.604,1.920,2.611]$.

These means are applicable, in particular, to calculate the complete elliptic
integrals of the first and second kind. Let us use the parameter $0\leq k\leq
1$.

The complete elliptic integral of the first kind $K(k)$ is defined as

$K(k)=\int_{0}^{1}\frac{dt}{\sqrt{(1-t^{2})(1-k^{2}t^{2})}}$

It can be computed in terms of the arithmetic-geometric mean:

$K(k)=\frac{\pi}{2\mathbf{AGM}(1,\sqrt{1-k^{2}})}$

On the other hand, for $k<1$, it can be computed in terms of the geometric-
harmonic mean:

$K(k)=\frac{\pi}{2}\mathbf{GHM}\left(1,\frac{1}{\sqrt{1-k^{2}}}\right)$

The complete elliptic integral of the second kind $E(k)$ is defined as

$E(k)=\int_{0}^{1}\sqrt{\frac{1-k^{2}t^{2}}{1-t^{2}}}dt$

It can be computed in terms of the modified arithmetic-geometric mean:

$E(k)=K(k)\mathbf{MAGM}(1,1-k^{2})$

The circumference of an ellipse is equal to

$2\pi\frac{\mathbf{MAGM}(a^{2},b^{2})}{\mathbf{AGM}(a,b)},$

where the semi-major and semi-minor axes are denoted $a$ and $b$.

On the other hand, $\pi$ can be expressed as

$\pi=\frac{\left(\mathbf{AGM}(1,\sqrt{2})\right)^{2}}{\mathbf{MAGM}(1,2)-1}$

So, to calculate $\pi$ one can run the commands

SPACE = R[]; FLOATPOS = 24; w = $\backslash$sqrt(2);
($\backslash$AGM(1,w))^2/($\backslash$MAGM(1,2)-1);

Every argument of functions $\mathbf{AGM}()$, $\mathbf{MAGM}()$, and
$\mathbf{GHM}()$ must be either a number or a variable, i.e., compound
expressions cannot be used as arguments.

Let a point mass be suspended from a pivot with a massless cord. The length of
the pendulum is denoted by $L$. It swings under gravitational acceleration
$g=9.80665m/s^{2}$. The maximum angle that the pendulum swings away from the
vertical, called the amplitude, is denoted by $\theta_{0}$. One can find the
period $T$ of the pendulum using the arithmetic-geometric mean

$T=\frac{2\pi}{\mathbf{AGM}(1,\cos(\theta_{0}/2))}\sqrt{\frac{L}{g}}$

If $L=1m$ and $\theta_{0}=120^{\circ}$, then $T=2.7546s$. To calculate the
period one can run the commands

SPACE = R64[]; FLOATPOS = 4; L = 1; g = 9.80665;

T =
$\backslash$value(2*$\backslash$pi*$\backslash$sqrt{L/g}/($\backslash$AGM(1,
0.5));

On the other hand, if $\theta_{0}$ is small, then the period is equal to
$2.0064s$.

Note that we use the $\mathbf{value}()$ command to evaluate $\pi$.

##  3 THE SYLVESTER MATRICES, THE RESULTANT, AND THE DISCRIMINANT

Let us consider two univariate polynomials $f(x)$ and $g(x)$, where
$\deg(f)=n$, $\deg(g)=m$, and $m\leq n$ hold. James Joseph Sylvester
introduced two matrices associated to $f(x)$ and $g(x)$. Please, refer to [3,
4]. More precisely, there are two different Sylvester matrices associated with
two univariate polynomials. Let us denote
$f(x)=f_{n}x^{n}+\cdots+f_{1}x+f_{0}$ and $g(x)=g_{m}x^{m}+\cdots+g+1x+g_{0}$.
The Sylvester matrix of the first kind was introduced in 1840 [21]. It is the
$(n+m)\times(n+m)$ matrix. Its determinant is called the resultant of $f$ and
$g$. For example, if $f=x^{3}+px+q$ and $g=3x^{2}+p$, then the Sylvester
matrix of the first kind is equal to

$\left(\begin{array}[]{ccccc}1&0&p&q&0\\\ 0&1&0&p&q\\\ 3&0&p&0&0\\\
0&3&0&p&0\\\ 0&0&3&0&p\\\ \end{array}\right)$

and its determinant equals $4p^{3}+27q^{2}$, i.e., it is the opposite of the
discriminant of $f$.

The Sylvester matrix of the second kind was introduced in 1853 as an
improvement of the Sturm theory [22]. It is the $(2n)\times(2n)$ matrix, where
$n\geq m$. The first and the second rows are

$\left(\begin{array}[]{ccccccccc}f_{n}&\cdots&f_{m+1}&f_{m}&\cdots&f_{0}&0&\cdots&0\\\
0&\cdots&0&g_{m}&\cdots&g_{0}&0&\cdots&0\\\ \end{array}\right)$

The next pair is the first pair, shifted one column to the right; the first
elements in the two rows are zero. The remaining rows are obtained the same
way as above. For example, if $f=x^{3}+px+q$ and $g=3x^{2}+p$, then the
Sylvester matrix of the second kind is equal to

$\left(\begin{array}[]{cccccc}1&0&p&q&0&0\\\ 0&3&0&p&0&0\\\ 0&1&0&p&q&0\\\
0&0&3&0&p&0\\\ 0&0&1&0&p&q\\\ 0&0&0&3&0&p\\\ \end{array}\right)$

Of course, if the resultant vanishes, then the determinant of the Sylvester
matrix of the second kind vanishes too.

The Sylvester matrix of the first kind can be calculated in MathPartner by a
ternary function called $\mathbf{sylvester}(\cdot,\cdot,0)$, where the third
argument is equal to zero. In the same way, the Sylvester matrix of the second
kind can be calculated in MathPartner by a ternary function called
$\mathbf{sylvester}(\cdot,\cdot,1)$, where the third argument is not equal to
zero. The first and the second arguments are univariate polynomials, for
example, $f(x)$ and $g(x)$. The variable must be the last one in the list of
variables. For example, if the polynomials over the ring of integers depend on
parameters $p$ and $q$, then the declaration in MathPartner can be SPACE =
Z[p,q,x].

The resultant of two univariate polynomials can be calculated as
$\mathbf{resultant}(f,g)$. The variable must be the last one in the list of
variables. For example, let us run

SPACE = Z[a, b, c, x]; f = a*x^2+b*x+c; g = 2*a*x+b; $\backslash$resultant(f,
g);

The output is equal to $4ca^{2}-b^{2}a$.

The discriminant of a univariate polynomial $f(x)=f_{d}x^{d}+\cdots+f_{0}$ is
equal to

$\mathbf{discriminant}(f)=\frac{(-1)^{d(d-1)/2}}{f_{d}}\mathbf{resultant}(f,f^{\prime})$

The discriminant can be calculated immediately. For example,

SPACE = Z[a, b, c, x]; f = a*x^2+b*x+c; $\backslash$discriminant(f);

The output is equal to $-4ca+b^{2}$. There exists another way to calculate the
discriminant of the univariate polynomial $x^{2}+bx+c$, where $b$ and $c$ are
parameters.

SPACE = Z[b, c, x]; f = x^2+b*x+c; -$\backslash$det($\backslash$sylvester(f,
$\backslash$D(f, x), 0));

The output is equal to $-4c+b^{2}$. Of course, $\mathbf{D}(f,x)$ calculates
the first derivative of $f$.

##  4 SYSTEMS OF ALGEBRAIC EQUATIONS

Let us show an application of the resultant of two univariate polynomials. For
this purpose, we consider a system of two polynomial equations in two
variables and eliminate a variable. Of course, variable elimination can be
done by computing a Gröbner basis. So, there exists another way to solve a
system of algebraic equations. Unfortunately, the Gröbner basis approach is
sometimes very complicated. Contrariwise, the approach based on the resultant
is often more effective. Let us consider the system

$\left\\{\begin{array}[]{rcl}x^{2}+y^{2}&=&1\\\ 2x^{2}+xy+y^{2}&=&1\\\
\end{array}\right.$

In this case, solutions to the system correspond to intersection points of the
circle and the ellipse. Let us consider two univariate polynomials depending
on one parameter $x$

$\begin{array}[]{rcl}f(y)&=&x^{2}+y^{2}-1\\\ g(y)&=&2x^{2}+xy+y^{2}-1\\\
\end{array}$

Its resultant is equal to $2y^{4}-3y^{2}+1$. On the other hand, the Gröbner
basis for the reverse lexicographical ordering consists of two polynomials
$x-2y^{3}+2y$ and $2y^{4}-3y^{2}+1$. The second polynomial is equal to the
resultant. Thus, every solution to the system satisfies the equation
$2y^{4}-3y^{2}+1=0$. So, one can eliminate this variable. There exist four
solutions to the equation $y_{1}=-1$, $y_{2}=-\sqrt{2}/2$, $y_{3}=\sqrt{2}/2$,
and $y_{4}=1$. The corresponding values of $x$ are $x_{1}=0$,
$x_{1}=\sqrt{2}/2$, $x_{1}=-\sqrt{2}/2$, and $x_{4}=0$. So, there are four
intersection points.

Next, let us show the corresponding program in MathPartner. The Gröbner basis
of a polynomial ideal can be obtained due to Bruno Buchberger. The algorithm
is implemented as $\mathbf{groebnerB}()$. The same basis can be calculated
using a matrix algorithm that is similar to the F4 algorithm. It is
implemented as $\mathbf{groebner}()$. The ordering is reverse lexicographical.
Note that functions should begin with the symbol $\backslash$.

SPACE = Z[y, x]; f = x^2+y^2 -1; g = 2*x^2+x*y+y^2 -1; $\backslash$groebner(f,
g);

The output consists of two polynomials $[x-2y^{3}+2y,2y^{4}-3y^{2}+1]$.
Another way to calculate the resultant is to run the commands

SPACE = Z[y, x]; f = x^2+y^2 -1; g = 2*x^2+x*y+y^2 -1;
$\backslash$det($\backslash$sylvester(f, g, 0));

The output consists of one univariate polynomial $2y^{4}-3y^{2}+1$. Of course,
one can run

SPACE = Z[y, x]; f = x^2+y^2 -1; g = 2*x^2+x*y+y^2 -1;
$\backslash$resultant(f, g);

To calculate roots of a polynomial one can run $\mathbf{solve}()$.

SPACE = Q[y]; $\backslash$solve(2*y^4 -3*y^2+1 = 0);

The output consists of four numbers expressed in radicals
$[-1,\sqrt{2}/2,((-1)\cdot\sqrt{2}/2),1]$.

Let us run the same command over the field of real numbers. We recommend using
option SPACE = R64[y]. It denotes the set of $64$-bit floating-point numbers
with $52$-digit mantissa, $11$-bit exponent, and one sign bit.

SPACE = R64[y]; $\backslash$solve(2*y^4 -3*y^2+1 = 0);

The output consists of four floating-point numbers $[1.00,-1.00,0.71,-0.71]$.

Of course, systems of linear algebraic equations can be solved with
$\mathbf{solve}()$. For example,

SPACE = Q[]; M = [[1, 2], [3, 1]]; b = [5, 5]; $\backslash$solve(M, b);

The output is equal to $[1,2]^{T}$. There exists another way

SPACE = Q[x, y]; $\backslash$solve([x+2*y = 5, 3*x+y = 5]);

The output is equal to $[1,2]$.

Moreover, one can solve a system of inequalities in one variable. For example,

SPACE = Q[x]; $\backslash$solve([x^2+4*x -5 > 0, x^2 -2*x -8 < 0]);

The output is equal to $(1,4)$. In the next example

SPACE = Q[x]; $\backslash$solve([x < 0, x > 2]);

The output is equal to the empty set $\emptyset$.

##  5 THE GREATEST COMMON DIVISOR OF TWO POLYNOMIALS

In this section we shall consider polynomials over either the field of
rational numbers or the ring of integers. The problem of calculating the
greatest common divisor of two polynomials is important for symbolic
computations, in particular, over a finite extension of the field of rational
numbers [5, 7, 10]. Unfortunately, the bit complexity of the Euclidean
algorithm is exponential. There exists a polynomial upper bound on the number
of arithmetic operations. But the size of a product of integers at
intermediate steps can be very large. For some discussion about the
computational complexity of powers of integers refer to [11].

A modified algorithm based on subresultant residues had been proposed by J.J.
Sylvester [22] and later improved by Walter Habicht [9] and Alkiviadis Akritas
[3].

The main result was obtained by Brown. He found a way to compute the
subresultant PRS without using matrix reduction. He proposed to modify the
Euclidean algorithm, reducing all coefficients by common factors so that they
coincide with the subresultant PRS [6]. This algorithm is applied in
MathPartner to compute the GCD of two polynomials. This approach was further
developed in work [4].

To calculate the greatest common divisor one can run $\mathbf{GCD}(f,g)$; for
example,

SPACE = Z[x]; $\backslash$GCD(9*x, 6*x+6);

The output is equal to $3$. To calculate Bézout coefficients one can run
$\mathbf{extendedGCD}(f,g)$. The least common multiple can be calculated too.

SPACE = Z[x]; $\backslash$LCM(9*x, 6*x+6);

The output is equal to $18x^{2}+18x$.

##  6 MATRIX FUNCTIONS

Today the list of matrix functions includes the transpose, adjugate,
conjugate, inverse, generalized inverse, and pseudo inverse of a matrix, the
matrix determinant, the kernel, the matrix echelon form, the characteristic
polynomial, the Bruhat decomposition, the triangular LSU decomposition, which
is an exact block recursive LU decomposition, the QR block recursive
decomposition, and the singular value decomposition. In addition, two block-
recursive functions have been implemented for calculating the Cholesky
decomposition of symmetric positive definite matrices: one function for sparse
matrices with the standard multiplication algorithm and another function for
dense matrices with multiplication according to the Winograd–Strassen
algorithm. The linear programming problems can be solved too.

For a given matrix $A$, the pseudo inverse of $A$ is a matrix $A^{-}$
satisfying both equalities $AA^{-}A=A$ and $A^{-}AA^{-}=A^{-}$. Furthermore,
the generalized inverse Moore–Penrose $A^{+}$ satisfies four equalities
$AA^{+}A=A$, $A^{+}AA^{+}=A^{+}$, $(A^{+}A)^{T}=A^{+}A$, and
$(AA^{+})^{T}=AA^{+}$. If $A$ is a square non-degenerate matrix, then three
inverses of $A$ coincide, i.e., $A^{-1}=A^{-}=A^{+}$. If a $n\times m$ matrix
$A$ can be decomposed as $A=BC$, where $B$ is a $n\times k$ matrix, $C$ is a
$k\times m$ matrix, and
$\mathrm{rank}(A)=\mathrm{rank}(B)=\mathrm{rank}(C)=k$, then
$A^{+}=C^{T}(CC^{T})^{-1}(B^{T}B)^{-1}B^{T}$. This idea was expressed by Vera
Nikolaevna Kublanovskaya [12]. About big matrices refer to [19].

For a given matrix $A$, one can calculate:

  * •

The transpose: $\mathbf{transpose}(A)$ or $A^{T}$;

  * •

The conjugate: $\mathbf{conjugate}(A)$ or $A^{\ast}$;

  * •

The matrix echelon form: $\mathbf{toEchelonForm}(A)$;

  * •

The kernel: $\mathbf{kernel}(A)$;

  * •

The rank: $\mathbf{rank}(A)$;

  * •

The determinant: $\mathbf{det}(A)$;

  * •

The inverse: $\mathbf{inverse}(A)$ or $A^{-1}$;

  * •

The adjugate: $\mathbf{adjoint}(A)$ or $A^{\star}$;

  * •

The generalized inverse Moore–Penrose: $\mathbf{genInverse}(A)$ or $A^{+}$;

  * •

The pseudo inverse: $\mathbf{pseudoInverse}(A)$;

  * •

The closure: $\mathbf{closure}(A)$ or $A^{\times}$. The closure of a matrix
$A$ is equal to the sum of matrices $I+A+A^{2}+A^{3}+\cdots$. For the
classical algebras it is equivalent to $(I-A)^{-1}$.

To calculate the characteristic polynomial of a matrix $A$, you should work
over the ring of polynomials in some new variable and run
$\mathbf{charPolynom}(A)$. For example, let us run the commands

SPACE=Z[x]; M=[[1, 2], [3, 5]]; f=$\backslash$charPolynom(M);

The output is equal to $f=x^{2}-6x-1$.

Let us take a closer look at some types of decomposition.

###  6.1 The Bruhat decomposition

To calculate the Bruhat decomposition of a matrix $A$ one can run
$\mathbf{BruhatDecomposition}(A)$. The result consists of three matrices
$[V,D,U]$, where both $V$ and $U$ are upper-triangular matrices, $D$ is a
permutation matrix multiplied by the inverse of the diagonal matrix [14]. If
all entries of the matrix $A$ are elements of commutative domain $R$, then all
entries of matrices $V$, $D^{-1}$, and $U$ belong to the same domain $R$. Let
us consider a $2\times 2$ matrix over $\mathbb{Z}$. For example,

$M=\left(\begin{array}[]{cc}-1&2\\\ 1&1\\\ \end{array}\right)$

Let us run the commands

SPACE = Z[]; M = [[-1, 2], [1, 1]]; $\backslash$BruhatDecomposition(M);

The output consists of three matrices

$\left[\left(\begin{array}[]{cc}3&-1\\\ 0&1\\\
\end{array}\right),\left(\begin{array}[]{cc}0&1/3\\\ 1&0\\\
\end{array}\right),\left(\begin{array}[]{cc}1&1\\\ 0&3\\\
\end{array}\right)\right]$

An entry of the middle matrix $D$ is not integer, but the inverse matrix has
integer entries.

$D^{-1}=\left(\begin{array}[]{cc}0&1\\\ 3&0\\\ \end{array}\right)$

###  6.2 The LSU decomposition

The LSU decomposition of a matrix $A$ can be calculated by means of the
command $\mathbf{LSU}(A)$. The result consists of three matrices $[L,S,U]$,
where $L$ is a lower-triangular matrix, $U$ is an upper-triangular matrix, $S$
is a permutation matrix multiplied by the inverse of a diagonal matrix. If all
entries of the matrix $A$ belong to a commutative domain $R$, then all entries
of matrices $L$, $S^{-1}$, and $U$ belong to the same domain $R$, refer to
[18]. Let us consider an example, where $M$ is a $2\times 2$ matrix.

$M=\left(\begin{array}[]{cc}1&2\\\ 3&1\\\ \end{array}\right)$

Let us run the commands

SPACE = Z[]; M = [[1, 2], [3, 1]]; $\backslash$LSU(M);

The output consists of three matrices

$\left[\left(\begin{array}[]{cc}1&0\\\ 3&-5\\\
\end{array}\right),\left(\begin{array}[]{cc}1&0\\\ 0&-1/5\\\
\end{array}\right),\left(\begin{array}[]{cc}1&2\\\ 0&-5\\\
\end{array}\right)\right]$

Both first and third matrices are triangular matrices over $\mathbb{Z}$. The
middle matrix $S$ has a rational entry, but the inverse matrix is defined over
$\mathbb{Z}$.

$S^{-1}=\left(\begin{array}[]{cc}1&0\\\ 0&-5\\\ \end{array}\right)$

To calculate the LSU decomposition of $A$ together with decomposition of the
pseudo inverse $A^{-}=WSM$, one can run the command $\mathbf{LSUWMdet}(A)$.
The result consists of five matrices and determinant of the largest non-
degenerate corner block $[L,S,U,W,M,\det]$, where $L$ and $U$ are lower and
upper triangular matrices, $S$ is a truncated weighted permutation matrix,
$SM$ and $WS$ are lower and upper triangular matrices. Moreover, $A=LSU$ and
$A^{-}=WSM$. If entries of the matrix $A$ belong to a commutative domain, then
all matrices, except for $S$, also belong to this domain. Let us run the
commands

SPACE = Z[]; M = [[1, 2], [3, 1]]; $\backslash$LSUWMdet(M);

The output consists of

$\left[\left(\begin{array}[]{cc}1&0\\\ 3&-5\\\
\end{array}\right),\left(\begin{array}[]{cc}1&0\\\ 0&-1/5\\\
\end{array}\right),\left(\begin{array}[]{cc}1&2\\\ 0&-5\\\
\end{array}\right),\left(\begin{array}[]{cc}-5&10\\\ 0&-5\\\
\end{array}\right),\left(\begin{array}[]{cc}-5&0\\\ 15&-5\\\
\end{array}\right),[[-5]]\right]$

Of course, three of these matrices coincide with three matrices in the
previous example.

Next, let us consider a matrix over the ring $\mathbb{Z}[x,y]$

$M=\left(\begin{array}[]{cc}y&x\\\ x&y\\\ \end{array}\right)$

and run the commands

SPACE = Z[x, y]; M = [[y, x], [x, y]]; $\backslash$LSU(M);

The output consists of three matrices

$\left[\left(\begin{array}[]{cc}y&0\\\ x&y^{2}-x^{2}\\\
\end{array}\right),\left(\begin{array}[]{cc}1/y&0\\\
0&1/\left(y^{3}-yx^{2}\right)\\\
\end{array}\right),\left(\begin{array}[]{cc}y&x\\\ 0&y^{2}-x^{2}\\\
\end{array}\right)\right]$

Entries of the middle matrix $S$ are rational functions. All entries of the
matrices $L$, $S^{-1}$, and $U$ are polynomials over $\mathbb{Z}$.

###  6.3 The QR block recursive decomposition

Let us consider a $2^{k}\times 2^{k}$ matrix $A$ over the field of reals. The
QR decomposition of $A$ can be calculated by means of the command
$\mathbf{QR}(A)$. Note that if the order is not equal to $2^{k}$ for any
integer $k$, then the algorithm does not work because it is based on block
recursion [17]. Let us consider an example, where $M$ is a $2\times 2$ matrix.

$M=\left(\begin{array}[]{cc}1&2\\\ 3&1\\\ \end{array}\right)$

Let us run the commands

SPACE = R64[]; M = [[1, 2], [3, 1]]; $\backslash$QR(M);

The output consists of two matrices

$\left[\left(\begin{array}[]{cc}0.32&-0.95\\\ 0.95&0.32\\\
\end{array}\right),\left(\begin{array}[]{cc}3.16&1.58\\\ 0&-1.58\\\
\end{array}\right)\right]$

The first matrix is orthogonal. The second matrix is upper-triangular. Their
product is equal to the initial matrix $M$.

###  6.4 The singular value decomposition

To calculate the singular value decomposition (SVD) of a matrix $A$, one can
run $\mathbf{SVD}(A)$. As a result, three matrices $[U,D,V]$ will be
calculated. The matrices $U$ and $V$ are unitary, the matrix $D$ is diagonal,
and $A=UDV$ holds. Let us consider an example, where $M$ is a $2\times 2$
matrix.

$M=\left(\begin{array}[]{cc}2&3\\\ 1&0\\\ \end{array}\right)$

Let us run the commands

SPACE = R64[]; FLOATPOS = 3; M = [[2, 3], [1, 0]]; $\backslash$SVD(M);

The output is equal to

$\left[\left(\begin{array}[]{cc}-0.987&-0.16\\\ -0.16&0.987\\\
\end{array}\right),\left(\begin{array}[]{cc}3.65&0\\\ 0&0.822\\\
\end{array}\right),\left(\begin{array}[]{cc}-0.585&-0.811\\\ 0.811&-0.585\\\
\end{array}\right)\right]$

###  6.5 The Cholesky decomposition

In general, the Cholesky decomposition is a decomposition of a Hermitian
positive-definite matrix into the product of a lower-triangular matrix and its
conjugate transpose, which is useful for efficient numerical solutions. It was
discovered by André-Louis Cholesky for real symmetric matrices [8]. And we
also suppose that matrices are real. So, every real symmetric positive-
definite matrix is equal to the product $LL^{T}$, where $L$ is a lower-
triangular matrix.

The Cholesky decomposition can be calculated for a symmetric and positive
definite matrix $A$ by means of the command $\mathbf{cholesky}(A)$. The result
consists of two lower triangular matrices $L$ and $S$ such that $A=LL^{T}$ and
$SL=I$. Let us consider an example, where $M$ is a $2\times 2$ matrix.

$M=\left(\begin{array}[]{cc}3&2\\\ 2&4\\\ \end{array}\right)$

Let us run the commands

SPACE = R64[]; FLOATPOS = 2; M = [[3, 2], [2, 4]]; $\backslash$cholesky(M);

The output is equal to

$\left[\left(\begin{array}[]{cc}1.73&0\\\ 1.15&1.63\\\
\end{array}\right),\left(\begin{array}[]{cc}0.58&0\\\ -0.41&0.61\\\
\end{array}\right)\right]$

For large dense matrices, whose size is greater than or equal to $128\times
128$, one can use a fast algorithm $\mathbf{cholesky}(A,1)$ that uses
multiplication of blocks by the Winograd–Strassen algorithm.

##  7 MODULAR ARITHMETIC

The current version of the MathPartner service supports operations over a
finite field $\mathbb{Z}/p\mathbb{Z}$, where $p$ is a prime number. One should
use either SPACE = Zp[] or SPACE = Zp32[]. The prime number $p$ is equal to
the constant MOD or MOD32, respectively. In the second case, $p$ satisfies the
inequality $p<2^{31}$. The default value is $268435399$. For example, working
over the field $\mathbb{Z}/5\mathbb{Z}$ one can run

SPACE = Zp32[x]; MOD32 = 5; $\backslash$GCD(x+2,x-3);

The output is equal to $x-3$ because $-3\equiv 2\pmod{5}$. On the other hand,
the same example over $\mathbb{Z}/7\mathbb{Z}$ leads to another answer.

SPACE = Zp32[x]; MOD32 = 7; $\backslash$GCD(x+2,x-3);

The output is equal to $1$.

All functions using only rational operations on input data can be calculated
over finite fields. In particular, for two polynomials over
$\mathbb{Z}/p\mathbb{Z}$ one can calculate the greatest common divisor
$\mathbf{GCD}()$ as well as Sylvester matrix $\mathbf{sylvester}()$. One can
calculate the Gröbner basis of an ideal in a polynomial ring using either
$\mathbf{groebner}()$ or $\mathbf{groebnerB}()$. One can also calculate the
determinant $\mathbf{det}()$, echelon form $\mathbf{toEchelonForm}()$,
characteristic polynomial $\mathbf{charPolynom}()$, Bruhat decomposition
$\mathbf{BruhatDecomposition}()$, LDU decomposition $\mathbf{LDU}()$, and
$\mathbf{LDUWMdet}()$ of a matrix.

##  8 CONCLUSION

Now the MathPartner service has become even better and allows us to solve new
problems in geometry and physics. In particular, new functions allow to
calculate the period of a simple pendulum as well as the circumference of an
ellipse in terms of the arithmetic-geometric and the modified arithmetic-
geometric means. The resultant of two univariate polynomials is a basic tool
of computer algebra because it allows solving systems of polynomial equations.
Matrix functions are also widely used to solve applied problems.

The reader is recommended to calculate examples of the considered quantities
using the MathPartner service. These exercises will help you remember and
understand the computer algebra methods better. On the other hand, new
algorithms can be implemented by the user through the branch and loop
operators. Moreover, the MathPartner service opens up the possibility of
distance learning.

## References

  * [1] Adlaj S. ‘‘An eloquent formula for the perimeter of an ellipse,’’ Notices of the American Mathematical Society, 2012, vol. 59, no. 8, pp. 1094–1099. DOI:10.1090/noti879
  * [2] Adlaj S. ‘‘An Arithmetic-Geometric Mean of a Third Kind!’’ In: England M., Koepf W., Sadykov T., Seiler W.M., Vorozhtsov E.V. (eds) Computer Algebra in Scientific Computing. 21st International Workshop, CASC 2019, Moscow, Russia, August 26–30, 2019, Proceedings, Lecture Notes in Computer Science, 2019, vol. 11661, pp. 37–56.
  * [3] Akritas A.G. Elements of Computer Algebra with Applications. John Wiley and Sons, NY, 1989.
  * [4] Akritas A.G., Malaschonok G.I., Vigklas P.S. ‘‘Subresultant polynomial remainder sequences obtained by polynomial divisions in Q[x] or in Z[x],’’ Serdica Journal of Computing, 2016, vol. 10, no. 3–4, pp. 197–217.
  * [5] Alaev P.E., Selivanov V.L. ‘‘Fields of algebraic numbers computable in polynomial time. I.’’ Algebra and Logic, 2020, vol. 58, no. 6, pp. 447–469. DOI: 10.1007/s10469-020-09565-0
  * [6] Brown W.S. ‘‘The Subresultant PRS Algorithm,’’ ACM Transactions on Mathematical Software, 1978, vol.4, no.3, pp. 237–-249  DOI: 10.1145/355791.355795
  * [7] Dolgov D.A. ‘‘Polynomial greatest common divisor as a solution of system of linear equations,’’ Lobachevskii Journal of Mathematics, 2018, vol. 39, no. 7, pp. 985–991. DOI: 10.1134/S1995080218070090
  * [8] Grcar J.F. ‘‘How ordinary elimination became Gaussian elimination,’’ Historia Mathematica, 2011, vol. 38, pp. 163–218. DOI:10.1016/j.hm.2010.06.003
  * [9] Habicht W. ‘‘Eine Verallgemeinerung des Sturmschen Wurzelzählverfahrens,’’ Commentarii Mathematici Helvetici, 1948, vol. 21, pp. 99–116. DOI:10.1007/BF02568028
  * [10] van der Hoeven J., Lecerf G.‘‘Fast computation of generic bivariate resultants,’’ Journal of Complexity, 2021, vol. 62, article 101499. DOI:10.1016/j.jco.2020.101499
  * [11] Kotochigov A.M., Suchkov A.I. ‘‘A method for reducing iteration in algorithms for building minimal additive chains,’’ Computer Tools in Education, 2020, no. 1, pp. 5–18 (in Russian). DOI: 10.32603/2071-2340-2020-1-5-18
  * [12] Kublanovskaya V.N. ‘‘Evaluation of a generalized inverse matrix and projector,’’ USSR Computational Mathematics and Mathematical Physics, 1966, vol. 6, no. 2, pp. 179–188. DOI: 10.1016/0041-5553(66)90064-4
  * [13] Lamarche F., Leroy C. ‘‘Evaluation of the volume of intersection of a sphere with a cylinder by elliptic integrals,’’ Computer Physics Communications, 1990, vol. 59, no. 2, pp. 359–369.
  * [14] Malaschonok G. ‘‘Generalized Bruhat decomposition in commutative domains.’’ In: Gerdt V.P., Koepf W., Mayr E.M., Vorozhtsov E.V. (eds) Computer Algebra in Scientific Computing. 15th International Workshop, CASC 2013, Berlin, Germany, September 9–13, 2013. Proceedings. Lecture Notes in Computer Science, 2013, vol. 8136, pp. 231–242.
  * [15] Malaschonok G.I. ‘‘Application of the MathPartner service in education,’’ Computer Tools in Education, 2017, no. 3, pp. 29–37 (in Russian).
  * [16] Malaschonok G.I. ‘‘MathPartner computer algebra,’’ Programming and Computer Software, 2017, vol. 43, no. 2, pp. 112–118. DOI:10.1134/S0361768817020086
  * [17] Malaschonok G. ‘‘Recursive matrix algorithms, distributed dynamic control, scaling, stability,’’ In: Proc. of 12th Int. Conf. on Comp. Sci. and Information Technologies (CSIT-2019). September 23–27, 2019, Yerevan, pp. 175–178.
  * [18] Malaschonok G.I. ‘‘LDU-factorization,’’ ArXiv e-print 2011.04108, 16 p. https://arxiv.org/abs/2011.04108
  * [19] Malaschonok G., Tchaikovsky I. ‘‘About big matrix inversion.’’ In: Abramov S.A., Sevastyanov L.A. (eds) Computer algebra: 4th International Conference Materials. Moscow, June 28–29, 2021. MAKS Press, Moscow, 2021, pp. 81–84. DOI: 10.29003/m2019.978-5-317-06623-9
  * [20] Mariani N.J., Mazza G.D., Martinez O.M., Barreto G.F. ‘‘Evaluation of radial voidage profiles in packed beds of low-aspect ratios,’’ The Canadian Journal of Chemical Engineering, 2000, vol. 78, no. 6, pp. 1133–1137. DOI:10.1002/cjce.5450780614
  * [21] Sylvester J.J. ‘‘A method of determining by mere inspection the derivatives from two equations of any degree,’’ Philosophical Magazine, 1840, vol. 16, pp. 132–135.
  * [22] Sylvester J.J. ‘‘On the theory of syzygetic relations of two rational integral functions, comprising an application to the theory of Sturm’s functions, and that of the greatest algebraical common measure,’’ Philosophical Transactions, 1853, vol. 143, pp. 407–548.
  * [23] Xu B.-X., Gao Y., Wang M.-Z. ‘‘Particle packing and the mean theory,’’ Physics Letters A, 2013, vol. 377, no. 3–4, pp. 145–147. DOI:10.1016/j.physleta.2012.11.022