<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
<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>
# Floating point error propagation in polynomial multiplication with Fast-Fourier Transform
- *Simple question*
Step1:  Examples
Let consider the polynomials $P(X) = 1 + X + X^3$ and $Q(X) = X^2 + X^5$, of degrees $n=3$ and $m=5$
Step2:  Their product is $(PQ)(X)$, of degree $n+m=8$
Step3:  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), it will be equal to $(PQ)(X)$.
The fit can be obtained, for instance, by Lagrange interpolation, which is not so efficient but easy to implement.
Here, I will simply use the numpy.polyfit function.
Step4:  Naive interpolation values
Let us consider the points $\lambda_i = i$, $i=0,\dots,n+m$.
Step5:  Ok, at least it seems to work!
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$.
The largest "should be 0" value (i.e., $\simeq 0$) value was
Step6:  Chebyshev nodes as interpolation values
Let us consider the points $\lambda_k = \cos\left(\frac{2k-1}{2d} \pi\right)$, $k=1,\dots,1+d=n+m+1$.
These are called the Chebyshev nodes.
Step7:  Ok, at least it seems to work!
But we saw that even with very small degrees ($n=3, m=5$), floating point errors were not so small on these points
Step8:  Benchmark
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
Step9:  Of course, our implementation is slower.
But on small polynomials, not so slower.
What about larger polynomials?
Step10:  On a larger example
<ASSISTANT_TASK:>
Python Code:
import numpy as np
np.version.full_version
Explanation: Table of Contents
<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>
# Floating point error propagation in polynomial multiplication with Fast-Fourier Transform
- *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?
In this [Jupyter notebook](https://www.jupyter.org/), [I](http://perso.crans.org/besson/) will try to understand this phenomena.
----
## Requirements
- 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.
End of explanation
P = [      1, 0, 1, 1]
n = len(P) - 1
P, n
Q = [1, 0, 0, 1, 0, 0]
m = len(Q) - 1
Q, m
Explanation: Examples
Let consider the polynomials $P(X) = 1 + X + X^3$ and $Q(X) = X^2 + X^5$, of degrees $n=3$ and $m=5$:
End of explanation
PQ = polymul(P, Q)
d = len(PQ) - 1
PQ, d
Explanation: Their product is $(PQ)(X)$, of degree $n+m=8$:
$$ \begin{align}(PQ)(X) &= (1 + X + X^3) (X^2 + X^5) \
&= X^2 + X^5 + X^3 + X^7 + X^5 + X^8 \
&= X^2 + X^3 + 2 X^5 + X^7 + X^8
\end{align} $$
End of explanation
assert d == n + m
Explanation: 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), it will be equal to $(PQ)(X)$.
The fit can be obtained, for instance, by Lagrange interpolation, which is not so efficient but easy to implement.
Here, I will simply use the numpy.polyfit function.
End of explanation
lambdas = np.arange(0, d + 1)
lambdas
values_P = np.polyval(P, lambdas)
values_P
values_Q = np.polyval(Q, lambdas)
values_Q
values_PQ = values_P * values_Q
values_PQ
PQ_sampled = np.polyfit(lambdas, values_PQ, d)
PQ_sampled
PQ
np.asarray(np.round(PQ_sampled), dtype=int)
Explanation: Naive interpolation values
Let us consider the points $\lambda_i = i$, $i=0,\dots,n+m$.
End of explanation
np.max(np.abs(PQ_sampled)[np.abs(PQ_sampled) < 0.9])
Explanation: Ok, at least it seems to work!
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$.
The largest "should be 0" value (i.e., $\simeq 0$) value was:
End of explanation
lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d))
lambdas
values_P = np.polyval(P, lambdas)
values_P
values_Q = np.polyval(Q, lambdas)
values_Q
values_PQ = values_P * values_Q
values_PQ
PQ_sampled2 = np.polyfit(lambdas, values_PQ, d)
PQ_sampled2
PQ
np.asarray(np.round(PQ_sampled2), dtype=int)
Explanation: Chebyshev nodes as interpolation values
Let us consider the points $\lambda_k = \cos\left(\frac{2k-1}{2d} \pi\right)$, $k=1,\dots,1+d=n+m+1$.
These are called the Chebyshev nodes.
End of explanation
np.max(np.abs(PQ_sampled2)[np.abs(PQ_sampled2) < 0.9])
Explanation: Ok, at least it seems to work!
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:
End of explanation
def mypolymul(P, Q):
    n = len(P) - 1
    m = len(Q) - 1
    d = n + m
    lambdas = np.cos(np.pi * (2 * np.arange(1, 2 + d) - 1) / (2 * d))
    values_P = np.polyval(P, lambdas)
    values_Q = np.polyval(Q, lambdas)
    values_PQ = values_P * values_Q
    PQ_sampled = np.polyfit(lambdas, values_PQ, d)
    # return PQ_sampled
    return np.asarray(np.round(PQ_sampled), dtype=int)
np.polymul(P, Q)
mypolymul(P, Q)
import warnings
warnings.simplefilter('ignore', np.RankWarning)
%timeit np.polymul(P, Q)
%timeit mypolymul(P, Q)
Explanation: Benchmark
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:
End of explanation
def random_polynomial(d=10, maxcoef=1):
    return np.random.randint(low=-maxcoef, high=maxcoef+1, size=d+1)
P = random_polynomial()
Q = random_polynomial()
P, Q
%timeit np.polymul(P, Q)
np.polymul(P, Q)
%timeit mypolymul(P, Q)
mypolymul(P, Q)
assert np.all(np.polymul(P, Q) == mypolymul(P, Q))
Explanation: Of course, our implementation is slower.
But on small polynomials, not so slower.
What about larger polynomials?
End of explanation
d = 100
maxcoef = 1
%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
%timeit mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef)
assert np.all(np.polymul(P, Q) == mypolymul(P, Q))
d = 10
maxcoef = 3
%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
%timeit mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef)
assert np.all(np.polymul(P, Q) == mypolymul(P, Q))
d = 10
maxcoef = 50
%timeit np.polymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
%timeit mypolymul(random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef))
P, Q = random_polynomial(d=d, maxcoef=maxcoef), random_polynomial(d=d, maxcoef=maxcoef)
assert np.all(np.polymul(P, Q) == mypolymul(P, Q))
Explanation: On a larger example:
End of explanation
<END_TASK>