{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Lecture 30 of MCS 320, we do linear algebra."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Systems of Linear Equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us define a random 3-by-3 system over the rational numbers, with numerators and denominators bounded by 100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[  -5/4 -57/79  13/25]\n",
       "[-11/84 -37/31 -34/93]\n",
       "[ -11/9     46  -1/18]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = random_matrix(QQ, 3, num_bound=100, den_bound=100)\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The matrix A is the coefficient matrix of our random linear system.  The right hand side vector is random vector with the same bounds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-22/43, -87/41, -62/27)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = random_vector(QQ, 3, num_bound=100, den_bound=100)\n",
    "b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The backslash operator (not to confuse with the division operator) is a convenient operator to solve the linear system defined by the coefficient matrix ``A`` and the right hand side vector ``b``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5005572229574/2065559623287, 252208486789/12393357739722, 3352882185150/688519874429)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = A\\b\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The residual vector is the difference between ``b`` and the product of ``A`` with the solution ``x``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0, 0, 0)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = b - A*x\n",
    "r"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We solved a random system over the rational numbers.  Observe the expression swell in the numbers of the solution ``x``.  We started out with numbers not longer than two decimal places and ended with numbers of more than six decimal places.  And that only for 3-dimensional problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can define matrices via a formula.  For example, the $(i,j)$-th element of the Hilbert matrix is $1/(i+j)$ for $i$ and $j$ starting at 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "h = lambda i,j: 1/((i+1)+(j+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "H = Matrix(QQ, 4, 4, h)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1/2 1/3 1/4 1/5]\n",
       "[1/3 1/4 1/5 1/6]\n",
       "[1/4 1/5 1/6 1/7]\n",
       "[1/5 1/6 1/7 1/8]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Matrix Factorizations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The four important matrix factorizations are\n",
    "\n",
    "1. The LU factorization.\n",
    "\n",
    "2. The QR factorization.\n",
    "\n",
    "3. The Eigenvalue Decomposition.\n",
    "\n",
    "4. The Singular Value Decomposition."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We work on the random 3-dimensional matrix ``A`` from above to illustrate those four factorizations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 the LU factorization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "P, L, U = A.LU()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1 0 0]\n",
       "[0 0 1]\n",
       "[0 1 0]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[              1               0               0]\n",
       "[          44/45               1               0]\n",
       "[         11/105 -143739/6005041               1]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[                -5/4               -57/79                13/25]\n",
       "[                   0           55346/1185             -141/250]\n",
       "[                   0                    0 -390538783/900756150]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[  -5/4 -57/79  13/25]\n",
       "[ -11/9     46  -1/18]\n",
       "[-11/84 -37/31 -34/93]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P*A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[  -5/4 -57/79  13/25]\n",
       "[ -11/9     46  -1/18]\n",
       "[-11/84 -37/31 -34/93]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L*U"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 the QR decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the QR decomposition, we must convert the matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(-5/4, -57/79, 13/25), (-11/84, -37/31, -34/93), (-11/9, 46, -1/18)]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.rows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[-1.25, -0.7215189873417721, 0.52],\n",
       " [-0.13095238095238096, -1.1935483870967742, -0.3655913978494624],\n",
       " [-1.2222222222222223, 46.0, -0.05555555555555555]]"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nbrs = [[CDF(n) for n in row] for row in A.rows()]\n",
    "nbrs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[               -1.25  -0.7215189873417721                 0.52]\n",
       "[-0.13095238095238096  -1.1935483870967742  -0.3655913978494624]\n",
       "[ -1.2222222222222223                 46.0 -0.05555555555555555]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C = Matrix(CDF, nbrs)\n",
    "C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "Q, R = C.QR()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[    1.0000000000000004  5.551115123125783e-17 1.8648277366750676e-16]\n",
       "[ 5.551115123125783e-17     1.0000000000000002 -7.632783294297951e-17]\n",
       "[1.8648277366750676e-16 -7.632783294297951e-17                    1.0]"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q*Q.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[    1.75313310577689   -31.46596530702512  -0.3047251230768725]\n",
       "[                 0.0    33.58330202196513  -0.3597878656704642]\n",
       "[                 0.0                  0.0 -0.42992880924251553]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[ -1.2500000000000002    -0.72151898734176   0.5200000000000001]\n",
       "[-0.13095238095238096    -1.19354838709677  -0.3655913978494625]\n",
       "[ -1.2222222222222223                 46.0 -0.05555555555555543]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Q*R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[               -1.25  -0.7215189873417721                 0.52]\n",
       "[-0.13095238095238096  -1.1935483870967742  -0.3655913978494624]\n",
       "[ -1.2222222222222223                 46.0 -0.05555555555555555]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 the Eigenvalue Decomposition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[-1.449367248669314?, -0.524868346991509? - 4.145964620413456?*I, -0.524868346991509? + 4.145964620413456?*I]"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "L = A.eigenvalues()\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The eigenvalue problem is equivalent to finding the roots of the characteristic polynomial and for this, the transition to algebraic numbers happens automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "D, V = A.right_eigenmatrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.449367248669314? & 0 & 0 \\\\\n",
       "0 & -0.524868346991509? - 4.145964620413456? \\sqrt{-1} & 0 \\\\\n",
       "0 & 0 & -0.524868346991509? + 4.145964620413456? \\sqrt{-1}\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.449367248669314? & 0 & 0 \\\\\n",
       "0 & -0.524868346991509? - 4.145964620413456? \\sqrt{-1} & 0 \\\\\n",
       "0 & 0 & -0.524868346991509? + 4.145964620413456? \\sqrt{-1}\n",
       "\\end{array}\\right)$$"
      ],
      "text/plain": [
       "[                       -1.449367248669314?                                          0                                          0]\n",
       "[                                         0 -0.524868346991509? - 4.145964620413456?*I                                          0]\n",
       "[                                         0                                          0 -0.524868346991509? + 4.145964620413456?*I]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "1 & 1 & 1 \\\\\n",
       "0.03664640772965092? & -0.691282990049901? + 0.04152266990464062? \\sqrt{-1} & -0.691282990049901? - 0.04152266990464062? \\sqrt{-1} \\\\\n",
       "-0.3325503262971190? & 0.4353035578097981? - 7.915394664754103? \\sqrt{-1} & 0.4353035578097981? + 7.915394664754103? \\sqrt{-1}\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "1 & 1 & 1 \\\\\n",
       "0.03664640772965092? & -0.691282990049901? + 0.04152266990464062? \\sqrt{-1} & -0.691282990049901? - 0.04152266990464062? \\sqrt{-1} \\\\\n",
       "-0.3325503262971190? & 0.4353035578097981? - 7.915394664754103? \\sqrt{-1} & 0.4353035578097981? + 7.915394664754103? \\sqrt{-1}\n",
       "\\end{array}\\right)$$"
      ],
      "text/plain": [
       "[                                           1                                            1                                            1]\n",
       "[                        0.03664640772965092? -0.691282990049901? + 0.04152266990464062?*I -0.691282990049901? - 0.04152266990464062?*I]\n",
       "[                        -0.3325503262971190?   0.4353035578097981? - 7.915394664754103?*I   0.4353035578097981? + 7.915394664754103?*I]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(V)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.449367248669314? & -0.524868346991509? - 4.145964620413456? \\sqrt{-1} & -0.524868346991509? + 4.145964620413456? \\sqrt{-1} \\\\\n",
       "-0.05311410314473801? & 0.5349840806605853? + 2.844240884324994? \\sqrt{-1} & 0.5349840806605853? - 2.844240884324994? \\sqrt{-1} \\\\\n",
       "0.4819875514693377? & -33.04542329550710? + 2.349786963655363? \\sqrt{-1} & -33.04542329550710? - 2.349786963655363? \\sqrt{-1}\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.449367248669314? & -0.524868346991509? - 4.145964620413456? \\sqrt{-1} & -0.524868346991509? + 4.145964620413456? \\sqrt{-1} \\\\\n",
       "-0.05311410314473801? & 0.5349840806605853? + 2.844240884324994? \\sqrt{-1} & 0.5349840806605853? - 2.844240884324994? \\sqrt{-1} \\\\\n",
       "0.4819875514693377? & -33.04542329550710? + 2.349786963655363? \\sqrt{-1} & -33.04542329550710? - 2.349786963655363? \\sqrt{-1}\n",
       "\\end{array}\\right)$$"
      ],
      "text/plain": [
       "[                       -1.449367248669314? -0.524868346991509? - 4.145964620413456?*I -0.524868346991509? + 4.145964620413456?*I]\n",
       "[                     -0.05311410314473801? 0.5349840806605853? + 2.844240884324994?*I 0.5349840806605853? - 2.844240884324994?*I]\n",
       "[                       0.4819875514693377? -33.04542329550710? + 2.349786963655363?*I -33.04542329550710? - 2.349786963655363?*I]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(V*D)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<html><script type=\"math/tex; mode=display\">\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.449367248669314? & -0.524868346991509? - 4.145964620413456? \\sqrt{-1} & -0.524868346991509? + 4.145964620413456? \\sqrt{-1} \\\\\n",
       "-0.05311410314473801? & 0.5349840806605853? + 2.844240884324994? \\sqrt{-1} & 0.5349840806605853? - 2.844240884324994? \\sqrt{-1} \\\\\n",
       "0.4819875514693377? & -33.04542329550710? + 2.349786963655363? \\sqrt{-1} & -33.04542329550710? - 2.349786963655363? \\sqrt{-1}\n",
       "\\end{array}\\right)</script></html>"
      ],
      "text/latex": [
       "$$\\newcommand{\\Bold}[1]{\\mathbf{#1}}\\left(\\begin{array}{rrr}\n",
       "-1.449367248669314? & -0.524868346991509? - 4.145964620413456? \\sqrt{-1} & -0.524868346991509? + 4.145964620413456? \\sqrt{-1} \\\\\n",
       "-0.05311410314473801? & 0.5349840806605853? + 2.844240884324994? \\sqrt{-1} & 0.5349840806605853? - 2.844240884324994? \\sqrt{-1} \\\\\n",
       "0.4819875514693377? & -33.04542329550710? + 2.349786963655363? \\sqrt{-1} & -33.04542329550710? - 2.349786963655363? \\sqrt{-1}\n",
       "\\end{array}\\right)$$"
      ],
      "text/plain": [
       "[                       -1.449367248669314? -0.524868346991509? - 4.145964620413456?*I -0.524868346991509? + 4.145964620413456?*I]\n",
       "[                     -0.05311410314473801? 0.5349840806605853? + 2.844240884324994?*I 0.5349840806605853? - 2.844240884324994?*I]\n",
       "[                       0.4819875514693377? -33.04542329550710? + 2.349786963655363?*I -33.04542329550710? - 2.349786963655363?*I]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "show(A*V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that $V D = A V$ or equivalently: $D = V^{-1} A V$, or $A = V D V^{-1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.4 Singular Value Decomposition"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To compute the Singular Value Decompsition, we must again use the matrix of ``CDF`` type, and work over the complex double field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, V = C.SVD()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both ``U`` and ``V`` are orthogonal matrices, that is: their transpose is the inverse matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[     1.0000000000000002   1.700029006457271e-16  1.7580338226852454e-16]\n",
       "[  1.700029006457271e-16      1.0000000000000002 -2.0469737016526324e-16]\n",
       "[ 1.7580338226852454e-16 -2.0469737016526324e-16                     1.0]"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U*U.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[     0.9999999999999991 -2.2551405187698492e-17   5.551115123125783e-17]\n",
       "[-2.2551405187698492e-17                     1.0                     0.0]\n",
       "[  5.551115123125783e-17                     0.0      0.9999999999999996]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V*V.transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[46.036789227343746                0.0                0.0]\n",
       "[               0.0 1.3707314283740828                0.0]\n",
       "[               0.0                0.0 0.4011228818407111]"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[      46.03678922734373   4.163336342344337e-17 -2.0816681711721685e-16]\n",
       "[   8.35909032964599e-15       1.370731428374083 -3.3306690738754696e-16]\n",
       "[  1.710491557421756e-15  2.7755575615628914e-17     0.40112288184071115]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U.transpose()*C*V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Matrices over Rings"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hermite normal form is computed with a unimodular matrix.  We make a random integer matrix with numbers that are 3 decimal places long, that is: in the range 100 and 999."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[503 842 787]\n",
       "[599 828 414]\n",
       "[857 948 500]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = random_matrix(ZZ, 3, x=100, y=999)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[       1        0  5995442]\n",
       "[       0        2 12585079]\n",
       "[       0        0 27081514]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "H, U = A.hermite_form(transformation=True)\n",
    "H"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[  15690  -27092    9727]\n",
       "[  32935  -56869   20418]\n",
       "[  70872 -122375   43937]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[       1        0  5995442]\n",
       "[       0        2 12585079]\n",
       "[       0        0 27081514]"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U*A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not all integer matrices are invertible, but unimodular matrices are those with inverses with integer coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "D, U, V = A.smith_form()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[       1        0        0]\n",
       "[       0        1        0]\n",
       "[       0        0 54163028]"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[       1        0        0]\n",
       "[       0        1        0]\n",
       "[       0        0 54163028]"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "U*A*V"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Resultants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resultant is the determinant of the Sylvester matrix, used to eliminate variables in polynomial systems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1/2*x^2 + 2*x*y + 1/2*y^2 + 1/3*x, 5*y^2 - 3*x + y - 3)"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P.<x,y> = PolynomialRing(QQ)\n",
    "p = P.random_element(degree=2, terms=5)\n",
    "q = P.random_element(degree=2, terms=5)\n",
    "(p, q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resultant of the system defined by ``p`` and ``q`` with respect to ``x`` is the condition on ``y`` for the system to have solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-25/2*y^4 + 25*y^3 + 30*y^2 - 14*y - 15/2"
      ]
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p.resultant(q,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-25/2*y^4 + 25*y^3 + 30*y^2 - 14*y - 15/2"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "q.resultant(p,x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[           -3 5*y^2 + y - 3             0]\n",
       "[            0            -3 5*y^2 + y - 3]\n",
       "[         -1/2     2*y + 1/3       1/2*y^2]"
      ]
     },
     "execution_count": 43,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = q.sylvester_matrix(p, x)\n",
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-25/2*y^4 + 25*y^3 + 30*y^2 - 14*y - 15/2"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "det(S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the roots of this polynomial will give all y-coordinates of the solutions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.3",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}