8
votes

I want to change the rounding mode for floating point operations in MATLAB. According to IEEE 754-2008, there are 5 strategies for rounding:

  • round to nearest, ties to even
  • round to nearest, ties away from zero
  • round toward zero
  • round up (toward positive infinity)
  • round down (toward negative infinity)

Does MATLAB supports these 5 strategies? How to change the rounding mode for floating point operations in MATLAB?

2
I don't think this is possible, but I might be wrong about that. The documentation doesn't even mention which of the modes is used.Cris Luengo
@CrisLuengo Yes, it's very strange that as a numerical computing software, MATLAB doesn't tell users about the rounding modes it supports.HYF
MATLAB does not support these 5 strategies. MATLAB simply use the strategy which is set in your processor (FPU). There are ways to change it but you'll have to use something closer to the metal, like C code (don't know if you can do it through mex). Inspiration here and here.Hoki
@Hoki Thanks for your reply, I have already learned to change rounding modes in C language. Do you know any numerical computing softwares that support different rounding modes?HYF
What exactly do you mean by the rounding mode for floating-point operations? Which operations? Can you show an example?nekomatic

2 Answers

7
votes

Answer

Kind of. There is an undocumented feature('setround') function call that you can use to get or set the rounding mode used by Matlab.

So, it can be done, but you shouldn’t do it. :)

WARNING: This is an undocumented, unsupported feature! Use at your own peril!

This feature('setround') supports 4 of the 5 IEEE-754 rounding modes: there’s only one “nearest” mode, and I don't know if it’s “ties to even” or “ties away from zero”.

Supported modes:

  • feature('setround') – Get current rounding mode
  • feature('setround', 0.5) – Round toward nearest (don’t know if it’s ties to even or away from zero)
  • feature('setround', Inf) – Round up (towards +Inf)
  • feature('setround', 0) – Round toward zero
  • feature('setround', -Inf) – Round down (towards -Inf)

Note on testing: The IEEE-754 rounding mode does not affect round() and its relatives. Rather, it governs how arithmetic operations behave around the limits of floating-point precision.

Demonstration

%ROUNDINGEXAMPLE Demonstrates IEEE-754 Rounding Mode control
%
% This uses a completely undocumented and unsupported feature!
% Not for production use!

%% Setup
clear; clc

n = 2000;
X = ones(n)*1E-30; % matrix with n^2 elements
defaultRoundingMode = feature('setround'); % store default rounding mode

%%
feature('setround',0.5);
r1 = prettyPrint('Nearest', sum(X(:)));
%{
  sign   exponent                       mantissa
     0 01110110001 0011010101111100001010011001101001110101010000011110
     | \_________/ \__________________________________________________/
     |      |             ______________________|___________________________
     |      |            /                                                  \
(-1)^0 2^( 945 - 1023) 1.0011010101111100001010011001101001110101010000011110 = 4e-24
%}

%%
feature('setround',-Inf);
r2 = prettyPrint('To -Infinity', sum(X(:)));
%{
  sign   exponent                       mantissa
     0 01110110001 0011010101111100001010011001101001011100000111000110
     | \_________/ \__________________________________________________/
     |      |             ______________________|___________________________
     |      |            /                                                  \
(-1)^0 2^( 945 - 1023) 1.0011010101111100001010011001101001011100000111000110 = 4e-24
%}

%%
feature('setround',Inf);
r3 = prettyPrint('To Infinity', sum(X(:)));
%{
  sign   exponent                       mantissa
     0 01110110001 0011010101111100001010011001101010100011101100100001
     | \_________/ \__________________________________________________/
     |      |             ______________________|___________________________
     |      |            /                                                  \
(-1)^0 2^( 945 - 1023) 1.0011010101111100001010011001101010100011101100100001 = 4e-24
%}

%%
feature('setround',0);
r4 = prettyPrint('To zero', sum(X(:)));
%{
  sign   exponent                       mantissa
     0 01110110001 0011010101111100001010011001101001011100000111000110
     | \_________/ \__________________________________________________/
     |      |             ______________________|___________________________
     |      |            /                                                  \
(-1)^0 2^( 945 - 1023) 1.0011010101111100001010011001101001011100000111000110 = 4e-24
%}

%%
feature('setround',defaultRoundingMode);
r5 = prettyPrint('No accumulated roundoff error', 4e-24);
%{
  sign   exponent                       mantissa
     0 01110110001 0011010101111100001010011001101010001000111010100111
     | \_________/ \__________________________________________________/
     |      |             ______________________|___________________________
     |      |            /                                                  \
(-1)^0 2^( 945 - 1023) 1.0011010101111100001010011001101010001000111010100111 = 4e-24
%}

%% Helper function
function r = prettyPrint(s, r)
    fprintf('%s:\n%65.60f\n\n', s, r); 
end

I get:

Nearest:
   0.000000000000000000000003999999999966490758963870373537264729

To -Infinity:
   0.000000000000000000000003999999999789077070014108839608005726

To Infinity:
   0.000000000000000000000004000000000118618095059505975310731249

To zero:
   0.000000000000000000000003999999999789077070014108839608005726

No accumulated roundoff error:
   0.000000000000000000000003999999999999999694801998206811298525

Acknowledgments

Thanks to Ryan Klots at MathWorks Technical Support for setting me straight on this and providing the nice demo code!

0
votes

@HYF: I found that feature('setround', 0.5) leads to round to even. I checked the following: a=1+2^-52 This means that mantissa looks like so: 1.0...01 where the last 1 is still in the mantissa. The leading 1 is not stored in IIIE754 format. ( I checked that 1+2^-52 == 1 but not 1+2^-53 == 1)

Then I computed b = a + 2^-53. Without rounding it is 1.0...01|1 where the last digit is to be rounded away. I found the following true: b==1+2^-51 We have b == 1.0...010.

We have several rounding modes submodes of round to nearest: This can be round to inf, round away from 0 or round to even.

Next we check -b==-1-2^-51 to be true which excludes round to inf still allowing round away from 0 or round to even.

Then I checked 1.5==1.5+2^-53.
Of course 1.5 = 1.10...0 binary and and 1.5+2^-53 = 1.10...0|1 without rounding, the last digit to be rounded away. Rounding away from 0 would be 1.10...01 and rounding to even is 1.10...0. So the latter is the case.