\documentclass{article}
\usepackage[latin1]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage[francais]{babel}
\usepackage{xspace}
\newcommand{\bs}{\symbol{92}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\N}{\mathbb{N}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}
\newtheorem{theorem}{Theorem}
\newtheorem{defn}{Definition}

\begin{document}

\title{Approximation des fonctions}
\author{Préparation agrégation, option C}
\date{05/06}

\maketitle

Ce texte présente quelques méthodes d'approximation de fonctions
qui servent en particulier à calculer les fonctions classiques
en utilisant des fonctions plus simples (addition, multiplication,
division, évaluation de polynômes).

\section{La méthode de Newton}
Pour calculer $F(a)$, on résoud une équation du type $f(x)=0$ dont
$x=F(a)$ est solution, et telle que le calcul de $g$ soit possible
(ou plus efficace). Par exemple pour calculer la racine $n$-ième
d'un réel positif $a$, on prend $f(x)=x^n-a$. 

On utilise souvent des fonctions dont la convexité est connue,
pour pouvoir appliquer le 
\begin{theorem} \index{Newton}
Soit $f$ de classe $C^2$ telle que
 $f(r)=0, f'(r)>0$ et $f'{'} \geq 0$ sur $[r,b]$. Alors
pour tout $u_0 \in [r,b]$ la suite de la méthode de Newton
\[ u_{n+1} = u_n -\frac{f(u_n)}{f'(u_n)},  \]
est définie, décroissante, minorée par $r$ et converge vers 
$r$. 
\end{theorem}
Le théorème des accroissements finis permet alors de calculer la
précision de la racine approchée trouvée en fonction de l'image
de cette racine~:
\[ 0 \leq u_n -r \leq \frac{f(u_n)}{f'(r)} \]
Typiquement le nombre de chiffres significatifs corrects est
multipli\'e par 2 \`a chaque it\'eration.
C'est par exemple cette m\'ethode qui permet de calculer 
la racine carr\'ee de $a>1$~:
\[ u_{n+1} = \frac{1}{2}(u_n+\frac{a}{u_n})\]

Le calcul de la valeur initiale $u_0$ de la suite a une grande
importance pour le nombre d'itérations à effectuer, on utilise
diverses méthodes, par exemple une des méthodes présentées ci-dessous,
la méthode de Newton pouvant servir à améliorer la précision (en particulier
dans le cadre des flottants multiprécision). 


\section{Les polynomes}
\subsection{Les séries entières.}
A l'intérieur du rayon de convergence d'une fonction analytique,
on peut utiliser le développement en séries entières de la fonction.
On calcule une majoration à priori du reste sur la zone considérée,
puis on cherche le rang à partir duquel ce reste est plus petit que
la précision souhaitée. Il y a trois grandes méthodes de majoration~:
le reste de Taylor (lorsque la dérivée $n$-ième de la fonction
est facile à calculer et majorer), les séries alternées ou
l'utilisation de la majoration suivante~:
s'il existe un point $x_0$ tel que 
$|a_n x_0^n|$ est born\'e par $M$, alors pour $|x|<|x_0|$ on 
peut majorer le reste de la s\'erie au rang $n$ par 
\[ |R_n| \leq M \frac{ |\frac{x}{x_0}|^n} {1-|\frac{x}{x_0}|} 
\]
Le calcul de $f$ à la précision souhaitée
se ramène ainsi à une évaluation de polynome.

Certaines séries se pretent bien à cette méthode, par exemples
les fonctions transcendantes directes (exp, sin, cos), d'autres
moins bien (atan, ln) car la convergence est un peu plus lente
(pour calculer $n$ bits de mantisse,
il faut $O(n)$ termes du d\'eveloppement au lieu de $O(n/\ln(n))$
termes).
Les syst\`emes de calcul en flottant multipr\'ecision utilisent
les s\'eries enti\`eres (apr\'es r\'eduction d'argument) pour 
l'exponentielle ou le calcul de $\ln(2)$. 
Richard Brent a montr\'e en 1976 qu'on pouvait \'evaluer les
fonctions usuelles en $O(M(n)\ln(n))$ op\'erations, o\`u $M(n)$
d\'esigne la complexit\'e de la multiplication des entiers de taille
$n$ (par exemple $M(n)=O(n \ln(n) \ln(\ln(n)))$ pour 
la m\'ethode de Sch\"onhage-Strassen). Cela fait intervenir 

\subsection{Interpolation polynomiale (Lagrange, ...)}
Cette fois, on approche $f$ par un polynôme qui est égal à $f$
en certains points (on peut aussi imposer que les dérivées de $f$
et $P$ soient égales en certains points). C'est la méthode d'interpolation
de Lagrange, pour laquelle on a une majoration de la différence
entre $f$ et $P$~:
Soit $f$ une fonction $n+1$ fois d\'erivable sur un intervalle $I=[a,b]$
de $\R$, $x_0,...,x_n$ des r\'eels distincts de $I$. 
Soit $P$ le polynome de Lagrange donn\'e par les $x_j$ et $y_j=f(x_j)$.
Pour tout r\'eel $x \in I$,
il existe un r\'eel $\xi_x \in [a,b]$ (qui d\'epend de $x$) tel
que~:
\[ f(x)-P(x) = \frac{f^{[n+1]}(\xi_x)}{(n+1)!} \prod_{j=0}^n(x-x_j) \]
Le choix des points $x_i$ (à $n$ fixé) rendant la majoration
ci-dessus optimale dépend de l'utilisation que l'on veut faire de $P$.
On peut prendre des points équirépartis sur l'intervalle, puis
s'en servir pour calculer une valeur approchée d'intégrale, ce sont
les méthodes de Newton-Cotes. Ce choix n'est pas optimal pour
l'ordre de la méthode de quadrature, si on laisse le choix des $x_i$
libres, on montre qu'il existe un choix rendant 
la méthode d'intégration exacte pour tous les polynômes de degré plus
petit que $2n+1$ (méthodes de Gauss).

On peut aussi souhaiter qu'une norme issue d'un produit
scalaire (par exemple la norme $L^2$ 
avec un poids) de $f-P$ soit minimale, ce qui revient à projeter
$f$ sur un sous-espace engendré par des polynômes de degré plus
petit qu'un degré fixé, il s'agit
alors d'en calculer une base (de polynomes orthogonaux). Le même
principe permet de construire les séries de Fourier en prenant
des polynômes trigonométriques.

\section{Les approximants de Padé}
Pour approcher $f$ on utilise une fraction rationnelle (dont le degré du
numérateur et du dénominateur sont bornés) à la place de polynomes.
Le calcul du numérateur et du dénominateur fait intervenir
l'algorithme de Bézout. 

Plus précisément, on cherche une fraction ayant même développement
de Taylor que $f$ en un point donné, 0 pour fixer les idées, à un
ordre fixé. On doit donc résoudre
\begin{equation} \label{eq:pade}
 P=\frac{N}{D} +O(x^{p+1})
\end{equation}
où $P$ est donné de degré inférieur ou égal à $p$, 
$N$ et $D$ sont les inconnues
de degré inférieur ou égal à $n$ et $d$, avec $n+d=p$ et $D$ premier
avec $x$. $D$ admet donc un inverse modulo $x^p$, si on multiplie
par $D$, on obtient~:
\[ D P = N \pmod{x^{p+1}}\]
On effectue l'algorithme de Bézout pour $A=x^{p+1}$ et $B=P$~:
\[ A U_k + B V_k=R_k, \quad U_0=1=V_1, U_1=0=V_0, R_0=x^{p+1}, R_1=P \]
et on s'arrête au premier rang $k$ tel que 
le degré de $R_k$ soit inférieur ou égal à $n$. On montre
ensuite que~:
\[ \mbox{deg}(V_{k})+\mbox{deg}(R_{k-1}) = ... =  
\mbox{deg}(V_1)+\mbox{deg}(R_0 ) = \mbox{deg}(A)=p+1
\]
ce qui entraine que le degré de $V_k$ est inférieur ou égal à $p+1-(n+1)=d$.
On obtient ainsi une solution de $D P = N \pmod{x^{p+1}}$ ($D=V_k$ et
$N=R_k$). Si $V_k$ est premier avec $x$, on conclut. On peut montrer
que si $V_k$ et $R_k$ ne sont pas premiers entre eux, alors
(\ref{eq:pade}) n'admet pas de solutions avec ces conditions de degr\'e.

Par exemple pour $e^x$ en prenant $p=6$ et $n=3$, on trouve
$(x^{3}+12\* x^{2}+60\* x+120)/(-x^{3}+12\* x^{2}-60\* x+120)$.

On peut montrer qu'il existe un $\varepsilon>0$ tel que l'erreur
absolue entre $f$ et son approximant de Pad\'e sur 
$[-\varepsilon,\varepsilon]$ est atteint aux bornes de l'intervalle,
et on peut calculer explicitement une valeur de $\varepsilon$.

Cette méthode est tr\`es utile en pratique, par exemple pour le calcul de
la fonction exponentielle en pr\'ecision machine, car elle 
nécessite d'une part moins
d'opérations (à cause de la symétrie numérateur/dénominateur)
et fournit de plus une approximation de meilleure qualité que
le développement de Taylor de $e^x$ (à l'ordre 6 ici).

On peut aussi approcher une fonction par une fraction rationnelle
interpolant la fonction en certains points (comme pour les polyn\^omes
de Lagrange mais avec une fraction au lieu d'un polyn\^ome). Le
probl\`eme \`a r\'esoudre est analogue mais modulo un produit de
$x-x_i$ au lieu de $x^{t+1}$. On parle d'approximant de Cauchy.

\section{Méthodes mixtes}
Pour le calcul des fonctions usuelles, on utilise la plupart du temps
une combinaison des méthodes ci-dessus et des propriétés de ces fonctions.

Par exemple, pour calculer le logarithme d'un flottant 
en base 2, on calcule une fois pour toutes $ln(2)$ à la 
précision souhaitée, puis on utilise la représentation
$a=2^e(1+m)$, donc $\ln(a)=e\ln(2)+\ln(1+m)$ avec $m \in[0,1[$.
On peut alors utiliser le développement en séries de $1+m$ mais
si $m$ est proche de 1, la convergence sera très lente (controlée
par la nature de série alternée). On peut alors décider de calculer
une approximation en prenant le développement à l'ordre disons 9,
puis on effectue quelques itérations de Newton de l'équation
$e^x-a=0$ en prenant comme valeur initiale cette valeur (par excès)
de $\ln(a)$. Cette méthode pourrait en particulier être utilisée
pour calculer le logarithme multi-précision à partir d'une
valeur approchée par excès en simple précision (en r\'ealit\'e,
on utilise la moyenne arithm\'etico-g\'eom\'etrique qui est utilis\'ee,
elle fait intervenir un encadrement d'une int\'egrale elliptique
par $ln(x)$ en utilisant les valeurs de $\ln(2)$ et $\pi$).

On peut aussi
utiliser des méthodes ad hoc pour calculer la valeur
initiale, par exemple pour la fonction
arctangeante sur un microprocesseur calculant en base 10, on
stocke les valeurs de $\arctan(10^{-k})$ pour les premières
valeurs de $k$. Soit $a$ l'angle dont on veut calculer
l'arctangente, on se ramène d'abord à $a$ dans l'intervalle $[0,1]$.
Soit $\alpha=\arctan(a)$. 
On pose alors $x=a$ et $y=1$. Soit $\beta=\arctan(10^{-1})$.
On calcule 
\[ \tan(\alpha-\beta)=\frac{x/y-1/10}{1+1/10 \* x/y}
= \frac{x-y/10}{y+x/10}\]
Si ce nombre est positif, on pose $x'=x-y/10$ et $y'=y+x/10$,
et on se ramène au calcul de $\arctan(x'/y')$. Si ce nombre est négatif,
cela signifie que le développement décimal de $\arctan(x/y)$ a comme
première décimale 0. On calcule ainsi la première décimale
du développement de $\arctan(a)$ en effectuant uniquement
des additions, des soustractions et des décalages (une division
par 10 étant un décalage). On continue alors de la même
manière pour les décimales suivantes qui sont dans la table
des arctangentes. On termine en effectuant la méthode de Newton.

\section{Suggestions de développement}
\begin{itemize}
%\item Les méthodes de quadrature (Newton-Cotes, Gauss, ...).
\item Complexit\'e des op\'erations de base sur les flottants 
multipr\'ecision, cas de la racine carr\'ee, de l'exponentielle
en utilisant le d\'eveloppement en s\'eries.
\item Approximation polynomiale par intervalles avec recollement
$C^1$, $C^2$, etc.
\item Calcul des fonctions usuelles (discussion sur les méthodes,
leur efficacité...)
\item Interpolation de Lagrange, calcul effectif (différences divisées),
applications (par exemple pour des algorithmes efficaces: polynome
caractéristique, PGCD à plusieurs variables,...)
%\item Polynomes orthogonaux.
\item Calcul d'approximants de Padé ou de Cauchy.
\item Comparaison entre différentes méthodes d'approximation pour
une même fonction (par exemple le logarithme nepérien, la fonction
sinus), dans l'optique du calcul en simple et multi-précision.
\end{itemize}

\end{document}

Pad\'e, suite
On avait $s(x)$ avec deg$(s)<=2t-1$,
il s'agissait de voir comment la solution $u,v,r$ calculee par Bezout
\begin{equation}   \label{eq:loc}
 u(x)  x^{2t}+ v(x) s(x) = r(x) 
\end{equation}
avec arr\^et pr\'ematur\'e au 1er reste $r(x)$ de degr\'e $<=t-1$ 
correspondait \`a l'\'equation
\[   l(x) s(x) = w(x) mod x^{2t} \]
avec deg$(l)<=t$ et deg$(w)<=t-1$

On a vu que deg$(v)<=t$.
On commence par factoriser la puissance de $x$ de degr\'e maximal $p$ dans
$v(x)$, et on simplifie (\ref{eq:loc}) par $x^p$. 
Quitte \`a changer $v$ et $r$, on se
ramene donc \`a:
\[   u(x)  x^{2t-p}+ v(x) s(x) = r(x) \]
avec $v(x)$ premier avec $x$, deg$(v)<= t-p$ et deg$(r)<= t-1-p$.
On simplifie ensuite par le pgcd de $v(x)$ et de $r(x)$
(qui divise $u(x)$ car premier avec $x$ puisqu'on a d\'ej\`a trait\'e les
puissances de $x$).
On a donc, quitte \`a changer de notation $u,v,r$ tels que
\[  u(x)  x^{2t-p}+ v(x) s(x) = r(x) \]
avec $v$ et $r$ premiers entre eux, $v$ premier avec $x$,
deg$(v)<=t-p$ et deg$(r)<=t-1-p$ (N.B.: $p=0$ en general)

On observe que $l(x)$ est premier avec $x$ ($0$ n'est pas racine de $l$).
On raisonne modulo $x^{2t-p}$, $l$ et $v$ sont inversibles modulo $x^{2t-p}$,
donc 
\[ s(x) = w(x)*inv(l) \pmod{ x^{2t-p}},
\quad s(x) = r(x)*inv(v) \pmod {x^{2t-p}} \]
donc 
\[ w(x)*inv(l)=r(x)*inv(v) \pmod{x^{2t-p}} 
\quad \Rightarrow \quad w(x)*v(x)=r(x)*l(x) \pmod {x^{2t-p}} \]
donc $w(x)*v(x)=r(x)*l(x)$ \`a cause des majorations de degres

D'autre part par construction $w(x)$ est premier avec $l(x)$ (car chacun
des facteurs de $l(x)$ divise tous les \'el\'ements de la somme d\'efinissant
$w(x)$ sauf un), donc $l(x)$ divise $v(x)$, et comme $v(x)$ est premier
avec $r(x)$, on en d\'eduit que $v(x)=C l(x)$ o\`u $C$ est une constante non
nulle, puis $r(x) = C w(x)$.

Bezout donne donc (apr\`es simplifications du couple $v(x), r(x)$ par
son pgcd) le polynome localisateur \`a une constante pr\`es (donc les
racines et les positions des erreurs), et on peut calculer
les valeurs des erreurs avec $v$ et $r$ car la constante $C$ se simplifie.