bugGNU Octave - Bugs: bug #55682, round(X,N) and round(X,N,type)

 
 

bug #55682: round(X,N) and round(X,N,type)

Submitter:  Joshua Gay <josh>
Submitted:  Sat 09 Feb 2019 01:30:37 PM UTC
   
 
Category:  Octave Function Severity:  1 - Wish
Priority:  3 - Low Item Group:  Feature Request
Status:  Confirmed Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 10 Jun 2024 06:45:25 PM UTC, comment #13: 

In response to comment #12, things are actually a bit easier.  There is a distinction between pure math and numerical computation.  And that distinction can be enormously subtle and hard to navigate.  But, round() is a function provided by Matlab and Octave seeks to be compatible with Matlab.  So, we actually have an easy way to answer the question of what "round (0.15, 1)" should return which is that it should return the value that the reference implementation Matlab returns.  In this case, that value is 0.2.  Similarly, Both Octave and Matlab rely on the C++ standard library for some reference functions such as realpow().


10 .^ -5 == 1e-5


returns false in Matlab, but true in Octave on my machine.  So, it does depend on the particular version of C++ libraries and probably compilation options as well since you report that this returned false for you.

Rik <rik5>
Group administrator
Mon 10 Jun 2024 05:26:02 PM UTC, comment #12: 

There are some subtleties to the definition of round(X, N) relating to the fact that many of the theoretical decimal rounding targets are not exactly representable in binary floating point.  Consider:

x = 0.15;   % has actual value 0.1499999999999999944488848768742172978818416595458984375
round(x, 1) % what should this return?

The theoretical rounding targets of interest are 0.1 and 0.2, and their corresponding nearest binary64 values are:

0.1000000000000000055511151231257827021181583404541015625
0.200000000000000011102230246251565404236316680908203125

Here, x would ideally round down regardless of whether we consider the theoretical or the binary64 rounding targets.  Unfortunately the implementation contemplated in earlier comments seems likely to round up in this case.  I am not sure whether there are inputs for which the theoretical and the binary64 rounding targets would imply different results.

Another consideration is that Octave does not specify the accuracy of the exponentiation operator or the realpow function.  They do not seem to be correctly rounded on my system.  For instance, `10 .^ -5 == 1e-5` is false.  Perhaps things are different on the C++ side; I am not familiar with that.

Underflow and overflow could also cause problems.

Even in cases where the theoretical rounding targets are exactly representable, setting aside issues of underflow and overflow, and assuming exponentiation is correctly rounded, there is still potential for trouble:  Computing `round(12345678.9, -5)` as `round(12345678.9 * 1e-5) / 1e-5` yields 12299999.99999999813735485076904296875 instead of 123e5.  The following might be helpful if exponentiation can be relied upon to return an exact result when it is representable:

   if (n > 0)
      retval = round(x .* 10 .^ n) ./ (10 .^ n);
   elseif (n < 0)
      retval = round(x ./ 10 .^ -n) .* (10 .^ -n);
   else
      retval = round(x);
   endif

I'm not certain that this never makes things worse though.

Anonymous
Sat 23 Dec 2023 09:41:58 PM UTC, comment #11: 

Came across this when looking to see if a round(X,N) report had been found. Didn't realize how deep this issue was.

In any case. Just updating with the detail that since the last work on this report MATLAB also made the following tweaks to the function in 2022a:


R2022a: Control tiebreak behavior
Specify how to break ties by using the TieBreaker name-value argument. For example, round(X,TieBreaker="tozero") rounds ties towards zero.

R2022a: round returns consistent results for ties
Starting in R2022a, the round function always rounds ties away from zero to the nearest multiple of 10—N with larger magnitude by default. For example:

X = 1.015:5.015;
N = 2;
Y = round(1.015:5.015,2)
Y =
    1.0200    2.0200    3.0200    4.0200    5.0200
In previous releases, the round function sometimes returned inconsistent results for ties by default. In the previous example, for instance, the second and third elements were rounded towards zero to 2.01 and 3.01, respectively.


Nicholas Jankowski <nrjank>
Group Member
Sat 25 May 2019 11:59:41 PM UTC, comment #10: 

Sigh

It is very clean that way, but I had hoped to limit the re-writing of the API.  This will require propagating a two-input round form up through libinterp and back down in to liboctave.  Presumably we should also maintain a legacy zero-input form, or arrange for defaults in the function prototype for the 2-input form such that it can be called with no arguments.

Can't use ceil() in place of round because it always rounds in the direction of positive infinity.


ceil (-1.5) == -1
round (-1.5) == -2



Rik <rik5>
Group administrator
Sat 25 May 2019 12:07:01 PM UTC, comment #9: 

I would try to limit the code in the Fround function to just doing argument parsing.  Something like:


{
  int nargin = args.length ();

  if (nargin < 1 || nargin > 3)
    print_usage ();

  if (nargin == 1)
    return ovl (args(0).ceil ());

  int n = args(1).xint_value ("N must be an integer value");

  enum octave::rounding_mode = octave::round_normal;
  if (nargin == 3)
    {
      std::string mode = xstring_value ("MODE must be a string");

      // Case insensitive compare?  Abbreviations?
      if (mode == "significant")
        mode = octave::round_significant;
      else if (mode == "decimals")
        mode = octave::round_decimals;
      else
        error ("round: MODE must be either \"decimals\" or \"significant\"");
    }

  return ovl (args(0).round (n, mode);
}


This change will require the definition of a new octave::rounding_mode enum and new two-argument round functions in octave_base_value and any class that needs them.  Using template functions at that point may help avoid some duplication.

Doing it this way eliminates the explicit checks on data types, which is something I would like to avoid.


John W. Eaton <jwe>
Group administrator
Fri 24 May 2019 05:24:21 PM UTC, comment #8: 

The function round is in mappers.cc.  I rewrote the first bit to


  if (nargin == 1)
    return ovl (args(0).round ());


which exactly preserves the functionality we have today where the round() function is a method on an octave_value object.

I then tried the following which works, but seems kind of slow because I am not doing an in-place map operation on the object x.


octave_value_list retval (1);

if (args(0).is_double_type ())
  {
    if (args(0).isreal ())
      {
        Matrix x = args(0).matrix_value ();
        x *= pow (10.0, N);
        x = x.map<double> (octave::math::round);
        x /= pow (10.0, N);
        retval(0) = octave_value (x);


I can work around this by writing my own static function in mappers.cc.


template <typename T>
void
do_map_round (Array<T>& x)
{
  for (octave_idx_type i = 0; i < x.numel (); i++)
    x.xelem (i) = octave::math::round (x.xelem (i));
}


followed by


octave_value_list retval (1);

if (args(0).is_double_type ())
  {
    if (args(0).isreal ())
      {
        NDArray x = args(0).array_value ();
        x *= pow (10.0, N);
        do_map_round (retval);
        x /= pow (10.0, N);
        retval(0) = octave_value (x);


Of course, I will have to write at least 4 if branches (double+real, double+complex, single+real, single+complex) and maybe another two branches for (sparse+real, sparse+complex).  Is there some other, more natural way, to tackle this?

A diff for the changes so far is attached.


Rik <rik5>
Group administrator
Fri 24 May 2019 05:23:15 PM UTC, comment #7: 

Seems like some stuff got truncated.

The function round is in mappers.cc.  I rewrote the first bit to


  if (nargin == 1)
    return ovl (args(0).round ());


which exactly preserves the functionality we have today where the round() function is a method on an octave_value object.

I then tried the following which works, but seems kind of slow because I am not doing an in-place map operation on the object x.


octave_value_list retval (1);

if (args(0).is_double_type ())
  {
    if (args(0).isreal ())
      {
        Matrix x = args(0).matrix_value ();
        x *= pow (10.0, N);
        x = x.map<double> (octave::math::round);
        x /= pow (10.0, N);
        retval(0) = octave_value (x);
verbatim-

I can work around this by writing my own static function in mappers.cc.

+verbatim+
template <typename T>
void
do_map_round (Array<T>& x)
{
  for (octave_idx_type i = 0; i < x.numel (); i++)
    x.xelem (i) = octave::math::round (x.xelem (i));
}


followed by


octave_value_list retval (1);

if (args(0).is_double_type ())
  {
    if (args(0).isreal ())
      {
        NDArray x = args(0).array_value ();
        x *= pow (10.0, N);
        do_map_round (retval);
        x /= pow (10.0, N);
        retval(0) = octave_value (x);


Of course, I will have to write at least 4 if branches (double+real, double+complex, single+real, single+complex) and maybe another two branches for (sparse+real, sparse+complex).  Is there some other, more natural way, to tackle this?

A diff for the changes so far is attached.


Rik <rik5>
Group administrator
Fri 24 May 2019 05:00:00 PM UTC, comment #6: 

Adding jwe to the CC list for some guidance on coding strategy.

This bug would be easy to resolve by writing a new m-file round.m and renaming the C++ function to something internal like _round_.  In Octave's languagage the operation to accomplish is just


Y=round(X*10^N)/10^N;


However, round is used so frequently by coders coming from Matlab that I think it is important that it be written in C++ for speed.  The input validation, written in an m-file, would be significant.

So I started down that route, but it has gotten ugly.

The function round is in mappers.cc.  I rewrote the first bit to


  if (nargin == 1)
    return ovl (args(0).round ());


which exactly preserves the functionality we have today where the round() function is a method on an octave_value object.

I then tried the following which works, but seems kind of slow because I am not doing an in-place map operation on the object x.


octave_value_list retval (1);

if (args(0).is_double_type ())
  {
    if (args(0).isreal ())
      {
        Matrix x = args(0).matrix_value ();
        x *= pow (10.0, N);
        x = x.map<double> (octave::math::round);
        x /= pow (10.0, N);
        retval(0) = octave_value (x);
verbatim-

I can work around this by writing my own static function in mappers.cc.

+verbatim+
template <typename T>
void
do_map_round (Array<T>& x)
{
  for (octave_idx_type i = 0; i < x.numel (); i++)
    x.xelem (i) = octave::math::round (x.xelem (i));
}


followed by


octave_value_list retval (1);

if (args(0).is_double_type ())
  {
    if (args(0).isreal ())
      {
        NDArray x = args(0).array_value ();
        x *= pow (10.0, N);
        do_map_round (retval);
        x /= pow (10.0, N);
        retval(0) = octave_value (x);


Of course, I will have to write at least 4 if branches (double+real, double+complex, single+real, single+complex) and maybe another two branches for (sparse+real, sparse+complex).  Is there some other, more natural way, to tackle this?

A diff for the changes so far is attached.



(file #46967)

Rik <rik5>
Group administrator
Tue 21 May 2019 05:14:09 PM UTC, comment #5: 

Documentation on Matlab's round function is available here (https://www.mathworks.com/help/matlab/ref/round.html).

They seem to implement the scheme suggested by Michael in comment #2.  For example, if you use the two-input form "round (x,n)" they require that x be a floating point class (single or double).  Similarly, they have an example of the possible confusion that can result when a floating point number looks like it should round in one direction, but doesn't because of extra precision.


format short
x = 112.05 - 110
x = 2.05
y = round (x, 1)
%% Expected value from above would be 2.1
y = 2.0
%% Check significant figures
format long
x
x = 2.049999999999997


So this wouldn't be that hard to implement.

We could alse decide whether we want to be nicer than Matlab and work with integer inputs when n != 0.  This would require converting to doubles, rounding, and then converting back to the original type.  Because of saturation, you could get some weird effects though.


round (uint8 (251), -2)
=> 255



Rik <rik5>
Group administrator
Mon 11 Feb 2019 08:36:09 PM UTC, comment #4: 

It's also compatible with Python's 2-argument round function.


Added in Matlab R2014b.

Mike Miller <mtmiller>
Group Member
Sun 10 Feb 2019 12:27:11 AM UTC, comment #3: 

My guess is that Matlab went this way so that they could grab users of Excel.  See the documentation for the Excel function round (https://support.office.com/en-us/article/ROUND-function-C018C5D8-40FB-4053-90B1-B3E7F61A213C).  It doesn't make it right, but they were catering to peopele who probably don't know the math behind the round function.

Rik <rik5>
Group administrator
Sat 09 Feb 2019 08:48:20 PM UTC, comment #2: 

Thank you.

The downstream bug I filed was to a matlab project that I would like to see compatible with octave. I suggestex to them they could do the same thing (and also said they could use . for multiplication and division.


Joshua Gay <josh>
Sat 09 Feb 2019 08:20:51 PM UTC, comment #1: 

What trivial problems make the guys at Mathworks feel the need to complexify their API... If you only need this to solve a given problem right now, just use


Y=round(X*10^N)/10^N;


But beware that these are then floating-point numbers, that is, if you get let's say 1.3 in this way and 2.7, it is not guaranteed that their sum is equal to a corresponding expression of 4.0 (different from proper rounding, which gives you integer-valued doubles that conform to exact arithmetics as long as you do not exceed bitmax()).

Michael Leitner <mleitner>
Sat 09 Feb 2019 01:30:37 PM UTC, original submission:  

Matlab supports three round functions.


Y = round(X)
Y = round(X,N)
Y = round(X,N,type)



It would be nice if octave could support these as well, although round(X,N) is what would help me the most right now.

I apologize if a feature request has already been filed against this -- I did spend several minutes searching.

Joshua Gay <josh>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #46967:  roundn.diff added by rik5 (3KiB - text/x-patch)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by nrjank (Posted a comment)
  • -email is unavailable- added by hardy
  • -email is unavailable- added by gyom
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by rik5 (Updated the item)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by josh (Submitted the item)
  • -email is unavailable- added by josh
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 11 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2020-12-12 rik5 Dependencies- bugs #59642 is dependent
    2020-11-25 hardy Carbon-Copy- Added hardy
    2019-09-26 gyom Carbon-Copy- Added gyom
    2019-05-24 rik5 Attached File- Added roundn.diff, #46967
        Carbon-Copy- Added jwe
    2019-02-11 mtmiller Severity3 - Normal 1 - Wish
        Priority5 - Normal 3 - Low
    2019-02-10 rik5 Item GroupMatlab Compatibility Feature Request
        StatusNone Confirmed
        Release4.4.0 dev
    2019-02-09 josh Carbon-Copy- Added -email is unavailable-

    Back to the top

    Powered by Savane 3.13-ec27.
    Corresponding source code