path: root/src/Lecture5/notebook/7-SageAlgebra.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This lecture's notes are in a different format: the presentations for the $\\LaTeX$ part were made with $\\LaTeX$, so this one is made with Sage, or rather with the [Jupyter Notebook](https://jupyter.org/).\n",
    "\n",
    "# The Jupyter Notebook\n",
    "**Reference:** [[1](https://jupyter.org/documentation)]\n",
    "\n",
    "The Jupyter Notebook is one of the default interfaces for SageMath, along with the command line interface. You can access it via web browser, but it is running locally on your device (notice the strange url: `http://localhost:8888/notebooks...`).\n",
    "\n",
    "You can create a new notebook by clicking on `New > SageMath 9.2`. You can also create a Python 3 notebook to write Python code.\n",
    "\n",
    "Jupyter saves and reads files in the `.ipynb` format. If you download the file for this lecture you can open it and follow the examples interactively.\n",
    "\n",
    "## Cells\n",
    "\n",
    "The notebook contains one or more *interactive cells* that you can run, like this one below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2/5"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Exercise: modify this cell to use the print() command\n",
    "2+2\n",
    "2/5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you are reading this from Jupyter rather than from the pdf file, you can edit the cell above and run it again. You can also add more cells by selecting `Insert` from the menu bar.\n",
    "\n",
    "Notice that only the last statement produces an output. You can force anything to be written as output with the `print()` command, which works like in Python. As an exercise, try to modify the cell above to provide more output!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markdown\n",
    "\n",
    "[Markdown](https://en.wikipedia.org/wiki/Markdown) is a simple markup language - think of LaTeX or html, but much simpler.\n",
    "You can add text to your notebook with Markdown cells by selecting `Cell > Cell Type > Markdown`.\n",
    "\n",
    "You can also include some LaTeX code in Markdown cells, with dollar signs $ or align environments:\n",
    "\n",
    "\\begin{align*}\n",
    "\\frac{(x+y)^2}{x+1} = \\frac{x^2+y^2}{x+1}\n",
    "\\end{align*}\n",
    "\n",
    "When you are done writing a Markdown cell, you can run it to see the well-formatted text. To edit the text again, double-click on the cell. Try doing it now to fix the formula above!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Symbolic expressions\n",
    "\n",
    "**Reference:** [[2](https://doc.sagemath.org/html/en/reference/calculus/sage/symbolic/expression.html)]\n",
    "\n",
    "Now, let's get started with Sage. One thing you might want to do is manipulating symbolic expressions, like the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x == -sqrt(6) - 1, x == sqrt(6) - 1]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f =  x^2 + 2*x - 5 == 0\n",
    "solve(f,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the single `=` is part of an assignment, as in Python: we are *assigning* to the variable `f` the value `x^2 + 2*x - 5 >= 0`, which in this case is an equation, so it contains the symbol `==`. Keep in mind the difference between the two!\n",
    "\n",
    "**Exercise:** change the code above to solve the corresponding inequality $x^2+2x-5\\geq 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mathematical variables\n",
    "\n",
    "Last time we saw what *variables* are in Python, and that they are a little bit different from the *Mathematical variables* that you use in Mathematics. In Sage, both concepts are present, but they are still distinct. For example in the cell above `f` is a variable in the sense of computer science, while `x` is a Mathematical variable.\n",
    "\n",
    "If you want to use Mathematical variables other than `x`, you first need to *declare* them with the `var()` command:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[y == -1/2*x - 1/2*sqrt(x^2 + 2*x + 9) - 1/2, y == -1/2*x + 1/2*sqrt(x^2 + 2*x + 9) - 1/2]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('y')\n",
    "solve(y^2 + (x+1)*y - 2 == 0, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try removing the first line in the cell above and see what error you get!\n",
    "\n",
    "Here is another example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[x == -1/2*a - 1/2*sqrt(a^2 - 4*b), x == -1/2*a + 1/2*sqrt(a^2 - 4*b)]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "var('a', 'b')\n",
    "f = x^2+a*x+b\n",
    "solve(f,x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some common constants are [already defined](https://doc.sagemath.org/html/en/reference/calculus/sage/symbolic/expression.html) in Sage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e^(pi*I)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will study symbolic expressions more in detail next time, in the context of calculus/analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic rings and fields\n",
    "\n",
    "**References:** [[3](https://doc.sagemath.org/html/en/reference/rings_standard/index.html)]\n",
    "[[4](https://doc.sagemath.org/html/en/reference/rings_numerical/index.html)]\n",
    "[[5](https://doc.sagemath.org/html/en/reference/finite_rings/index.html)]\n",
    "\n",
    "As you should know, a *field* is a Mathematical structure with two operations, addition and multiplication, which respect certain rules (distributivity, associativity, commutativity...). Some examples of fields are the Rational numbers $\\mathbb Q$, the Real numbers $\\mathbb R$ and the Complex numbers $\\mathbb C$, but there are many more. As you should also know, a *(commutative) ring* is like a field, except not all elements different from $0$ need have a multiplicative inverse. For example the integers $\\mathbb Z = \\{ \\dots, -1, 0, 1, 2, \\dots\\}$ are a ring, but not a field.\n",
    "\n",
    "These structures are already implemented in Sage. Some of the most common are listed in the following table:\n",
    "\n",
    "|Mathematical object|Math symbol|Sage name|\n",
    "|------------------:|:---------:|:--------|\n",
    "|Integers|$\\mathbb Z$|`ZZ`|\n",
    "|Rational numbers|$\\mathbb Q$|`QQ`|\n",
    "|Real numbers|$\\mathbb R$|`RR`|\n",
    "|Complex numbers|$\\mathbb C$|`CC`|\n",
    "|Integers modulo $n$|$\\mathbb Z/n\\mathbb Z$|`Integers(n)`|\n",
    "|Finite fields|$\\mathbb F_p$|GF(p)|\n",
    "|$\\dots$|$\\dots$|$\\dots$|"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you write a number or an expression, Sage will figure out where it \"lives\", choosing the most restrictive interpretation possible. For example `3` will be interpreted to be an integer, even if it is also a rational number, a real number and a complex number."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parents and coercion\n",
    "**Reference:** [[6](https://doc.sagemath.org/html/en/tutorial/tour_coercion.html)]\n",
    "\n",
    "You can check where an object \"lives\" with the `parent()` command. It works more or less like the Python command `type()`, but it gives a more Mathematically inclined answer. Check the reference link [6] above if you want more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Rational Field"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Edit this cell to find out the type of other objects that we used\n",
    "parent(3/5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sometimes Sage does not give you the best possible interpretation, so you can force something to be interpreted as living in a smaller ring as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Symbolic Ring\n",
      "Integer Ring\n"
     ]
    }
   ],
   "source": [
    "minus_one = e^(pi*I)\n",
    "minus_one_coerced = ZZ(e^(pi*I))  # coercion\n",
    "print(parent(minus_one))\n",
    "print(parent(minus_one_coerced))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Remark.** Notice that there is a fundamental difference between the rings `RR` and `CC` and all the others in the table above: the real and complex numbers are *approximated*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3\n",
      "3.00000000000000\n"
     ]
    }
   ],
   "source": [
    "print(QQ(3))\n",
    "print(RR(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also choose the precision of this approximation using the alternative name `RealField`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Real Field with 53 bits of precision\n",
      "Real Field with 1000 bits of precision\n"
     ]
    }
   ],
   "source": [
    "print(RR)\n",
    "print(RealField(prec=1000))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Polynomial rings\n",
    "\n",
    "**Reference:** [[7](https://doc.sagemath.org/html/en/reference/polynomial_rings/index.html)]\n",
    "\n",
    "If you want to work with polynomials over a certain ring it is better to use this specific construction, rather than the symbolic expressions introduced above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Multivariate Polynomial Ring in x, y, z over Real Field with 53 bits of precision"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polring.<x,y,z> = RR[]  # Alternative: polring.<x,y,z> = PolynomialRing(RR)\n",
    "polring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can use as many variables as you like, and you can replace `RR` with any ring. In the example above `polring` is just the name of the variable (in the computer science sense) associated with this polynomial ring.\n",
    "\n",
    "## Operations on polynomials\n",
    "\n",
    "The usual Mathematical operations are available on polynomial rings, including Euclidean division `//` and remainder `%`. There is also the single-slash division `/`, but the result may not be a polynomial anymore.\n",
    "\n",
    "**Exercise:** use the `parent()` command to find out what the quotient of two polynomials is.\n",
    "\n",
    "**Question:** what happens if you remove the first line in the cell below? What if we used the variable `y` instead of `x`?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x + 1\n",
      "-4\n",
      "(x^2 + 2*x - 3)/(x + 1)\n"
     ]
    }
   ],
   "source": [
    "polring.<x> = QQ[]\n",
    "p = x^2 + 2*x - 3   # Don't forget * for multiplication!\n",
    "q = p // (x+1)\n",
    "r = p %  (x+1)\n",
    "f = p /  (x+1)\n",
    "print(q)\n",
    "print(r)\n",
    "print(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can do more complex operations. Try out `roots()` and `factor` in the cell below.\n",
    "\n",
    "**Remark.** Notice how the result can change substantially if you change the base ring.\n",
    "\n",
    "**Remark.** [Factorizations](https://doc.sagemath.org/html/en/reference/structure/sage/structure/factorization.html) are a particular object in Sage. They are kinda like a list, but not really. You can get a list of pairs (factor, power) with `list(factor(f))`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(t + 1) * (t^2 - 3) * (t^2 + 1)\n",
      "[(-1, 1)]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(y + 1) * x"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "polring_onevar.<t> = QQ[]\n",
    "\n",
    "f = t^5 + t^4 - 2*t^3 - 2*t^2 - 3*t - 3\n",
    "print(factor(f))\n",
    "print(f.roots())  # Result: list of pairs (root,multiplicity)\n",
    "\n",
    "polring_manyvar.<x,y,z> = QQ[]\n",
    "factor(x*y+x)\n",
    "\n",
    "# The following line gives an error, because the polynomial\n",
    "# is understood to possibly have many variables:\n",
    "#(x^2-1).roots()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matrices and vectors\n",
    "\n",
    "**References:** [[8](https://doc.sagemath.org/html/en/reference/matrices/index.html)], but in particular the subections [[9](https://doc.sagemath.org/html/en/reference/matrices/sage/matrix/docs.html)] and [[10](https://doc.sagemath.org/html/en/reference/matrices/sage/matrix/matrix2.html)]\n",
    "\n",
    "In Sage you can easily manipulate matrices and vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[   1    2    3]\n",
      "[   0    0    1]\n",
      "[   4   -3 22/7] \n",
      "\n",
      "[1/2   0   0]\n",
      "[  7   0   0]\n",
      "[  1   1   1] \n",
      "\n",
      "(3/2, 21, 6) \n",
      "\n",
      "[  -7/2    -10   80/7]\n",
      "[    17     -4   15/7]\n",
      "[ 241/7  -18/7 869/49] \n",
      "\n",
      "Rank of A = 3\n",
      "Rank of B = 2\n"
     ]
    }
   ],
   "source": [
    "A = matrix([[1,2,3],[0,0,1],[4,-3,22/7]])\n",
    "B = matrix([[1/2,0,0],[7,0,0],[1,1,1]])\n",
    "v = vector([3,4,-1])\n",
    "\n",
    "print(A, \"\\n\")  # \\n just means \"newline\"\n",
    "print(B, \"\\n\")\n",
    "print(B*v, \"\\n\")\n",
    "print(A^2 + 2*B - A*B, \"\\n\")\n",
    "\n",
    "print(\"Rank of A =\", rank(A)) # You can also use A.rank()\n",
    "print(\"Rank of B =\", rank(B))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** in the cell above, compute the determinant, inverse and characteristic polynomial of the matrix `A`. *Hint: look at the reference [10] above (the functions are listed in alphabetic order).*\n",
    "\n",
    "As for polynomials, you can specify where a matrix or a vector lives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Full MatrixSpace of 2 by 2 dense matrices over Complex Field with 53 bits of precision"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = matrix(CC, [[0,1],[1,0]])\n",
    "parent(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also solve linear systems and compute eigenvalues and eigenvectors of a matrix\n",
    "\n",
    "**Warning.** In linear algebra there are distinct concepts of *left* and *right* eigenvalues (and eigenvector). The one you know is probably that of **right** eigen-{value,vector}, that is an element $\\lambda$ of the base field and a non-zero vector $\\mathbf v$ with $A\\mathbf v=\\lambda\\mathbf v$. The other concept corresponds to the equality $\\mathbf v^TA=\\lambda \\mathbf v$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.289916349448506, 0.0241596957873755)"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix(RR, [[sqrt(59),32],[-1/4,3]])\n",
    "v = vector(RR, [3,0])\n",
    "A.solve_right(v) # Solve Ax=v. Alternative: A \\ v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[\n",
       "(-0.3722813232690144?, Vector space of degree 2 and dimension 1 over Algebraic Field\n",
       "User basis matrix:\n",
       "[                   1 -0.6861406616345072?]),\n",
       "(5.372281323269015?, Vector space of degree 2 and dimension 1 over Algebraic Field\n",
       "User basis matrix:\n",
       "[                 1 2.186140661634508?])\n",
       "]"
      ]
     },
     "execution_count": 64,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = Matrix(QQ, [[1,2],[3,4]])\n",
    "A.eigenspaces_right() # Also: A.eigenvalues(), A.eigenvectors_right()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also extract a specific submatrix by selecting only some rows and columns, with a syntax similar to that of Python's lists. Check out more examples in the reference [9] above, and try them in the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-14   2   0  -1   1  -2  -1]\n",
      "[  0  -8   0   9  -2  11   1]\n",
      "[  0   3   1  -1   1   1 221]\n",
      "[ -1   2   1 -25 -10   4   0]\n",
      "[ -3   0   0   2  16  -1  -2]\n",
      "[  1  -3   3 -41   1   0   0]\n",
      "[ -2   1   0   0  -6   2  12] \n",
      "\n",
      "[ 0  9 -2]\n",
      "[ 1 -1  1] \n",
      "\n",
      "[-14   2   0  -1   1  -2  -1] \n",
      "\n",
      "[-14   2   0  -1   1]\n",
      "[  1  -3   3 -41   1]\n",
      "[  0   3   1  -1   1]\n"
     ]
    }
   ],
   "source": [
    "A = MatrixSpace(ZZ, 7).random_element()\n",
    "print(A, \"\\n\")\n",
    "print(A[1:3,2:5], \"\\n\") # Rows from 1 to 3, columns from 2 to 5\n",
    "print(A[0,0:], \"\\n\")    # First row, all columns\n",
    "print(A[[0,5,2],0:5])   # Rows 0, 5 and 2 (in this order) and columns 0 to 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** write a sage function that computes the determinant of an $n\\times n$ matrix $A=(a_{ij})$ using Laplace's rule by the first row, that is \n",
    "\\begin{align*}\n",
    "    \\operatorname{det}A = \\sum_{j=1}^n (-1)^ja_{0j}M_{0j}\n",
    "\\end{align*}\n",
    "where $M_{0j}$ is the determinant of the $(n-1)\\times(n-1)$ matrix obtained by removing the $0$-th row and the $j$-th column from $A$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_det(A):\n",
    "    if not A.is_square():\n",
    "        print(\"Error: matrix is not square\")\n",
    "    \n",
    "    n = A.nrows()  # size of the matrix\n",
    "    \n",
    "    # Continue from here!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Number Theory\n",
    "\n",
    "**Reference:** [[11](https://doc.sagemath.org/html/en/reference/rings_standard/sage/rings/integer.html)]\n",
    "\n",
    "Sage includes a large library of functions for computing with the integers, see the link above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3^2 * 3607 * 3803\n",
      "True\n",
      "True\n",
      "619703040\n",
      "9\n",
      "13548070123626141\n"
     ]
    }
   ],
   "source": [
    "n = 123456789\n",
    "m = 987654321\n",
    "p = 3607\n",
    "\n",
    "print(factor(n))\n",
    "print(is_prime(p))\n",
    "print(p.divides(n))\n",
    "print(euler_phi(m))\n",
    "print(gcd(n, m))\n",
    "print(lcm(n, m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Primes\n",
    "\n",
    "**Reference:** [[12](https://doc.sagemath.org/html/en/reference/sets/sage/sets/primes.html)]\n",
    "\n",
    "The set of prime numbers is called `Primes()`. It is like an infinite list: for example you can get the one-millionth prime number or you can use this list to create other lists. You can also check what the first prime number larger than a given number is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Set of all prime numbers: 2, 3, 5, 7, ...\n",
      "31 15485867\n",
      "47\n",
      "[79, 83, 89, 97]\n"
     ]
    }
   ],
   "source": [
    "PP = Primes()\n",
    "print(PP)\n",
    "print(PP[10], PP[10^6])\n",
    "print(PP.next(44))\n",
    "\n",
    "First_Thousand_Primes = PP[0:1000]\n",
    "print([p for p in First_Thousand_Primes if p < 100 and p > 75])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Chinese remainder theorem (CRT)\n",
    "\n",
    "We say that two integers $a$ and $b$ are *congruent* modulo another integer $n>0$ if they have the same remainder when divided by $n$. We denote this by $a\\equiv b\\pmod n$, or in Python/Sage syntax `a % n == b % n`.\n",
    "\n",
    "The Chinese remainder theorem states that if $a,b\\in\\mathbb Z$ and $n,m\\in \\mathbb Z_{>0}$ are such that $\\gcd(n,m)=1$ then the system of congruences\n",
    "\n",
    "\\begin{align*}\n",
    "\\begin{cases}\n",
    "    x \\equiv a \\pmod n\\\\\n",
    "    x \\equiv b \\pmod m\n",
    "\\end{cases}\n",
    "\\end{align*}\n",
    "\n",
    "has exactly one solution modulo $mn$. This means that there is one and only one number $x$ with $0\\leq x<mn$ such that $x\\equiv a\\pmod n$ and $x\\equiv b\\pmod m$.\n",
    "\n",
    "The procedure to find such a number is not too hard to describe (you might see it in an algebra or number theory course), but it can be a bit long. Luckily, Sage can do this for you:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "74306 2 798\n"
     ]
    }
   ],
   "source": [
    "a = 2\n",
    "b = -1\n",
    "n = 172\n",
    "m = 799\n",
    "\n",
    "if gcd(n,m) != 1:\n",
    "    print(\"The numbers are not comprime, I can't solve this!\")\n",
    "else:\n",
    "    x = crt(a, b, n, m)\n",
    "    print(x, x%n, x%m)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise.** There is a more general version of the Chinese remainder theorem which says that if $a_0, a_1, \\dots, a_k\\in\\mathbb Z$ and $n_0, n_2, \\dots, n_k\\in\\mathbb Z_{>0}$ are such that $\\gcd(n_i, n_j)=1$ for $i\\neq j$, then the system of congruences\n",
    "\n",
    "\\begin{align*}\n",
    "\\begin{cases}\n",
    "    x \\equiv a_0 \\pmod {n_0}\\\\\n",
    "    x \\equiv a_1 \\pmod {n_1}\\\\\n",
    "    \\dots \\\\\n",
    "    x \\equiv a_k \\pmod {n_k}\n",
    "\\end{cases}\n",
    "\\end{align*}\n",
    "\n",
    "has exactly one solution modulo $\\prod_{i=0}^kn_i$. Use the `crt()` function to find a solution to such a system.\n",
    "*Hint: start by running the command `help(crt)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "#help(crt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cryptography: RSA\n",
    "\n",
    "[Cryptography](https://en.wikipedia.org/wiki/Cryptography) is the discipline that studies methods to communicate secrets in such a way that any unauthorized listener would not be able to understand the message.\n",
    "\n",
    "A simple cryptographic protocol could be changing every letter of your text following a fixed scheme (or *cypher*), for example by turning every A into a B, every B into a C and so on. However this is not a very secure method, for many reasons. One of them is that at some point the people who want to communicate need to agree on what method to use, and anyone listening to that conversation would be able to decypher every subsequent conversation. A public-key cryptographic protocol solves this problem.\n",
    "\n",
    "## Public-key cryptography\n",
    "\n",
    "Public-key cryptographic protocols, such as RSA, work like this: there are two keys, a *private* key that is only known to person A (traditionally called Alice in every example), and a *public* key that does not need to be secret.\n",
    "\n",
    "The public key is used to *encrypt* the message (that is to \"lock\" it, or \"hyde\" it), but one needs the private key to *decrypt* it. Imagine having two keys for your door, but one can only be used to lock it, while the other only to open it.\n",
    "\n",
    "The message exchange works like this: suppose that person B (Bob) wants to send a secret message to Alice. Then Alice secretely generates a private and a public key and sends only the public one to Bob. Now Bob encrypts the message and sends it to Alice, who can use her private key to decrypt it. Even if Eve (short for *eavesdropper*, an unauthorized listener) listens to every message exchanged, she won't be able to decypher the secret: the private key has never left Alice's house!\n",
    "\n",
    "Notice that such a protocol is *asymmetric*: if Alice wanted to send a secret to Bob in reply, Bob would need to generate a pair of keys of his own.\n",
    "\n",
    "Let's see how we can do this in practice, using number theory!\n",
    "\n",
    "## RSA\n",
    "\n",
    "As many other cryptography protocols, RSA is based on a Mathematical process that is easy to do in one direction, but very hard to invert. In this case the hard process is integer factorization, that is decomposing an integer number as a product of primes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True True False\n"
     ]
    }
   ],
   "source": [
    "p = 100003100019100043100057100069\n",
    "q = 100144655312449572059845328443\n",
    "n = p*q\n",
    "print(is_prime(p), is_prime(q), is_prime(p*q))\n",
    "\n",
    "# Use the command below to see how long it takes\n",
    "#timeit(\"factor(n)\", number=1, repeat=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to generate the keys, Alice picks a number $n$ which is the product of two large primes $p$ and $q$ of more or less the same size. Finding such primes is relatively easy compared to factoring the number $n$ she obtained. Then she computes the Euler totient $\\varphi(n)=(p-1)(q-1)$ of $n$, which she can do because she knows that $n=pq$ - it would be impossible otherwise!\n",
    "\n",
    "Then Alice can compute two integers $(d,e)$ such that $de\\equiv 1\\pmod{\\varphi(n)}$. She will send the numbers $n$ and $d$ to Bob and keep $e$ secret. In this case the public key is the pair $(n,d)$, while $e$ is the private key.\n",
    "\n",
    "Of course, she does all of this using Sage!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(419199544978969, 235530823946467, 80799425863927)"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def two_large_primes():\n",
    "    p, q = 0, 0\n",
    "    # We make sure that they are different\n",
    "    while p == q:\n",
    "        p = Primes()[randint(10^6, 2*10^6)]\n",
    "        q = Primes()[randint(10^6, 2*10^6)]\n",
    "    return p, q\n",
    "\n",
    "def random_unit_mod(N):\n",
    "    R = Integers(N)\n",
    "    d = R(0)\n",
    "    # We make sure that it is invertible\n",
    "    while not d.is_unit():\n",
    "        d = R.random_element()\n",
    "    return d\n",
    "\n",
    "def Alice_generate_keys():\n",
    "    p, q = two_large_primes()\n",
    "    n = p*q\n",
    "    phi_n = (p-1)*(q-1) # euler_phi(n) is slow!\n",
    "    \n",
    "    d = random_unit_mod(phi_n)\n",
    "    e = d^-1\n",
    "    return n, d, e\n",
    "\n",
    "Alice_generate_keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, how does Bob encrypt his message? Let's say he wants to send to Alice the number $m$ with $1<m<n$ (In practice he would like to send her some text with emojis, or maybe a voice message; but for computers everything is a number, and there are different ways to translate any sort of information to a number. He just chooses one of the many standard methods that already exist, no cryptography is needed in this step. If the message $m$ is too long, he can split it up in some pieces and repeat the process multiple times.)\n",
    "\n",
    "Now he computes $m^d\\pmod n$ and sends it back to Alice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "149461597163501"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Bob_encrypt(m, n, d):\n",
    "    R = Integers(n)\n",
    "    return R(m)^d  # Assume that n is large enough\n",
    "    \n",
    "message = 42424242\n",
    "Bob_encrypt(message, 419199544978969, 235530823946467)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since $de\\equiv 1\\pmod{\\varphi(n)}$, it follows that $(m^d)^e\\equiv m\\pmod n$ (see [Wikipedia: Euler's theorem](https://en.wikipedia.org/wiki/Euler%27s_theorem)). So for Alice it is very easy to get back the original message:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "42424242"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def Alice_decrypt(m_encrypted, n, e):\n",
    "    R = Integers(n)\n",
    "    return R(m_encrypted)^e\n",
    "\n",
    "Alice_decrypt(149461597163501, 419199544978969, 80799425863927)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another assumption on which RSA relies is that even if one knows $M=m^e$ and $e$, extracting the $e$-th root of $M$ modulo $n$ (and thus obtaining $m$) is very hard. Currently the best known way to do this is by factorizing $n$ first, which is considered to be a very hard problem. However, there is no proof that faster algorithms can't be devised.\n",
    "\n",
    "Moreover, one day we will overcome the current technological difficulties and quantum computers will be available. Quantum computers are not just \"more powerful\" than classical hardware, but they work based on completely different logical foundations and they make the factorization problem much easier to solve: for example [Shor's algorithm](https://en.wikipedia.org/wiki/Shor%27s_algorithm) takes advantage of this different logic and can factorize numbers quickly, if run on a quantum computer.\n",
    "\n",
    "To this day the largest number factorized with a quantum computer is $21=3\\times 7$. Nonetheless, quantum-safe cryptography protocols (i.e. based on problems that are hard to solve also with quantum computers) have already been developed."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}