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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
|
<html lang="en">
<head>
<title>Setting up your IEEE environment - GNU Scientific Library -- Reference Manual</title>
<meta http-equiv="Content-Type" content="text/html">
<meta name="description" content="GNU Scientific Library -- Reference Manual">
<meta name="generator" content="makeinfo 4.11">
<link title="Top" rel="start" href="index.html#Top">
<link rel="up" href="IEEE-floating_002dpoint-arithmetic.html" title="IEEE floating-point arithmetic">
<link rel="prev" href="Representation-of-floating-point-numbers.html" title="Representation of floating point numbers">
<link rel="next" href="IEEE-References-and-Further-Reading.html" title="IEEE References and Further Reading">
<link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
<!--
Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 The GSL Team.
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.3 or
any later version published by the Free Software Foundation; with the
Invariant Sections being ``GNU General Public License'' and ``Free Software
Needs Free Documentation'', the Front-Cover text being ``A GNU Manual'',
and with the Back-Cover Text being (a) (see below). A copy of the
license is included in the section entitled ``GNU Free Documentation
License''.
(a) The Back-Cover Text is: ``You have the freedom to copy and modify this
GNU Manual.''-->
<meta http-equiv="Content-Style-Type" content="text/css">
<style type="text/css"><!--
pre.display { font-family:inherit }
pre.format { font-family:inherit }
pre.smalldisplay { font-family:inherit; font-size:smaller }
pre.smallformat { font-family:inherit; font-size:smaller }
pre.smallexample { font-size:smaller }
pre.smalllisp { font-size:smaller }
span.sc { font-variant:small-caps }
span.roman { font-family:serif; font-weight:normal; }
span.sansserif { font-family:sans-serif; font-weight:normal; }
--></style>
</head>
<body>
<div class="node">
<p>
<a name="Setting-up-your-IEEE-environment"></a>
Next: <a rel="next" accesskey="n" href="IEEE-References-and-Further-Reading.html">IEEE References and Further Reading</a>,
Previous: <a rel="previous" accesskey="p" href="Representation-of-floating-point-numbers.html">Representation of floating point numbers</a>,
Up: <a rel="up" accesskey="u" href="IEEE-floating_002dpoint-arithmetic.html">IEEE floating-point arithmetic</a>
<hr>
</div>
<h3 class="section">41.2 Setting up your IEEE environment</h3>
<p><a name="index-IEEE-exceptions-2658"></a><a name="index-precision_002c-IEEE-arithmetic-2659"></a><a name="index-rounding-mode-2660"></a><a name="index-arithmetic-exceptions-2661"></a><a name="index-exceptions_002c-IEEE-arithmetic-2662"></a><a name="index-division-by-zero_002c-IEEE-exceptions-2663"></a><a name="index-underflow_002c-IEEE-exceptions-2664"></a><a name="index-overflow_002c-IEEE-exceptions-2665"></a>The IEEE standard defines several <dfn>modes</dfn> for controlling the
behavior of floating point operations. These modes specify the important
properties of computer arithmetic: the direction used for rounding (e.g.
whether numbers should be rounded up, down or to the nearest number),
the rounding precision and how the program should handle arithmetic
exceptions, such as division by zero.
<p>Many of these features can now be controlled via standard functions such
as <code>fpsetround</code>, which should be used whenever they are available.
Unfortunately in the past there has been no universal API for
controlling their behavior—each system has had its own low-level way
of accessing them. To help you write portable programs GSL allows you
to specify modes in a platform-independent way using the environment
variable <code>GSL_IEEE_MODE</code>. The library then takes care of all the
necessary machine-specific initializations for you when you call the
function <code>gsl_ieee_env_setup</code>.
<div class="defun">
— Function: void <b>gsl_ieee_env_setup</b> ()<var><a name="index-gsl_005fieee_005fenv_005fsetup-2666"></a></var><br>
<blockquote><p><a name="index-g_t_0040env_007bGSL_005fIEEE_005fMODE_007d-2667"></a>This function reads the environment variable <code>GSL_IEEE_MODE</code> and
attempts to set up the corresponding specified IEEE modes. The
environment variable should be a list of keywords, separated by
commas, like this,
<pre class="display"> <code>GSL_IEEE_MODE</code> = "<var>keyword</var>,<var>keyword</var>,..."
</pre>
<p class="noindent">where <var>keyword</var> is one of the following mode-names,
<ul>
<li><code>single-precision</code>
<li><code>double-precision</code>
<li><code>extended-precision</code>
<li><code>round-to-nearest</code>
<li><code>round-down</code>
<li><code>round-up</code>
<li><code>round-to-zero</code>
<li><code>mask-all</code>
<li><code>mask-invalid</code>
<li><code>mask-denormalized</code>
<li><code>mask-division-by-zero</code>
<li><code>mask-overflow</code>
<li><code>mask-underflow</code>
<li><code>trap-inexact</code>
<li><code>trap-common</code>
</ul>
<p>If <code>GSL_IEEE_MODE</code> is empty or undefined then the function returns
immediately and no attempt is made to change the system's IEEE
mode. When the modes from <code>GSL_IEEE_MODE</code> are turned on the
function prints a short message showing the new settings to remind you
that the results of the program will be affected.
<p>If the requested modes are not supported by the platform being used then
the function calls the error handler and returns an error code of
<code>GSL_EUNSUP</code>.
<p>When options are specified using this method, the resulting mode is
based on a default setting of the highest available precision (double
precision or extended precision, depending on the platform) in
round-to-nearest mode, with all exceptions enabled apart from the
<span class="sc">inexact</span> exception. The <span class="sc">inexact</span> exception is generated
whenever rounding occurs, so it must generally be disabled in typical
scientific calculations. All other floating-point exceptions are
enabled by default, including underflows and the use of denormalized
numbers, for safety. They can be disabled with the individual
<code>mask-</code> settings or together using <code>mask-all</code>.
<p>The following adjusted combination of modes is convenient for many
purposes,
<pre class="example"> GSL_IEEE_MODE="double-precision,"\
"mask-underflow,"\
"mask-denormalized"
</pre>
<p class="noindent">This choice ignores any errors relating to small numbers (either
denormalized, or underflowing to zero) but traps overflows, division by
zero and invalid operations.
<p>Note that on the x86 series of processors this function sets both the
original x87 mode and the newer <span class="sc">mxcsr</span> mode, which controls SSE
floating-point operations. The SSE floating-point units do not have a
precision-control bit, and always work in double-precision. The
single-precision and extended-precision keywords have no effect in
this case.
</p></blockquote></div>
<p class="noindent">To demonstrate the effects of different rounding modes consider the
following program which computes e, the base of natural
logarithms, by summing a rapidly-decreasing series,
<pre class="example"> e = 1 + 1/2! + 1/3! + 1/4! + ...
= 2.71828182846...
</pre>
<pre class="example"><pre class="verbatim"> #include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_ieee_utils.h>
int
main (void)
{
double x = 1, oldsum = 0, sum = 0;
int i = 0;
gsl_ieee_env_setup (); /* read GSL_IEEE_MODE */
do
{
i++;
oldsum = sum;
sum += x;
x = x / i;
printf ("i=%2d sum=%.18f error=%g\n",
i, sum, sum - M_E);
if (i > 30)
break;
}
while (sum != oldsum);
return 0;
}
</pre></pre>
<p class="noindent">Here are the results of running the program in <code>round-to-nearest</code>
mode. This is the IEEE default so it isn't really necessary to specify
it here,
<pre class="example"> $ GSL_IEEE_MODE="round-to-nearest" ./a.out
i= 1 sum=1.000000000000000000 error=-1.71828
i= 2 sum=2.000000000000000000 error=-0.718282
....
i=18 sum=2.718281828459045535 error=4.44089e-16
i=19 sum=2.718281828459045535 error=4.44089e-16
</pre>
<p class="noindent">After nineteen terms the sum converges to within <!-- {$4 \times 10^{-16}$} -->
4 \times 10^-16 of the correct value.
If we now change the rounding mode to
<code>round-down</code> the final result is less accurate,
<pre class="example"> $ GSL_IEEE_MODE="round-down" ./a.out
i= 1 sum=1.000000000000000000 error=-1.71828
....
i=19 sum=2.718281828459041094 error=-3.9968e-15
</pre>
<p class="noindent">The result is about
<!-- {$4 \times 10^{-15}$} -->
4 \times 10^-15
below the correct value, an order of magnitude worse than the result
obtained in the <code>round-to-nearest</code> mode.
<p>If we change to rounding mode to <code>round-up</code> then the final result
is higher than the correct value (when we add each term to the sum the
final result is always rounded up, which increases the sum by at least
one tick until the added term underflows to zero). To avoid this
problem we would need to use a safer converge criterion, such as
<code>while (fabs(sum - oldsum) > epsilon)</code>, with a suitably chosen
value of epsilon.
<p>Finally we can see the effect of computing the sum using
single-precision rounding, in the default <code>round-to-nearest</code>
mode. In this case the program thinks it is still using double precision
numbers but the CPU rounds the result of each floating point operation
to single-precision accuracy. This simulates the effect of writing the
program using single-precision <code>float</code> variables instead of
<code>double</code> variables. The iteration stops after about half the number
of iterations and the final result is much less accurate,
<pre class="example"> $ GSL_IEEE_MODE="single-precision" ./a.out
....
i=12 sum=2.718281984329223633 error=1.5587e-07
</pre>
<p class="noindent">with an error of
<!-- {$O(10^{-7})$} -->
O(10^-7), which corresponds to single
precision accuracy (about 1 part in 10^7). Continuing the
iterations further does not decrease the error because all the
subsequent results are rounded to the same value.
<hr>The GNU Scientific Library - a free numerical library licensed under the GNU GPL<br>Back to the <a href="/software/gsl/">GNU Scientific Library Homepage</a></body></html>
|
