\documentclass[a4paper,12pt]{amsart}
\usepackage{kbexo}


%%%%%%%%%%%%%%%%%%%%%
%% Alice's packages
%%%%%%%%%%%%%%%%%%%%%
\usepackage{xcolor}

%%%%%%%%%%%%%%%%%%%
%% Alice's macros
%%%%%%%%%%%%%%%%%%%
\newcommand{\comment}[1]{{\color{gray}#1}}

\def\pow{\^{}\hskip0pt}

\course{Algèbre et calcul formel -- Agrégation}
\num{4}

\begin{document}
\exotitle[Multiplication de polynômes : Karatsuba, FFT]
\medskip

On propose ici d'écrire des procédures de multiplication de polynômes.
Chacune de ces procédures prendra comme entrée les deux polynômes sous forme
de listes et rendra le produit également sous forme de liste. Ces procédures
seront faciles à appliquer à des polynômes grâce aux transformations
liste-polynôme et polynôme-liste.


\exo[Gérer l'aléa]
Commençons par quelques commandes utiles pour gérer l'aléa dans Sage. 

\Q La plupart des ensembles de Sage ont une méthode {\tt .random\_element()}, qui permet de tirer des éléments aléatoires dans l'ensemble. Si l'ensemble est fini, l'élément est tiré uniformément dans l'ensemble. Sinon il est tiré selon une distribution choisie par Sage (vous pouvez obtenir plus d'info grâce à l'aide, en faisant  {\tt ZZ.random\_element?} par exemple).
Testez par exemple

\vspace{2mm}
\noindent{\tt ZZ.random\_element()}

\noindent{\tt ZZ.random\_element(a,b)} \comment{(l'élément tiré sera uniforme dans [a,b[)}

\noindent{\tt QQ['x'].random\_element()}

\noindent{\tt QQ['x'].random\_element(d)} \comment{(renvoie un polynôme aléatoire de degré $d$)}

\begin{verbatim}
Integers(7).random_element()
GF(32).random_element()
\end{verbatim}

\Q Lorsqu'on veut faire une démo et que certains bouts du code utilisent de l'aléa, on aime bien fixer l'aléa, pour que la démo fasse toujours la même chose. C'est faisable en sage en utilisant {\tt set\_random\_seed(42)} (vous pouvez remplacer $42$ par votre entier préféré). Attention, si vous évaluez plusieurs fois la même cellule, et que vous avez mis votre {\tt set\_random\_seed(42)} dans une cellule précédente, l'aléa va quand même changer. Par exemple

\begin{verbatim}
set_random_seed(42)
print(ZZ.random_element())
print(ZZ.random_element())
\end{verbatim}

Vous devriez obtenir $-31$, puis $5$. Normalement, n'importe qui utilisant Sage avec la même version que vous devrait obtenir les mêmes valeurs.

\exo[Karatsuba]
On désire calculer le produit de deux polynômes $P,Q\in R[X]$ de degrés
$<n$, où $R$ est un anneau commutatif et $n$ une puissance de $2$.
L'approche naïve a une complexité
algébrique de $O(n^2)$ opérations dans $R$. Une façon d'améliorer ce résultat
est la suivante.
Si $n=1$, $P$ et $Q$ sont des polynômes constants et on les
multiplie normalement. Sinon, on pose $m=n/2$.
On \'ecrit 
$$P=X^{m}P_1+P_0\quad \text{ and }\quad Q=X^{m}Q_1+Q_0,$$
o\`u $P_1$, $P_0$, $Q_1$ et $Q_0$ sont des polyn\^omes de
degr\'esstrictement inférieurs à $m$.  
On a alors

\begin{align*}
PQ & =  X^{n}P_1 Q_1+X^{m}(P_1 Q_0 + P_0 Q_1)+ P_0 Q_0\\
& =   X^{n}P_1 Q_1+X^{m}\left((P_1 + P_0)(Q_1 + Q_0) - P_1 Q_1 - P_0 Q_0\right) + P_0 Q_0,
\end{align*}

de telle sorte que nous avons juste à calculer trois produits
$$P_0Q_0,\quad P_1Q_1 \quad\text{et}\quad (P_0+P_1)(Q_0+Q_1)$$
de polynômes de degrés $<m$. On utilise cette idée de façon
récursive.

\Q Quelle est la complexité de cet algorithme ?

\Q Soit donc $n$ une puissance de $2$. Écrire une procédure récursive
\texttt{Karatsuba($P,Q,n$)} utilisant le principe rappelé ci-dessus et
renvoyant le polynôme $PQ$.

Dans cette procédure, $P$ et $Q$ seront rentrés comme des listes, et le
polynôme $PQ$ obtenu sera également une liste.

\Q Tester cette procédure sur les polynômes $a_0+a_1x$ et
$a_2+a_3x$ de $\mathbb Z[a_0,a_1,a_2,a_3][x]$.

\Q  Écrire ensuite une procédure {\tt MulK(A,P,Q)} qui, étant donnés deux
polynômes $P$ et $Q$ définis dans l'anneau de polynômes {\tt A} utilise
\texttt{Karatsuba} pour donner le produit $PQ$, lui aussi défini dans
{\tt A}.

\Q Tester cette procédure sur des polynômes symboliques de degré 3.

\Q La tester numériquement avec de gros polynômes pris au hasard.

\def\F{\mathcal{F}}
\exo\label{FFT} [FFT]
Soit $n = 2^k$ une puissance de $2$, avec $k>0$. 
Soit $\omega$ une racine primitive $n$-ième de l'unité,
par exemple $\omega=e^{2i\pi/n}$.
On définit un isomorphisme de $\CM$-algèbres $\F_\omega:\CM[X]/(X^n-1)\to \CM^n$
par
$$\F_{\omega}(R)=(R(1),R(\omega),\dots,R(\omega^{n-1})).$$
On identifiera une classe $\bar{R}\in\CM[X]/(X^n-1)$ avec son représentant
$R\in\CM[X]$ de degré $<n$, ainsi qu'avec le $n$-uplet $(R_0,\ldots,
R_{n-1})$ de ses coefficients, dans $\CM^n$.

On évalue $\F_{\omega}(R)$ en calculant récursivement deux $\F$ de degrés
$<m=n/2$ par le biais des formules
\begin{eqnarray*}
  R(\omega^p)&=&\sum_{j=0}^{m-1}R_{2j}\alpha^{jp}
  +\omega^p\sum_{j=0}^{m-1}R_{2j+1}\alpha^{jp}\\
  R(\omega^{p+m})&=&\sum_{j=0}^{m-1}R_{2j}\alpha^{jp}
  -\omega^p\sum_{j=0}^{m-1}R_{2j+1}\alpha^{jp},
\end{eqnarray*}
où $0\leq p<m$ et où $\alpha=\omega^2$.
\medskip

\Q Programmer l'algorithme récursif s'appuyant sur cette remarque. La
procédure \texttt{FFT} recevra en entrées $R$, $\omega$ et $n$, et retournera
$\F_{\omega}(R)$. Là encore, le polynôme $R$ sera défini par une liste. On
prendra garde à ne pas calculer $\omega^p$ à chaque étape de la boucle sur
$p$. 

\Q Quelle est la complexité de votre algorithme ? (En terme d'opérations dans le corps de base $K$).

\exo [Produit rapide de polynômes par FFT]
Ici encore $n=2^k$ avec $k>0$. Soient $P$ et $Q$ deux polynômes de $\CM[X]$
vérifiant deg$(PQ)<n$. On identifiera encore 
$P$ et $Q$ aux $n$-uplets $(P_0,\dots,P_{n-1})$ et $(Q_0,\dots,Q_{n-1})$
formés de leurs coefficients. On rappelle que l'on a alors
$$\F_{\omega}(PQ)=\F_{\omega}(P)\cdot \F_{\omega}(Q)\footnote{ici $(u_i)_{0\leq i<n }\cdot (v_i)_{0\leq i<n }=(u_iv_i)_{0\leq i<n }$.},$$
et que pour tout polynôme $R\in{\mathbb C}[X]$ de degré $<n$ on a
\medskip
$$\F_{\omega^{-1}}(\F_{\omega}(R))=\F_{\omega}(\F_{\omega^{-1}}(R))=nR.$$
Écrire une procédure prenant en arguments $P$, $Q$ et $n$ et retournant $PQ$,
procédure qui prendra bien sûr appui sur la procédure \texttt{FFT} de
l'exercice précédent.

\misctitle{Remarque} Pour s'assurer que  $\deg(PQ)<n$ on pourra imposer  à
$P$ et $Q$ d'être tous deux de degrés $<m=n/2$.


\exo[La fonction {\tt fft} de sage]
\Q Expérimenter les commandes suivantes.
\begin{verbatim}
   A=[RR(1) for i in range(8)]
   s=IndexedSequence(A,range(8))
   t=s.fft();t
   lt=t.list();lt
\end{verbatim}
  \Q Comparer les résultats donnés par cette fonction {\tt fft} avec
les résultats de l'exercice~\ref{FFT}.

\exo[FFT sur un corps fini]
Appliquer l'algorithme de l'exercice 2 à un corps $\mathbb F_p$, par exemple
$\mathbb F_{257}$ et $n=16$. Essayer aussi l'exercice 3 pour de tels
exemples.

\exo[Filtres]
Soit $f$ une fonction périodique de période $1$ de $\mathbb R$ dans $\mathbb
C$. Il suffit de connaître $f$ sur $[0,1]$ pour connaître $f$. Soient $N$ une
puissance de $2$ et $w=e^{2i\pi/N}$. Soit
$$
F=\left[f(0),f\left(\frac1N\right),\dots,f\left(\frac{N-1}N\right)\right].
$$
Pour tout $n$,
$$
F[n]=\frac1N\sum_{k=0}^{N-1}\mathcal F_{w^{-1}}(F)[k]w^{nk}.
$$
Plus simplement, si l'on note
$$
\hat F=\frac1N\mathcal F_{w^{-1}}(F),
$$
$$
F[n]=\sum_{k=0}^{N-1}\hat F[k]w^{nk}.
$$
Ainsi, la transformation $F\mapsto\hat F$ permet de passer de $F$ à la suite
des composantes fréquencielles de $F$.

\Q Étudions un exemple simple. Tous les calculs seront faits en flottants.
Soit
$$
f(t)=\sin(18\pi t)+3\cos(10\pi t).
$$
Soient $N=16$ et $w=e^{i\pi/8}$. On définit $F$ comme ci-dessus.
\qq En utilisant votre fonction pour la FFT, calculer $\hat F$.
\qq 
%Soient
%$$
%lF=[F(0),\dots,\F(15)]
%\hbox{\quad et\quad}
%lFc=[\hat F(0),\dots,\hat F(15)],
%$$
Représenter sur un graphique la ligne brisée définie par les points
$$\left[(0,F[0]),\dots,\left(\frac{15}{16},F[15]\right)\right]$$
(on pourra utiliser la fonction {\tt line}).

Représenter sur un autre graphique les points
$$\left[(0,|\hat F[0]|),\dots,\left(\frac{15}{16},|\hat F[15]|\right)\right]$$
en utilisant la fonction
{\tt point}. Commenter ce graphique.

\qq Calculer $\mathcal F_w(\hat F)$, et représenter sur un graphique la ligne
brisée définie par cette fonction. La comparer avec la ligne brisée
correspondant à $F$.

Pour mieux les comparer, on peut les placer sur le même graphique en
s'inspirant du modèle suivant, où $A$ et $B$ sont des listes, et où $x$ est
la liste des $\dfrac i{16}$.

\begin{verbatim}
  GrapheA = line([(x[i],A[i]) for i in range(16)],color='red')
  GrapheB = line([(x[i],B[i]) for i in range(16)])
  plot(GrapheA + GrapheB)
\end{verbatim}
Les deux premières lignes sont des affectations. La troisième trace les deux
lignes sur un même graphique. On pourrait aussi utiliser le fait que Sage affiche par défaut la dernière ligne de la cellule et ne pas mettre le {\tt plot} dans la dernière ligne. Par contre, si vous rajoutez une commandes après dans la cellule, votre graphique ne s'affichera plus...
Cela peut aussi s'écrire
\begin{verbatim}
  GrapheA = line(zip(x,A),color='red')
  GrapheB = line(zip(x,B))
  plot(GrapheA + GrapheB)
\end{verbatim}

\qq Introduisons un bruit dans le signal $F$.
\begin{verbatim}
  def h():
      return ZZ.random_element(-500,500)/1000
\end{verbatim}
définit une fonction aléatoire. Soit le bruit
\begin{verbatim}
  B=[h() for i in range(16)]
\end{verbatim}
Ajoutons ce bruit à $F$.
\begin{verbatim}
  FB=[F[i]+B[i] for i in range(16)]
\end{verbatim}
Faire un graphe dessinant la ligne définie par $FB$ et la comparer à celle de
$F$. On suppose qu'au lieu du signal $F$, on ait reçu le signal brouillé
$FB$. On va essayer de reconstituer $F$, en considérant que les fréquences
ajoutée par le bruit sont de faible amplitude.

\qq Calculer $\widehat{FB}$ et dessiner l'ensemble des points correspondants.
Comparer avec les points donnés par $\hat F$.

\qq\label{nettoyage} Soit $G$ définie par : $G[i]=\widehat{FB}[i]$ si
$\abs{\widehat{FB}[i]} > 1/4$ et $G[i]=0$ sinon. Calculer $\mathcal F_w(G)$ et
comparer avec $F$. Commenter.

\Q Refaire l'exercice avec $N=1024$ et $w=e^{2i\pi/1024}$.

\Q\label{polynome} Soit $f$ la fonction périodique de période $1$ définie
sur $[0,1[$ par
$$
f(t)=1920(t-1/6)(t^2-1/2)(t-1/2)(t-7/8)(t-1/10)(t-1)t.
$$
Refaire l'exercice précédent sur cette fonction avec $N=1024$ et
$w=e^{2i\pi/1024}$. Quand on nettoie le signal brouillé
(question~\ref{nettoyage}), on peut modifier le seuil à partir duquel on met
$\widehat{FB}[k]$ à $0$.

\Q On peut aussi choisir de filtrer les fréquences situées dans certains
intervalles. Reprendre la question~\ref{polynome} en mettant à $0$ les
$\widehat{FB}[k]$ tels que $k\in[[11,1013]]$ (cela revient à filtrer les hautes fréquences, responsables des oscillations serrées).

\Q Expérimenter sur la fonction définie sur $[0,1[$ par
$$
f(t)=1920(t-1/6)(t^2-1/2)(t-1/2)(t-7/8)(t-1/10)(t-1).
$$

\misctitle{Remarque} Ça ne marche pas bien sur cette fonction, avez-vous une explication ? \comment{(Indice: regardez les valeurs de $f(0)$ et $f(1)$)}
\end{document}
