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

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

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

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>
Project Member
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;
    error ('idivide: unrecognized rounding type "%s"', op);

Hopefully there's a better way to do it.

Arun Giridhar <arungiridhar>
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>
Project Member
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>
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>
Project Member
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?

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>
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 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]

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

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:


foo = 40000000000000000
bar =
  199999999  200000000  200000001

q1 =
  200000001  200000000  199999999

q2 =
  200000001  200000000  199999999

ans =
  1  1  1

Arun Giridhar <arungiridhar>


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

Attach Files:

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -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)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 2 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2021-10-09 mmuetzel StatusNone => Confirmed
        Operating SystemGNU/Linux => Any

    Back to the top

    Powered by Savane 3.7