%% Last modified: <Analyse_ieee.tex modifié par  Guillaume CONNAN le mardi 28 avril 2015 à 12h 45min 35s>


\chapterimage{girl_robot}
\chaptertext{Des avions qui s'écrasent, des  fusées qui s'écrasent, des missiles
  qui tuent les  amis, les indices boursiers totalement faux,  les résultats des
  élections totalement faux, les jeux qui plantent,...

La  méconnaissance  du  traitement  des  nombres sur  machine  (et  non  pas  le
traitement lui-même !) peut être tragique ou au mieux catastrophique.

Nous allons initier  notre exploration de ce monde mystérieux  que vous utilisez
chaque jour en plongeant dans la machine.}

\chapter{IEEE 754-2008}

\section{Les nombres à virgule flottante}


\subsection{Les problèmes}


\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
\begin{figure}
  \includegraphics[width=3cm]{Kahan}
\caption{W.\bsc{Kahan (1933 -)}}
\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5.2cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}

Wiliam \textsc{Kahan} est un mathématicien et informaticien canadien, lauréat du
prix \textsc{Turing}  en 1989 et principal  artisan de la norme  IEEE 754 publiée
pour la première fois en 1985.

Les problèmes  liés à  la manipulation  des VF (nombres  à virgule  flottante, ou
\textit{floating point (FP)} in english of speciality), sont nombreux.

Citons notamment le célèbre cas du missile américain \textit{Patriot} qui manqua
sa cible, un  missile \textit{Scud} irakien, durant la première  guerre du Golfe
en février 1991. Le \textit{Patriot} poursuit  sa cible en mesurant le temps mis
par les ondes radars pour revenir après rebond sur la cible. Or le système mesurait le temps en dizaines de seconde (drôle
d'idée) puis multipliait ce  temps par 0,1 pour obtenir le  temps en secondes et
utilisait un registre à 24 bits en virgule fixe.

Après cent  heures d'utilisation  (l'armée US  ne chôme pas),  un écart  de 0,34
seconde est apparu (pourquoi? cf en TD...). Un \textit{Scud} volant à 1.676 m/s,
le  \textit{Patriot} a  visé 500  m à  côté. Cette  erreur causa  la mort  de 28
soldats tandis que 98 étaient blessés.



Il y a eu beaucoup d'autres cas célèbres, même chez les plus puissants:

%\begin{figure}
\begin{center}
  \includegraphics[width=\linewidth]{google_float.jpg}
\end{center}
%\end{figure}

Bon, ne soyons pas de mauvaise foi:

\begin{haskellcode}
*Main> 5 - 4.9 - 0.1
-3.608224830031759e-16
\end{haskellcode}

Même des logiciels spécialisés, très chers (Maple \copyright), se trompent:

\maple

\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{mupad}
> evalf(sin(2^100));
                             .4491999480
> evalf(sin(2^100),20);
                        -.58645356896925826300
> evalf(sin(2^100),30);
                   .199885621653625738215132811525
> evalf(sin(2^100),40);
              -.8721836054182673097807197782134705593243
\end{minted}
            

Heureusement, le libre sauve la mise...Giac/XCAS:


\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{mupad}
1>> evalf(sin(2^100))
-0.872183605418
\end{minted}



Haskell

\begin{haskellcode}
*Main> sin(2^100)
-0.8721836054182673
\end{haskellcode}


Un exemple rigolo: les élections  au Schleswig-Holstein en 1992 (merci à
Paul \textsc{Zimmermann}).

Les partis ayant moins de 5\% des voix n’ont pas de siège.

Au dépouillement, le parti écologiste obtient exactement 5\%.

Après l’annonce officielle des résultats, on découvre que le
parti écologiste n’a en fait que 4.97\% des voix!

Le programme affichant les résultats arrondit à un seul chiffre
après la virgule, et donc a arrondi 4.97 à 5.0.

L'erreur corrigée, le siège est retiré aux verts, et attribué au
SPD, qui avec 45 sièges sur 89 obtient la majorité absolue.


Un autre exemple qui fait pitié (merci à William \textsc{Kahan}):



%\begin{figure}
\begin{center}
  \includegraphics[width=0.6\linewidth]{excel}
\end{center}
%\end{figure}


Un dernier qui vous touchera plus:



%\begin{figure}
\begin{center}
  \includegraphics[width=\linewidth]{jeu_erreur}
\end{center}
%\end{figure}


\medskip

Le problème des VF est apparu sur machine en 1941: Konrad \textsc{Zuse} a proposé son Z3 avec une
mantisse de 14 bits, un exposant de 7 bits, un bit de signe le tout en base 2.


En effet, les humains préfèrent compter  en base 10...et n'arrêtent pas de faire
des erreurs d'arrondis suite à leur manie (pourquoi? cf en TD...).
Certaines  machines calculent  donc  en  base 10  (les  calculettes, Maple\copyright,  les
machines dédiées à  la Finance,...). La norme IEEE 754-2008  introduit aussi des
calculs en base 10.



La machine pour l'instant est basée sur une logique à deux états et préfère donc
la base 2 (ou 8 ou 16).


\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
\begin{figure}
  \includegraphics[width=3cm]{brusentsov}
\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-6.5cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}


Certaines machines  ont pourtant  été conçues  selon une  logique à  trois états
(Vrai, Faux, j'en sais rien) comme le \selectlanguage{russian}
\foreignlanguage{russian}{Сетунь}
\selectlanguage{francais}
 soviétique mis au point en
1958  par  une  équipe dirigée  par  \foreignlanguage{russian}{Николай  Петрович
  Брусенцов} 
et  le   célèbre  mathématicien  \foreignlanguage{russian}{Сергей Львович Соболев}. 
L'histoire    de   cette    machine   est    édifiante:   la    production   des
\foreignlanguage{russian}{Сетунь}  fut confiée  à  l'\textit{usine des  machines
  mathématiques}      qui     en      fabriqua     très      peu     car      le
\foreignlanguage{russian}{Сетунь} était trop bon marché!  De plus il a fonctionné des
années durant sans aucune panne ni besoin de pièces de rechange ce qui n'est pas
très    bien    vu    dans    le    milieu    informatique    :-)    Enfin    le
\foreignlanguage{russian}{Сетунь}    réglait    le     problème    crucial    du
\textit{tiers-exclus}   et    obéissait   aux   canons   les    plus   exigeants
d'E.  \textsc{Dijkstra}.  Bref,  les  autorités soviétiques  ont  décidé  de  le
découper  en morceaux  et  de mettre  les  restes dans  une  décharge. De  zélés
informaticiens ont toutefois réussi à en  sauvé un exemplaire en le cachant dans
un grenier.



 Vous utilisez vous-même une logique ternaire en SQL pour gérer
\texttt{NULL}.
Les machines quantiques peuvent travailler avec des \textit{qtrits}. Vous verrez
ça un jour, jeunes padawans.



Autre problème: le profusion d'appellations qui sèment la confusion! En C (cf dinechin):

\begin{itemize}
\item \texttt{char}  est un entier codé  sur 8 bits. C'est  l'abréviation du nom
  commun \textit{character};
\item \texttt{int} est l'abréviation du  non commun \textit{integer} qui désigne
  un    nombre     entier    comme     en    mathématiques     mais    cependant
  $2147483647+1=-2147483648$...
\item \texttt{short} et  \texttt{long} sont des adjectifs  désignant des entiers
  de taille différentes;
\item \texttt{float} est un verbe;
\item \texttt{single}? Simple quoi?
\item \texttt{double}? Mais le double de quoi?
\item \texttt{long} \texttt{double}: deux adjectifs à suivre ?
\end{itemize}

Alors il faut clarifier tout ça!



\subsection{Quelle précision?}

En  physique,  on ne  peut  guère  aller en  deçà  de  grandeurs de  l'ordre  de
$10^{-15}$ alors pourquoi s'embêter?

\begin{haskellcode}
*Main> 3 * 0.1
0.30000000000000004
*Main> sum $ take 10000000 $ repeat 0.1
999999.9998389754
\end{haskellcode}


\textbf{Il faut  parfois aller bien  plus loin que} $\mathbf{10^{-15}}$  \textbf{pour avoir
  une précision de} $\mathbf{10^{-15}}$\textbf{: c'est bien là le problème.}

On cherchera avec attention l'exercice\vref{ex::der}

Pourtant, avec 50 bits  on peut coder la distance Terre-Lune  avec une erreur de
l'ordre de la taille d'une bactérie!


Qu'est-ce  qui se  passe? On  ne peut  plus faire  de calcul?  On ne  pourra pas
fabriquer de jeux vidéos pour dégommer plein de méchants?

Pas  de panique:  il faut  juste  avoir en  tête que  \textbf{LES NOMBRES  RÉELS
  N'EXISTENT PAS: TOUT CE QUE VOUS AVEZ VU AU LYCÉE N'EST QU'ILLUSION!}

Oubliez donc  les réels: les  nombres que  vous allez manipuler  sont différents
mais sont  tout aussi  maîtrisables. Il suffit  de connaître  leur arithmétique,
\textbf{l'arithmétique des nombres à virgule flottante} qui obéissent à la norme IEEE 754.

Des problèmes subsistent néanmoins. Par exemple, on sait très bien ce que donnera un calcul
selon un format  donné...le problème est qu'on ne sait  pas toujours quel format
est utilisé...



\section{La norme IEEE 754}

\subsection{En bref}

\textbf{Nous nous contenterons, sauf mention contraire, de travailler en base 2.
}
Vous en  avez parlé  en archi: un  VF est  représenté par un  bit de  signe, une
mantisse (ou significande en franglais) et un exposant:

$$v = (-1)^s\times m\times 2^E$$

Les formats  habituels sont  le \textit{binary32}  ou \textit{single}  avec $(\#
E,\# m)=(8,24)$ et le\textit{ binary64} ou \textit{double} avec  $(\#
E,\# m)=(11,53)$.

Pour simplifier  nos présentations,  on illustrera souvent  nos propos  avec des
représentations \textit{toy7} plus petites, comme par exemple  $(\#
E,\# m)=(3,4)$.

La mantisse doit  dans la mesure du  possible vérifier $1 \leqslant  m<2$ ce qui
minimise   l'exposant  et   apporte   plus  d'informations.   Par  exemple,   en
\textit{toy7}:

$$
\begin{aligned}
  v&=0,010{\text{\textcolor{red}{\small \textit{01}}}}&\times &2^0\\
  v&=0,100{\text{\textcolor{red}{\small \textit{1}}}}&\times &2^{-1}\\
  v&=1,001&\times &2^{-2}\\
\end{aligned}
$$

C'est la \textit{forme normale} d'un VF. De plus, cela se permet de se passer du
premier 1 qui est implicite.  On notera 

$$m=1,f$$


On peut également gagner de la place  en ne stockant pas le signe de l'exposant:
il suffit de le translater de $2^{(\# E) - 1}-1$...POURQUOI?

$$e=E+2^{(\# E) - 1}-1$$



\begin{exemple}
  Représentons $0,75_{10}$ en \textit{toy7}.

$$0,75_{10} = \frac{3_{10}}{2_{10}^{2}}=11_2\times 2^{-2}=1,100_2\times2^{-1}$$

L'exposant stocké $e$ vérifie avec le décalage: 

$$e - \left(2^{3-1}-1\right)=E=-1 \lequiv e=-1+3=2=10_2
$$

Ainsi:

\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    0&0&1&0&1&0&0\\
    \hline
  \end{tabular}
\end{center}

\end{exemple}


On doit  également s'occuper des exposants  extrêmes qui ne correspondent  pas à
des formes normales (cf TD).

\begin{NB}[Petit algo pour convertir un nombre de valeur absolue inférieure à 1 en base 2]
  
Soit $x=a_1\times 2^{-1}+a_2\times 2^{-2}+a_3\times 2^{-3}+\cdots+a_n\times 2^n$
l'écriture de $x$ tel que $|x|<1$ en base 2 avec les $a_i$ égaux à 0 ou 1.


Alors $2x  = a_1 + a_2\times  2^{-1}+a_3\times 2^{-2}+\cdots+a_n\times 2^{n-1}$:
$a_1$ est donc la partie entière de $2x$.

Ensuite  $2(2x-a_1)=a_2+a_3\times  2^{-1}+\cdots+a_n\times 2^{n-2}$:  $a_2$  est
alors la partie entière de $2(2x-a_1)$ und so weiter...

\end{NB}


\subsection{Les normaux, les sous-normaux et les paranormaux}

\setindent{-0.7cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3.7cm]{super_science22}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5.5cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}


Nous venons de voir que l'exposant E, qui est un entier signé, est stocké sur $\#
E$ bits sous la forme $e=E+2^{(\# E) - 1}-1$.

 Comme  $e$ est  un  entier naturel  stocké  sur  $\# E$  bits,  il varie  entre
 $e_{\min}=0$ et $e_{\max}=2^{\# E}-1$.

Par   exemple,   en   \textit{toy7},   $\#E  =   3$   donc   $e_{\min}=000$   et
$e_{\max}=111=1000-1$.

CEPENDANT les valeurs extrêmes de $E$ sont réservées pour des cas spéciaux.

 Cela  donne:

$$
\begin{aligned}
E_{\min}&=\left(e_{\min}-2^{(\# E) - 1}+1\right)+1=0-2^{3-1}+1+1=-2 \\
E_{\max}&=\left(e_{\max}-2^{(\# E) - 1}+1\right)-1=(2^3-1)-2^{3-1}+1-1=7-4+1-1=3
\end{aligned}
$$



En effet, des problèmes apparaissent avec les formes normales:
\begin{itemize}
\item on ne peut pas représenter zéro sous forme normale;
\item des nombres du  type $0.\underbrace{000...00}_{\# E}11\times 2^0=1.1\times
  2^{-(\#E +1)}$ ne sont
  pas représentables sous  forme normale car il faudrait un  exposant plus petit
  que $E_{\min}$.
\end{itemize}

La  norme IEEE  754  a  donc introduit  les  nombres \textit{sous-normaux}  pour
remplir  le vide  qui existerait  entre 0  et le  plus petit  nombre normal  (cf
exercice \vref{ex::scale}).

Pour  les signaler  et  les distinguer  des  normaux, la  norme  impose que  leur
exposant décalé $e$ soit 0: $e_{\min}$ est donc une valeur réservée.

Un cas  particulier: si  $e=e_{min}$ et  $f=0$, alors la  norme ne  considère ce
nombre  ni comme  un nombre  normal  ni comme  un nombre  sous-normal. C'est  la
représentation de zéro.


\begin{NB}
  Il y a deux zéros!  En \textit{toy7}:

\begin{multicols}{2}
$+0$:
\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    0&0&0&0&0&0&0\\
    \hline
  \end{tabular}
\end{center}

$-0$:
\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    1&0&0&0&0&0&0\\
    \hline
  \end{tabular}
\end{center}
\end{multicols}

\medskip

Des tests \texttt{x == 0} pourraient alors  échouer! La norme impose donc que les
deux zéros soient considérés comme égaux:

\begin{multicols}{2}
\begin{haskellcode}
*Main> let x = 1 / 1e500
*Main> x
0.0
*Main> x == 0
True
*Main> let y = -1 / 1e500
\end{haskellcode}

\begin{haskellcode}
*Main> y
-0.0
*Main> y == 0
True
*Main> x == y
True
\end{haskellcode}
\end{multicols}


\end{NB}


\setindent{-0.7cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=4.5cm]{buzz}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}


Reprenons  à  présent un  exemple  classique:  vous  voulez calculer  une  norme
  $N=  \sqrt{x^2+y^2}  $  avec  $x=1\times 2^3$  et  $y=1,1\times  2^3$  en
\textit{toy7}.

Les carrés de ces nombres posent problème: leur exposant est trop grand. On est dans un cas de
\textit{dépassement} ou d'\textit{overflow} en english de spécialité.

Une  première idée  serait de  remplacer $x^2$  et $y^2$  par le  plus grand  nombre
disponible $1,111\times 2^4$ soit:


\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    0&1&1&1&1&1&1\\
    \hline
  \end{tabular}
\end{center}


Alors $x^2$ puis $y^2$ puis $x^2+y^2$ seraient remplacés par $1,111\times 2^4$.

On aurait  donc $N=  \left(1,111\times 2^4\right)^{1/2}$ arrondi  à $1,010\times
2^2$ ce qui est dramatiquement  faux et pourraient envoyer des \textit{Patriots}
sur le toit du bâtiment B.

La norme introduit donc deux infinis codés en \textit{toy7}:

\begin{multicols}{2}
$+\infty$:
\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    0&1&1&1&0&0&0\\
    \hline
  \end{tabular}
\end{center}
  

$-\infty$:
\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    1&1&1&1&0&0&0\\
    \hline
  \end{tabular}
\end{center}
\end{multicols}


Soit, plus généralement, $e=e_{\max}$ et $f=0$.


On retrouve alors les limites habituelles vues au lycée (ah, quand même...):

\begin{multicols}{2}
\begin{haskellcode}
*Main> 1 + 1e500
Infinity
*Main> let x = 1e500
*Main> x
Infinity
*Main> 1 + x
Infinity
\end{haskellcode}

\begin{haskellcode}

*Main> x^3
Infinity
*Main> 1 / x
0.0
*Main> 4 - x
-Infinity
\end{haskellcode}
\end{multicols}
Nous verrons bien d'autres exemples en exercice.



\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3cm]{nan}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-5.5cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}


Vous avez aussi parlé au lycée des \textbf{FORMES INDÉTERMINÉES} lors de l'étude
des limites. Par exemple, que dire de:

$$
\displaystyle\lim_{x\to +\infty}x^2 -x
$$

Vous êtes muets? Demandons à la machine:


\begin{haskellcode}
*Main> x^2 - x
NaN
\end{haskellcode}



Une manière exotique de régler le problème est d'étudier l'algorithme suivant:

\begin{center}
  \href{http://www.food.com/recipe/madhur-jaffreys-naan-bread-446809}{http://www.food.com/recipe/madhur-jaffreys-naan-bread-446809}
\end{center}


Sinon, il  suffit de savoir  que la  Hitroizeu introduit le  \texttt{NaN}
(comme dans  \textit{Not a  Number}) dans  les cas  où le  calcul demandé  est \og
indéfini\fg{}.


\begin{multicols}{2}
  \begin{haskellcode}
*Main> let x = 1e500
*Main> x - x
NaN
*Main> x / x
NaN
*Main> sqrt(-1)
NaN
*Main> 0 / 0
NaN

  \end{haskellcode}
  


\begin{haskellcode}
-- Bizarre ?
*Main> x - x
NaN
*Main> x - x == x - x
False
*Main> x^2 - x
NaN
*Main> x * (x - 1)
Infinity
\end{haskellcode}
\end{multicols}


Un  \texttt{NaN}   n'est  pas  comparable.   Trier  un  tableau   comportant  un
\texttt{NaN} ne pose donc pas de problème:

\begin{haskellcode}
*Main> let x = 0/0
*Main> x
NaN
*Main> x > 3
False
*Main> x < 3
False
*Main> x == 3
False
*Main> x == x
False
\end{haskellcode}



Il faut  penser à  un problème  célèbre dont  fut cette  fois victime  la marine
étasunienne.

\begin{aparte}
L'affaire est beaucoup moins dramatique mais plus croustillante que celle des \textit{Patriots}. La Marine US
s'équipe de PC montés  avec...Windows NT 4.0 et en profite  pour réduire le 10\%
le nombre de marins  à bord du USS Yorktown. Mais, comme le  dira un des experts
civils après la panne:
\og
{\itshape 
Using Windows NT, which is known to have some failure modes, on a warship is
similar to hoping that luck will be in our favor
}
\fg{}

Que s'est-il passé le 21 septembre 1997?  Un opérateur a tapé par erreur un zéro sur son clavier et
tout le navire a été bloqué en pleine mer :-)

\og{\itshape
  Your \$2.95 calculator, for example, gives you a zero when you
try to divide a number by zero, and does not stop executing the next set of instructions.
It seems that the computers on the Yorktown were not designed to tolerate such a simple
failure.
}\fg{}

On espère que les missiles atomiques sont mieux programmés...

\end{aparte}



Revenons  sur cette  anecdote et  transplantons  là dans  un environnement  plus
pacifique.


Vous cherchez les zéros de, disons, comme ça, au hasard, $ax^2+bx+c$. Comme vous
avez  de  vagues souvenirs  de  votre  jeunesse, vous  programmez:  retourne-moi
$-b\pm\frac{\sqrt{Δ}}{2a}$ avec $Δ=b^2-4ac$.

On suppose que toutes les valeurs sont proches de zéro et les mesures difficiles
à ce niveau de précision. 

\begin{multicols}{2}
Un premier essai donne:
\begin{haskellcode}
*Main> let a = 1e-6
*Main> let b = 1e-4
*Main> let c = 1e-10
*Main> let d = b^2 - 4*a*c
*Main> d
9.9999996e-9
*Main> -b - sqrt(d)/(2*a)
-50.000099
*Main> sqrt(d)
9.999999799999999e-5
\end{haskellcode}


Un deuxième essai, très proche, donne:
\begin{haskellcode}
Main> let a = 1e-6
*Main> let b = 1e-10
*Main> let c = 1e-4
*Main> let d = b^2 - 4*a*c
*Main> d
-3.9999999999e-10
*Main> -b - sqrt(d)/(2*a)
NaN
*Main> sqrt(d)
NaN
\end{haskellcode}
\end{multicols}

Ouf! Le \textit{Yorktown} ne coulera pas!

Ainsi, on  n'a plus  besoin de  surcharger le  programme de  conditionnelles. Le
calcul continue.


Reprenons nos formes indéterminées:

\begin{haskellcode}
*Main> let x = 1e500
*Main> x^2 - x
NaN
\end{haskellcode}


En effet,  $\infty -  \infty$ est  une forme  indéterminée. Que  faisiez-vous au
lycée pour lever cette indétermination?


\begin{haskellcode}
*Main> x*(x - 1)
Infinity
\end{haskellcode}

C'est réglé...
Autre curiosité:

\begin{haskellcode}
*Main> x - x
NaN
*Main> x - x == x - x
False
\end{haskellcode}

On peut débusquer un \texttt{NaN} en effectuant le test \texttt{x /= x}.


Selon la  Hitroizeu, un \texttt{NaN}  est représenté avec l'exposant  maximum et
une mantisse non nulle: il y a donc toute une famille de \texttt{NaN}!
En voici un exemple en \textit{toy7}:

\begin{center}
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    $s$ (1  bit) &  \multicolumn{3}{c||}{$e$ (3  bits)} &  \multicolumn{3}{c|}{$f$ (3
      bits)}\\
    \hline
    1&1&1&1&0&1&0\\
    \hline
  \end{tabular}
\end{center}

\subsection{Aparté:  de l'importance  pour  un programmeur  de  bien penser  son
  message d'erreur}


\setindent{-0.85cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=4cm]{airbus}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-3.7cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}


Encore un drame...Le  1\er juin 2009, le vol Air  France 447 finit tragiquement:
les 228 personnes à bord meurent englouties au large des côtes du Brésil.

Beaucoup de bruit  a été fait autour  de cet accident, beaucoup de  gens ont été
blâmés, en particulier les copilotes.

Voici le  problème: pour économiser  du carburant, les avions  volent maintenant
très haut et en  limite de décrochage. Il faut un  guidage informatique très fin
pour le  permettre. Ce soir-là,  l'avion est pris dans  une tempête et  le froid
extrême obstrue les capteurs qui ne peuvent plus mesurer la vitesse du vent.

Le  pilote  automatique envoie  alors  comme  message  à ses  collègues  humains
des messages trop vagues:

\textit{Les messages de panne successivement affichés sur l’ECAM n’ont pas permis à
l’équipage de faire un diagnostic rapide et efficace de la situation dans laquelle
l’avion se trouvait, en particulier de l’obstruction des sondes Pitot. Il n’a jamais été
en mesure de faire le lien entre les messages qui sont apparus et la procédure à
appliquer, alors que la lecture de l’ECAM et des messages doit faciliter l’analyse de
la situation et permettre d’organiser le traitement des pannes. Plusieurs systèmes
avaient  pourtant identifié  l’origine du  problème  mais n’ont  généré que  des
messages} (extrait du  rapport final du Bureau d'Enquêtes et  d'Analyses pour la
sécurité de l'aviation civile).

Donner les conséquences et pas les causes! Sans
indications, les pilotes perdent de  précieuses secondes et ne peuvent effectuer
la man\oe uvre habituelle permettant de contrer un décrochage.

\subsection{Successeur d'un VF}\label{sec::succ}

Il n'est  pas possible de définir  un successeur d'un  réel mais un VF,  lui, en
possède un.

Écrivons un flottant avec sa mantisse entière en posant $$n=\# f$$


\begin{remarque}
  Dans la littérature, on considère souvent $p=n+1$.
\end{remarque}


Par exemple, pour \textit{toy7}, $n=3$, pour \textit{binary32}, $n=23$, pour \textit{binary64}, $n=52$.

$$v=M(v)\times 2^{E(v)-n}$$

$M(v)$ étant un entier, il admet un  successeur. Si l'on reste dans les limites des
VF normaux du système:

$$\operatorname{succ} (v) = (M(v)+1)\times 2^{E(v)-n} = v +  2^{E(v)-n}$$

On note

$${\rm ulp}(v)=2^{E(v)-n}$$


\begin{NB}[ULP]

\end{NB}



Par exemple, en \textit{binary64}, soit $x$ un réel et $y=RN(x)$.

On a $y=(-1)^s\times 1,b_1b_2...b_{52}\times 2^E$.

Alors ${\rm ulp}(y)=0,000...001\times 2^E=2^{E-52}$


\subsection{Reconnaissance des flottants}

La machine traite une chaîne de bits \begin{tabular}[h]{|c||c||c|}
    \hline
    $s$&$e$&$f$\\
    \hline
  \end{tabular} ainsi:

\begin{center}

  \begin{tabular}[h!]{llcl}

\hline
$e=0...0$ & $f=0$ & $\rightarrow$ & $v=(-1)^s\times 0.0$ \\
$e=0...0$ & $f\neq  0$ & $\rightarrow$ & $v=(-1)^s\times  0.f \times 2^{1-E_{\max}}$\\
$e=1...1$ & $f=  0$ & $\rightarrow$ & $v=(-1)^s\times \infty$\\
$e=1...1$ & $f\neq 0$ & $\rightarrow$ & $\text{NaN}$\\
$0...0<e<1...1$ & $f\neq  0$ & $\rightarrow$ & $v=(-1)^s\times  1.f \times 2^{e-E_{\max}}$\\
\hline

  \end{tabular}

\end{center}

\subsection{Tableau récapitulatif}
On notera $ε_m=2^{-n}$ l'\textit{epsilon} de  la machine, c'est-à-dire le successeur de
1.

On notera $λ$ le plus petit VF normal positif.

On notera $μ$ le plus petit VF sous-normal positif.

On notera $Ω$ le plus grand VF normal.

\begin{recherche}
  Quellle relation existe-t-il entre $μ$, $ε_m$ et $λ$?
\end{recherche}


\setindent{-10mm}


  \begin{tabular}{|*{9}{c|}}
\hline
    Format & $\# E$ & $\# f$ & $E_{\min}$ & $E_{\max}$ & $ε_m$ & $λ$ & $μ$ & $Ω$ \\
\hline
\textit{toy7} & 3 & 3 & $-2$& $+3$ & $2^{-3}=1/8$ & $2^{-2}=1/4$ & $2^{-5}=1/32$ & $1,111\times 2^3=15$ \\
\textit{binary32} & 8 & 23 & $-126$ &$+127$ &$2^{-23}\approx 1,2\times 10^{-7}$ &\textcolor{white}{$2^{-23}\approx 1,2\times 10^{-7}$} & \textcolor{white}{$2^{-23}\approx 1,2\times 10^{-7}$}& \textcolor{white}{$2^{-23}\approx 1,2\times 10^{-7}$}\\
\textit{binary64} & 11 & 52 & $-1022$& & & & & \\
\hline
  \end{tabular}\label{tab::reca}

\setindent{\indentationcours}

\section{Algèbre des nombres VF}


\subsection{L'ensemble des VF}

Nous allons  à présent travailler  dans notre  ensemble de flottants  normaux et
sous-normaux. Nous noterons par exemple $\bbv_{32}$ les normaux et sous-normaux de
\textit{binary32} et nous noterons $\overline{\bbv_{32}}$ ce même
ensemble enrichi des deux infinis.

Dans le cas général, nous noterons tout simplement $\bbv_b$ voire $\bbv$.

\begin{NB}
Ainsi $\overline{\bbv_b}$ est un ensemble FINI.
\end{NB}

On  va  munir cet  ensemble  d'opérations  usuelles  et étudier  les  structures
algébriques produites.

\subsection{Comparaison}
\begin{recherche}
Il  est  très  simple  et  rapide  de  comparer  deux  VF:  comment  la  machine
procède-t-elle? Quel est l'avantage de ce stockage des VF?
\end{recherche}

\subsection{Addition}

On décompose l'action en trois étapes:

\begin{enumerate}
\item on commence par ramener les deux nombres au même exposant, en l'occurrence
  le plus grand des deux;
\item on ajoute les deux mantisses \textit{complètes} en tenant compte du signe;
\item on renormalise le nombre obtenu.
\end{enumerate}



Par exemple, en \textit{toy7}, additionnons 1,1 et 0,0111 ou plutôt $1,1\times
2^{0}$  et  $1,11\times  2^{-2}$  ou  plutôt  (en  \textit{toy7},  le  décalage
d'exposant est de 3):

\medskip

 1,1:
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    0&0&1&1&1&0&0\\
    \hline
  \end{tabular}
  \hfill 0,0011:
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
   \hline
    0&0&0&1&1&1&0\\
    \hline
  \end{tabular}

\medskip

On réécrit 0,00101 avec  le même exposant que 1,1 mais  attention à la précision
de la mantisse!

1,100

0,001\textcolor{blue}{\itshape 11}

$----$

1,101\textcolor{blue}{\itshape 11}


On ne peut pas garder les deux derniers 1 du plus petit nombre car on n'a pas la place de
les stocker!

Que faire? Il faut arrondir en choisissant selon quatre méthodes usuelles:


\begin{definition}[Les arrondis]
\begin{enumerate}
\item l'arrondi au plus proche (RN) qui  arrondit au VF...le plus proche. En cas
  d'égalité,  on choisit  la valeur  paire  (donc qui  se  termine par  un 0  en
  binaire);
\item l'arrondi  vers 0  (RZ) qui  arrondit à  la valeur  de plus  petite valeur
  absolue: c'est la troncature;
\item l'arrondi vers $+\infty$ (RU) qui  arrondit à la valeur supérieure la plus
  petite;
\item l'arrondi vers $-\infty$ (RD) qui arrondit à la valeur inférieure la plus grande.
\end{enumerate}
Le mode d'arrondi par défaut est le premier.
\end{definition}



Posons $x=1,10111$  (qui n'est  pas un VF!):  il est compris  entre les  deux VF
$1,101$ et $1,110$.

$$
1,101\textcolor{blue}{ \mathit{00}} \leqslant \underbrace{1,101\textcolor{blue}{ \mathit{11}}}_x \leqslant 1,110\textcolor{blue}{ \mathit{00}}
$$

$$
\underbrace{1,100+ \textcolor{blue}{\mathit{1,00}}ulp(x)}_v \leqslant  \underbrace{1,100+ \textcolor{blue}{ \mathit{1,11}}ulp(x)}_{x} \leqslant \underbrace{1,100+ \textcolor{blue}{\mathit{10,00}}ulp(x)}_{\operatorname{succ}(v)}
$$

$|x-v|=0,11ulp(x) $,  $\ |x- \operatorname{succ}(v)|=0,01ulp(x)$  donc $RN(x)=\operatorname{succ}(v)$




%\begin{figure}
\begin{center}
  \includegraphics{DrapeauFloat.6}
\end{center}
%\end{figure}



$$RN(1,1 + 0,001) = 1,110$$

Il  y a  énormément de  choses à  dire sur  ces problèmes  d'arrondis. Nous  leur
consacrerons une section.

Il y a beaucoup de choses à dire sur les sommes itérées. Nous leur
consacrerons une section.



%asso


\begin{recherche}
  Est-ce que l'addition des VF est associative? 
  
  Est-on sûr de l'ordre dans le quel un compilateur calcule $a+b+c+d$?
\end{recherche}


\subsection{Multiplication}


On suit un ordre logique:

\begin{enumerate}
\item on \og xore\fg{} les bits de signe;
\item on additionne les exposants réels et on décale ou plutôt on additionne les
  exposants décalés et on retire la valeur d'un décalage;
\item on multiplie les mantisses;
\item on normalise.
\end{enumerate}

Effectuons par exemple le produit de 101,1 par -10,01, i.e. :

 101,1:
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    0&1&0&1&0&1&1\\
    \hline
  \end{tabular}
  \hfill -10,01:
  \begin{tabular}[h]{|c||c|c|c||c|c|c|}
   \hline
    1&1&0&0&0&0&1\\
    \hline
  \end{tabular}

  \begin{enumerate}
  \item 0 xor 1 donne 1: le bit de signe est 1;
  \item $101 +  100 - 11 = 1001 -  11 = 110 $: l'exposant décalé  est $110$ donc
    l'exposant est 3;
  \item $1,011\times 1,001=1,011 + 1,011\times 0,001=1,011 + 0,001011=1,100011$;
\item le produit est donc $1,100011\times 2^3$. Le passage à la forme normale va
  arrondir le produit. On a 
$$1,100<x<1,101$$

avec   $|x-1,100|= 0,011ulp(x)$   et  $|x-1,101|= 0,101ulp(x)$   donc   $$RN(101,1\times
(-10,01))=-1,100\times 2^3$$

i.e. \begin{tabular}[h]{|c||c|c|c||c|c|c|}
    \hline
    1&1&1&0&1&0&0\\
    \hline
  \end{tabular}
  \end{enumerate}

Ici encore, le VF ne correspond pas au réel...

Il est temps de prendre le problème à bras le corps, mais avant:


\begin{recherche}
  Est-ce que la multiplication des VF est associative? 

Est-ce que la multiplication des VF est distributive sur l'addition?
\end{recherche}

% distrib
\section{Réels, arrondis et flottants}

\textit{Dans la suite de nos aventures, nous ne considérerons que le mode d'arrondi par défaut
sauf mention contraire.
}



\subsection{Un problème, dix problèmes, mille problèmes}

Nous commençons à bien connaître le problème:

\begin{haskellcode}
*Main> 3*0.1 == 0.3
False
\end{haskellcode}

On pourrait alors se dire qu'il ne faut jamais comparer deux VF $x$ et $y$ avec \texttt{==}
mais plutôt  les considérer égaux  s'ils vérifient  $|x-y|<d$ avec $d$  un seuil
d'erreur. Mais résoudre ce problème ainsi va en créer de nombreux autres:

\begin{itemize}
\item comment fixer $d$?
\item si $x$ et $y$ ont des signes différents?
\item $x\equiv y \lequiv |x-y|<d$ définit-elle une relation d'équivalence?
\item ...
\end{itemize}

\begin{recherche}
  Répondez à la troisième question...
\end{recherche}

Il faut penser différemment!

\begin{DANGER}
Beaucoup de propriétés  des réels ne sont plus vraies  dans $\bbv$. Par exemple,
le programme  suivant, que nous  allons étudier,  n'a aucun intérêt  dans $\bbr$
mais est primordial dans $\bbv$:

\begin{haskellcode}
f1 a =
  if (a + 1.0) - a /= 1.0 then a  
  else f1 (2.0 * a)
       
f2 a b =
  if (a + b) - a == b then b
  else f2 a (b + 1.0)
       
b = f2 (f1 1) 1
\end{haskellcode}
\end{DANGER}

Mais comme souvent  en mathématique et en informatique,  \textbf{un petit dessin
  aide à visualiser une situation complexe}.



\begin{DANGER}
  Attention: tout VF est  en fait un rationnel mais tout  rationnel n'est pas un
  VF! Par exemple...
\end{DANGER}

La seule certitude est qu'un VF se représente correctement...lui-même.

Mais  je vous  rappelle  que ce  que  vous  tapez sur  votre  clavier n'est  pas
forcément un VF si vous tapez en base 10.

En particulier, si $x$ est un VF, $x=RN(x)$.

La norme assure de plus que pour toute opération arithmétique standard, 

$$RN(x) \top RN(y) == RN(x\top y)$$


Il faut aussi penser que rien n'est gratuit:

\begin{NB}
  S'il  y a  un \textit{epsilon}  dans votre  calcul, c'est  à vous  de vous  en
  occuper et de le  contrôler: la norme vous donne les  moyens de contrôler mais
  ce n'est pas la machine qui fait le boulot...c'est VOUS!

  N'oubliez pas  que vous avec  un avantage sur  la machine: vous  connaissez la
  valeur réelle du réel que vous manipulez ou faites manipuler.
\end{NB}



Voici une petite réflexion encore plus subtile...


\begin{DANGER}
 $0,001\times 2^0$  ce n'est pas
tout à fait pareil que sa forme normalisée $1,000\times 2^{-3}$.


Les 3 zéros rajoutés à droite n'ont sûrement aucune valeur!

Par  exemple si  on  calcule $11,111  -  11,110$,  mais sur  6  bits, on  aurait
peut-être un autre résultat qui dépend de bits ici inconnus...

Il y a donc plusieurs échelle de certitude!

\end{DANGER}

Quand des nombres  sont très proches, les soustraire peut  engendrer ce genre de
phénomène  qu'on appelle  \textit{annulation catastrophique}:  l'\textit{erreur relative}
est bien  plus grande que l'\textit{erreur  absolue} (nous reparlerons de  ces notions un
peu plus loin).)


On pourra s'intéresser à l' \vref{ex::delta} pour méditer sur le danger
relatif de  certaines annulations (cancellations in  angliche): toute annulation
n'est pas catastrophique!



\begin{DANGER}
  Les machines calculent vite,  c'est un bien mais aussi un  mal: si on effectue
  une  somme avec  une toute  petite erreur,  une grande  somme peut  nous faire
  perdre tous nos bits  en une seconde car la somme de  beaucoup de très petites
  erreurs peut créer une grosse erreur...
\end{DANGER}

Mais le plus grand danger vient de:

\begin{quote}
  {\itshape
    95 \% of the folks out there are completely clueless about floating-point
}
\begin{flushright}
  James \textsc{Gosling} (M. Java) - 28 février 1998
\end{flushright}
\end{quote}



\begin{center}
  \includegraphics[width=0.5\textwidth]{ingenieur}
\end{center}



\subsection{Maîtriser les erreurs}


Tout  d'abord, une  déception: il  n'existe  pas de  traitement automatique  des
erreurs sur les VF.

Cependant un  espoir: il  existe des  moyens de traiter  les problèmes  selon le
contexte...mais qui nécessitent un bagage minimum en mathématique ;-)


D'abord un peu de vocabulaire. Il faudra bien distinguer:

\begin{description}
\item[Précision]  (\textit{precision}): c'est  le nombre  de bits  utilisés pour
  représenter un nombre. La précision concerne  donc le format utilisé pour écrire ou
  stocker  ou   arrondir  un  nombre.  Par   exemple  \texttt{3},  \texttt{3.0},
  \texttt{3.0000},  \texttt{3.0e0} n'ont  pas la  même précision,  précision qui
  dépend en plus du langage.
\item[Exactitude] (\textit{accuracy}): c'est ce qui  relie un nombre au contexte
  dans lequel il est employé. Connaître la distance qui sépare votre domicile de
  votre lieu de travail  au mètre près peut être considéré  comme très proche du
  résultat exact mais connaître votre
  taille  au mètre  près  l'est moins...L'exactitude  est liée  à  la mesure  de
  l'erreur.
\end{description}


Nous définirons deux types de mesure d'erreur:

\begin{definition}[Erreurs]
  Soit $x$ un nombre et $\widehat{x}$ le nombre qui le représente. On distingue:
  \begin{description}
  \item[l'erreur absolue] $|x - \widehat{x}|$
  \item[l'erreur relative] $η=\frac{x-\widehat{x}}{x}$ alors $\widehat{x}=x(1+η)$
  \end{description}
commise en prenant $\widehat{x}$ à la place de $x$.



\end{definition}


On  prend   souvent  comme  approximation   de  l'erreur  relative   le  rapport
  $|\log(x) - \log(\widehat{x})|$.   Nous
justifierons cette approximation lors de l'étude du calcul différentiel.

Pour reprendre un exemple  de W. \cite{kahanJava},
$3,177777777777777$  est une  approximation plutôt  précise (16  décimales) mais
inexacte (2 décimales) de $\pi$.



Comme nous  allons le  voir dans de  nombreux exemples  (e.g. \vref{ex::delta}),
analyser  les erreurs  avec  un peu  de mathématique  permet  de contrôler  son
programme.


On peut mesurer  les erreurs à deux niveaux. Lorsqu'on  étudie un algorithme qui
donne un résultat $r$  en fonction d'une donnée $d$, on  peut considérer que $r$
est une fonction de $d$:$r=f(d)$.

On peut analyser l'erreur commise sur $r$ ou l'erreur commise sur $d$.

\begin{definition}
  L'erreur  \textit{aval} (\textit{forward  error}) est  la différence  entre le
  résultat mesuré $\widehat{r}$ et le résultat théorique $r$.

L'erreur  \textit{amont}, ou  erreur \textit{inverse}  (\textit{backward error})
est le plus petit $|Δ d|$ tel que $f(d + Δ d ) = \widehat{r}$.
\end{definition}














\begin{quote}
  {\itshape
   Feel nervous, but feel in control.
It’s not dark magic, it’s science.
}
\begin{flushright}
  Florent \textsc{de Dinechin}
\end{flushright}
\end{quote}

C'est ce peu de mathématique que nous allons découvrir dans ce cours...

Mais un point non négligeable que nous n'aborderons pas (nous l'évoquerons) est
comment le processeur, l'OS, le langage (mal)traitent la norme IEEE 754...

\begin{NB}
  La norme ne vous transforme pas en  Harry \textsc{Potter}: il ne suffit pas de
  se dire qu'il y a une norme et que  tout va bien se passer en appuyant sur son
  clavier-baguette magique!

  La norme vous donne les moyens de  maîtriser les forces du mal \textbf{mais ça
    peut demander un sacré effort de réflexion (c'est-à-dire de mathématique :-) ) }.
\end{NB}


%\begin{figure}
\begin{center}
  \includegraphics[width=\linewidth]{AS400}
\end{center}
%\end{figure}


\section{Que la force de l'erreur soit avec vous}

\subsection{Atelier Padawan \# 1: majorer l'erreur}

\subsubsection{Précision machine}

\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3.4cm]{padawan}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-4cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}


Soit $x$ un nombre. On peut l'encadrer par deux VF successifs:

$$M\times 2^{E-n}\leqslant  x < (M+1)\times 2^{E-n}=M\times 2^{E-n} +2^{E-n}$$

avec M un entier, $n$ la longueur de la pseudo-mantisse.

L'approximation $\widehat{x}$ est donc la plus  proche de ces deux bornes de cet
intervalle de longueur $2^{E-n}$. Avec l'arrondi par défaut:

$$|x-\widehat{x}|\leqslant \frac{1}{2}2^{E-n}$$

donc l'\textbf{erreur relative} peut être majorée

$$\left|\frac{x-\widehat{x}}{x}\right|\leqslant
\left|\frac{x-\widehat{x}}{m\times 2^E}\right|\leqslant
\frac{1}{2}\frac{2^{E-n}}{m\times 2^E}=\frac{1}{2}\frac{2^{-n}}{m}=\frac{1}{2}\frac{ε_M}{m}$$

avec $m$ la mantisse normalisée.

Si $\widehat{x}$ n'est pas un VF sous-normal, alors $1 \leqslant m <2$ donc


$$\left|\frac{x-\widehat{x}}{x}\right|\leqslant\frac{1}{2}ε_M$$

L'erreur relative  commise est donc  majorée par la  moitié du successeur  de 1,
i.e. l'\textit{epsilon} de la machine.

Si $\widehat{x}$ est sous-normal, alors l'\textbf{erreur absolue} est majorée:

$$|x-\widehat{x}|\leqslant \frac{1}{2}2^{E_{\min}-1-n}$$


Ainsi, dans tous les cas:



\begin{NB}[Précision]

$\widehat{x}=x(1+η_1)+η_2$
avec:

\begin{itemize}
\item si  $\widehat{x}$ est normal $|η_1|\leqslant  \frac{1}{2}ε_M$ est l'erreur
  relative et $η_2=0$
\item si $\widehat{x}$ est sous-normal $η_1=0$ et $|η_2|\leqslant  \frac{1}{2}2^{E_{\min}-1-n} $
\end{itemize}


On appelle  $u=\frac{1}{2}ε_M$ l'unité d'arrondi (\textit{unit  roundoff}) ou la
précision machine.

\end{NB}

Par exemple,  la précision  machine en \textit{binary64}  vaut $2^{-52}/2\approx
1.1\times 10^{-16}$.

\begin{DANGER}[Précision $\neq$ exactitude!]

$R((a\cdot    b)\cdot   c)=R(a\cdot    b)\cdot   R(c)=a\cdot    b(1+η_{ab})\cdot
c(1+η_c)=a\cdot b\cdot c(1+η_{ab})(1+η_c)\approx a\cdot b\cdot c(1+η_{ab}+η_c)$
  
La précision est toujours $u$ mais l'exactitude est d'environ $2u$...

\begin{center}
\textbf{L'exactitude n'est pas assurée par la précison!}
\end{center}
\end{DANGER}


\begin{remarque}
  On  préfère étudier  l'erreur  relative  car elle  ne  dépend  pas du  facteur
  d'échelle:  si  $x$ et  $\widehat{x}$  sont  multipliés  par un  facteur  $k$,
  l'erreur absolue aussi mais l'erreur relative reste inchangée.
\end{remarque}


Oui, bon, OK, mais cette majoration est parfois très pessimiste et ne permet pas
d'aller très loin.

\subsubsection{Éliminiation}

L'élimination (\textit{cancellation} in english of speciality) est le phénomène qui apparaît lorsqu'on soustrait deux nombres
très proches.

Avec    les   notations    habituelles,   $\widehat{a}=a(1+η_a)$,
$\widehat{b}=b(1+η_b)$, $x=a-b$ et $\widehat{x}==\widehat{a-b}\widehat{a}-\widehat{b}$.
Alors 

$$
\left|\frac{x-\widehat{x}}{x}\right|=\left|\frac{-a η_a+b η_b}{a-b}\right| \leqslant \max(|η_a|,|η_b|)\frac{|a|+|b|}{|a-b|}
$$

On n'obtient ici qu'une majoration, très pessimiste d'ailleurs.

Cependant on constate que l'erreur est  d'autant plus grande que $|a-b| \ll |a|+|b|$ et de plus cette
erreur augmente les erreurs déjà présentes dans $\widehat{a}$ et $\widehat{b}$.

Ainsi, l'erreur  causée par l'élimination  remettent au premier plan  de petites
erreurs faites précédemment.

Parfois  cette erreur  n'est pas  forcément  dramatique. D'abord  parce que  les
données initiales sont peut-être exactes et donc aucune erreur n'est amplifiée.
Ensuite, l'erreur causée peut être \og écrasée\fg{} par une valeur plus grande.
Nous étudierons ce phénomène par exemple dans l'\vref{ex::delta}.

\subsection{Atelier Padawan \# 2: quitter la réalité (un peu difficile...)}


\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3.4cm]{vador2}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-4.2cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}



Dans l'atelier précédent,  nous avons travaillé avec des  outils \og réels\fg{}:
majorations, minorations avec  des opérations définies sur $\bbr$  et nous avons
oublié  qu'en fait  nous travaillons  dans $\bbv$:  il est  temps d'enlever  nos
scaphandres pour être plus à l'aise...


\textit{Exceptionnellement dans cette section nous travaillerons en base $β$ quelconque}.

\subsubsection{Trouver la base (algorithme de \textsc{Malcolm-Gentleman})}\label{algo::base}

Observez ce  qui se passe  et expliquez pourquoi cela  se passe sachant  que sur
JavaScript, il n'y a qu'un seul type numérique: \textit{binary64}...
\js
 \begin{multicols}{3}
\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{javascript}
Rhino 1.7 release 3 

js> Math.pow(2, 53) 
9007199254740992
js> Math.pow(2, 53) + 1
9007199254740992
js> Math.pow(2, 53) + 2
9007199254740994
\end{minted}


\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{javascript}
js> Math.pow(2, 53) + 3
9007199254740996
js> Math.pow(2, 53) + 4
9007199254740996
js> Math.pow(2, 53) + 5
9007199254740996
js> Math.pow(2, 54) 
18014398509481984
\end{minted}

\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{javascript}
js> Math.pow(2, 54) + 2
18014398509481984
js> Math.pow(2, 54) + 3
18014398509481988
js> Math.pow(2, 55) 
36028797018963970
js> Math.pow(2, 55) + 4
36028797018963970
\end{minted}
\end{multicols}


\haskell


Et maintenant, voici un algorithme en impératif qui, si on travaillait dans $\bbr$ ne
terminerait  pas mais  que  vous  allez pouvoir  comprendre  après notre  petite
observation sur JavaScript::

\begin{algo}[numbered]
\ALGO{de \textsc{Malcolm-Gentleman}}
\STATE{$a$\recoit{}1.0, $b$\recoit{}1.0}
\WHILE{$R(R(a+1.0) - a) == 1.0$}
\STATE{$a$\recoit{} $R(2\times a)$}
\ENDWHILE
\WHILE{$R(R(a+b) - a) \neq b$}
\STATE{$b$\recoit{} $b+1$}
\ENDWHILE
\RETURN{b}
\end{algo}


\begin{haskellcode}
base = boucle2 (boucle1 1) 1
  where 
    boucle1 a =
      if (a + 1.0) - a /= 1.0 then a  
      else boucle1 (2.0 * a)
    boucle2 a b =
      if (a + b) - a == b then b
      else boucle2 a (b + 1.0)
\end{haskellcode}

Regardons la première boucle.
On y construit une suite $(a_i)$ vérifiant pour tout naturel $i$ $a_{i+1}=2a_i$
avec $a_0=1.0$. Une récurrence immédiate démontre que $a_i=2^i$ pour tout entier
naturel $i$.


...Mouais,  c'est  exact  SI  l'on travaille  dans  $\bbr$...Mais  nous
travaillons dans $\bbv_β$!

\begin{quote}
  {\itshape Dans quel espace tu travailles tu me diras et  si ta
    proposition est vraie je saurai}

  \begin{flushright}
    Mathemator - 44 ap. GC
  \end{flushright}
\end{quote}


En fait notre espace de travail est $\bbv_β$ avec $β$ notre base (qui
est  donc  un  entier  supérieur  ou  égal   à  2)  dont  la  mantisse  (pas  la
pseudo-mantisse...) est de longueur $p$.


Les lignes 4 à 6 construisent une suite $a_{i}=R(2a_{i-1})$ avec $a_0=1.0$ mais est-on sûr que
l'on a toujours $a_i=2^i$ dans ces conditions?

On démontre  d'abord par  récurrence que $a_i=R(2^i)=2^i$  tant que  $i$ vérifie
$2^i  \leqslant β^p-1$:  $2^i=2\times 2^{i-1}$  est un  entier qui  s'écrit avec
moins de $p$  chiffres de la base $β$  et en deçà de $β^p-1$, on  peut passer d'un
flottant à son successeur par des pas inférieurs à 1.

Alors $R(a_i+1.0)=a_i+1.0$: pas d'arrondi. 

Puis $R(R(a_i+1.0) -a_i)=R(a_i+1.0-a_i)=1.0$: la condition de la première boucle
est donc vérifiée tant que $2^i <β^p$. Que se passe-t-il après?


 


Dès que $β^p-1$ est atteint, les choses changent.

Considérons la première itération $i_s$ telle que $2^{i_s}\geqslant β^p$.

On   a   $$a_{i_s}=R(2\times  a_{i_s-1})=R(2\times   2^{i_s-1})<R(2\times   β^p)
\leqslant R(β\times β^{p})=β^{p+1}$$

car d'une part  $ a_{i_s-1}=2^{i_s-1}$, d'autre part  $a_{i_s-1} \leqslant β^p-1
<β^p$ et enfin $2 \leqslant β$.

Ainsi $$β^p \leqslant a_{i_s} < β^{p+1}$$


L'exposant de la forme normale de $a_{i_s}$  est donc $p$ et d'après le résultat
suivant démontré à la section \vref{sec::succ}, le successeur  d'un VF $v$
vérifie (en notant $n=p-1$ la longueur de la pseudo-mantisse):

$$\operatorname{succ}(v)=v+β^{E(v)-n}=v+β^{E(v)-p+1}$$

Ici                                                                        ${\rm
  succ}(a_{i_s})=a_{i_s}+β^{E(a_{i_s}-p+1)}=a_{i_s}+β^{p-p+1}=a_{i_s}+β$.

On  en  déduit   que  $a_{i_s}+1.0$  est  entre  $a_{i_s}$   et  son  successeur
$a_{i_s}+β$.

Donc,  selon  l'arrondi,  $R(a_{i_s}+1.0)$  vaut  $a_{i_s}$  ou  $a_{i_s}+β$  et
finalement
$R(R(a_{i_s}+1.0)-a_{i_s})$ vaut  0 ou $β$ mais  en  aucun cas  1.0: on sort
donc de la boucle.

\begin{center}
\textit{On  en  déduit que  \texttt{boucle1  1}  vaut  $a_{i_s}$ donc  que  $β^p
  \leqslant \mathtt{boucle1} \ 1 < β^{p+1}$}
\end{center}



Notons $a= \mathtt{boucle1} \ 1$. Son successeur est $a+β$.

Passons à la seconde boucle.

\medskip

Tant que $b<β$...mais je vous laisse conclure la démonstration :-)



\begin{NB}[Moralité]

Cette   démonstration   est   plutôt   technique  et   difficile   en   première
lecture. Cependant elle est pleine d'enseignements:

\begin{itemize}
\item elle nous
montre bien que nous changeons de monde: un test du style 
\lstinline|if (a + 1.0) - a /= 1.0| paraîtrait bien hors de propos si nous raisonnions avec des
réels;
\item nous avons vu que l'on pouvait raisonner rigoureusement sur les VF;
\item nous pouvons  malgré tout retenir la trame intuitive  de la démonstration:
  tant que  notre nombre entier  est représentable  avec $p$ chiffres,  on reste
  exact. Les problèmes  arrivent lorsqu'on n'a plus assez de  place pour stocker
  tous les chiffres: il y a de la perte d'information...
\item Ce genre de raisonnement se  retrouve très souvent pour travailler sur les
  erreurs commises: il faudra donc en retenir la substantifique moelle...
\end{itemize}

\end{NB}


N'oubliez pas que certains langages ne travaillent pas en base 2, par exemple Maple\copyright:

\maple


\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{mupad}
> a := 1.0:
> b := 1.0:
> while evalf(evalf(a + 1.0) - a) - 1.0 = 0.0 do
      a := 2.0 * a;
  end_while:
> while evalf(evalf(a + b) - a) - b <> 0.0 do
      b := b + 1.0
  end_while:
> b;
                                 10.0
\end{minted}
                               
\haskell


\subsubsection{Nouvelles opérations}

Pour alléger les  notations et bien comprendre que nous  changeons de monde mais
que  la structure  des  calculs sur  les flottants  reste  algébrique (mais  des
propriétés intéressantes comme l'associativité disparaissent), nous noterons:

\begin{notation}
  Les opérations arithmétiques  de base sont bien définies sur  $\bbv$ donc nous
  noterons à présent, si l'arrondi n'a pas besoin d'être précisé:

  \begin{itemize}
  \item $R(x+y)$ sous la forme $x\oplus y$;
  \item $R(x-y)$ sous la forme $x\ominus y$;
  \item $R(x\times y)$ sous la forme $x\otimes y$;
  \item $R(x/y)$ sous la forme $x \oslash y$.
  \end{itemize}

Ainsi on notera par exemple $a\oplus b = a + b + {\rm err}(a\oplus b)$.

\end{notation}




\subsubsection{Trouver la précision}\label{algo::prec}

Encore un algorithme de M.A. \textsc{Malcolm} proposé en 1972.

Cette fois, vous vous occuperez de la démonstration tout(e) seul(e)...


\begin{haskellcode}
precision = toPrec 1.0 0
  where 
    toPrec a i = 
      if (a + 1.0) - a /= 1.0 then i 
      else toPrec (base * a) (i + 1)
\end{haskellcode}




\subsection{Atelier Padawan \# 3: somme de flottants proches}


\subsubsection{Quelques lemmes utiles}


\setindent{-0.5cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3.8cm]{vador3}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-3.8cm}
\begin{minipage}[t]{\linewidth}
\end{minipage}



Un lemme  en mathématique est  un résultat intermédiaire  qui sert à  prouver un
résultat  plus  important.  Un  lemme  peut  être  cependant  très  compliqué  à
démontrer: la démonstration du \textit{lemme fondamental} a valu à Bao Châu \textsc{Ngô}
d'obtenir en 2010 la médaille \textsc{Fields}...


Fini de jouer...Les démonstrations des algorithmes précédents sont intéressantes 
mais les résultats  obtenus sont anecdotiques: on sait bien  dans quelle base on
travaille et avec quelle précision.


Passons à présent  à des résultats beaucoup plus utiles  maintenant que nous sommes
au parfum des démonstrations sur les flottants. 


Un  des grands  dangers de  la manipulation  des flottants  est d'effectuer  des
soustractions  car   on  peut  perdre  beaucoup   d'information  par  annulation
(\textit{cancellation}).

Pour avoir un meilleur contrôle de la chose, nous allons étudier un lemme  énoncé par
Pat  H. \cite{sterbenz}  (1927 -  2008)...mais  pour le  démontrer, nous  aurons
besoin de quelques sous-lemmes ;-)


\begin{DANGER}[Convention]
  Pour  éviter de  déduire des  lemmes suivants  des théorèmes  totalement faux,
  n'oubliez pas que \textbf{DANS CE QUI SUIT X ET Y SONT DES NOMBRES À VIRGULE FLOTTANTE}!
\end{DANGER}


\begin{lemme}[Majoration de l'erreur d'une somme]\label{lemme::majSom}
Posons $x\oplus y =  x + y + \operatorname{err}(x\oplus y)$. Alors, s'il  n'y a pas de
dépassement de capacité, 

$$\left|\operatorname{err}(x\oplus y)\right| \leqslant \min  \left(|x|,|y|\right) 
$$

On a bien sûr un résultat analogue pour la différence

  
\end{lemme}


On a $x\oplus y = x + y + \operatorname{err}(x\oplus y)$ (on est dans $\bbr$...)

De plus $x = x + y + (-y)$ (on est toujours dans $\bbr$...)


$x$ et  $x\oplus y$  sont donc deux  VF distants respectivement  de $|y|$  et de
$\left|\operatorname{err}(x\oplus y)\right|$ de $x+y$ (que de de, de deux et de 2: pffff....)


Mais la norme  nous assure que c'est $x\oplus  y$ le VF le plus  proche de $x+y$
donc  l'erreur  commise  en  le  choisissant  est  moindre  que  celle  faite  en
choisissant $-y$.

On en déduit que $\left|\operatorname{err}(x\oplus y)\right| \leqslant |y|$.

On montre de même que  $\left|\operatorname{err}(x\oplus y)\right| \leqslant |x|$.


\begin{remarque}
  De l'importance de disposer avec la IEEE 754 de la \textbf{meilleure approximation}!
\end{remarque}



\begin{corollaire}[M\o ller, Knuth, Dekker]
  L'erreur commise $\left|\operatorname{err}(x\oplus y)\right|$ peut être exprimée
  exactement sur $p$ bits.
\end{corollaire}

Ce  n'est  qu'un  corollaire  d'un  sous-lemme,  mais  c'est  un  résultat  très
important!

Supposons, sans perdre de généralité, que $|x| \geqslant |y|$.


Comme $x$ et $y$ sont des VF, le  plus petit bit significatif de $\operatorname{err}(x\oplus y)$ est
au moins de magnitude celle de ulp($y$).

De plus  $\left|\operatorname{err}(x\oplus y)\right| \leqslant |y|$,  donc la mantisse
entière de $\operatorname{err}(x\oplus y)$ a une longueur inférieure à $p$ bits.

Finalement   le  mantisse   entière  de   $\left|\operatorname{err}(x\oplus  y)\right|
\leqslant |y|$ est exactement exprimable sur $p$ bits.



\begin{lemme}\label{lemme::plusExact}
  Supposons que $|x+y|\leqslant \min(|x|,|y|)$, alors $x\oplus y = x + y$.

On obtient un résultat analogue pour la soustraction.
\end{lemme}


La démonstration va ressembler à la précédente.

Supposons, sans perdre de généralité, que $|x| \geqslant |y|$.


 Le  plus petit bit significatif de $x+y$ est
au moins de magnitude celle de ulp($y$).


Or $|x+y| \leqslant |y|$ donc la mantisse
entière de $x+y$  a une longueur inférieure à $p$ bits.

Finalement le mantisse entière de $x+y$ est exactement exprimable sur $p$ bits.


\begin{DANGER}[N'oubliez pas la convention!]\label{danger::convention}
  Nous avons démontré nos lemmes en considérant des VF $x$ et $y$.

  Est-ce que le résultat suivant contredit notre dernier lemme?

  \begin{haskellcode}
*Main> 1 - 0.9
9.999999999999998e-2
  \end{haskellcode}

\end{DANGER}



\subsubsection{Lemme de Sterbenz}


 

\begin{lemme}[Lemme de \textsc{Sterbenz} (1973)]\label{sterbenz}
  Soit $(x,y)\in \bbv^2$  vérifiant $\dfrac{x}{2} \leqslant y  \leqslant 2x$. On
  a $$x\ominus y=x-y$$
La différence de deux VF suffisamment proches est donc exacte.
\end{lemme}


Une petite  démonstration?... Allez, le  lemme précédent va sûrement  être utile
pour avoir un petit schéma de démonstration qui tient sur une ligne :-)

Pour  une démonstration  plus  complète,  on pourra  par  exemple  se référer  à
\cite{MullerEtAl2010} page 123.


On peut se restreindre à des nombres positifs.

On a deux cas:

\begin{enumerate}
\item  $\mathbf{x <  y}$:  alors  $x <  y  \leqslant  2x$ donc  $0
  < y-x \leqslant x \leqslant y$;
\item $\mathbf{x  \geqslant y}$: alors  $ \frac{x}{2} \leqslant y  \leqslant x$
  donc $ -\frac{x}{2} \leqslant y-x \leqslant 0$ et par suite $0 \leqslant x - y
  \leqslant \frac{x}{2} \leqslant y \leqslant x$.
\end{enumerate}

Dans tous les cas $|x-y| \leqslant \min(|x|,|y|)$ et on peut utiliser le lemme 
\vref{lemme::plusExact}.




\begin{remarques}\label{rem::}
  Concrètement, à quoi correspond cette condition $\dfrac{x}{2} \leqslant y \leqslant 2x$?

Demandez-vous ce  que l'on peut dire  de l'écart maximum entre  les exposants de
$x$ et $y$.


On admettra que \textbf{le résultat reste vrai même en cas de sous-capacité}.


\end{remarques}





\subsection{Atelier Padawan \# 4: somme compensée de flottants quelconques} 



\setindent{0cm}
\begin{flushleft}
\begin{minipage}[t]{3cm}
%\begin{figure}
  \includegraphics[width=3.3cm]{vador}
%\caption{\foreignlanguage{russian}{Николай Петрович Брусенцов}(1925 -)}
%\end{figure}
\end{minipage}
\end{flushleft}
\setindent{\indentationcours}
\vspace{-4.8cm}
\begin{minipage}[t]{\linewidth}


\end{minipage}


Les lemmes précédents  ne permettent de s'assurer de l'exactitude  de la somme
que dans des cas restreints.

Que  se passe-t-il  quand on  ne peut  pas s'assurer  à chaque  instant que  les
opérandes sont suffisamment proches?


\subsubsection{Fast2Sum}



Voici une  première version  officiellement due à  T.J. \textsc{Dekker}  en 1971
mais qui  avait déjà  été introduite par  l'incontournable W.  \textsc{Kahan} en
1965.






\begin{theoreme}[Fast2Sum - \textsc{Dekker} \& \textsc{Kahan}]\label{th::Fast2Sum}
On considère  deux VF $x$  et $y$ tels que  $|x| \geqslant |y|$  et l'algorithme
suivant:

\begin{algo}[numbered]
\STATE{$s$ \recoit{} $x\oplus y$ }
\STATE{$y_v$\recoit{} $s\ominus x$}
\STATE{$d$ \recoit{} $y \ominus y_v$}
\RETURN{$(s,d)$}
\end{algo}
  
Alors $x+y=s+d$ avec $s=x\oplus y$ et $d={\rm err}(x\oplus y)$. De plus $s$ et $d$ ne se chevauchent pas.

\end{theoreme}





ce qui donne dans notre langage favori:


\begin{haskellcode}
fastTwoSum x y  
    | abs x >= abs y = 
      let s  = x + y   in
      let yv = s - x  in 
      let d  = y - yv in 
      (s, d)
    | otherwise =  fastTwoSum y x
\end{haskellcode}



La ligne 1 donne $s=x + y + \operatorname{err}(x \oplus y)$.




De plus, comme $|x| \geqslant |y|$, c'est $x$ \og qui gagne\fg{} donc $s$ et $x$
ont le même signe. 

La ligne 2  calcule un $y$ virtuel  ($y_v$) qui vaut $s \ominus  x$: on aimerait
bien que ça fasse $s-x$ c'est-à-dire $y + \operatorname{err}(x \oplus y)$. 

Enfin la  ligne 3 calcule  $y\ominus y_v$: on aimerait  bien aussi que  ça fasse
$y-y_v$  c'est-à-dire $-  \operatorname{err}(x  \oplus y)$  comme  ça on  récupérerait
exactement l'erreur commise. 

En plus, comme $d=- \operatorname{err}(x \oplus  y)=$, on a $|d| \leqslant \frac{1}{2}
\operatorname{ulp}(s)$ donc $s$ et $d$ ne se chevauchent pas dans leur somme.


\smallskip

On aimerait  bien, on  aimerait bien...En  fait tout  s'arrange car  nous allons
pouvoir utiliser le lemme de \textsc{Sterbenz}.

Il suffit de distinguer deux cas:

\begin{itemize}
\item  \textbf{si} $\mathbf{x}$  \textbf{et} $\mathbf{y}$  \textbf{sont de  même
    signe    OU    si}    $\mathbf{|y|    \leqslant    \frac{|x|}{2}}$:    alors
  $\frac{x}{2}\leqslant s \leqslant 2x$ et on peut appliquer le lemme;
\item \textbf{sinon}:  on a $x$ et  $y$ de signes opposés  ET $y> \frac{|x|}{2}$
  alors si $y$  est négatif $x/2 < -y  < x$ et sinon $-x/2<y<-x$.  Dans les deux
  cas, d'après  le lemme  de \textsc{Sterbenz}, $s$  est calculée  exactement et
  alors $y_v=y$.
\end{itemize}


On illustrera cet algorithme par un exemple simple en TD dans l'\vref{exo::fast2sum}


On illustrera encore mieux le problème par cette... illustration:





\newcolumntype{C}[1]{>{\centering}b{\widthof{#1}}}
\newcommand\MyBox[2][XXX]{%
  \fbox{\tabular[b]{C{#1}} #2 \endtabular}}


$x$\hspace{1em} \MyBox[mmmmmmmmm]{$x_1$}
\MyBox[mmmm]{$x_2$}

$+y$ \hspace{10em}\MyBox[mmmm]{$y_1$}
\MyBox[mmmmmmmmm]{$y_2$}

$=s$   \hspace{0.15em}\MyBox[mmmmmmmmm]{$x_1$}
\MyBox[mmmm]{$x_2+y_1$} \ $y_2$ est perdu

$-x$ \hspace{10.15em}\MyBox[mmmm]{$y_1$}
\MyBox[mmmmmmmmm]{0}

$-y$\hspace{15.75em} \MyBox[mmmmmmmmm]{$-y_2$}
\MyBox[mmmm]{0} \ on récupère l'erreur

Ainsi, $x+y=s+d$ se lit:



$x$\hspace{1em} \MyBox[mmmmmmmmm]{$x_1$}
\MyBox[mmmm]{$x_2$}

$+y$ \hspace{10em}\MyBox[mmmm]{$y_1$}
\MyBox[mmmmmmmmm]{$y_2$}

$=s$   \hspace{0.15em}\MyBox[mmmmmmmmm]{$x_1$}
\MyBox[mmmm]{$x_2+y_2$} \ $y_2$ est perdu

$+d$\hspace{15.75em} \MyBox[mmmmmmmmm]{$y_2$}
\MyBox[mmmm]{0} \ on récupère l'erreur




\begin{DANGER}
$s$ et $d$ ne se chevauchent pas: $s\oplus d=s$!
  
\end{DANGER}






\subsubsection{2Sum}

Il  s'avère  que sur  les  machines  modernes,  la comparaison  préliminaire  du
Fast2Sum peut s'avérer plus coûteuse que les trois additions supplémentaires qui
vont suivre.

De plus, le Fast2Sum ne fonctionne qu'en base 2 ou 3 alors que l'algorithme suivant
fonctionne dans toutes les bases.



Mais  il faut  encore  rentrer dans  plus de  détails  comme l'éventualité  d'un
parallélisme (lignes 2 et 3), etc...vous verrez ça en Master...



\begin{theoreme}[2Sum - \textsc{Knuth} ]\label{th::2Sum}
On considère  deux VF $x$  et $y$   et l'algorithme
suivant:

\begin{algo}[numbered]
\STATE{$s$ \recoit{} $x\oplus y$ }
\STATE{$y_v$\recoit{} $s\ominus x$, $x_v$\recoit{} $s\ominus y$}
\STATE{$e_y$\recoit{} $y\ominus y_v$, $e_x$\recoit{} $x\ominus x_v$}
\STATE{$d$ \recoit{} $e_x \oplus e_y$}
\RETURN{$(s,d)$}
\end{algo}
  
Alors $x+y=s+d$ avec $s=x\oplus y$ et $d={\rm err}(x\oplus y)$. De plus $s$ et $d$ ne se chevauchent pas.

\end{theoreme}





\subsubsection{Somme compensée d'un nombre quelconque de flottants}


Voici  qu'entre à  nouveau dans  l'arène W.  \textsc{Kahan} avec  son algorithme
\textit{compensated summation} de 1965 (un peu en avance sur son temps...).

Le \og \textit{one liner}\fg{} en Haskell:\label{algo::sommeKahan}

\begin{haskellcode}
sommeKahan liste =
  let (s,e) = foldl (\ (s,c) x -> twoSum s (x + c)) (0,0) liste in s
\end{haskellcode}

Sauriez-vous le traduire en impératif avec une boucle \texttt{pour}? Il vaudrait
mieux pour répondre aux questions de l'\vref{exo::sum}.

Est-ce optimum? En 1972, Michèle \textsc{Pichat}\citep{Pichat1972CdS} mais aussi
Arnold \textsc{Neumaier} (mais l'article est en allemand) ont l'idée de sommer à
part les erreurs  et de soustraire le tout  à la fin. Cet algorithme  (dit de la
\textit{sommation en cascade}) est repris en
2008 par Siegfried M. \textsc{Rump}, Takeshi \textsc{Ogita}
et   Shin'Ichi   \textsc{Oishi}\citep{Rump:2008:AFS:1461600.1461604}   mais   en
introduisant le fastSum...faites de même!






 Le schéma suivant vous inspirera sûrement:

\begin{center}
\tikzstyle{int}=[draw, fill=blue!20, minimum size=2em]
\tikzstyle{init} = [pin edge={<-,thin,black}]
\tikzstyle{fin} = [pin edge={->,thin,black}, pin position = below]
\tikzstyle{sort} = [pin edge={->,thin,black}, pin position = right]
\tikzstyle{entre} = [pin edge={<-,thin,black}, pin position = left]

\begin{tikzpicture}[node distance=2.9cm,auto,>=latex']
  \node [int, pin={[init]$x_1$}, pin={[fin]$e_1$},pin={[sort]$s_1$}] (a) {\texttt{twoSum}}; %bloc 1
  \node (b) [left of=a,node distance=1.5cm, coordinate] {a}; % noeud d'entrée
  \node [int,pin={[entre]}, pin={[init]$x_2$}, pin={[fin]$e_2$},pin={[sort]$s_2$}] (c) [right of=a] {\texttt{twoSum}}; % bloc 2
  \node [int,pin={[entre]}, pin={[init]$x_3$}, pin={[fin]$e_3$},pin={[sort]$s_3$}] (d) [right of=c] {\texttt{twoSum}}; % bloc 3
  \node [int,pin={[entre]}, pin={[init]$x_4$}, pin={[fin]$e_4$}] (e) [right of=d] {\texttt{twoSum}}; % bloc 4
  \node [coordinate] (end) [right of=e, node distance=1.5cm]{}; % noeud fin
  \draw[->] (b)node[left] {$x_0$} --(a) ; % entrée
  \draw[->] (e) --(end)node[right] {$s$} ; % fin
 \end{tikzpicture}

\end{center}


Voyons tout de suite un exemple. On veut effectuer la somme en \textit{binary64}
de $x_0=2^{53}-1$, $x_1=2^{53}$ et $x_2=-(2^{54}-2)$.

Que vaut la somme exacte (sans machine car la machine se trompe...)?


Voyons maintenant le schéma du fonctionnement sur machine:

\begin{center}
\tikzstyle{int}=[draw, fill=blue!20, minimum size=2em]
\tikzstyle{init} = [pin edge={<-,thin,black}]
\tikzstyle{fin} = [pin edge={->,thin,black}, pin position = below]
\tikzstyle{sort} = [pin edge={->,thin,black}, pin position = right]
\tikzstyle{entre} = [pin edge={<-,thin,black}, pin position = left]

\begin{tikzpicture}[node distance=3.1cm,auto,>=latex']
  \node [int, pin={[init]$2^{53}$}, pin={[fin]$-1$},pin={[sort]$2^{54}$}] (a) {\texttt{twoSum}}; %bloc 1
  \node (b) [left of=a,node distance=1.5cm, coordinate] {a}; % noeud d'entrée
  \node [int,pin={[entre]}, pin={[init]$-(2^{54}-2)$}, pin={[fin]$0$}] (c) [right of=a] {\texttt{twoSum}}; % bloc 2
  \node [coordinate] (end) [right of=c, node distance=1.5cm]{}; % noeud fin
  \draw[->] (b)node[left] {$2^{53}-1$} --(a) ; % entrée
  \draw[->] (c) --(end)node[right] {$2$} ; % fin
 \end{tikzpicture}
\end{center}

\begin{recherche}\label{exo::2sum}
  Expliquez les résultats trouvés. Que va faire l'algorithme de somme compensée?
  Et celui de somme en cascade?
\end{recherche}




On  note  $u=\frac{ε_M}{2}$ et  $γ_n=\frac{n  u}{1-nu}$.  Alors nos  trois  amis
\cite{Rump:2008:AFS:1461600.1461604} ont démontré le théorème suivant:

\begin{theoreme}
 En                  appliquant                 l'algorithme                  de
 \textsc{Pichat}-\textsc{Rump}-\textsc{Ogita}-\textsc{Oishi}   à   une   famille
 $(x_i)_{1 \leqslant i \leqslant  n}$ de flottants et si $nu  <1$ alors, même en
 cas de sous-capacité mais sans sur-capacité, le résultat $s$ retourné par
 l'algorithme vérifie:

$$
\left|s-      \displaystyle\sum_{i=1}^{n}x_i\right|\leqslant       u      \left|
E  \displaystyle\sum_{i=1}^{n}x_i\right|+γ_{n-1}^2 \displaystyle\sum_{ i=1}^{n}\left|x_i\right|
$$



\end{theoreme}

Ça fait peur, hein?


\section{Et la multiplication? Et la division? Et le sinus? Et...} 

Pour la multiplication, on procède à peu près de la même manière. On peut
gagner du temps si  le processeur utilise la FMA (Fused  Multiple Add) mais ceci
est une autre histoire...


Et la division?  On cherche l'inverse $x$  de $a$: on est donc  amené à résoudre
l'équation $\frac{1}{x}-a=0$ d'inconnue $a$ qui n'est pas linéaire.

Et la racine carrée? Là aussi, on quitte le linéaire.

Il nous  faut donc explorer  des outils un  peu plus sophistiqués...ce  que nous
allons faire dans le chapitre suivant.

% Il nous faut aller bien plus loin pour explorer autre chose que des sommes!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% EXERCICES
%
%%%%%


\vfill

\eject













\modeexercice

\section{EXERCICES}



\begin{exercice}
  Comment est codé  2 en tant qu'entier  32 bits? En tant  que binary32? Comment
  est-il codé en mémoire sous forme hexadécimale?
\end{exercice}



\begin{exercice}
Donnez au format \textit{binary32} et \textit{binary64} les représentations
mémoires hexadécimales des nombres 125,25 puis 0,375 puis $-0,375$ puis 0,44.

\end{exercice}

\begin{exercice}
Donnez la représentation mémoire hexadécimale au format \textit{binary32} de $2^{-130}$.
\end{exercice}

\begin{exercice}


Comment  traduiriez-vous  par une  égalité  qu'un  nombre $x\in  [0,1[$  s'écrit
$b_1b_2...b_n$ dans une base $β$?


Soit l'algorithme:

\begin{algo}
\FUNC{multSuccessives}{\pfarg{x}{réel} \pfarg{base}{entier naturel}}{liste
  d'entiers naturels}
\STATE{nb\recoit{}x}
\STATE{dec = []}
\WHILE{nb $\neq$ 0}
\STATE{bi\recoit{} $\mathtt{nb\times base}$}
\STATE{dec \recoit{} dec ++ liste($\lfloor\mathtt{ bi} \rfloor$)}
\STATE{nb\recoit{} $\{\mathtt{bi}\}$} 
\ENDWHILE
\RETURN{dec}
\end{algo}

On notera $\{x\}$ la partie fractionnaire de $x$.

Que pensez-vous de sa précondition? De sa postcondition? de sa terminaison? De sa correction? De sa complexité? 
\end{exercice}


\begin{exercice}
  On suppose qu'on dispose d'une fonction \mintinline{java}{FloatToIntBits} 
qui renvoie l'écriture machine en base 2  sur 32 bits d'un réel et des fonctions
bit à bit habituelles.

Donnez des  fonctions renvoyant le signe,  l'exposant, la mantisse d'un  réel 32
bits.


Pour vous faire plaisir, vous pourrez coder tout ça en Java...Par exemple, voici
une fonction qui donne la  représentation mémoire hexadécimale d'un floattant 32
bits:

\begin{javacode}
static String memoire(float reel){
    int bits;
    bits = Float.floatToIntBits(reel);
    return (Integer.toHexString(bits));
}
\end{javacode}


\end{exercice}




\begin{exercice}\label{ex::der}
\begin{multicols}{2} 
  Que pensez-vous de ça :

  \begin{haskellcode}
der f h x =
  (f(x + h) - f(x)) / h
  
cosApp = der sin 0.0000001
  
sinApp = \x -> - (der cosApp 0.0000001 x)

cosApp2 = der sinApp 0.00000001 
  \end{haskellcode}


Qui donne:

\begin{haskellcode}
*Main> sinApp (pi/6)
0.5051514762044462




*Main> cosApp2 (pi/2)
1110223.0246251565
\end{haskellcode}
\end{multicols}
\end{exercice}


%


\begin{exercice}
  Pourquoi \lstinline+3 * 0.1 /= 0.3+? 
Commentez:

\minihaskell

\begin{haskellcode}
Prelude Numeric> floatToDigits 2 0.1
([1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1 0,1],-3)
Prelude Numeric> floatToDigits 2 (2*0.1)
([1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1],-2)
Prelude Numeric> floatToDigits 2 (3*0.1)
([1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1],-1)
Prelude Numeric> floatToDigits 2 0.3
([1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],-1)

Prelude Numeric> floatToDigits 16 0.1
([1,9,9,9,9,9,9,9,9,9,9,9,9,10],0)
Prelude Numeric> floatToDigits 16 (3*0.1)
([4,12,12,12,12,12,12,12,12,12,12,12,13],0)
Prelude Numeric> floatToDigits 16 0.3
([4,12,12,12,12,12,12,12,12,12,12,12,12,12],0)
\end{haskellcode}

\haskell

\end{exercice}



%

\begin{exercice}

  Pourquoi le \textit{Patriot} a raté son \textit{Scud} de 600m?

Représentez 1/10 avec une mantisse de 24 bits et tronquez le résultat.

Calculez l'erreur en seconde puis  après 100 heures d'utilisation. Sachant qu'un
Scud vole à $\SI{1676}{m.s^{-1}}$, quelle est environ l'erreur commise en mètres?

Et si on avait arrondi au plus proche?


\end{exercice}

% 

\begin{exercice}
  Représentez 1,32 sur  \textit{toy16}: 1 bit de signe, 4  bits d'exposant et 11
  bits pour $f$.
\end{exercice}


%

\begin{exercice}\label{ex::scale}
  En utilisant vos carreaux, représentez les VF normaux de \textit{toy7}.
  Faites de même avec les sous-normaux.
\end{exercice}

%

\begin{exercice}
Quel est le \og PGCD\fg{} de tous les \textit{binary64}?
\end{exercice}

%

\begin{exercice}
  Représentez  $1,0011$ et  $-1,0011$ en  \textit{toy7} selon  les quatre  modes
  d'arrondi usuels.

  Représentez également $0,00011$ en \textit{toy7} puis $0,0000001$.
\end{exercice}

%

\begin{exercice}
  Calculez $ulp(1)$ en \textit{binary32}, \textit{binary64} et \textit{toy7}.
\end{exercice}



%


\begin{exercice}
 %  En \textit{toy7}, \textit{binary32} et \textit{binary64}, quels sont:
% \begin{multicols}{2} 
%  \begin{itemize}
%   \item la longueur de l'exposant;
%   \item la longueur de la pseudo-mantisse;
%   \item l'exposant minimum;
%   \item l'exposant maximum;
%   \item le plus petit VF normal positif;
%   \item le plus petit VF sous-normal positif;
%   \item l'$ε$ de la machine.
%   \end{itemize}
% \end{multicols}

Complétez  le tableau \vpageref{tab::reca}.


\end{exercice}


%

\begin{exercice}
  
Le recours au  dixième peut être mortel  sur machine. Il vaut  mieux diviser par
deux mais est-ce pour autant sans danger?

  Effectuez la division répétée de 1 par 2 en \textit{toy7}. Que cela vous inspire-t-il?



\end{exercice}


%


\haskell

%

\begin{exercice}
\begin{enumerate}

\item   Comment expliquer les résultats suivants:

  \begin{haskellcode}
*Main> let f(x) = x^2 / sqrt(x^3 + 1)
*Main> f(1e100)
1.0e50
*Main> f(1e150)
0.0
*Main> f(1e200)
NaN
  \end{haskellcode}

\item  On  pose  \lstinline+let  f(x)  =  exp(x)/x+.  Prévoyez  ce  que  donnera
  \lstinline+f(1e500)+,        \lstinline+f(1e100)+,       \lstinline+f(1e-100)+,
  \lstinline+f(0)+.





\end{enumerate}



\end{exercice}




% 

\begin{exercice}
  Petit jeu  à 2:  effectuez toute  sorte de  multiplications et  d'additions en
  \textit{toy7} en comparant avec votre voisin(ne).
\end{exercice}


%

\begin{exercice}
  En utilisant \lstinline+floatToDigits+  de la bibliothèque \lstinline+Numeric+
  écrivez  une  fonction  qui  calcule  la  somme  de  deux  flottants  puis  la
  multiplication de deux flottants en n'oubliant pas les cas particuliers...

\end{exercice}




%

\begin{exercice}
  Soit $a$ et $b$ deux VF tels que $|a| \geqslant |b|$. 
Démontrez  que $a+b$
  est un multiple de $2^{E_b-p}$. 

Donnez alors une nouvelle démonstration du lemme \vref{lemme::majSom}
\end{exercice}









%

\begin{exercice}
  Terminez la démonstration de la section \vref{algo::base} et démontrez le programme de la section \vref{algo::prec}
\end{exercice}

%


\begin{exercice}
  Répondez à la question de l'encadré \textit{Danger} \vpageref{danger::convention}
\end{exercice}


%

\begin{exercice}[Somme exacte?]
  

Voici un petit code Haskell (ça vous rappelle un TP en
Archi-Prog avec Monsieur F.) : expliquez ce curieux \texttt{([8],0)}.
\minihaskell

\begin{haskellcode}
*Main> let ss n = sum [0.1 | k <- [1..n]]
*Main> import Numeric
*Main Numeric> [floatToDigits 16 (ss j) | j <- [5,10..20]]
([8],0),
([15,15,15,15,15,15,15,15,15,15,15,15,15,8],0),
([1,8,0,0,0,0,0,0,0,0,0,0,0,1],1),
([2,0,0,0,0,0,0,0,0,0,0,0,0,2],1)]
\end{haskellcode}

Mais attention! L'hexadécimal, ce n'est pas le binaire:

\begin{haskellcode}
*Main Numeric> [floatToDigits 2 (ss n) | n <- [1..10]]
[
([1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1],-3),
([1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1],-2),
([1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1],-1),
([1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,1],-1),
([1],0),
([1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],0),
([1,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],0),
([1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1],0),
([1,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1,0,0,1,1],0),
([1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1],0)
] 
\end{haskellcode}
\haskell

\end{exercice}



%

\begin{exercice}\label{exo::fast2sum}
Appliquez  l'algorithme  Fast2Sum  aux  deux   cas  suivants  en  posant  chaque
opération sachant qu'on travaille avec un mantisse de 4 bits:
\begin{multicols}{2}
  \begin{enumerate}
  \item $x=1111\times 2^2$ et $y=1001$;
  \item $x=1001\times 2^1$ et $y=-1011$.
  \end{enumerate}
\end{multicols}
\end{exercice}





%


\begin{exercice}
  Écrivez  une  fonction \texttt{twoSum}  en  Haskell  qui calcule  les  nombres
  \texttt{s} et \texttt{d} de l'algorithme du théorème \vref{th::2Sum}.
\end{exercice}


%

\begin{exercice}\label{exo::sum}
\begin{enumerate}
\item   Traduisez  la  fonction   \texttt{sommeKahan}  \vpageref{algo::sommeKahan}  en
  impératif en pseudo-spaghetti puis par exemple en C ou (et) en Java.
\item Répondez aux questions de l'encadré de la section \vref{exo::2sum}.
\item On voudrait calculer $ \displaystyle\sum_{k=1}^{\num{100000}}RN\left(\frac{1}{k}\right)$.
    \begin{enumerate}
    \item Utilisez une somme récursive standard.
    \item Utilisez \texttt{sommeKahan}.
    \item        Voici       les        premières       bonnes        décimales:
      $\num{12,0901461298634280}$. Estimez  les erreurs commises en  les exprimant
      en ulps.
    \item Et $ \displaystyle\sum_{k=1}^{\num{100000}}\frac{1}{k}$ dans tout ça?
    \end{enumerate}
\end{enumerate}

On pourra utiliser le squelette suivant:

\haskell

\begin{haskellcode}
import qualified Data.List as L -- pour avoir reverse et sort

-- twoSum : somme compensée de 2 flottants sans branchement
twoSum x y =
  let s       = x + y            in
  let [xv,yv] = [s - y , s - x ] in
  let [ex,ey] = [x - xv, y - yv] in
  (s , ex + ey)
      
-- somme basique (+) des éléments d'une liste classée dans l'ordre croissant
sommeC liste =

-- somme basique (+) des éléments d'une liste classée dans l'ordre décroissant
sommeD liste =

-- somme compensée (Kahan) : on rajoute les raisins perdus après chaque cep 
-- secoué
sommeK liste = 
  let (sn,en) = foldl (\ (sk, ek) xk' -> twoSum sk (xk' + ek)) (0,0) liste in sn

-- somme en cascade (Pichat) : on collecte les raisins perdus à part et on 
-- les rajoute à la fin
sommeP liste =
\end{haskellcode}

\end{exercice}


%

\begin{exercice}\label{ex::delta}
  J'ose espérer que vous savez résoudre  une équation du type $ax^2+bx+c=0$ avec
  $a\neq 0$...

Les  racines, si  elles  existent,  sont données  par  une  formule bien  connue
dépendant de $a$, $b$ et $c$:

$$r=\frac{-b\pm     \sqrt{b^2-4ac}}{2a}=\frac{-b\pm      \sqrt{Δ}}{2a}     \quad
\text{avec}\quad Δ = b^2-4ac$$





Où peuvent se cacher d'éventuelles  annulations catastrophiques? Étudiez ces cas
avec attention, voyez si vous pouvez éviter les éliminations catastrophiques en
réécrivant   les   formules   un   peu   dans   l'esprit   de   la   \og   levée
d'indétermination\fg{}.

\begin{enumerate}
\item Que  se passe-t-il lorsque  $b^2 \gg |4ac|$?  En quoi la  formule $\frac{-
    \left(b+\operatorname{signe}(b)\sqrt{Δ}\right)}{2a} $ peut aider? 
\item  Que se  passe-t-il  lorsque  $b^2\approx 4ac$?  Peut-on  y remédier?  Que
  peut-on dire de $Δ$ par rapport à $b^2$? Y a-t-il élimination catastrophique?
\item  Que  se  passe-t-il  dans   le  cas  de  l'équation  $10^{200}x^2-3\times
  10^{200}x+2\times 10^{200}=0$?
\item Et dans le cas de $10^{-200}x^2-3x+2\times 10^{200}=0$?
\item Moralité?


\end{enumerate}



\end{exercice}





%

\begin{exercice}[Bits de garde, d'arrondi, collant]
  Soustrayez $1,101\times 2^{-2}$ puis $1,111\times 2^{-2}$ à $1\times 2^0$ en \textit{toy7} en gardant
  tous les bits, même ceux qui dépassent  le format de la mantisse puis refaites
  le calcul en faisant disparaître les bits qui dépassent de la mantisse.

Pourquoi parle-ton de \textit{bit de garde}?

En fait, il y a toujours trois bits à droite de la mantisse: G le bit de garde (\textit{guard}),
R le bit d'arrondi (\textit{rounding}) et S le bit collant (\textit{sticky}).

Les deux  premiers sont comme  des paniers  qui recueillent les  bits translatés
(\textit{shifted}). Le dernier  est un panier absorbant  qui est un ou  inclusif des bits
translatés au-delà de R.
\end{exercice}






%


\begin{exercice}[Calcul de exp 1 et puissances de 10]
On veut calculer une approximation de  e ($= \exp(1)$) comme limite au voisinage
de l'infini de $ \left(1+\frac{1}{n}\right)^n $.

  \begin{enumerate}
  \item Rappel:  $f$ est dérivable au voisinage de $x$ si, et seulement si, $
    \displaystyle\lim_{t\to  x} \frac{f(t)-f(x)}{t-x}\in\bbr$.  Dans  ce cas  on
    note $f'(x)$ cette valeur.

  \item Utiliser ce rappel pour calculer $
    \displaystyle\lim_{t\to 0}\frac{\ln(1+t)}{t}$.

  \item Déduisez-en $
    \displaystyle\lim_{n\to +\infty} n \ln \left(1+\frac{1}{n}\right) $ puis $
    \displaystyle\lim_{n\to +\infty}\left(1+\frac{1}{n}\right)^n$.

\item On note \lstinline|expn n = (1 + 1/n)**n|. Comparez et commentez:

  \begin{haskellcode}
*Main> [abs (exp 1 - (expn (10.0**k)))    | k <- [0..20]]
*Main> [abs (exp 1 - (expn (2.0**(3*k)))) | k <- [0..20]]
  \end{haskellcode}


  \end{enumerate}


\end{exercice}

%


\begin{exercice}[0,1 en base 2]
  Nous démontrerons dans  le chapitre suivant que, pour tout  réel $q$ vérifiant
  $|q|<1$ on a $$ \displaystyle\sum_{k=k_0}^{+\infty}q^k=\frac{q^{k_0}}{1-q}$$

  \begin{enumerate}
  \item       Démontrez       que      $       \displaystyle\sum_{i=1}^{+\infty}
    \left(2^{-4i}+2^{-4i-1}\right)=\frac{1}{10}$.
  \item Déduisez-en le développement de 1/10 en base 2.

  \item           Montrez          que,           en          \textit{binary32},
    $\frac{\widehat{x}-x}{x}=-\frac{ε_M}{8}$.
  \end{enumerate}
\end{exercice}

% 

\begin{exercice}[Max IEEE]

On peut calculer le maximum de deux réels ainsi:

\begin{algo}[noframe]
\FUNC{max}{\pfarg{$x$,$y$}{flottants}}{flottant}
\IF{$x>y$}
\STATE{$\max = x$}
\ELSE
\STATE{$\max = y$}
\ENDIF
\end{algo}

Est-ce adapté à la recherche du maximum de deux flottants en arithmétique IEEE?
  
\end{exercice}



%

\begin{exercice}[Options GCC et double arrondi]

\ce

Voici un petit code C pour effectuer une multiplication \citep{MullerEtAl2010}:

\begin{ccode}
#include <stdio.h>

int
main(void)
{
  double a = 1848874847.0;
  double b = 19954562207.0;
  double c;  c = a * b;
  printf("c = %20.19e\n", c);
  return 0;
}
\end{ccode}
 

et  voici  les  résultats  obtenus  en utilisant  deux  options  de  compilation
différentes:


\shell

\begin{itemize}
\item Avec \lstinline+mfpmath=387+:



\begin{bashcode}
$ gcc -mfpmath=387 multi.c -o multi_387
$ ./multi_387
c = 3.6893488147419103232e+19
\end{bashcode}

\item Avec \lstinline+mfpmath=sse+:

\begin{bashcode}
$ gcc -mfpmath=sse multi.c -o multi_sse
$ ./multi_sse
c = 3.6893488147419111424e+19
\end{bashcode}
\end{itemize}


Que se passe-t-il sachant qu'on trouve dans \lstinline+man gcc+ les lignes suivantes:


\begin{bashcode}
 -mfpmath=unit
           Generate floating-point arithmetic for selected unit unit.  The choices for unit are:

           387 Use the standard 387 floating-point coprocessor present on the majority of chips and emulated otherwise.  Code compiled with this option runs almost everywhere.  The temporary results are computed in 80-bit precision  instead of the precision specified by the type, resulting in slightly different results compared to most of other chips.  This is the default choice for i386 compiler.

           sse Use scalar floating-point instructions present in the SSE instruction set.  This instruction set is supported by Pentium3 and newer chips, in the AMD line by Athlon-4, Athlon-xp and Athlon-mp chips.  The earlier version of SSE instruction set supports only single-precision arithmetic, thus the double and extended-precision arithmetic are still done using 387.  A later version, present only in Pentium4 and the future AMDx86-64 chips, supports double-precision arithmetic too.
\end{bashcode}

et qu'en précision infinie on trouve:

\haskell

\begin{haskellcode}
Prelude>  1848874847 * 19954562207
36893488147419107329
Prelude> 36893488147419107329 - 36893488147419111424
-4095
Prelude>  36893488147419107329 -  36893488147419103232
4097
\end{haskellcode}


et que la représentation binaire exacte du produit est:

$$
\overbrace{\underbrace{10000000000000000000000000000000000000000000000000000}_{53
    \text{ bits}}10000000000}^{64 \text{ bits}}01
$$


et que le calcul effectué avec l'arrondi au plus près en \textit{binary64} est

\begin{haskellcode}
*Main Numeric> floatToDigits 2 (1848874847 * 19954562207)
([1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
    0,0,1],66)
\end{haskellcode}


On pourra méditer sur le schéma suivant:



%\begin{figure}
\begin{center}
  \includegraphics{DrapeauFloat.7}
\end{center}
%\end{figure}


\end{exercice}


%

\begin{exercice}[Double arrondi]

Nous   montrerons    dans   le   chapitre   suivant    que   $\sqrt{1-u}\approx
1-\frac{1}{2}u-\frac{1}{8}u^2$ pour de petites valeurs de $u$.


Le \textit{double étendu} est  un format sur 80 bits avec  15 bits d'exposant et
64 bits de pseudo-mantisse.

  Quel est le  résultat de $\sqrt{1-2^{-53}}$ en \textit{binary64}  avec ou sans
  double arrondi depuis un format étendu?

On  pourra   écrire  le  développement  en   binaire  sur  110  bits   de  $1  -
\frac{1}{2}u-\frac{1}{8}u^2$ et arrondir selon le mode  387 ou selon le mode SSE
comme cela a été illustré dans l'exercice précédent..




\end{exercice}



%

\begin{exercice}[MPFR]

La bibliothèque  MPFR qui  est une  extension GNU  depuis 2009  a été  fondée au
départ en 1998 par Guillaume \textsc{Hanrot}, Jean-Michel \textsc{Muller}, Joris
van der \textsc{Hoeven} et Paul \textsc{Zimmermann}.

Elle  a d'abord  été conçue  pour  permettre d'obtenir  des résultats  corrects,
portables et  reproductible, parfois au  dépends de  la vitesse mais  elle reste
très efficace.

Par exemple:

\ce

\begin{minted}[linenos,breaklines,breakindent=1cm,frame=lines,framesep=2mm,fontsize=\small,fontfamily=txtt]{c}
#include <stdio.h>
#include <mpfr.h>

int main (void) 
{
  mpfr_t un , dix , trois, zero_1, zero_3, trois_fois_zero_1, diff ; /* Déclarations simultanées */
  int test;

  mpfr_set_default_prec (53) ; /* Fixer la precision par defaut */
  mpfr_set_emin ( -1073) ; /* Fixer l' exposant emin ( en realite , emin - precision +2) */
  mpfr_set_emax (1024) ; /* Fixer l' exposant emax ( en realite , emax +1) */
  mpfr_inits (zero_1,zero_3, trois_fois_zero_1, (mpfr_ptr) 0); /* Précision par defaut 53 bits */  
  mpfr_inits2 (128, un, trois, dix , diff,  (mpfr_ptr) 0) ; /* Précision de 128 bits exactement */
  mpfr_set_d (un , 1.0 , MPFR_RNDN ); /* MPFR_RND Division en arrondi au plus près */
  mpfr_set_d (trois , 3.0 , MPFR_RNDN );
  mpfr_set_d (dix , 10. , MPFR_RNDN );

  mpfr_div (zero_1 , un , dix , MPFR_RNDN );
  mpfr_div (zero_3 , trois , dix , MPFR_RNDN );
  mpfr_mul (trois_fois_zero_1, trois, zero_1, MPFR_RNDN  );
  mpfr_sub (diff, trois_fois_zero_1,  zero_3, MPFR_RNDN  );

  test = mpfr_equal_p (trois_fois_zero_1, zero_3);
  
  mpfr_printf (" Avec %3d bits de précision , 1/10     ~ %RNe \n" , mpfr_get_prec (zero_1) , zero_1);
  mpfr_printf (" Avec %3d bits de précision , 1/10     ~ %.38RNe \n" , mpfr_get_prec (zero_1) , zero_1);
  mpfr_printf (" Avec %3d bits de précision , 1/10     ~ %.38Rg \n" , mpfr_get_prec (zero_1) , zero_1);
  mpfr_printf (" Avec %3d bits de précision , 3/10     ~ %.38RNe \n" , mpfr_get_prec (zero_3) , zero_3);
  mpfr_printf (" Avec %3d bits de précision , 3 * 1/10 ~ %.38RNe \n" , mpfr_get_prec (trois_fois_zero_1) , trois_fois_zero_1);
  printf (" 3 * 1/10 == 3/10 ? (0 si faux, /= 0 si vrai) :  %u \n" , test);
  mpfr_printf (" Avec %3d bits de précision , la différence 3*1/10 - 0.3 vaut  ~ %.38RNe \n" , mpfr_get_prec (diff) , diff);

  mpfr_clears ( un , dix , trois, zero_1, zero_3, trois_fois_zero_1, diff, (mpfr_ptr) 0);
  return 0;
}
\end{minted}


renvoie

\shell

\begin{bashcode}
$ gcc -o mpfr_essai mpfr_essai.c -lmpfr -lgmp
$ ./mpfr_essai

 Avec  53 bits de precision , 1/10     ~ 1.0000000000000001e-01 
 Avec  53 bits de precision , 1/10     ~ 1.00000000000000005551115123125782702118e-01 
 Avec  53 bits de precision , 1/10     ~ 0.10000000000000000555111512312578270212 
 Avec  53 bits de precision , 3/10     ~ 2.99999999999999988897769753748434595764e-01 
 Avec  53 bits de precision , 3 * 1/10 ~ 3.00000000000000044408920985006261616945e-01 
 3 * 1/10 == 3/10 ? (0 si faux, /= 0 si vrai) :  0 
 Avec 128 bits de precision , la différence 3*1/10 - 0.3 vaut  ~ 5.55111512312578270211815834045410156250e-17 
\end{bashcode}


Créez alors une fonction qui calcule la somme des inverses des entiers de 1 à $\num{100000}$ en
prenant la précision en argument.

Sur  Sage,  on peut  utiliser  la  bibliothèque MPFR  à  l'aide  de la  fonction
\texttt{RealField(p)} avec \texttt{p} la précision:

\sagec

\begin{pythoncode}
sage:  RealField(128)(RealField(53)(1)/ RealField(53)(10))
0.10000000000000000555111512312578270212
sage:  RealField(128)(RealField(53)(3)/ RealField(53)(10))
0.29999999999999998889776975374843459576
sage:  RealField(128)(RealField(53)(3) *  RealField(53)(0.1))
0.30000000000000004440892098500626161695
\end{pythoncode}


Il  existe aussi  sur Haskell  une  bibliothèque \verb+hmpfr+  qui nécessite  de
recompiler \verb+ghc+ avec une version alternative de \verb+gmp+.
\haskell






\end{exercice}


\begin{exercice}[6 points]
Répondez aux questions suivantes en justifiant vos réponses:
  \begin{enumerate}
  \item Est-ce que l'addition des VF est associative?
  \item Est-ce que la multiplication des VF est associative?
  \item Est-ce que la multiplication est distributive sur l'addition? 
  \end{enumerate}
\end{exercice}


%4

\begin{exercice}[6 points]
\begin{enumerate}

\item   Comment expliquer les résultats suivants:

\begin{multicols}{2}
  \begin{lstlisting}
*Main> let f(x) = x^2 / sqrt(x^3 + 1)

*Main> f(1e100)
1.0e50
\end{lstlisting}
  \begin{lstlisting}
*Main> f(1e150)
0.0
*Main> f(1e200)
NaN
  \end{lstlisting}
\end{multicols}


\item On pose \lstinline|f(x) = (x**2 - x) / (2 * (sqrt (x**4 - x**2)))|.  Prévoyez ce que donnera
  \lstinline+f(1e500)+,        \lstinline+f(1e100)+,       \lstinline+f(1e-100)+,
  \lstinline+f(0)+.


\item Comment améliorer les calculs faits par la machine?


\end{enumerate}



\end{exercice}

\pagebreak

%5

\begin{exercice}[4 points]
On considère la fonction suivante:



\begin{lstlisting}
p = toP 1.0 0
  where 
    toP a cpt = 
      if (a + 1.0) - a /= 1.0 then cpt 
      else toP (2 * a) (cpt + 1)
\end{lstlisting}

Imaginons que Haskell travaille uniquement en  \textit{toy7} (1.0 est donc un VF
codé sur 7 bits). Que vaut alors p? Justifiez votre réponse.


\end{exercice}


\begin{exercice}
  
Exemple dû à Jean-Michel Muller  (in \cite{MullerEtAl2010})


M. X a récemment été à sa banque (Chaotic Bank Society), pour connaître les nouvelles offres proposées aux
meilleurs clients. Son banquier lui propose l'offre suivante: \og vous déposez tout d'abord e - 1 euros sur votre compte
(où e = 2.7182818... est la base du logarithme népérien). La première année, nous prenons 1 euro sur votre compte de
frais de gestion. Par contre, la deuxième année est plus intéressante pour vous, car nous multiplions votre capital restant
par 2 et prenons 1 euro de frais de gestion. La troisième année est encore plus intéressante, car nous multiplions votre
capital par 3 et prenons 1 euro de frais de gestion. Et ainsi de suite : la n-ième année, nous multiplions votre capital par n
et prenons 1 euro de frais de gestion. Intéressant, non? \fg{} Pour prendre sa décision, M. X décide de demander l'aide
d'un informaticien et se retourne vers vous.

\end{exercice}


\begin{exercice}[Bac]

Le sujet du Bac S des Centres Étrangers 2012 fait étudier la suite:

$$
I_{n+2}= \frac{1}{2} {\rm e} -  \frac{n+1}{2} I_n, \qquad I_1= \frac{1}{2} ({\rm
  e} -1)
$$
  


Calculez $I_{59}$ à la machine...



\end{exercice}

%

% \begin{exercice}
% On considère la fonction $f$ définie sur $\bbv$ par $f(x)=(((x-0,5)+x)-0,5)+x$.

% Démontrez  que   $f(x)$  ne  s'annule  jamais   sur  $\bbv$,  que  ce   soit  en
% \textit{binary32} ou en \textit{binary64}.

% \textit{Joker}: utilisez le lemme de \textsc{Sterbenz} trois fois.
% \end{exercice}

\modenormal


%\nocite{*}




%%% Local Variables: 
%%% TeX-master: "PolyAnalyse14"
%%% End: 

