{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"toc": "true"
},
"source": [
"# Table of Contents\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Floating point error propagation in polynomial multiplication with Fast-Fourier Transform\n",
"\n",
"- *Simple question*: when using the [FFT](https://en.wikipedia.org/wiki/Multiplication_algorithm#Fourier_transform_methods) (or any Fourier transform methods) for multiplying two polynomials, how to deal with floating point error propagation?\n",
"\n",
"In this [Jupyter notebook](https://www.jupyter.org/), [I](http://perso.crans.org/besson/) will try to understand this phenomena."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Requirements\n",
"- The [`numpy.polymul`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polymul.html#numpy.polymul) and [`numpy.polyfit`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html#numpy.polyfit) functions."
]
},
{
"cell_type": "code",
"execution_count": 96,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'1.12.1'"
]
},
"execution_count": 96,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"import numpy as np\n",
"np.version.full_version"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Examples\n",
"\n",
"Let consider the polynomials $P(X) = 1 + X + X^3$ and $Q(X) = X^2 + X^5$, of degrees $n=3$ and $m=5$:"
]
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"([1, 0, 1, 1], 3)"
]
},
"execution_count": 116,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = [ 1, 0, 1, 1]\n",
"n = len(P) - 1\n",
"P, n"
]
},
{
"cell_type": "code",
"execution_count": 117,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"([1, 0, 0, 1, 0, 0], 5)"
]
},
"execution_count": 117,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Q = [1, 0, 0, 1, 0, 0]\n",
"m = len(Q) - 1\n",
"Q, m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Their product is $(PQ)(X)$, of degree $n+m=8$:\n",
"$$ \\begin{align*}(PQ)(X) &= (1 + X + X^3) (X^2 + X^5) \\\\\n",
"&= X^2 + X^5 + X^3 + X^7 + X^5 + X^8 \\\\\n",
"&= X^2 + X^3 + 2 X^5 + X^7 + X^8\n",
"\\end{align*} $$"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"(array([1, 0, 1, 2, 0, 1, 1, 0, 0]), 8)"
]
},
"execution_count": 36,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PQ = polymul(P, Q)\n",
"d = len(PQ) - 1\n",
"PQ, d"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we evaluate both $P(X)$ and $Q(X)$, on $n+m$ different points, $\\lambda_i \\in \\mathbb{R}$ or $\\in\\mathbb{N}$, then we can *fit* a polynomial of degree $n+m = \\delta(PQ)$ on these sampling points, and by uniqueness (thanks to [the Fundamental Theorem of Algebra](https://en.wikipedia.org/wiki/Fundamental_theorem_of_algebra#Corollaries)), it will be equal to $(PQ)(X)$.\n",
"\n",
"The fit can be obtained, for instance, by Lagrange interpolation, which is not so efficient but easy to implement.\n",
"Here, I will simply use the [`numpy.polyfit`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html#numpy.polyfit) function."
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"assert d == n + m"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Naive interpolation values\n",
"\n",
"Let us consider the points $\\lambda_i = i$, $i=0,\\dots,n+m$."
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0, 1, 2, 3, 4, 5, 6, 7, 8])"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lambdas = np.arange(0, d + 1)\n",
"lambdas"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1, 3, 11, 31, 69, 131, 223, 351, 521])"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"values_P = np.polyval(P, lambdas)\n",
"values_P"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0, 2, 36, 252, 1040, 3150, 7812, 16856, 32832])"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"values_Q = np.polyval(Q, lambdas)\n",
"values_Q"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0, 6, 396, 7812, 71760, 412650,\n",
" 1742076, 5916456, 17105472])"
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"values_PQ = values_P * values_Q\n",
"values_PQ"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1.00000000e+00, -5.06843863e-12, 1.00000000e+00,\n",
" 2.00000000e+00, 2.43215915e-09, 9.99999993e-01,\n",
" 1.00000001e+00, -9.77014560e-09, 5.62310258e-09])"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PQ_sampled = np.polyfit(lambdas, values_PQ, d)\n",
"PQ_sampled"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 0, 1, 2, 0, 1, 1, 0, 0])"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PQ"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 0, 1, 2, 0, 1, 1, 0, 0])"
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.asarray(np.round(PQ_sampled), dtype=int)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ok, at least it seems to work!\n",
"\n",
"But we saw that even with very small degrees ($n=3, m=5$), floating point errors were not so small on these wrongly chosen points $\\lambda_i = i$, $i=0,\\dots,n+m$.\n",
"The largest \"should be 0\" value (i.e., $\\simeq 0$) value was:"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"9.7701456003317433e-09"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.max(np.abs(PQ_sampled)[np.abs(PQ_sampled) < 0.9])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Chebyshev nodes as interpolation values\n",
"\n",
"Let us consider the points $\\lambda_k = \\cos\\left(\\frac{2k-1}{2d} \\pi\\right)$, $k=1,\\dots,1+d=n+m+1$.\n",
"These are called the [Chebyshev nodes](https://en.wikipedia.org/wiki/Chebyshev_nodes)."
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.98078528, 0.83146961, 0.55557023, 0.19509032, -0.19509032,\n",
" -0.55557023, -0.83146961, -0.98078528, -0.98078528])"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d))\n",
"lambdas"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 2.92424164, 2.40629924, 1.72705159, 1.20251551, 0.79748449,\n",
" 0.27294841, -0.40629924, -0.92424164, -0.92424164])"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"values_P = np.polyval(P, lambdas)\n",
"values_P"
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 1.86948796, 1.08874542, 0.36158742, 0.03834284, 0.03777763,\n",
" 0.25572914, 0.29393801, 0.05439157, 0.05439157])"
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"values_Q = np.polyval(Q, lambdas)\n",
"values_Q"
]
},
{
"cell_type": "code",
"execution_count": 88,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 5.46683454, 2.61984727, 0.62448014, 0.04610786, 0.03012707,\n",
" 0.06980086, -0.11942679, -0.05027096, -0.05027096])"
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"values_PQ = values_P * values_Q\n",
"values_PQ"
]
},
{
"cell_type": "code",
"execution_count": 91,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:1: RankWarning: Polyfit may be poorly conditioned\n",
" \"\"\"Entry point for launching an IPython kernel.\n"
]
},
{
"data": {
"text/plain": [
"array([ 1.21880049e+00, 1.55720018e-14, 5.62399019e-01,\n",
" 2.00000000e+00, 2.73500613e-01, 1.00000000e+00,\n",
" 9.45299877e-01, -2.11823147e-16, 1.70937883e-03])"
]
},
"execution_count": 91,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PQ_sampled2 = np.polyfit(lambdas, values_PQ, d)\n",
"PQ_sampled2"
]
},
{
"cell_type": "code",
"execution_count": 92,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 0, 1, 2, 0, 1, 1, 0, 0])"
]
},
"execution_count": 92,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PQ"
]
},
{
"cell_type": "code",
"execution_count": 93,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 0, 1, 2, 0, 1, 1, 0, 0])"
]
},
"execution_count": 93,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.asarray(np.round(PQ_sampled2), dtype=int)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ok, at least it seems to work!\n",
"\n",
"But we saw that even with very small degrees ($n=3, m=5$), floating point errors were not so small on these points: the largest \"should be 0\" value (i.e., $\\simeq 0$) value was:"
]
},
{
"cell_type": "code",
"execution_count": 94,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"text/plain": [
"0.56239901918986446"
]
},
"execution_count": 94,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.max(np.abs(PQ_sampled2)[np.abs(PQ_sampled2) < 0.9])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Benchmark\n",
"\n",
"Stupidly, let us check if our naive implementation of $(P, Q) \\mapsto PQ$ by evaluation-and-interpolation is more or less efficient than `numpy.polyfit`:"
]
},
{
"cell_type": "code",
"execution_count": 123,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def mypolymul(P, Q):\n",
" n = len(P) - 1\n",
" m = len(Q) - 1\n",
" d = n + m\n",
" lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d))\n",
" values_P = np.polyval(P, lambdas)\n",
" values_Q = np.polyval(Q, lambdas)\n",
" values_PQ = values_P * values_Q\n",
" PQ_sampled = np.polyfit(lambdas, values_PQ, d)\n",
" # return PQ_sampled\n",
" return np.asarray(np.round(PQ_sampled), dtype=int)"
]
},
{
"cell_type": "code",
"execution_count": 124,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 0, 1, 2, 0, 1, 1, 0, 0])"
]
},
"execution_count": 124,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.polymul(P, Q)"
]
},
{
"cell_type": "code",
"execution_count": 125,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([1, 0, 1, 2, 0, 1, 1, 0, 0])"
]
},
"execution_count": 125,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"mypolymul(P, Q)"
]
},
{
"cell_type": "code",
"execution_count": 121,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"33 µs ± 995 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
"181 µs ± 7.68 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
]
}
],
"source": [
"import warnings\n",
"warnings.simplefilter('ignore', np.RankWarning)\n",
"\n",
"%timeit np.polymul(P, Q)\n",
"%timeit mypolymul(P, Q)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Of course, our implementation is slower.\n",
"But on small polynomials, not so slower.\n",
"\n",
"What about larger polynomials?"
]
},
{
"cell_type": "code",
"execution_count": 126,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def random_polynomial(d=10, maxcoef=1):\n",
" return np.random.randint(low=-maxcoef, high=maxcoef+1, size=d+1)"
]
},
{
"cell_type": "code",
"execution_count": 134,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([ 1, 1, -1, -1, 0, 1, 1, 0, 1, 1, 0]),\n",
" array([-1, 0, -1, 0, 0, -1, 1, 0, -1, 1, 1]))"
]
},
"execution_count": 134,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"27.4 µs ± 1.67 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
]
},
{
"data": {
"text/plain": [
"array([-1, -1, 0, 0, 1, -1, -1, 1, -3, -2, 1, 0, -1, -3, 0, 3, 0,\n",
" 0, 2, 1, 0])"
]
},
"execution_count": 134,
"metadata": {},
"output_type": "execute_result"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"350 µs ± 7.41 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
]
},
{
"data": {
"text/plain": [
"array([-1, -1, 0, 0, 1, -1, -1, 1, -3, -2, 1, 0, -1, -3, 0, 3, 0,\n",
" 0, 2, 1, 0])"
]
},
"execution_count": 134,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"P = random_polynomial()\n",
"Q = random_polynomial()\n",
"P, Q\n",
"%timeit np.polymul(P, Q)\n",
"np.polymul(P, Q)\n",
"%timeit mypolymul(P, Q)\n",
"mypolymul(P, Q)\n",
"assert np.all(np.polymul(P, Q) == mypolymul(P, Q))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"On a larger example:"
]
},
{
"cell_type": "code",
"execution_count": 138,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"48.6 µs ± 1.36 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
"11.9 ms ± 1.7 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/lib/python3.5/dist-packages/ipykernel_launcher.py:7: DeprecationWarning: elementwise == comparison failed; this will raise an error in the future.\n",
" import sys\n"
]
},
{
"ename": "AssertionError",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrandom_polynomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxcoef\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxcoef\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom_polynomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxcoef\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxcoef\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpolymul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mmypolymul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m: "
]
}
],
"source": [
"d = 100\n",
"maxcoef = 1\n",
"%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))\n",
"%timeit mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))\n",
"\n",
"P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef)\n",
"assert np.all(np.polymul(P, Q) == mypolymul(P, Q))"
]
},
{
"cell_type": "code",
"execution_count": 139,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"38.9 µs ± 2.57 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
"369 µs ± 15.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
]
},
{
"ename": "AssertionError",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrandom_polynomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxcoef\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxcoef\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom_polynomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxcoef\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxcoef\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpolymul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mmypolymul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m: "
]
}
],
"source": [
"d = 10\n",
"maxcoef = 3\n",
"%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))\n",
"%timeit mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))\n",
"\n",
"P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef)\n",
"assert np.all(np.polymul(P, Q) == mypolymul(P, Q))"
]
},
{
"cell_type": "code",
"execution_count": 140,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"36 µs ± 1.21 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
"452 µs ± 62.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
]
},
{
"ename": "AssertionError",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrandom_polynomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxcoef\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxcoef\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom_polynomial\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0md\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmaxcoef\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmaxcoef\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mall\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpolymul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0mmypolymul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mQ\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m: "
]
}
],
"source": [
"d = 10\n",
"maxcoef = 50\n",
"%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))\n",
"%timeit mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))\n",
"\n",
"P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef)\n",
"assert np.all(np.polymul(P, Q) == mypolymul(P, Q))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our method is slower.\n",
"And wrong.\n",
"\n",
"That's sad."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"## Conclusion\n",
"\n",
"Implementing naively the multiplication of 1D polynomials with evaluation-and-interpolation does not work.\n",
"\n",
"- It is slower that the FFT based method (available in any numerical computation environment), e.g., `numpy.polymul` in Python with NumPy.\n",
"- And it does not work. Booum!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"----\n",
"> Thanks for reading!\n",
"\n",
"> See [this repo on GitHub](https://github.com/Naereen/notebooks/) for more notebooks, or [on nbviewer.jupyter.org](https://nbviewer.jupyter.org/github/Naereen/notebooks/).\n",
"\n",
"> That's it for this demo! See you, folks!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.5.3"
},
"notify_time": "5",
"toc": {
"colors": {
"hover_highlight": "#DAA520",
"running_highlight": "#FF0000",
"selected_highlight": "#FFD700"
},
"moveMenuLeft": true,
"nav_menu": {
"height": "172px",
"width": "251px"
},
"navigate_menu": true,
"number_sections": true,
"sideBar": true,
"threshold": 4,
"toc_cell": true,
"toc_position": {
"height": "480px",
"left": "0px",
"right": "1068px",
"top": "142px",
"width": "212px"
},
"toc_section_display": "block",
"toc_window_display": true
}
},
"nbformat": 4,
"nbformat_minor": 2
}