{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"# Basic Algebraic Operations\n",
"\n",
"This worksheet provides an introduction to representations of basic mathematical objects and how arithmetic can be performed on them.\n",
"\n",
"The fundamental objects we'll consider are integers and polynomials.\n",
"\n",
"The fundamental arithmetic operations we'll consider are addition/subtraction, multiplication, and division."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The “Integer” Datatype"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Programming languages such as C allow declaring \"int\" datatypes.\n",
"\n",
"But \"int\" is not actually a general integer type. On machines where \"int\"s are represented by 4 bytes (32-bit) they can only represent integers in the range $-2{,}147{,}483{,}648$ to $2{,}147{,}483{,}647$.\n",
"On 64-bit machines ints can only represent integers up to around $\\pm9.2$ quintillion.\n",
"\n",
"Just as in mathematics, we would like our computations to have no such limitation."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"123^45"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"But wait a minute: the CPU in the computer that we're running on is limited to performing computations on integers of a specific \"word size\" (typically 64 bits, i.e., 8 bytes).\n",
"\n",
"So how is it possible to perform computations on integers of _arbitrary_ size?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Representation of Integers"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"The main idea is to split an integer of arbitrary size into \"word size\" chunks. The base-$2^{64}$ representation of an integer (assuming a word size of size 64 bits) is an appropriate way to achieve this.\n",
"\n",
"Any integer $a$ can be written in the form\n",
"\n",
"$$ a = (-1)^s \\sum_i a_i \\, 2^{64 i} \\notag $$\n",
"\n",
"where $a_i$ are the \"base-$2^{64}$ digits\" of $a$ (so that $0 \\leq a_i < 2^{64}$) and $s$ is 0 or 1, depending on if $a$ is positive or negative.\n",
"\n",
"In other words, arbitrary integers can be represented by a list $[s, a_0, a_1, \\dotsc, a_{n-1}]$. Here $n$ is the \"length\" of an integer and is given by $\\lfloor\\log_{2^{64}}(a)\\rfloor+1$. This expression is a bit unwiedly—to make our lives easier we often analyze things \"asymptotically\". In this case, it is sufficient to know that $a$ is represented by a list of length $O(\\log(a))$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Let's see an example of this representation:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = 123^45; show(a)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"A = a.digits(base=2^64); show(A)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"b = add(A[i]*2^(64*i) for i in range(len(A))); show(b)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(a==b)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Representation of Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Polynomials have similar properties to integers in many (but not all) aspects.\n",
"\n",
"A polynomial in the variable $x$ is of the form $f_0 + f_1 x + \\dotsb + f_n x^n$ where the $f_i$ are the **coefficients** of the polynomial and $n$ is the **degree** of the polynomial (the final coefficient $f_n$ is assumed to not be zero).\n",
"\n",
"If the coefficients are from $\\mathbb Z$ (integers) then we write $f\\in\\mathbb Z[x]$. Other possible domains for the coefficients include $\\mathbb Q$ (rational numbers), $\\mathbb R$ (real numbers), or $\\mathbb C$ (complex numbers). More generally, any algebraic structure that supports addition and multiplication can be used (such structures are called rings)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f = 1 + 2*x + 3*x^2\n",
"show(f)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"A polynomial can be represented as a list of its coefficients $[f_0, f_1, \\dotsc, f_n]$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f.list()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is not the only possible representation, but it suffices for now."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Specifying the coefficient ring in Sage\n",
"\n",
"The polynomial $f=1+2x+3x^2$ defined above is an expression in terms of $x$. By default Sage treats `x` as a _symbolic variable_ in what Sage calls the _symbolic ring_.\n",
"\n",
"Sage also allows you to control the ring of coefficients to use when constructing and performing operations on polynomials. The following line enables defining polynomials in $x$ over the integers:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"R. = ZZ['x']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This line actually defines two things:\n",
"\n",
"* The Sage variable `R` to be $\\mathbb Z[x]$ (the ring of polynomials in $x$)\n",
"* The Sage variable `x` to be the symbolic variable `'x'`, the indeterminant of `R`\n",
"\n",
"The above definition can also be replaced with `R. = ZZ[]` because if an explicit symbolic variable is not given Sage will match the indeterminant symbolic variable with the name given to the indeterminant. The syntax `R. = PolynomialRing(ZZ)` also works.\n",
"\n",
"Note that for some operations the underlying polynomial ring will be important.\n",
"For example, when $x$ is the indeterminant of $R=\\mathbb{Z}[x]$ the polynomial $x^2+1$ does not factor—but if $R=\\mathbb{C}[x]$ then the polynomial $x^2+1$ does factor as $(x-i)(x+i)$.\n",
"\n",
"After defining the polynomial ring as above, `R` can be used as a function which creates polynomials in the ring with a given list of coefficients:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"R([1,2,3])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Addition of Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Addition is one of the simplest arithmetic operations that can be performed on polynomials.\n",
"\n",
"Suppose $f = \\sum_{i=0}^n f_i \\, x^i$ and $g = \\sum_{i=0}^n g_i \\, x^i$ and let $h = f + g$ be their sum. We have $h = \\sum_{i=0}^n (f_i + g_i) \\, x^i$.\n",
"\n",
"In other words, we can add their coefficients pointwise. This leads to a simple algorithm for addition:\n",
"\n",
" for i from 0 to n do\n",
" h[i] = f[i] + g[i]\n",
" \n",
"This uses $n+1 = O(n)$ coefficient additions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Addition of Integers"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"This is similar to the polynomial case but already we'll see that the integer arithmetic case is a bit more difficult. For simplicity, we'll assume the numbers are positive and of the same length.\n",
"\n",
"Suppose that $a = \\sum_{i=0}^n a_i \\, 2^{64 i}$ and $b = \\sum_{i=0}^n b_i \\, 2^{64 i}$ and let $c = a + b$ be their sum. We have that $c = \\sum_{i=0}^n (a_i + b_i) \\, 2^{64 i}$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"### Is this a valid final answer?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"The summation is correct, but it may be the case that $a_i + b_i \\geq 2^{64}$ and therefore that $a_i+b_i$ would not be a valid base-$2^{64}$ digit."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Addition of Integers II"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"First, we can represent $a_0 + b_0 = c_0 + \\gamma_0\\cdot2^{64}$ where $0\\leq c_0<2^{64}$ and $\\gamma_0$ is either $0$ or $1$.\n",
"\n",
"Second, we can represent $a_1+b_1+\\gamma_0 = c_1 + \\gamma_1\\cdot2^{64}$ where $0\\leq c<2^{64}$ and $\\gamma_1$ is either $0$ or $1$. This also works for each remaining coefficient.\n",
"\n",
"This results in the following addition algorithm:\n",
"\n",
" gamma[-1] = 0\n",
" for i from 0 to n+1 do\n",
" gamma[i] = 0\n",
" c[i] = a[i] + b[i] + gamma[i-1]\n",
" if c[i] >= 2^64 then\n",
" gamma[i] = 1\n",
" c[i] = c[i] - 2^64"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"### What is the asymptotic complexity in word size additions?"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"This uses $O(n)$ word size additions. There are $O(n)$ loop iterations and a constant number of additions on each iteration."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Relationship Between Integers and Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Note the similarity between the representation $a = \\sum_{i=0}^n a_i \\, 2^{64 i}$ for (positive) integers and $f = \\sum_{i=0}^n f_i x^i$ for polynomials.\n",
"\n",
"Even though integers are perhaps more \"familiar\" we tend to prioritize the polynomial case in this course—this makes life easier!\n",
"\n",
"Often algorithms that operate on polynomials can be adapted (with a bit of massaging) to work with integers. Not always, however: some algorithms only work on polynomials."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Multiplication of Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Suppose $f = \\sum_{i=0}^n f_i \\, x^i$ and $g = \\sum_{i=0}^m g_i \\, x^i$ and let $h = f \\cdot g$ be their product. Then $h = \\sum_{i=0}^n (f_i\\cdot g) x^i$ and this implies that the product can be found by multiplying $g$ by the coefficients of $f$ and then summing the results together.\n",
"\n",
"This is similar to the integer multiplication algorithm you may have learned in elementary school.\n",
"\n",
" for i from 0 to n do\n",
" h[i] = f[i]*g*x^i\n",
" h = add(h[i] for i in range(n))"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"### Example"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"f = 1 + 2*x + 3*x^2\n",
"F = f.list() \n",
"g = 4 + 5*x + 6*x^2\n",
"\n",
"show(f.list())\n",
"show(g.list())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p = [0, 0, 0]\n",
"for i in range(3):\n",
" p[i] = expand(F[i]*g*x^i)\n",
" show(p[i].list())"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"h = add(p[i] for i in [0..2]); show(h)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"show(h==f*g)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Cost in coefficient multiplications and additions"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"First, line 2 is run $n+1$ times.\n",
"\n",
"Line 2 requires $m+1$ multiplications each time it runs. Note that the multiplication by $x^i$ does not require any coefficient multiplications. Instead, multiplying $f_i\\cdot g$ by $x^i$ adds $i$ zeros to the front of its list of coefficients, i.e.,\n",
"\n",
"$$ [f_i\\cdot g_0, f_i\\cdot g_1, \\dotsc, f_i\\cdot g_m] \\cdot x^i = [\\underbrace{0,0,\\dotsc,0}_{\\text{$i$ zeros}},f_i\\cdot g_0, f_i\\cdot g_1, \\dotsc, f_i\\cdot g_m] . $$\n",
"\n",
"Thus, line 2 requires $(n+1)\\cdot(m+1)=O(n\\cdot m)$ coefficient multiplications.\n",
"\n",
"Line 3 adds together $n+1$ polynomials of maximum degree $n+m$. Thus, this uses $O(n(n+m))$ cofficient additions.\n",
"\n",
"If you are careful and only add together coefficients of the polynomials that are nonzero this actually requires exactly $n\\cdot m$ coefficient additions (assuming all $n+1$ coefficients of $f$ and all $m+1$ coefficients of $g$ are nonzero)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multiplication of integers"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"A similar algorithm works for multiplication of integers.\n",
"\n",
"It requires $O(n\\cdot m)$ word size additions and multiplications in order to multiply an integer of length $n$ with an integer of length $m$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Division of Polynomials"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Given two polynomials $a$ and $b$ the _division with remainder_ problem is to find $q$ and $r$ (quotient and remainder) which satisfy\n",
"\n",
"$$ a = q\\cdot b + r \\quad\\, \\text{where} \\quad\\, \\deg(r) < \\deg(b) . \\notag $$\n",
"\n",
"Caveat: For such a representation to exist the leading coefficient of $b$ must be invertible. Over $\\mathbb Z$ this means the leading coefficient must be $\\pm1$. Over $\\mathbb Q$ or $\\mathbb R$ any nonzero coefficient works."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Example"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Suppose $b=x^2+2x+3$ and $a=4x^4+5x^3+6x^2+7x+8$.\n",
"\n",
"By hand and some algebraic manipulation you can work out an expression $a=q\\cdot b+r$:\n",
"\n",
"\\begin{align*}\n",
"a &= 0\\cdot b + 4x^4 + 5x^3 + 6x^2 + 7x + 8 \\\\\n",
"&= (4x^2)\\cdot b - 3x^3 - 6x^2 + 7x + 8 \\\\\n",
"&= (4x^2-3x)\\cdot b + 16 x + 8\n",
"\\end{align*}\n",
"\n",
"So $q=4x^2-3x$ and $r=16 x + 8$."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Sage provides a `quo_rem` function to compute a quotient and remainder simultaneously."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"a = 4*x^4 + 5*x^3 + 6*x^2 + 7*x + 8\n",
"b = x^2 + 2*x + 3\n",
"q, r = a.quo_rem(b)\n",
"show([q,r])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Suppose $a$ has degree $n$, $b$ has degree $m$, and the leading coefficient of $b$ is 1 (for simplicity).\n",
"\n",
" r = a\n",
" for i from n-m, n-m-1, ..., 0 do\n",
" if deg(r) = m+i then\n",
" q[i] = leading coefficient of r\n",
" # Set the leading coefficient of r to zero\n",
" r = r - q[i]*x^i*b\n",
" else\n",
" q[i] = 0\n",
" q = add(q[i]*x^i for i in range(n-m+1))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Sage implementation for polynomial division"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"outputs": [],
"source": [
"# Return the quotient and remainder of a divided by b in R = Z[x]\n",
"def division(a, b):\n",
" r = a\n",
" n = a.degree(a)\n",
" m = b.degree(b)\n",
" q = [0 for i in [0..n-m]]\n",
" for i in range(n-m,-1,-1):\n",
" if r.degree() == m+i:\n",
" q[i] = r.leading_coefficient()\n",
" r = r - q[i]*x^i*b\n",
" return R([q[i] for i in [0..n-m]]), r\n",
"q, r = division(a, b)\n",
"show([q, r])\n",
"show(a==q*b+r)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Analysis in coefficient multiplications and additions"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"The line that updates $r$ requires $O(m)$ multiplications and additions each time it is run, which is at most $n-m+1$ times.\n",
"\n",
"These are the only coefficient additions and multiplications! In total the algorithm uses $O((n-m+1)\\cdot m)$ coefficient operations."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Integer case"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"The same basic algorithm for division works in the integer case but the details are somewhat more complicated.\n",
"\n",
"It uses $O((n-m+1)\\cdot m)$ word operations to divide an integer of length $n$ by an integer of length $m$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Take-aways"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"Arbitrary-precision integers and polynomials have many similarities.\n",
"\n",
"| | Polynomials | Integers |\n",
"|--------------------|:----------------------:|:---------------:|\n",
"| **size** | degree | length |\n",
"| **algorithm cost** | coefficient operations | word operations |\n",
"\n",
"If $a$ and $b$ are polynomials/integers of degree/length $n$ and $m$ then the cost of each algebraic operation is:\n",
"\n",
"* Addition/subtraction: $O(n + m)$\n",
"* Multiplication: $O(n\\cdot m)$\n",
"* Division with remainder: $O((n - m + 1)\\cdot m)$\n",
"\n",
"Generally, this means that addition/subtraction run in **linear** time in the size of the input while multiplication/division run in **quadratic** time in the size of the input."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"An important and interesting problem is now raised: **Can we do better?**\n",
"\n",
"In the case of addition/subtraction: the classical algorithms are already optimal.\n",
"\n",
"In the case of multiplication/division: **yes**, we can do better!\n",
"\n",
"Seeing how to solve these and related problems quickly is a driving motivation of this course.\n",
"\n",
"A huge amount of work has gone into solving these problems quicker, both theoretically and in practice."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "-"
}
},
"source": [
"For example, the GMP library is a C library for arbitrary precision integer arithmetic used by numerous applications such as public-key cryptography and compilers for various programming languages.\n",
"They have run extensive benchmarks of various algorithms for multiplication:\n",
"\n",
"![Multiplication benchmark results](The-Fastest-Arithmetic-in-the-World-new.png)\n",
"\n",
"The horizontal axis represents the length of the first operand and the vertical axis represents the length of the second operand (both on a log scale). The colours denote the algorithm that gives the best performance in each case."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "SageMath 9.5",
"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.4"
},
"latex_envs": {
"LaTeX_envs_menu_present": true,
"autoclose": false,
"autocomplete": false,
"bibliofile": "biblio.bib",
"cite_by": "apalike",
"current_citInitial": 1,
"eqLabelWithNumbers": true,
"eqNumInitial": 1,
"hotkeys": {
"equation": "Ctrl-E",
"itemize": "Ctrl-I"
},
"labels_anchors": true,
"latex_user_defs": false,
"report_style_numbering": false,
"user_envs_cfg": false
},
"latex_metadata": {
"author": "Curtis Bright",
"date": "September 14, 2022",
"title": "Computational Mathematics: Handout 02"
}
},
"nbformat": 4,
"nbformat_minor": 4
}