bugGNU Octave - Bugs: bug #61319, idivide fails to distinguish...

 
 

bug #61319: idivide fails to distinguish between ceil and floor for inputs close to perfect squares

Submitter:  Arun Giridhar <arungiridhar>
Submitted:  Sat 09 Oct 2021 12:55:58 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 07 Mar 2022 11:55:56 AM UTC, comment #26: 

IIUC, this is working correctly now.

Closing as fixed.

If someone would like to look into speeding up the current implementation, we can re-open or use a new report.

Markus Mützel <mmuetzel>
Group administrator
Sun 27 Feb 2022 12:56:15 PM UTC, comment #25: 

No problem to mention me on the list of contributors (that's my real name), but I don't particularly ask for it (it's a modest contribution). I don't really know what this list is actually, so it's up to you :)

Gaël Bonithon <tamaranch>
Sun 27 Feb 2022 12:39:59 PM UTC, comment #24: 

Thanks again for your contribution.

I slightly changed the BISTs (to limit the scope of the shared variables and add a bug number) and pushed your change to the stable branch here:
https://hg.savannah.gnu.org/hgweb/octave/rev/290e7e3f859f

I noticed too late that you aren't in the list of contributors yet. Would you like to be added to that list? Can I use the name "Gaël Bonithon" for that?

Marking as ready for test.

Markus Mützel <mmuetzel>
Group administrator
Sun 27 Feb 2022 11:21:54 AM UTC, comment #23: 

Here is the patch with some tests from this thread and some comments to immediately have in mind a compact formulation.

(file #52939)

Gaël Bonithon <tamaranch>
Sat 26 Feb 2022 04:32:19 PM UTC, comment #22: 

Imho, the version from comment #18 is fine. If someone would like to tackle this in C++, this can be done in follow-up changes.
But we should probably add test cases to avoid regressing if/when this is done.

Could you please add some BISTs? Preferably ones that are failing with the current implementation, but are fixed with the changes in your patch from comment #18.

Markus Mützel <mmuetzel>
Group administrator
Wed 23 Feb 2022 07:07:42 PM UTC, comment #21: 

It actually works for `idivide(a,b)` i.e. `idivide(a,b,"fix")`, but not for `idivide(a,b,"ceil")`, you are right.

Never mind, back to #18. Is this version acceptable, or should we consider going "lower than `.m`"?

Gaël Bonithon <tamaranch>
Wed 23 Feb 2022 06:51:42 PM UTC, comment #20: 

It would be nice if it worked like that (and you could drop the assignment of x and replace the x ./ y by 1 ./ y, which would perhaps be even faster). But it won't work for


a=int64(4e16)+int64(1);
b=int64(2e8);
idivide(a,b)


isn't it?

Michael Leitner <mleitner>
Wed 23 Feb 2022 05:59:42 PM UTC, comment #19: 

Another approach: normalization. The extra cost in computation time compared to the current version is only 25% for `y` vector and 50% for `y` scalar (independently of the rounding rule this time), and the resulting code is simpler.

(file #52915)

Gaël Bonithon <tamaranch>
Tue 22 Feb 2022 12:18:35 PM UTC, comment #18: 

Regarding the patch at the `.m` level, to integrate the modification proposed by @mmuetzel in #12, I think that unfortunately we can't escape a distinction between scalar and vector for `x` and `y`, which makes the code much bigger. But I think that any more compact form pays off in computation time.

With this, we go from 100% increase in computation time (first version of the patch) to about 60%, compared to the current version.

Note that the worst case is in fact when `y` is a scalar and `x` is a vector. This is already noticeable with the current version, but only for "round", which is the only case that uses integer arithmetic.

(file #52902)

Gaël Bonithon <tamaranch>
Tue 22 Feb 2022 09:20:28 AM UTC, comment #17: 

Yes, there is overhead, but what is shown in bug #62006 is not directly applicable here: the comparison should not be between an empty loop and a loop with an empty .oct file function, but between a loop with a trivial built-in function and a loop with a trivial .oct file function. Because we would replace the z = x ./ y; (which I think to be exactly the same as z=rdivide(x,y), with the same overhead, without having tested), by this .oct file function call. Thus, is there more overhead for.oct file functions than for built-ins?

And in my naive view (never been closer to the internals of octave than .oct files) it would not mean diving into the interpreter: we would not need a new operator (in the sense of a new symbol that the interpreter has to parse), but just a new built-in function named somewhat like idivide_nat, that does native (right) integer division, exactly the same as rdivide (specialized to integers and to element-wise operation), but only without the post-processing. And (to repeat, in my naive view) adding a new built-in function should not be so hard, but definitely has also disadvantages (function look-up for any function would become slower, as the list to be searched becomes longer).

Michael Leitner <mleitner>
Tue 22 Feb 2022 07:39:24 AM UTC, comment #16: 

IIUC, `rdivide` (i.e. the `/` operator) is implemented on the "class level" of the integer classes.
The "most-straightforward" implementation for a new "C++ divide" would probably be an .oct file. Rik recently opened bug #62006 about a surprisingly large overhead for calling those functions.
Alternatively, it might be possible to implement a new "C++ divide" operator similarly to how the "rdivide" operator is implemented currently. But that would mean to dive pretty deep into the interpreter. (It might still be useful to have that "C++ divide" operator in other occasions as well.)
Anyway, it's hard to predict how performant which approach would be compared to the other short of implementing them and timing their performance.

It would be ok for me to keep the implementation on the .m file level for the moment. But if someone would like to dive into this and try different implementations, that could be a nice project to familiarize themselves with the .oct file interface and/or the interpreter.

Markus Mützel <mmuetzel>
Group administrator
Mon 21 Feb 2022 07:36:27 PM UTC, comment #15: 

Is there more function call overhead? z=x./y and z=rdivide(x,y) have the same speed (and result, I hope). rdivide is a built-in function, and we would need a private .oct file that does exactly what your code snippet does, but without computing w and without the if. Do .oct file functions have more overhead than built-ins? The only thing that is complicated would be to cast x or y to the appropriate integer representation if one of them is floating-point, which will be happening above the code snippet, otherwise it would be quite a trivial function. But in any case it should be possible to just duplicate the code and also make it to a a built-in function. Then you would do the normal division for the "round" case, otherwise precompute z by the new compiled function, where "fix" would then return immediately, and the other two do the post-processing.

While we are at it, one could make the whole idivde a compiled function. But is it worth it? I think not.


Michael Leitner <mleitner>
Mon 21 Feb 2022 07:17:54 PM UTC, comment #14: 

Integer division for Octave integer types doesn't follow the same rules as in C++ for compatibility with Matlab.

Afaict, it is implemented here:
https://hg.savannah.gnu.org/hgweb/octave/file/9a95ccd6c417/liboctave/util/oct-inttypes.h#l491

  // Division with rounding to nearest.  Note that / and % are
  // probably computed by a single instruction.

  static T div (T x, T y)
  {
    if (y != 0)
      {
        T z = x / y;
        T w = x % y;

        if (w >= y-w)
          z += 1;

        return z;
      }
    else
      return x ? octave_int_base<T>::max_val () : 0;
  }


You are right that we'll essentially need to revert that post-processing. But afaict, there is no way to "just" call the C++ division from .m code. (Please, correct me if that is wrong.)

It would probably be possible to unpack the underlying rep in an .oct file function, apply the C++ integer division to that and repack into an Octave array. But with the function call overhead, I don't know if that would be any faster...

Markus Mützel <mmuetzel>
Group administrator
Mon 21 Feb 2022 07:06:15 PM UTC, comment #13: 

Markus, you are missing the x(y_sel) instead of x in the "floor" case. Apart from that, it looks correct.

Yes, integer division is slower than floating-point division. But note that C99 and later specifies that integer division rounds towards zero, that is, what here seems to be the "fix" option. C++ (what octave is written in) seems to be the same. However, octave integer division per default does what the "round" mode does, thus there has to be some post-processing to be done on the result of x/y (on the C++ level). What I do not get now is why it is not even slower. On the other hand, there should be a way to speed up the "fix" case: just use C++ integer division without the post-processing. Or am I missing the point somewhere?

Michael Leitner <mleitner>
Mon 21 Feb 2022 06:40:20 PM UTC, comment #12: 

Thanks! I haven't tested thoroughly. But I can also see the slowdown here.

I believe that integer arithmetic is always a bit slower than floating point arithmetic. So, most probably we'll need to live with some slow-down if we want to fix this bug.

Nevertheless, it might be possible to get rid of some multiplications and calls to the `sign` and `abs` function.

With the attached modifications to your patch the function is a bit faster for the "ceil" and "floor" cases (albeit not as fast as we were with the floating point division).

Maybe something similar is possible for the "fix" case.

I haven't put much thought into this admittedly. So, please check if this still does the right thing...

(file #52898)

Markus Mützel <mmuetzel>
Group administrator
Mon 21 Feb 2022 06:02:38 PM UTC, comment #11: 

A quick profiling shows that the patch is about twice as slow as the current version in the worst cases ("ceil" and "floor"). The difference is a bit smaller for "fix" (no call to `abs()`), and there is of course no difference for "round".

Gaël Bonithon <tamaranch>
Sun 20 Feb 2022 08:02:00 PM UTC, comment #10: 

Hello,

Here is a proposal for a fix. I followed the option proposed by @mleitner, which seemed to me the most suitable.

Since the `typ` variable was no longer useful, I also revised and simplified the type test for `x` and `y`. This is not quite equivalent to what was done before for the `char` and `logical` types. If one wants to preserve the possibility of giving two `char` or two `logical` as input, a little extra work is needed.



(file #52895)

Gaël Bonithon <tamaranch>
Sat 09 Oct 2021 04:38:53 PM UTC, comment #9: 

@mleitner: Sorry. I took too long to write my first comment and didn't refresh the page. I completely missed your comment.

That approach sounds reasonable.

Markus Mützel <mmuetzel>
Group administrator
Sat 09 Oct 2021 04:32:11 PM UTC, comment #8: 

@mmuetzel: I think that whether the inputs are larger than flintmax or not is not the only point -- of course, if you have an input that cannot be represented as double and you cast it to double, it will be complicated to get it right, but the initial case of foo = 40000000000000000 is exactly representable (it is only a product of a large number of fives, which is representable, and the twos in the factorization end up in the exponent). No, the point is that in this case the double division gives as result (the double representation of) an integer, which does not change with round, ceil or floor. So using double division also does not work due to this point, and I can suspect that this can happen also for inputs smaller than flintmax.

And I really do think that it is better to test for c*b==a rather than subtracting the result of mod, if only because an integer multiplication should be much faster than an integer division (or the modulo operation). I would think that something along the lines of c+=(c*b>a)*sign(x) should work (for each mode a different expression), where x is one of a, b or c (I would have to work it out on paper).

Michael Leitner <mleitner>
Sat 09 Oct 2021 04:19:29 PM UTC, comment #7: 

This is where I am now after some testing. The tricky ones are the cases near intmin, intmax, and the unsigned cases. I think it might be better to split signed and unsigned types separately.


  if (strcmp (op, "fix"))
    ii = (y > 0) & (x > 0) | (y < 0) & (x < 0);      ## ii == same signs for x and y
    z(ii) = (x(ii) - mod (x(ii), y(ii))) ./ y(ii);   ## same as floor for positive
    z(~ii) = (-x(~ii) - mod (-x(~ii), y(~ii))) ./ y(~ii);   ## opposite sign for negative
    ## FIXME: unsigned types will break if user is mixing unsigned and signed.
  elseif (strcmp (op, "round"))
    z = x ./ y;
  elseif (strcmp (op, "floor"))
    z = (x - mod (x, y)) ./ y;   ## FIXME: underflow near intmin?
  elseif (strcmp (op, "ceil"))
    r = mod (x, y);
    z = (x - r) ./ y;            ## FIXME: underflow near intmin?
    z(r>0) += 1;
  else
    error ('idivide: unrecognized rounding type "%s"', op);
  endif


Hopefully there's a better way to do it.

Arun Giridhar <arungiridhar>
Group Member
Sat 09 Oct 2021 03:58:00 PM UTC, comment #6: 

Does that expression work if `foo` and `bar` are of opposite sign?

I might be misunderstanding. But that expression looks like it might be better suited for the "fix" case...

Markus Mützel <mmuetzel>
Group administrator
Sat 09 Oct 2021 03:36:32 PM UTC, comment #5: 

You're right -- calling floor and ceil separately won't do any good here because they're already integers. I wasn't thinking.

Good idea about rem and mod. Since they are written in C++, perhaps they can access native-mode integer division without rounding. Given two numbers foo and bar, this seems to work reliably for large n:


r = rem (foo,bar);
q = (foo - r) ./ bar;   # always floor

assert(all(q .* bar + r == foo)),


What should the fix and ceil expressions be?

Arun Giridhar <arungiridhar>
Group Member
Sat 09 Oct 2021 03:08:41 PM UTC, comment #4: 

IIUC, for integers smaller flintmax, converting to double and using `fix`, `floor`, or `ceil` probably works fine.
For larger integers, this does no longer work in general.

I don't understand what you mean by rounding the numerator or denominator separately. Aren't they already integer for integer input? What would rounding do in that case?

IIUC, the result of integer division is always rounded to the closest integer. (At least that is what the docstring of `idivide` says.)
Maybe we could use `rem` or `mod` to detect whether that was rounded up or down and correct the result accordingly?

Markus Mützel <mmuetzel>
Group administrator
Sat 09 Oct 2021 02:50:39 PM UTC, comment #3: 

c = a/b causes implicit rounding for integers in Octave, not native mode integer division. That was the reason for idivide if we need the floor or other mode. But here idivide itself is rounding incorrectly for numbers above flintmax, so testing c = a/b would need a workaround too.

How about this:


c = a/b;     ## causes integer rounding
if (c*b > a && mode is floor)
  c -= 1;
elseif (c*b < a && mode is ceil)
  c += 1;
etc for the other modes


But this seems flaky too. Not sure how to access native C++ integer division in an Octave function without going through double?

Anonymous
Sat 09 Oct 2021 02:04:27 PM UTC, comment #2: 

Why going through floating point: because it is easier. But as you have found, it does not work if the results are so large and in exact arithmetic so close to an integer that in double precision they are exact that integer, because then no rounding does happen.

I would say that there is no Matlab compatibility to break: there is a documentation, which agrees between Matlab and octave, so octave should follow its documentation (which at the moment it doesn't do), and if this is resolved and the results don't agree, then it is a Matlab bug (because it would mean that it doesn't follow its own documentation).

As to why not integer division: yes, I also think that this would be the way to do it. One would do the integer division according to its semantics (this is compiled code, and hopefully there is no bug there) c=a/b, then test whether c*b==a (in the native multiplication of the integer class), and if this is not true, add or subtract 1 to c (or do not change it at all) according to the chosen rounding mode. This would be not less performant, and should work as far as I can tell.

Rounding numerator and denominator separately would not work, I think.

Michael Leitner <mleitner>
Sat 09 Oct 2021 01:07:36 PM UTC, comment #1: 

The lines from scripts/general/idivide.m which cause the incorrect rounding are here:


    106   if (strcmp (op, "fix"))
    107     z = cast (fix (double (x) ./ double (y)), typ);
    108   elseif (strcmp (op, "round"))
    109     z = x ./ y;
    110   elseif (strcmp (op, "floor"))
    111     z = cast (floor (double (x) ./ double (y)), typ);
    112   elseif (strcmp (op, "ceil"))
    113     z = cast (ceil (double (x) ./ double (y)), typ);
    114   else
    115     error ('idivide: unrecognized rounding type "%s"', op);
    116   endif


Should it be done so that the numerator is rounded up while the denominator is rounded down? Will it break Matlab compatibility if this is changed?

Also, why are we going through floating point for integer division?

Arun Giridhar <arungiridhar>
Group Member
Sat 09 Oct 2021 12:55:58 PM UTC, original submission:  

The idivide function gives the same results for ceil and floor options for certain inputs close to perfect squares, even though that violates the arithmetic bounds for ceil and floor:

Test:

% test idivide
inputs = [  uint64(4e16) + (-4:+4), ...
            uint64(9e16) + (-8:+8), ...
            uint64(299792458) ^ 2 + (-12 : +4), ...
            uint64(16e18) + (-1024 : +1024), ...
         ];

function testfun (foo)
  bar = uint64 (round(sqrt(foo))) + (-1:1);
  q1 = idivide (foo, bar, "ceil");
  q2 = idivide (foo, bar, "floor");
  q1 == q2     # should be [0 0 0] or [0 1 0]
endfunction

for n = inputs(:)'
  testfun (n); # incorrectly returns [1 1 1]
end


The last line from the test should return [0 1 0] if the input is a perfect square, or [0 0 0] if it is not. Instead, for certain inputs close to perfect squares, we get this impossible result:

Result:

foo = 40000000000000000
bar =
  199999999  200000000  200000001

q1 =
  200000001  200000000  199999999

q2 =
  200000001  200000000  199999999

ans =
  1  1  1


Arun Giridhar <arungiridhar>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #52939:  bug61319-v5.patch added by tamaranch (5KiB - text/x-patch)
file #52915:  bug61319-v4.patch added by tamaranch (1KiB - text/x-patch)
file #52902:  bug61319-v3.patch added by tamaranch (3KiB - text/x-patch)
file #52898:  bug61319-v2.patch added by mmuetzel (2KiB - text/plain)
file #52895:  bug61319.patch added by tamaranch (2KiB - text/x-patch)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by tamaranch (Updated the item)
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by arungiridhar (Submitted the item)
  •  

    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
    2022-03-07 mmuetzel StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-02-27 mmuetzel StatusPatch Submitted Ready For Test
    2022-02-27 tamaranch Attached File- Added bug61319-v5.patch, #52939
    2022-02-26 mmuetzel StatusConfirmed Patch Submitted
    2022-02-23 tamaranch Attached File- Added bug61319-v4.patch, #52915
    2022-02-22 tamaranch Attached File- Added bug61319-v3.patch, #52902
    2022-02-21 mmuetzel Attached File- Added bug61319-v2.patch, #52898
    2022-02-20 tamaranch Attached File- Added bug61319.patch, #52895
    2021-10-09 mmuetzel StatusNone Confirmed
        Operating SystemGNU/Linux Any

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code