summaryrefslogtreecommitdiffstats log msg author committer range
blob: bea98038851a5b5d27ea3c0d476e3360fa9c6e1d (plain)
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 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 !