File: rounding.htm

package info (click to toggle)
boost 1.33.1-10
  • links: PTS
  • area: main
  • in suites: etch, etch-m68k
  • size: 100,948 kB
  • ctags: 145,103
  • sloc: cpp: 573,492; xml: 49,055; python: 15,626; ansic: 13,588; sh: 2,099; yacc: 858; makefile: 660; perl: 427; lex: 111; csh: 6
file content (567 lines) | stat: -rw-r--r-- 26,090 bytes parent folder | download | duplicates (2)
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
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
    "http://www.w3.org/TR/html4/loose.dtd">
<html>
<head>
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <link rel="stylesheet" type="text/css" href="../../../../boost.css">
  <title>Rounding Policies</title>
</head>

<body lang="en">
<h1>Rounding Policies</h1>

<p>In order to be as general as possible, the library uses a class to compute
all the necessary functions rounded upward or downward. This class is the
first parameter of <code>policies</code>, it is also the type named
<code>rounding</code> in the policy definition of <code>interval</code>.</p>

<p>By default, it is <code>interval_lib::rounded_math&lt;T&gt;</code>. The
class <code>interval_lib::rounded_math</code> is already specialized for the
standard floating types (<code>float</code> , <code>double</code> and
<code>long double</code>). So if the base type of your intervals is not one
of these, a good solution would probably be to provide a specialization of
this class. But if the default specialization of
<code>rounded_math&lt;T&gt;</code> for <code>float</code>,
<code>double</code>, or <code>long double</code> is not what you seek, or you
do not want to specialize <code>interval_lib::rounded_math&lt;T&gt;</code>
(say because you prefer to work in your own namespace) you can also define
your own rounding policy and pass it directly to
<code>interval_lib::policies</code>.</p>

<h2>Requirements</h2>

<p>Here comes what the class is supposed to provide. The domains are written
next to their respective functions (as you can see, the functions do not have
to worry about invalid values, but they have to handle infinite
arguments).</p>
<pre>/* Rounding requirements */
struct rounding {
  // defaut constructor, destructor
  rounding();
  ~rounding();
  // mathematical operations
  T add_down(T, T); // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T add_up  (T, T); // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T sub_down(T, T); // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T sub_up  (T, T); // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T mul_down(T, T); // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T mul_up  (T, T); // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T div_down(T, T); // [-&#x221e;;+&#x221e;]([-&#x221e;;+&#x221e;]-{0})
  T div_up  (T, T); // [-&#x221e;;+&#x221e;]([-&#x221e;;+&#x221e;]-{0})
  T sqrt_down(T);   // ]0;+&#x221e;]
  T sqrt_up  (T);   // ]0;+&#x221e;]
  T exp_down(T);    // [-&#x221e;;+&#x221e;]
  T exp_up  (T);    // [-&#x221e;;+&#x221e;]
  T log_down(T);    // ]0;+&#x221e;]
  T log_up  (T);    // ]0;+&#x221e;]
  T cos_down(T);    // [0;2&#x3c0;]
  T cos_up  (T);    // [0;2&#x3c0;]
  T tan_down(T);    // ]-&#x3c0;/2;&#x3c0;/2[
  T tan_up  (T);    // ]-&#x3c0;/2;&#x3c0;/2[
  T asin_down(T);   // [-1;1]
  T asin_up  (T);   // [-1;1]
  T acos_down(T);   // [-1;1]
  T acos_up  (T);   // [-1;1]
  T atan_down(T);   // [-&#x221e;;+&#x221e;]
  T atan_up  (T);   // [-&#x221e;;+&#x221e;]
  T sinh_down(T);   // [-&#x221e;;+&#x221e;]
  T sinh_up  (T);   // [-&#x221e;;+&#x221e;]
  T cosh_down(T);   // [-&#x221e;;+&#x221e;]
  T cosh_up  (T);   // [-&#x221e;;+&#x221e;]
  T tanh_down(T);   // [-&#x221e;;+&#x221e;]
  T tanh_up  (T);   // [-&#x221e;;+&#x221e;]
  T asinh_down(T);  // [-&#x221e;;+&#x221e;]
  T asinh_up  (T);  // [-&#x221e;;+&#x221e;]
  T acosh_down(T);  // [1;+&#x221e;]
  T acosh_up  (T);  // [1;+&#x221e;]
  T atanh_down(T);  // [-1;1]
  T atanh_up  (T);  // [-1;1] 
  T median(T, T);   // [-&#x221e;;+&#x221e;][-&#x221e;;+&#x221e;]
  T int_down(T);    // [-&#x221e;;+&#x221e;]
  T int_up  (T);    // [-&#x221e;;+&#x221e;]
  // conversion functions
  T conv_down(U);
  T conv_up  (U);
  // unprotected rounding class
  typedef ... unprotected_rounding;
};</pre>

<p>The constructor and destructor of the rounding class have a very important
semantic requirement: they are responsible for setting and resetting the
rounding modes of the computation on T. For instance, if T is a standard
floating point type and floating point computation is performed according to
the Standard IEEE 754, the constructor can save the current rounding state,
each <code>_up</code> (resp. <code>_down</code>) function will round up
(resp. down), and the destructor will restore the saved rounding state.
Indeed this is the behavior of the default rounding policy.</p>

<p>The meaning of all the mathematical functions up until
<code>atanh_up</code> is clear: each function returns number representable in
the type <code>T</code> which is a lower bound (for<code> _down</code>) or
upper bound (for <code>_up</code>) on the true mathematical result of the
corresponding function. The function <code>median</code> computes the average
of its two arguments rounded to its nearest representable number. The
functions <code>int_down</code> and <code>int_up</code> compute the nearest
integer smaller or bigger than their argument. Finally,
<code>conv_down</code> and <code>conv_up</code> are responsible of the
conversions of values of other types to the base number type: the first one
must round down the value and the second one must round it up.</p>

<p>The type <code>unprotected_rounding</code> allows to remove all controls.
For reasons why one might to do this, see the <a
href="#Protection">protection</a> paragraph below.</p>

<h2>Overview of the provided classes</h2>

<p>A lot of classes are provided. The classes are organized by level. At the
bottom is the class <code>rounding_control</code>. At the next level come
<code>rounded_arith_exact</code>, <code>rounded_arith_std</code> and
<code>rounded_arith_opp</code>. Then there are
<code>rounded_transc_dummy</code>, <code>rounded_transc_exact</code>,
<code>rounded_transc_std</code> and <code>rounded_transc_opp</code>. And
finally are <code>save_state</code> and <code>save_state_nothing</code>. Each
of these classes provide a set of members that are required by the classes of
the next level. For example, a <code>rounded_transc_...</code> class needs
the members of a <code>rounded_arith_...</code> class.</p>

<p>When they exist in two versions <code>_std</code> and <code>_opp</code>,
the first one does switch the rounding mode each time, and the second one
tries to keep it oriented toward plus infinity. The main purpose of the
<code>_opp</code> version is to speed up the computations through the use of
the "opposite trick" (see the <a href="#perf">performance notes</a>). This
version requires the rounding mode to be upward before entering any
computation functions of the class. It guarantees that the rounding mode will
still be upward at the exit of the functions.</p>

<p>Please note that it is really a very bad idea to mix the <code>_opp</code>
version with the <code>_std</code> since they do not have compatible
properties.</p>

<p>There is a third version named <code>_exact</code> which computes the
functions without changing the rounding mode. It is an "exact" version
because it is intended for a base type that produces exact results.</p>

<p>The last version is the <code>_dummy</code> version. It does not do any
computations but still produces compatible results.</p>

<p>Please note that it is possible to use the "exact" version for an inexact
base type, e.g. <code>float</code> or <code>double</code>. In that case, the
inclusion property is no longer guaranteed, but this can be useful to speed
up the computation when the inclusion property is not desired strictly. For
instance, in computer graphics, a small error due to floating-point roundoff
is acceptable as long as an approximate version of the inclusion property
holds.</p>

<p>Here comes what each class defines. Later, when they will be described
more thoroughly, these members will not be repeated. Please come back here in
order to see them. Inheritance is also used to avoid repetitions.</p>
<pre>template &lt;class T&gt;
struct rounding_control
{
  typedef ... rounding_mode;
  void set_rounding_mode(rounding_mode);
  void get_rounding_mode(rounding_mode&amp;);
  void downward ();
  void upward   ();
  void to_nearest();
  T to_int(T);
  T force_rounding(T);
};

template &lt;class T, class Rounding&gt;
struct rounded_arith_... : Rounding
{
  void init();
  T add_down(T, T);
  T add_up  (T, T);
  T sub_down(T, T);
  T sub_up  (T, T);
  T mul_down(T, T);
  T mul_up  (T, T);
  T div_down(T, T);
  T div_up  (T, T);
  T sqrt_down(T);
  T sqrt_up  (T);
  T median(T, T);
  T int_down(T);
  T int_up  (T);
};

template &lt;class T, class Rounding&gt;
struct rounded_transc_... : Rounding
{
  T exp_down(T);
  T exp_up  (T);
  T log_down(T);
  T log_up  (T);
  T cos_down(T);
  T cos_up  (T);
  T tan_down(T);
  T tan_up  (T);
  T asin_down(T);
  T asin_up  (T);
  T acos_down(T);
  T acos_up  (T);
  T atan_down(T);
  T atan_up  (T);
  T sinh_down(T);
  T sinh_up  (T);
  T cosh_down(T);
  T cosh_up  (T);
  T tanh_down(T);
  T tanh_up  (T);
  T asinh_down(T);
  T asinh_up  (T);
  T acosh_down(T);
  T acosh_up  (T);
  T atanh_down(T);
  T atanh_up  (T);
};

template &lt;class Rounding&gt;
struct save_state_... : Rounding
{
  save_state_...();
  ~save_state_...();
  typedef ... unprotected_rounding;
};</pre>

<h2>Synopsis.</h2>
<pre>namespace boost {
namespace numeric {
namespace interval_lib {

<span style="color: #FF0000">/* basic rounding control */</span>
template &lt;class T&gt;  struct rounding_control;

<span style="color: #FF0000">/* arithmetic functions rounding */</span>
template &lt;class T, class Rounding = rounding_control&lt;T&gt; &gt; struct rounded_arith_exact;
template &lt;class T, class Rounding = rounding_control&lt;T&gt; &gt; struct rounded_arith_std;
template &lt;class T, class Rounding = rounding_control&lt;T&gt; &gt; struct rounded_arith_opp;

<span style="color: #FF0000">/* transcendental functions rounding */</span>
template &lt;class T, class Rounding&gt; struct rounded_transc_dummy;
template &lt;class T, class Rounding = rounded_arith_exact&lt;T&gt; &gt; struct rounded_transc_exact;
template &lt;class T, class Rounding = rounded_arith_std&lt;T&gt; &gt; struct rounded_transc_std;
template &lt;class T, class Rounding = rounded_arith_opp&lt;T&gt; &gt; struct rounded_transc_opp;

<span style="color: #FF0000">/* rounding-state-saving classes */</span>
template &lt;class Rounding&gt; struct save_state;
template &lt;class Rounding&gt; struct save_state_nothing;

<span style="color: #FF0000">/* default policy for type T */</span>
template &lt;class T&gt;  struct rounded_math;
template &lt;&gt;  struct rounded_math&lt;float&gt;;
template &lt;&gt;  struct rounded_math&lt;double&gt;;

<span style="color: #FF0000">/* some metaprogramming to convert a protected to unprotected rounding */</span>
template &lt;class I&gt; struct unprotect;

} // namespace interval_lib
} // namespace numeric
} // namespace boost</pre>

<h2>Description of the provided classes</h2>

<p>We now describe each class in the order they appear in the definition of a
rounding policy (this outermost-to-innermost order is the reverse order from
the synopsis).</p>

<h3 id="Protection">Protection control</h3>

<p>Protection refers to the fact that the interval operations will be
surrounded by rounding mode controls. Unprotecting a class means to remove
all the rounding controls. Each rounding policy provides a type
<code>unprotected_rounding</code>. The required type
<code>unprotected_rounding</code> gives another rounding class that enables
to work when nested inside rounding. For example, the first three lines below
should all produce the same result (because the first operation is the
rounding constructor, and the last is its destructor, which take care of
setting the rounding modes); and the last line is allowed to have an
undefined behavior (since no rounding constructor or destructor is ever
called).</p>
<pre>T c; { rounding rnd; c = rnd.add_down(a, b); }
T c; { rounding rnd1; { rounding rnd2; c = rnd2.add_down(a, b); } }
T c; { rounding rnd1; { rounding::unprotected_rounding rnd2; c = rnd2.add_down(a, b); } }
T d; { rounding::unprotected_rounding rnd; d = rnd.add_down(a, b); }</pre>

<p>Naturally <code>rounding::unprotected_rounding</code> may simply be
<code>rounding</code> itself. But it can improve performance if it is a
simplified version with empty constructor and destructor. In order to avoid
undefined behaviors, in the library, an object of type
<code>rounding::unprotected_rounding</code> is guaranteed to be created only
when an object of type <code>rounding</code> is already alive. See the <a
href="#perf">performance notes</a> for some additional details.</p>

<p>The support library defines a metaprogramming class template
<code>unprotect</code> which takes an interval type <code>I</code> and
returns an interval type <code>unprotect&lt;I&gt;::type</code>  where the
rounding policy has been unprotected. Some  information about the types:
<code>interval&lt;T, interval_lib::policies&lt;Rounding, _&gt;
&gt;::traits_type::rounding</code> <b>is</b> the same type as
<code>Rounding</code>, and <code>unprotect&lt;interval&lt;T,
interval_lib::policies&lt;Rounding, _&gt; &gt; &gt;::type</code> <b>is</b>
the same type as <code>interval&lt;T,
interval_lib::policies&lt;Rounding::unprotected, _&gt; &gt;</code>.</p>

<h3>State saving</h3>

<p>First comes <code>save_state</code>. This class is responsible for saving
the current rounding mode and calling init in its constructor, and for
restoring the saved rounding mode in its destructor. This class also defines
the <code>unprotected_rounding</code> type.</p>

<p>If the rounding mode does not require any state-saving or initialization,
<code>save_state_nothing</code> can be used instead of
<code>save_state</code>.</p>

<h3>Transcendental functions</h3>

<p>The classes <code>rounded_transc_exact</code>,
<code>rounded_transc_std</code> and <code>rounded_transc_opp</code> expect
the std namespace to provide the functions exp log cos tan acos asin atan
cosh sinh tanh acosh asinh atanh. For the <code>_std</code> and
<code>_opp</code> versions, all these functions should respect the current
rounding mode fixed by a call to downward or upward.</p>

<p><strong>Please note:</strong> Unfortunately, the latter is rarely the
case. It is the reason why a class <code>rounded_transc_dummy</code>  is
provided which does not depend on the functions from the std namespace. There
is no magic, however. The functions of <code>rounded_transc_dummy</code> do
not compute anything. They only return valid values. For example,
<code>cos_down</code> always returns -1. In this way, we do verify the
inclusion property for the default implementation, even if this has strictly
no value for the user. In order to have useful values, another policy should
be used explicitely, which will most likely lead to a violation of the
inclusion property. In this way, we ensure that the violation is clearly
pointed out to the user who then knows what he stands against. This class
could have been used as the default transcendental rounding class, but it was
decided it would be better for the compilation to fail due to missing
declarations rather than succeed thanks to valid but unusable functions.</p>

<h3>Basic arithmetic functions</h3>

<p>The classes <code>rounded_arith_std</code> and
<code>rounded_arith_opp</code> expect the operators + - * / and the function
<code>std::sqrt</code> to respect the current rounding mode.</p>

<p>The class <code>rounded_arith_exact</code> requires
<code>std::floor</code> and <code>std::ceil</code> to be defined since it can
not rely on <code>to_int</code>.</p>

<h3>Rounding control</h3>

<p>The functions defined by each of the previous classes did not need any
explanation. For example, the behavior of <code>add_down</code> is to compute
the sum of two numbers rounded downward. For <code>rounding_control</code>,
the situation is a bit more complex.</p>

<p>The basic function is <code>force_rounding</code> which returns its
argument correctly rounded accordingly to the current rounding mode if it was
not already the case. This function is necessary to handle delayed rounding.
Indeed, depending on the way the computations are done, the intermediate
results may be internaly stored in a more precise format and it can lead to a
wrong rounding. So the function enforces the rounding. <a
href="#extreg">Here</a> is an example of what happens when the rounding is
not enforced.</p>

<p>The function <code>get_rounding_mode</code> returns the current rounding
mode, <code>set_rounding_mode</code> sets the rounding mode back to a
previous value returned by <code>get_rounding_mode</code>.
<code>downward</code>, <code>upward</code> and <code>to_nearest</code> sets
the rounding mode in one of the three directions. This rounding mode should
be global to all the functions that use the type <code>T</code>. For example,
after a call to <code>downward</code>, <code>force_rounding(x+y)</code> is
expected to return the sum rounded toward -&#x221e;.</p>

<p>The function <code>to_int</code> computes the nearest integer accordingly
to the current rounding mode.</p>

<p>The non-specialized version of <code>rounding_control</code> does not do
anything. The functions for the rounding mode are empty, and
<code>to_int</code> and <code>force_rounding</code> are identity functions.
The <code>pi_</code> constant functions return suitable integers (for
example, <code>pi_up</code> returns <code>T(4)</code>).</p>

<p>The class template <code>rounding_control</code> is specialized for
<code>float</code>, <code>double</code> and <code>long double</code> in order
to best use the floating point unit of the computer.</p>

<h2>Template class <tt>rounded_math</tt></h2>

<p>The default policy (aka <code>rounded_math&lt;T&gt;</code>) is simply
defined as:</p>
<pre>template &lt;class T&gt; struct rounded_math&lt;T&gt; : save_state_nothing&lt;rounded_arith_exact&lt;T&gt; &gt; {};</pre>

<p>and the specializations for <code>float</code>, <code>double</code> and
<code>long double</code> use <code>rounded_arith_opp</code>, as in:</p>
<pre>template &lt;&gt; struct rounded_math&lt;float&gt;       : save_state&lt;rounded_arith_opp&lt;float&gt; &gt;       {};
template &lt;&gt; struct rounded_math&lt;double&gt;      : save_state&lt;rounded_arith_opp&lt;double&gt; &gt;      {};
template &lt;&gt; struct rounded_math&lt;long double&gt; : save_state&lt;rounded_arith_opp&lt;long double&gt; &gt; {};</pre>

<h2 id="perf">Performance Issues</h2>

<p>This paragraph deals mostly with the performance of the library with
intervals using the floating-point unit (FPU) of the computer. Let's consider
the sum of [<i>a</i>,<i>b</i>] and [<i>c</i>,<i>d</i>] as an example. The
result is [<code>down</code>(<i>a</i>+<i>c</i>),
<code>up</code>(<i>b</i>+<i>d</i>)], where <code>down</code> and
<code>up</code> indicate the rounding mode needed.</p>

<h3>Rounding Mode Switch</h3>

<p>If the FPU is able to use a different rounding mode for each operation,
there is no problem. For example, it's the case for the Alpha processor: each
floating-point instruction can specify a different rounding mode. However,
the IEEE-754 Standard does not require such a behavior. So most of the FPUs
only provide some instructions to set the rounding mode for all subsequent
operations. And generally, these instructions need to flush the pipeline of
the FPU.</p>

<p>In this situation, the time needed to sum [<i>a</i>,<i>b</i>] and
[<i>c</i>,<i>d</i>] is far worse than the time needed to calculate
<i>a</i>+<i>b</i> and <i>c</i>+<i>d</i> since the two additions cannot be
parallelized. Consequently, the objective is to diminish the number of
rounding mode switches.</p>

<p>If this library is not used to provide exact computations, but only for
pair arithmetic, the solution is quite simple: do not use rounding. In that
case, doing the sum [<i>a</i>,<i>b</i>] and [<i>c</i>,<i>d</i>] will be as
fast as computing <i>a</i>+<i>b</i> and <i>c</i>+<i>d</i>. Everything is
perfect.</p>

<p>However, if exact computations are required, such a solution is totally
unthinkable. So, are we penniless? No, there is still a trick available.
Indeed, down(<i>a</i>+<i>c</i>) = -up(-<i>a</i>-<i>c</i>) if the unary minus
is an exact operation. It is now possible to calculate the whole sum with the
same rounding mode. Generally, the cost of the mode switching is worse than
the cost of the sign changes.</p>

<p>Let's recapitulate. Before, when doing an addition, there were three
rounding mode switches (down, up and restore). Now, with this little trick,
there are only two switches (up and restore). It is better, but still a
bottleneck when many operations are nested. Indeed, the generated code for
[<i>a</i>,<i>b</i>] + [<i>c</i>,<i>d</i>] + [<i>e</i>,<i>f</i>] will probably
look like:</p>
<pre>up();
t1 = -(-a - c);
t2 = b + d;
restore();
up();
x = -(-t1 - e);
y = t2 + f;
restore();</pre>

<p>If you think it is possible to do much better, you are right. For example,
this is better (and probably optimal):</p>
<pre>up();
x = -(-a - c - e);
y = b + d + f;
restore();</pre>

<p>Such a code will be generated by a compiler if the computations are made
without initialization and restoration of the rounding mode. However, it
would be far too easy if there were no drawback: because the rounding mode is
not restored in the meantime, operations on floating-point numbers must be
prohibited. This method can only be used if all the operations are operations
on intervals (or operations between an interval and a floating point
number).</p>

<h4 id="perfexp">Example</h4>

<p>Here is an example of the Horner scheme to compute the value of a polynom.
The rounding mode switches are disabled for the whole computation.</p>
<pre>// I is an interval class, the polynom is a simple array
template&lt;class I&gt;
I horner(const I&amp; x, const I p[], int n) {

  // save and initialize the rounding mode
  typename I::traits_type::rounding rnd;

  // define the unprotected version of the interval type
  typedef typename boost::numeric::interval_lib::unprotect&lt;I&gt;::type R;
  
  const R&amp; a = x;
  R y = p[n - 1];
  for(int i = n - 2; i &gt;= 0; i--) {
    y = y * a + (const R&amp;)(p[i]);
  }
  return y;

  // restore the rounding mode with the destruction of rnd
}</pre>

<p>Please note that a rounding object is especially created in order to
compensate for the protection loss. Each interval of type I is converted in
an interval of type R before any operations. If this conversion is not done,
the result is still correct, but the interest of this whole optimization has
disappeared. Whenever possible, it is good to convert to <code>const
R&amp;</code> instead of <code>R</code>: indeed, the function could already
be called inside an unprotection block so the types <code>R</code> and
<code>I</code> would be the same interval, no need for a conversion.</p>

<h4>Uninteresting remark</h4>

<p>It was said at the beginning that the Alpha processors can use a specific
rounding mode for each operation. However, due to the instruction format, the
rounding toward plus infinity is not available. Only the rounding toward
minus infinity can be used. So the trick using the change of sign becomes
essential, but there is no need to save and restore the rounding mode on both
sides of an operation.</p>

<h3 id="extreg">Extended Registers</h3>

<p>There is another problem besides the cost of the rounding mode switch.
Some FPUs use extended registers (for example, float computations will be
done with double registers, or double computations with long double
registers). Consequently, many problems can arise.</p>

<p>The first one is due to to the extended precision of the mantissa. The
rounding is also done on this extended precision. And consequently, we still
have down(<i>a</i>+<i>b</i>) = -up(-<i>a</i>-<i>b</i>) in the extended
registers. But back to the standard precision, we now have
down(<i>a</i>+<i>b</i>) &lt; -up(-<i>a</i>-<i>b</i>) instead of an equality.
A solution could be not to use this method. But there still are other
problems, with the comparisons between numbers for example.</p>

<p>Naturally, there is also a problem with the extended precision of the
exponent. To illustrate this problem, let <i>m</i> be the biggest number
before +<i>inf</i>. If we calculate 2*[<i>m</i>,<i>m</i>], the answer should
be [<i>m</i>,<i>inf</i>]. But due to the extended registers, the FPU will
first store [<i>2m</i>,<i>2m</i>] and then convert it to
[<i>inf</i>,<i>inf</i>] at the end of the calculus (when the rounding mode is
toward +<i>inf</i>). So the answer is no more accurate.</p>

<p>There is only one solution: to force the FPU to convert the extended
values back to standard precision after each operation. Some FPUs provide an
instruction able to do this conversion (for example the PowerPC processors).
But for the FPUs that do not provide it (the x86 processors), the only
solution is to write the values to memory and read them back. Such an
operation is obviously very expensive.</p>

<h2>Some Examples</h2>

<p>Here come several cases:</p>
<ul>
  <li>if you need precise computations with the <code>float</code> or
    <code>double</code> types, use the default
    <code>rounded_math&lt;T&gt;</code>;</li>
  <li>for fast wide intervals without any rounding nor precision, use
    <code>save_state_nothing&lt;rounded_transc_exact&lt;T&gt;
  &gt;</code>;</li>
  <li>for an exact type (like int or rational with a little help for infinite
    and NaN values) without support for transcendental functions, the
    solution could be <code>save_state_nothing&lt;rounded_transc_dummy&lt;T,
    rounded_arith_exact&lt;T&gt; &gt; &gt;</code> or directly
    <code>save_state_nothing&lt;rounded_arith_exact&lt;T&gt; &gt;</code>;</li>
  <li>if it is a more complex case than the previous ones, the best thing is
    probably to directly define your own policy.</li>
</ul>
<hr>

<p>Revised: 2004-02-17<br>
Copyright (c) Guillaume Melquiond, Sylvain Pion, Herv Brnnimann, 2002.
Polytechnic University.<br>
Copyright (c) Guillaume Melquiond, 2004. ENS Lyon.</p>
</body>
</html>