Diff of /src/fesetprec.m [000000] .. [aea747]  Maximize  Restore

  Switch to unified view

a b/src/fesetprec.m
1
## Copyright (C) 2008 Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
2
##
3
## This program is free software; you can redistribute it and/or modify
4
## it under the terms of the GNU General Public License as published by
5
## the Free Software Foundation; either version 3 of the License, or
6
## (at your option) any later version.
7
##
8
## This program is distributed in the hope that it will be useful,
9
## but WITHOUT ANY WARRANTY; without even the implied warranty of
10
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
11
## GNU General Public License for more details.
12
##
13
## You should have received a copy of the GNU General Public License
14
## along with this program; if not, see <http://www.gnu.org/licenses/>.
15
##
16
17
## -*- texinfo -*-
18
## @deftypefn {Function File} {} fesetprec (@var{mode})
19
## Change the precision of the floating point arithmetics. At present, @code{fesetprec} may work @strong{only} on systems with x86 or x86_64 processors. 
20
##
21
## Possible values of @var{mode} are 
22
## @table @samp
23
## @item 1 
24
## @itemx 24
25
## Sets the x87 FPU precision mode to single precision (24-bit significand)
26
## @item 2 
27
## @itemx 53
28
## Sets the x87 FPU precision mode to double precision (53-bit significand)
29
## @item 3 
30
## @itemx 64
31
## Sets the x87 FPU precision mode to double extended precision (64-bit significand)
32
## @end table
33
## 
34
## When successful, @code{fesetprec} returns 0. It is always recommended to verify the result experimentally by calling the function @code{fetestenv}, which tests the rounding mode and the precision of the floating point arithmetics. 
35
##
36
## Modern processors, including x86 and x86_64 architectures, support changing certain properties of their floating point arithmetics. Here, we focus @strong{exclusively} on x86 and x86_64.
37
## These processors presently feature two types of floating point processing units: the older x87 FPU and the newer SSE unit. Only x87 FPU allows changing the precision, so @code{fesetprec} applies only to those instructions which are executed by x87 FPU. @code{fesetprec} changes the precision of the four basic arithmetic operations and of the square root on x87 FPU. 
38
## 
39
## Since floating point operations on x86 or x86_64 may be executed either by x87 FPU or SSE units, @code{fesetprec} may affect all, some, or no Octave's operations. The rule of the thumb, however, is that 64-bit systems usually do their floating point on SSE by default, when 32-bit systems use x87 FPU. You can always check how your system reacts to changes made by @code{fesetprec} by verifying them with @code{fetestenv}.
40
## 
41
## It is important to know that all expressions and their intermediate results processed by Octave's interpreter are cast to temporary @code{double} variables and thus it is impossible to execute an interpreted expression in double extended precision. On the other hand, those libraries called by Octave which rely on x87 FPU will usually perform the computations using x87 FPU 80-bit registers, i.e. in double extended precision. This may lead to some subtle differences in how the precision mode affects the computations. Let's take a look at the simplest example:
42
## 
43
## @example
44
## macheps = builtin('eps');
45
## fesetprec(3);
46
## x = [1, macheps/2, -1]*[1;1;1];
47
## y = 1 + macheps/2 - 1;
48
## @end example 
49
## 
50
## The built-in function @code{eps} returns the machine epsilon in @emph{double precision}. Mathematically the expressions for @code{x} and @code{y} are equivalent, but in reality they may be not. The expression for @code{x} will call a BLAS subroutine and on a 32-bit system it will most likely be executed on x87 FPU using 64-bit double extended precision - so @code{x} will be greater than zero. On the other hand, the expression for @code{y} will be executed through the interpreter and these calculations will store all intermediate results in double precision, so that @code{y} will be equal to zero. In contrast to 32-bit systems, many 64-bit systems will do all floating point on SSE - and thus both @code{x} and @code{y} will then be equal to zero.
51
##
52
## It is also known that some numerical libraries set their own precision modes or execute floating point operations on different units (say, they may use only the SSE unit when Octave uses exclusively the x87 FPU...). Therefore calls to such numerical libraries may also cancel changes made previously by @code{fesetprec}:
53
##
54
## @example
55
## fesetprec(1); fesetround(0); fetestenv
56
## @print{} Machine epsilon: 1.19209e-07
57
## @print{} Rounding mode:   0
58
## hilb(2)\ones(2,1); fetestenv
59
## @print{} Machine epsilon: 2.22045e-16
60
## @print{} Rounding mode:   0.5
61
## @end example
62
##
63
## (the above test has been run on an x86 32-bit system using Octave 3.0.2 package provided along with the Fedora 9 Linux distribution).
64
##
65
## Changes issued by @code{fesetprec} @strong{are} reflected by Octave's function @code{eps}: calling @code{fesetprec} replaces the built-in function @code{eps} with a dispatch @code{fetestenv}. In order to perform a more detailed check of actual parameters of the floating point arithmetics, use @code{fetestenv} directly.
66
## 
67
## A possible application of @code{fesetprec}, following the idea by W. Kahan, is an experimental (in)stability analysis. Running some code with various settings of the floating point arithmetics may reveal some instability of the numerical algorithm being used in the code or reveal possible ill conditioning of the very problem being solved.
68
## 
69
## Literature
70
## 
71
## @enumerate 
72
## @item W. Kahan, "How Futile are Mindless Assessments of Roundoff in Floating-Point Computation?", available on the Web: @indicateurl{http://www.cs.berkeley.edu/~wkahan/Mindless.pdf}.
73
## @item W. Kahan, "Why I can Debug some Numerical Programs You can't", available on the Web: @indicateurl{http://www.cs.berkeley.edu/~wkahan/Stnfrd50.pdf}.
74
## @item Intel 64 and IA-32 Architectures Software Developer's Manual, Volume 1: Basic Architecture, May 2007.
75
## @end enumerate
76
## @seealso{fetestenv, fesetround, system_dependent, eps, fe_system_dependent}
77
## @end deftypefn
78
79
## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
80
##
81
82
## Revision history: 
83
##
84
##    2008-10-10, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
85
##        Initial release
86
##
87
88
function err = fesetprec(mode)
89
if (nargin ~= 1)
90
  print_usage();
91
end
92
93
systemtype = computer();
94
95
# Presently, we support either x86_64...
96
supported = ~isempty(findstr('x86_64', systemtype));
97
98
# ...or x86
99
#     (Actually we only support i486 and above, 
100
#     but in many cases the default compiler settings are for i386.
101
#     We also do believe i286 and older are dead or at the museum.)
102
supported |= strcmp("i86", systemtype([1,3,4]));
103
104
err = -2;
105
if (supported && exist('system_dependent'))
106
  switch mode
107
  case {1, 24}
108
      err = fe_system_dependent('precision','single');
109
  case {2, 53}
110
      err = fe_system_dependent('precision','double');
111
  case {3, 64}
112
      err = fe_system_dependent('precision','double ext');
113
  otherwise
114
      warning('Cannot set the precision mode to %g', mode);
115
      usage('fesetprec (MODE), with MODE in {1, 2, 3}');
116
  end
117
  if (err == 0)
118
      % replace original eps function with fetestenv
119
      dispatch('eps', 'any'); % restore the built-in version
120
      dispatch('eps', 'fetestenv', 'any'); % obscure by fetestenv
121
  else 
122
      warning('Problems setting precision mode %g.', mode);
123
      warning('Use fetestenv() to verify.');
124
  end
125
else
126
  warning('This system is not supported by fesetprec');
127
end
128
end
129
%!demo
130
%!  fesetprec(1);
131
%!  fetestenv();
132
%!  fesetprec(3);
133
%!  fetestenv();

Get latest updates about Open Source Projects, Conferences and News.

Sign up for the SourceForge newsletter:





No, thanks