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.00.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.00.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 compilerdependent.
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 1714 + 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 !
