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

  Switch to unified view

a b/src/fetestenv.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} {} [@var{macheps}, @var{round}, @var{machepsb}, @var{roundb}] = fetestenv ()
19
## Check some properties of floating point arithmetics: 
20
## the rounding mode and the machine epsilon, as seen by
21
## Octave's interpreted code (@var{round} and @var{macheps}, respectively) and, if requested, as seen by BLAS (@var{roundb} and @var{machepsb}, respectively). 
22
##
23
## The function performs a variant of the standard 
24
## @display
25
## 1 + @var{macheps} ~= 1
26
## @end display
27
## test in order to check the actual machine epsilon @var{macheps} 
28
## and the present rounding mode @var{round} used in all interpreted 
29
## expressions in Octave.
30
##
31
## As Octave interpreted code cannot achieve more than double precision,
32
## and some BLAS implementations may impose their internal precision and rounding modes,
33
## the user may request to check the actual rounding and precision used in BLAS calls.
34
## To this end, we perform a machine epsilon/rounding mode check of the form 
35
## @display
36
## [1, @var{machepsb}, -1]*[1; 1; 1] ~= 0.
37
## @end display
38
## This Octave code results in a call to a BLAS routine, which in turn may use different floating point arithmetic settings than
39
## available to Octave. 
40
## Rounding mode and machine epsilon (both as seen by BLAS) are returned 
41
## in @var{roundb} and @var{machepsb}, respectively.
42
##
43
## See the manual pages of @command{fesetprec} for more details and subtleties.
44
##
45
## Note that the results of this test may or may not apply to library functions called by Octave. The simplest way to verify if the function you use in Octave is affected by the change of the rounding mode or the precision, is to run your function on the same data with different settings. Small differences between various modes indicate these modes do affect the behaviour of your function. No differences mean they probably do not apply. Big differences most likely indicate either an instability of your function or ill conditioning of your problem.
46
##
47
##
48
## @seealso{fesetround, fesetprec, system_dependent, eps, fe_system_dependent}
49
## @end deftypefn
50
51
## Author: Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
52
##
53
54
## Revision history: 
55
##
56
##    2008-10-10, Piotr Krzyzanowski <piotr.krzyzanowski@mimuw.edu.pl>
57
##        Initial release
58
##
59
60
function [macheps, round, machepsb, roundb] = fetestenv()
61
if (nargin != 0)
62
  warning ("Ignoring extra arguments");
63
end
64
65
# first run the standard "1 + eps == 1" test, modified to accommodate 
66
# various rounding modes
67
x     = [  1,  1, -1, -1];
68
xeps  = [  1, -1,  1, -1];
69
do 
70
  xeps = xeps ./ 2;
71
  stop = [(x + xeps > x)([1,3]), (x + xeps < x)([2,4])];
72
until ( ~prod(stop) )
73
74
# machine eps is thus a twice as big number
75
macheps = xeps(1)*2;
76
77
# use an even smaller number to reveal rounding direction
78
xeps = xeps ./ 2;
79
stop = [(x + xeps > x)([1,3]), (x + xeps < x)([2,4])];
80
81
round = sense_round(stop);
82
83
if(nargout != 2)
84
# Check macheps and rounding as in BLAS' routine
85
# (Inspired by some discussion on the Octave's mailing list)
86
87
  x   = [  1,  1, -1, -1];
88
  xeps    = [  1, -1,  1, -1];
89
  do 
90
      xeps = xeps ./ 2;
91
      stop = [([1,1,1]*[x; xeps; -x]>0)([1,3]), ...
92
          ([1,1,1]*[x; xeps; -x]<0)([2,4])];
93
  until ( ~prod(stop) )
94
95
  # machine eps is thus a twice as big number
96
  machepsb = xeps(1)*2;
97
98
  # use an even smaller number to reveal rounding direction
99
  xeps = xeps ./ 2;
100
  stop = [([1,1,1]*[x; xeps; -x]>0)([1,3]), ...
101
      ([1,1,1]*[x; xeps; -x]<0)([2,4])];
102
103
  roundb = sense_round(stop);
104
end
105
106
if(nargout == 0)
107
  printf('Interpreter:\n');
108
  printf('\tMachine epsilon: %g\n', macheps); 
109
  printf('\tRounding mode:   %g\n', round);   
110
  printf('BLAS:\n');
111
  printf('\tMachine epsilon: %g\n', machepsb);    
112
  printf('\tRounding mode:   %g\n', roundb);  
113
end
114
115
endfunction
116
117
function rnd = sense_round(stop)
118
# tries to sniff the rounding mode from the sign checks 
119
# of the "1+meps-1" type expressions provided in 'stop' argument
120
switch stop
121
case [0, 0, 0, 0]
122
  # all additions/subtractions collapsed to +1 or -1; we are in the
123
  # round-to-nearest mode
124
  rnd = 0.5;
125
case [1, 1, 0, 0]
126
  # round-up mode
127
  rnd = +Inf;
128
case [0, 0, 1, 1]
129
  # round-down mode
130
  rnd = -Inf;
131
case [0, 1, 1, 0]
132
  # round-to-zero mode
133
  rnd = 0;
134
otherwise
135
  # am I dreaming?
136
  rnd = NaN;
137
  error('Unknown rounding mode detected.');
138
end
139
endfunction
140
%!demo
141
%!  fesetround(+Inf); fesetprec(1);
142
%!  fetestenv();
143
%!  fesetround(0.5); fesetprec(3);
144
%!demo
145
%!  fesetround(0); fesetprec(3);
146
%!  fetestenv();
147
%!  fesetround(0.5); fesetprec(3);
148
%!test 
149
%!  if( fesetround(+Inf) == 0) 
150
%!    [eps, round] = fetestenv();
151
%!    assert(round, Inf);
152
%!    fesetround(0.5);
153
%!  endif
154
%!test 
155
%!  if( fesetround(0) == 0) 
156
%!    [eps, round] = fetestenv();
157
%!    assert(round, 0);
158
%!    fesetround(0.5);
159
%!  endif