Added some -not-so-funny- examples

9 files changed, 520 insertions, 23 deletions

diff --git a/scilab_doc/scilabisnotnaive/introduction.tex b/scilab_doc/scilabisnotnaive/introduction.tex
index b697405..db7affb 100644
--- a/scilab_doc/scilabisnotnaive/introduction.tex
+++ b/scilab_doc/scilabisnotnaive/introduction.tex
@@ -5,12 +5,42 @@ While most mathematic books deals with exact formulas,
 Scilab uses algorithms which are specifically designed for 
 computers. 
 
+As a practical example of the problem considered in this
+document, consider the following experiments. The following is an example of 
+a Scilab 5.1 session, where we compute 0.1 by two ways.
+
+\begin{verbatim}
+-->format(25)
+-->0.1
+ ans  =
+    0.1000000000000000055511  
+-->1.0-0.9
+ ans  =
+    0.0999999999999999777955  
+\end{verbatim}
+
+I guess that for a person who has never heard of these problems,
+this experiment may be a shock. To get things clearer, let's 
+check that the sinus function is approximated.
+
+\begin{verbatim}
+-->format(25)
+-->sin(0.0)
+ ans  =
+    0.  
+-->sin(%pi)
+ ans  =
+    0.0000000000000001224647  
+\end{verbatim}
+
 The difficulty is generated by the fact that, while 
 the mathematics treat with \emph{real} numbers, the 
 computer deals with their \emph{floating point representations}.
 This is the difference between the 
 \emph{naive}, mathematical, approach, and the \emph{numerical},
 floating-point aware, implementation.
+(The detailed explanations of the previous examples are presented  
+in the appendix of this document.)
 
 In this article, we will show examples of these problems by 
 using the following theoric and experimental approach.
@@ -36,18 +66,11 @@ digits in the computed value $x_c$. For example, if the relative
 error $e_r=10^{-6}$, then the number of significant digits is 6.
 
 When the expected value is zero, the relative error cannot 
-be computed, and we then use the absolute error 
+be computed, and we then use the absolute error instead 
 \begin{eqnarray}
 e_a=|x_c-x_e|.
 \end{eqnarray}
 
-A central reference on this subject is the article 
-by Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic", 
-\cite{WhatEveryComputerScientist}.
-If one focuses on numerical algorithms, the "Numerical Recipes" \cite{NumericalRecipes},
-is another good sources of solutions for that problem.
-The work of Kahan is also central in this domain, for example \cite{Kahan2004}.
-
 Before getting into the details, it is important to 
 know that real variables in the Scilab language are stored in 
 \emph{double precision} variables. Since Scilab is 
@@ -56,3 +79,5 @@ variables are stored with 64 bits precision.
 As we shall see later, this has a strong influence on the 
 results.
 
+
+
diff --git a/scilab_doc/scilabisnotnaive/numericalderivative.tex b/scilab_doc/scilabisnotnaive/numericalderivative.tex
index c95a3b3..7ce0768 100644
--- a/scilab_doc/scilabisnotnaive/numericalderivative.tex
+++ b/scilab_doc/scilabisnotnaive/numericalderivative.tex
@@ -12,12 +12,6 @@ In the third part, we analyse
 why and how floating point numbers must be taken into account when the 
 numerical derivatives are to compute.
 
-A reference for numerical derivates 
-is \cite{AbramowitzStegun1972}, chapter 25. "Numerical Interpolation, 
-Differentiation and Integration" (p. 875).
-The webpage \cite{schimdtnd} and the book \cite{NumericalRecipes} give
-results about the rounding errors.
-
 \subsection{Theory}
 
 The basic result is the Taylor formula with one variable \cite{dixmier}
@@ -307,7 +301,7 @@ a function of one variable is presented in figure \ref{robust-numericalderivativ
 \label{robust-numericalderivative}
 \end{figure}
 
-\subsection{One step further}
+\subsection{One more step}
 
 In this section, we analyse the behaviour of \emph{derivative}
 when the point $x$ is either large $x \rightarrow \infty$, 
@@ -426,3 +420,12 @@ large values of $x$, but produces a slightly sub-optimal step size when
 $x$ is near 1. The behaviour near zero is the same, i.e. both commands 
 produce wrong results when $x \rightarrow 0$ and $x\neq 0$.
 
+\subsection{References}
+
+A reference for numerical derivates 
+is \cite{AbramowitzStegun1972}, chapter 25. "Numerical Interpolation, 
+Differentiation and Integration" (p. 875).
+The webpage \cite{schimdtnd} and the book \cite{NumericalRecipes} give
+results about the rounding errors.
+
+
diff --git a/scilab_doc/scilabisnotnaive/pythagoreansum.tex b/scilab_doc/scilabisnotnaive/pythagoreansum.tex
new file mode 100644
index 0000000..5833c69
--- /dev/null
+++ b/scilab_doc/scilabisnotnaive/pythagoreansum.tex
@@ -0,0 +1,180 @@
+\section{The Pythagorean sum}
+
+In this section, we analyse the computation of the Pythagorean sum,
+which is used in two different computations, that is the norm of a complex
+number and the 2-norm of a vector of real values.
+
+In the first part, we briefly present the mathematical formulas for these 
+two computations.
+We then present the naïve algorithm based on these mathematical formulas. 
+In the second part, we make some experiments in Scilab and compare our
+naïve algorithm with the \emph{abs} and \emph{norm} Scilab primitives.
+In the third part, we analyse 
+why and how floating point numbers must be taken into account when the 
+Pythagorean sum is to compute.
+
+\subsection{Theory}
+
+\subsection{Experiments}
+
+% TODO : compare both abs and norm.
+
+\lstset{language=Scilab}
+\lstset{numbers=left}
+\lstset{basicstyle=\footnotesize}
+\lstset{keywordstyle=\bfseries}
+\begin{lstlisting}
+// Straitforward implementation
+function mn2 = mynorm2(a,b)
+  mn2 = sqrt(a^2+b^2)
+endfunction
+// With scaling
+function mn2 = mypythag1(a,b)
+  if (a==0.0) then
+    mn2 = abs(b);
+  elseif (b==0.0) then
+    mn2 = abs(a);
+  else
+    if (abs(b)>abs(a)) then
+      r = a/b;
+      t = abs(b);
+    else
+      r = b/a;
+      t = abs(a);
+    end
+    mn2 = t * sqrt(1 + r^2);
+  end
+endfunction
+// With Moler & Morrison's
+// At most 7 iterations are required.
+function mn2 = mypythag2(a,b)
+  p = max(abs(a),abs(b))
+  q = min(abs(a),abs(b))
+  //index = 0
+  while (q<>0.0)
+    //index = index + 1
+    //mprintf("index = %d, p = %e, q = %e\n",index,p,q)
+    r = (q/p)^2
+    s = r/(4+r)
+    p = p + 2*s*p
+    q = s * q
+  end
+  mn2 = p
+endfunction
+function compare(x)
+  mprintf("Re(x)=%e, Im(x)=%e\n",real(x),imag(x));
+  p = abs(x);
+  mprintf("%20s : %e\n","Scilab",p);
+  p = mynorm2(real(x),imag(x));
+  mprintf("%20s : %e\n","Naive",p);
+  p = mypythag1(real(x),imag(x));
+  mprintf("%20s : %e\n","Scaling",p);
+  p = mypythag2(real(x),imag(x));
+  mprintf("%20s : %e\n","Moler & Morrison",p);
+endfunction
+// Test #1 : all is fine
+x = 1 + 1 * %i;
+compare(x);
+// Test #2 : more difficult when x is large
+x = 1.e200 + 1 * %i;
+compare(x);
+// Test #3 : more difficult when x is small
+x = 1.e-200 + 1.e-200 * %i;
+compare(x);
+\end{lstlisting}
+
+\begin{verbatim}
+***************************************
+Example #1 : simple computation with Scilab 5.1
+x(1)=1.000000e+000, x(2)=1.000000e+000
+              Scilab : 1.414214e+000
+               Naive : 1.414214e+000
+             Scaling : 1.414214e+000
+    Moler & Morrison : 1.414214e+000
+***************************************
+Example #2 : with large numbers ?
+              Scilab : Inf
+               Naive : Inf
+             Scaling : 1.000000e+200
+    Moler & Morrison : 1.000000e+200
+***************************************
+Example #3 : with small numbers ?
+x(1)=1.000000e-200, x(2)=1.000000e-200
+              Scilab : 0.000000e+000
+               Naive : 0.000000e+000
+             Scaling : 1.414214e-200
+    Moler & Morrison : 1.414214e-200
+***************************************
+> Conclusion : Scilab is naive !
+Octave 3.0.3
+***************************************
+octave-3.0.3.exe:29> compare(x);
+***************************************
+x(1)=1.000000e+000, x(2)=1.000000e+000
+              Octave : 1.414214e+000
+               Naive : 1.414214e+000
+             Scaling : 1.414214e+000
+    Moler & Morrison : 1.414214e+000
+***************************************
+x(1)=1.000000e+200, x(2)=1.000000e+000
+              Octave : 1.000000e+200
+               Naive : Inf
+             Scaling : 1.000000e+200
+    Moler & Morrison : 1.000000e+200
+***************************************
+octave-3.0.3.exe:33> compare(x)
+x(1)=1.000000e-200, x(2)=1.000000e-200
+              Octave : 1.414214e-200
+               Naive : 0.000000e+000
+             Scaling : 1.414214e-200
+    Moler & Morrison : 1.414214e-200
+***************************************
+> Conclusion : Octave is not naive !
+
+With complex numbers.
+***************************************
+
+Re(x)=1.000000e+000, Im(x)=1.000000e+000
+              Scilab : 1.414214e+000
+               Naive : 1.414214e+000
+             Scaling : 1.414214e+000
+    Moler & Morrison : 1.414214e+000
+***************************************
+Re(x)=1.000000e+200, Im(x)=1.000000e+000
+              Scilab : 1.000000e+200
+               Naive : Inf
+             Scaling : 1.000000e+200
+    Moler & Morrison : 1.000000e+200
+***************************************
+Re(x)=1.000000e-200, Im(x)=1.000000e-200
+              Scilab : 1.414214e-200
+               Naive : 0.000000e+000
+             Scaling : 1.414214e-200
+    Moler & Morrison : 1.414214e-200
+***************************************
+> Conclusion : Scilab is not naive !
+\end{verbatim}
+
+\subsection{Explanations}
+
+\subsection{References}
+
+The paper by Moler and Morrisson 1983 \cite{journals/ibmrd/MolerM83} gives an 
+algorithm to compute the Pythagorean sum $a\oplus b = \sqrt{a^2 + b^2}$
+without computing their squares or their square roots. Their algorithm is based on a cubically
+convergent sequence.
+The BLAS linear algebra suite of routines \cite{900236} includes the SNRM2, DNRM2
+and SCNRM2 routines which conpute the euclidian norm of a vector.
+These routines are based on Blue \cite{355771} and Cody \cite{Cody:1971:SEF}.
+In his 1978 paper \cite{355771}, James Blue gives an algorithm to compute the 
+Euclidian norm of a n-vector $\|x\| = \sqrt{\sum_{i=1,n}x_i^2}$. 
+The exceptionnal values of the \emph{hypot} operator are defined as the 
+Pythagorean sum in the IEEE 754 standard \cite{P754:2008:ISF,ieee754-1985}.
+The \emph{\_\_ieee754\_hypot(x,y)} C function is implemented in the 
+Fdlibm software library \cite{fdlibm} developed by Sun Microsystems and 
+available at netlib. This library is used by Matlab \cite{matlab-hypot}
+and its \emph{hypot} command.
+
+
+
+
diff --git a/scilab_doc/scilabisnotnaive/quadratic.tex b/scilab_doc/scilabisnotnaive/quadratic.tex
index 122e80d..1a888f9 100644
--- a/scilab_doc/scilabisnotnaive/quadratic.tex
+++ b/scilab_doc/scilabisnotnaive/quadratic.tex
@@ -229,12 +229,6 @@ IEEE overflow exception and set the result as \emph{Inf}.
 
 \subsection{Explanations}
 
-The technical report by G. Forsythe \cite{Forsythe1966} is 
-especially interesting on that subject. The paper by Goldberg 
-\cite{WhatEveryComputerScientist} is also a good reference for the 
-quadratic equation. One can also consult the experiments performed by 
-Nievergelt in \cite{Nievergelt2003}.
-
 The following tricks are extracted from the 
 \emph{quad} routine of the \emph{RPOLY} algorithm by
 Jenkins \cite{Jenkins1975}. This algorithm is used by Scilab in the 
@@ -498,7 +492,20 @@ should limit the possible overflows since $|b'/c| < 1$.
 These normalization tricks are similar to the one used by Smith in the 
 algorithm for the division of complex numbers \cite{Smith1962}.
 
-
-
+\subsection{References}
+
+The 1966 technical report by G. Forsythe \cite{Forsythe1966} 
+presents the floating point system and the possible large error 
+in using mathematical algorithms blindly. An accurate way of solving 
+a quadratic is outlined. A few general remarks are made about 
+computational mathematics. The 1991 paper by Goldberg 
+\cite{WhatEveryComputerScientist} is a general presentation of the floating
+point system and its consequences. It begins with background on floating point 
+representation and rounding errors, continues with a discussion
+of the IEEE floating point standard and concludes with examples of how
+computer system builders can better support floating point. The section
+1.4, "Cancellation" specificaly consider the computation of the roots
+of a quadratic equation.
+One can also consult the experiments performed by Nievergelt in \cite{Nievergelt2003}.
 
 
diff --git a/scilab_doc/scilabisnotnaive/references.tex b/scilab_doc/scilabisnotnaive/references.tex
new file mode 100644
index 0000000..8f7463e
--- /dev/null
+++ b/scilab_doc/scilabisnotnaive/references.tex
@@ -0,0 +1,10 @@
+\section{References}
+
+A central reference on this subject is the article 
+by Goldberg, "What Every Computer Scientist Should Know About Floating-Point Arithmetic", 
+\cite{WhatEveryComputerScientist}.
+If one focuses on numerical algorithms, the "Numerical Recipes" \cite{NumericalRecipes},
+is another good sources of solutions for that problem.
+The work of Kahan is also central in this domain, for example \cite{Kahan2004}.
+
+
diff --git a/scilab_doc/scilabisnotnaive/scilabisnotnaive.bib b/scilab_doc/scilabisnotnaive/scilabisnotnaive.bib
index 6de801c..36384a0 100644
--- a/scilab_doc/scilabisnotnaive/scilabisnotnaive.bib
+++ b/scilab_doc/scilabisnotnaive/scilabisnotnaive.bib
@@ -241,3 +241,119 @@ publisher= {Third Edition, Addison Wesley, Reading, MA}}
   bibsource = {DBLP, http://dblp.uni-trier.de}
 }
 
+ @article{journals/ibmrd/MolerM83,
+title = {Replacing Square Roots by Pythagorean Sums.},
+author = {Cleve B. Moler and Donald Morrison},
+journal = {IBM Journal of Research and Development},
+number = {6},
+pages = {577-581},
+url = {http://dblp.uni-trier.de/db/journals/ibmrd/ibmrd27.html#MolerM83},
+volume = {27},
+year = {1983},
+description = {dblp},
+date = {2002-01-03},
+keywords = {dblp }
+} 
+
+@techreport{900236,
+ author = {Lawson,, C. L. and Hanson,, R. J. and Kincaid,, D. R. and Krogh,, F. T.},
+ title = {Basic Linear Algebra Subprograms for FORTRAN Usage},
+ year = {1977},
+ source = {http://www.ncstrl.org:8900/ncstrl/servlet/search?formname=detail\&id=oai%3Ancstrlh%3Autexas_cs%3AUTEXAS_CS%2F%2FCS-TR-77-72},
+ publisher = {University of Texas at Austin},
+ address = {Austin, TX, USA},
+ }
+
+@article{355771,
+ author = {Blue,, James L.},
+ title = {A Portable Fortran Program to Find the Euclidean Norm of a Vector},
+ journal = {ACM Trans. Math. Softw.},
+ volume = {4},
+ number = {1},
+ year = {1978},
+ issn = {0098-3500},
+ pages = {15--23},
+ doi = {http://doi.acm.org/10.1145/355769.355771},
+ publisher = {ACM},
+ address = {New York, NY, USA},
+ }
+@InCollection{Cody:1971:SEF,
+  author =       "W. J. Cody",
+  title =        "Software for the Elementary Functions",
+  crossref =     "Rice:1971:MS",
+  pages =        "171--186",
+  year =         "1971",
+  bibdate =      "Thu Sep 15 18:56:47 1994",
+  bibsource =    "garbo.uwasa.fi:/pc/doc-soft/fpbiblio.txt",
+  acknowledgement = ack-nj,
+}
+
+@techreport{ieee754-1985,
+        author = {Stevenson, David  },
+        booktitle = {ANSI/IEEE Std 754-1985},
+        citeulike-article-id = {1677576},
+        journal = {ANSI/IEEE Std 754-1985},
+        keywords = {floating-point, para08},
+        month = {August},
+        posted-at = {2008-04-21 12:47:40},
+        priority = {1},
+        title = {IEEE standard for binary floating-point arithmetic},
+        url = {http://ieeexplore.ieee.org/xpls/abs\_all.jsp?arnumber=30711},
+        year = {1985}
+}
+
+@Book{P754:2008:ISF,
+  author =       "{IEEE Task P754}",
+  title =        "{IEEE 754-2008, Standard for Floating-Point
+                 Arithmetic}",
+  publisher =    pub-IEEE-STD,
+  address =      pub-IEEE-STD:adr,
+  pages =        "58",
+  day =          "29",
+  month =        aug,
+  year =         "2008",
+  DOI =          "http://doi.ieeecomputersociety.org/10.1109/IEEESTD.2008.4610935",
+  ISBN =         "0-7381-5753-8 (paper), 0-7381-5752-X (electronic)",
+  ISBN-13 =      "978-0-7381-5753-5 (paper), 978-0-7381-5752-8
+                 (electronic)",
+  bibdate =      "Thu Sep 25 09:50:30 2008",
+  URL =          "http://ieeexplore.ieee.org/servlet/opac?punumber=4610933;
+                  http://en.wikipedia.org/wiki/IEEE_754-2008",
+  abstract =     "This standard specifies interchange and arithmetic
+                 formats and methods for binary and decimal
+                 floating-point arithmetic in computer programming
+                 environments. This standard specifies exception
+                 conditions and their default handling. An
+                 implementation of a floating-point system conforming to
+                 this standard may be realized entirely in software,
+                 entirely in hardware, or in any combination of software
+                 and hardware. For operations specified in the normative
+                 part of this standard, numerical results and exceptions
+                 are uniquely determined by the values of the input
+                 data, sequence of operations, and destination formats,
+                 all under user control.",
+  acknowledgement = ack-nhfb,
+}
+
+@article{matlab-hypot,
+note={\url{http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/ref/hypot.html&http://www.google.fr/search?q=hypot&ie=utf-8&oe=utf-8&aq=t&rls=org.mozilla:fr:official&client=firefox-a}},
+title={Matlab - hypot : square root of sum of squares},
+}
+
+@article{fdlibm,
+title = {A Freely Distributable C Math Library},
+author = {Sun Microsystems, Inc.},
+year = {1993},
+note = {\url{http://www.netlib.org/fdlibm}},
+}
+
+@book{261217,
+ author = {Muller,, Jean-Michel},
+ title = {Elementary functions: algorithms and implementation},
+ year = {1997},
+ isbn = {0-8176-3990-X},
+ publisher = {Birkhauser Boston, Inc.},
+ address = {Secaucus, NJ, USA},
+ }
+
+
diff --git a/scilab_doc/scilabisnotnaive/scilabisnotnaive.pdf b/scilab_doc/scilabisnotnaive/scilabisnotnaive.pdf
index c7e1f0d..7ef8eaa 100644
--- a/scilab_doc/scilabisnotnaive/scilabisnotnaive.pdf
+++ b/scilab_doc/scilabisnotnaive/scilabisnotnaive.pdf
diff --git a/scilab_doc/scilabisnotnaive/scilabisnotnaive.tex b/scilab_doc/scilabisnotnaive/scilabisnotnaive.tex
index 0c35b37..aced19f 100644
--- a/scilab_doc/scilabisnotnaive/scilabisnotnaive.tex
+++ b/scilab_doc/scilabisnotnaive/scilabisnotnaive.tex
@@ -122,10 +122,17 @@ is not sufficiently accurate.
 
 \include{complexdivision}
 
+
 \include{conclusion}
 
 \clearpage
 
+\appendix
+
+\include{simpleexperiments}
+
+%\include{pythagoreansum}
+
 %% Bibliography
 
 
diff --git a/scilab_doc/scilabisnotnaive/simpleexperiments.tex b/scilab_doc/scilabisnotnaive/simpleexperiments.tex
new file mode 100644
index 0000000..bea9803
--- /dev/null
+++ b/scilab_doc/scilabisnotnaive/simpleexperiments.tex
@@ -0,0 +1,149 @@
+\section{Simple experiments}
+
+In this section, we analyse the examples given in the introduction of this 
+article.
+
+\subsection{Why $0.1$ is rounded}
+
+In this section, we present a brief explanation for the 
+following Scilab session.
+
+\begin{verbatim}
+-->format(25)
+-->x1=0.1
+ x1  =
+    0.1000000000000000055511  
+-->x2 = 1.0-0.9
+ x2  =
+    0.0999999999999999777955  
+-->x1==x2
+ ans  =
+  F  
+\end{verbatim}
+
+In fact, only the 17 first digits $0.100000000000000005$ are 
+significant and the last digits are a artifact of Scilab's 
+displaying system.
+
+The number $0.1$ can be represented as the normalized number 
+$1.0 \times 10^{-1}$. But the binary floating point representation
+of $0.1$ is approximately \cite{WhatEveryComputerScientist} 
+$1.100110011001100110011001... \times 2^{-4}$. As you see, the decimal
+representation is made of a finite number of digits while the 
+binary representation is made of an infinite sequence of 
+digits. Because Scilab computations are based on double precision numbers
+and because that numbers only have 64 bits to represent the number, 
+some \emph{rounding} must be performed.
+
+In our example, it happens that $0.1$ falls between two 
+different binary floating point numbers. After rounding, 
+the binary floating point number is associated with the decimal 
+representation "0.100000000000000005", that is "rounding up" 
+in this case. On the other side, $0.9$ is also not representable 
+as an exact binary floating point number (but 1.0 is exactly represented). 
+It happens that, after the substraction "1.0-0.9", the decimal representation of the 
+result is "0.09999999999999997", which is different from the rounded 
+value of $0.1$.
+
+\subsection{Why $sin(\pi)$ is rounded}
+
+In this section, we present a brief explanation of the following 
+Scilab 5.1 session, where the function sinus is applied to the 
+number $\pi$.
+
+\begin{verbatim}
+-->format(25)
+-->sin(0.0)
+ ans  =
+    0.  
+-->sin(%pi)
+ ans  =
+    0.0000000000000001224647  
+\end{verbatim}
+
+Two kinds of approximations are associated with the previous 
+result
+\begin{itemize}
+\item $\pi=3.1415926...$ is approximated by Scilab 
+as the value returned by $4*atan(1.0)$,
+\item the $sin$ function is approximated by a polynomial.
+\end{itemize}
+
+This article is too short to make a complete presentation 
+of the computation of elementary functions. The interested 
+reader may consider the direct analysis of the Fdlibm library
+as very instructive \cite{fdlibm}.
+The "Elementary Functions" book by Muller \cite{261217}
+is a complete reference on this subject.
+
+In Scilab, the "sin" function is directly performed by a 
+fortran source code (sci\_f\_sin.f) and no additionnal 
+algorithm is performed directly by Scilab.
+At the compiler level, though, the "sin" function is 
+provided by a library which is compiler-dependent. 
+The main structure of the algorithm which computes 
+"in" is probably the following 
+
+\begin{itemize}
+\item scale the input $x$ so that in lies in a restricted
+interval, 
+\item use a polynomial approximation of the local 
+behaviour of "sin" in the neighbourhood of 0, with a guaranteed
+precision.
+\end{itemize}
+
+In the Fdlibm library for example, the scaling interval is 
+$[-\pi/4,\pi/4]$. 
+The polynomial approximation of the sine function has the general form
+
+\begin{eqnarray}
+sin(x) &\approx& x + a_3x^3 + \ldots + a_{2n+1} x^{2n+1}\\
+&\approx & x + x^3 p(x^2)
+\end{eqnarray}
+
+In the Fdlibm library, 6 terms are used.
+
+For the inverse tan "atan" function, which is 
+used to compute an approximated value of $\pi$, the process is the same.
+All these operations are guaranteed with some precision.
+For example, suppose that the functions are guaranteed with 14 significant
+digits. That means that 17-14 + 1 = 3 digits may be rounded in the process.
+In our current example, the value of $sin(\pi)$ is approximated 
+with 17 digits after the point as "0.00000000000000012". That means that
+2 digits have been rounded. 
+
+\subsection{One more step}
+
+In fact, it is possible to reduce the number of 
+significant digits of the sine function to as low as 0 significant digits.
+The mathematical theory is $sin(2^n \pi) = 0$, but that is not true with
+floating point numbers. In the following Scilab session, we 
+
+\begin{verbatim}
+-->for i = 1:5
+-->k=10*i;
+-->n = 2^k;
+-->sin(n*%pi)
+-->end
+ ans  =
+  - 0.0000000000001254038322  
+ ans  =
+  - 0.0000000001284135242063  
+ ans  =
+  - 0.0000001314954487872237  
+ ans  =
+  - 0.0001346513391512239052  
+ ans  =
+  - 0.1374464882277985633419  
+\end{verbatim}
+
+For $sin(2^{50})$, all significant digits are lost. This computation
+may sound \emph{extreme}, but it must be noticed that it is inside the 
+double precision possibility, since $2^{50} \approx 3.10^{15} \ll 10^{308}$.
+The solution may be to use multiple precision numbers, such as in the 
+Gnu Multiple Precision system.
+
+If you know a better algorithm, based on double precision only, 
+which allows to compute accurately such kind of values, the Scilab 
+team will surely be interested to hear from you !
+