summaryrefslogtreecommitdiffstats
path: root/scilab_doc/scilabisnotnaive/simpleexperiments.tex
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 !