Floating_point_error_propagation_in_polynomial_multiplication_with_Fast-Fourier_Transform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": "true"
   },
   "source": [
    "# Table of Contents\n",
    " <p><div class=\"lev1 toc-item\"><a href=\"#Floating-point-error-propagation-in-polynomial-multiplication-with-Fast-Fourier-Transform\" data-toc-modified-id=\"Floating-point-error-propagation-in-polynomial-multiplication-with-Fast-Fourier-Transform-1\"><span class=\"toc-item-num\">1&nbsp;&nbsp;</span>Floating point error propagation in polynomial multiplication with Fast-Fourier Transform</a></div><div class=\"lev2 toc-item\"><a href=\"#Requirements\" data-toc-modified-id=\"Requirements-11\"><span class=\"toc-item-num\">1.1&nbsp;&nbsp;</span>Requirements</a></div><div class=\"lev2 toc-item\"><a href=\"#Examples\" data-toc-modified-id=\"Examples-12\"><span class=\"toc-item-num\">1.2&nbsp;&nbsp;</span>Examples</a></div><div class=\"lev3 toc-item\"><a href=\"#Naive-interpolation-values\" data-toc-modified-id=\"Naive-interpolation-values-121\"><span class=\"toc-item-num\">1.2.1&nbsp;&nbsp;</span>Naive interpolation values</a></div><div class=\"lev3 toc-item\"><a href=\"#Chebyshev-nodes-as-interpolation-values\" data-toc-modified-id=\"Chebyshev-nodes-as-interpolation-values-122\"><span class=\"toc-item-num\">1.2.2&nbsp;&nbsp;</span>Chebyshev nodes as interpolation values</a></div><div class=\"lev3 toc-item\"><a href=\"#Benchmark\" data-toc-modified-id=\"Benchmark-123\"><span class=\"toc-item-num\">1.2.3&nbsp;&nbsp;</span>Benchmark</a></div><div class=\"lev2 toc-item\"><a href=\"#Conclusion\" data-toc-modified-id=\"Conclusion-13\"><span class=\"toc-item-num\">1.3&nbsp;&nbsp;</span>Conclusion</a></div>"
   ]
  },
  {
   "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<ipython-input-138-c9d7506c5aaa>\u001b[0m in \u001b[0;36m<module>\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<ipython-input-139-49289ccd0af4>\u001b[0m in \u001b[0;36m<module>\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<ipython-input-140-f5f8ae306925>\u001b[0m in \u001b[0;36m<module>\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
}