bugGNU Octave - Bugs: bug #35611, Bug in bincoeff with large n and k

 
 

bug #35611: Bug in bincoeff with large n and k

Submitted by:  None
Submitted on:  Thu 23 Feb 2012 03:28:28 PM UTC  
 
Category: NoneSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: Wont FixAssigned to: None
Originator Name: Marcus MeyerOriginator Email: -unavailable-
Open/Closed: ClosedRelease: dev
Operating System: Microsoft Windows

Add a New Comment(Rich Markup)
   

You are not logged in

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

 

Fri 24 Feb 2012 04:31:33 PM UTC, comment #4:

Marco,

This is just what I suspected. The underlying routine for gammaln probably differs by a few eps depending on the backend library. For your system, the result is slightly under the true result and rounding brings it up to the correct value. For the bug reporter's system the result is slightly higher than the true result and rounding causes an error.

To recap, there isn't any consistency that Octave can exploit and so this bug can't be "fixed". On the other hand, there is a clear workaround which is to use the appropriate function, nchoosek for integers, and bincoeff for real values.

Rik <rik5>
Project Administrator
Fri 24 Feb 2012 08:16:51 AM UTC, comment #3:

Dear all,

just to inform that with 3.6.0 in linux, I get, for bincoeff, 1.73103094564398e+13 before rounding, 1.73103094564400e+13 after rounding and for nchoosek 1.73103094564400e+13. The same is true for 3.2.4.
This time I'm lucky that the results are right, but should I think that my system is different from all the others?

Cheers,

Marco

Marco Caliari <caliari>
Project Member
Thu 23 Feb 2012 05:58:19 PM UTC, comment #2:

Check the documentation for nchoosek and bincoeff. They compute the same mathematical value, but the underlying algorithms are very different. In particular, nchoosek is optimized for integer inputs and gives the correct answer in both Octave and Matlab. bincoeff is a much broader function and can handle real inputs as well as negative inputs.

Quoting from the documentation for bincoeff:

Quoting from the documentation for nchoosek:

nchoosek (100,10) = 1.7e13 which is representable exactly as an integer in doubele precision arithmetic (you can check bitmax() which is 9e15 for doubles). Thus, nchoosek is really the correct function to use for this example.

The line to take a look at is 98 of bincoeff.m which I reproduce below.

Before the rounding, the value of b is 17310309456440.8. In this case, if we used truncation, fix(), rather than round we would get the correct answer. However, I'm note sure this is a general result. I believe round was chosen so that 50% of the time we would guess low and 50% of the time we would guess high and on average the algorithm would be more accurate. It would be valuable if you could prove any relation that might exist between the input and the fractional part of the bincoeff algorithm. For example, is it always the case for integers less than bitmax that the result should use fix() rather than round()?

Rik <rik5>
Project Administrator
Thu 23 Feb 2012 03:39:08 PM UTC, comment #1:

The code for bincoeff isn't exceptionally complicated and it's fairly self-contained:

http://hg.savannah.gnu.org/hgweb/octave/file/d17237256856/scripts/miscellaneous/bincoeff.m#l55

Do you want to attempt a fix yourself? There's an obvious rounding error.

Jordi GutiƩrrez Hermoso <jordigh>
Project Administrator
Thu 23 Feb 2012 03:28:28 PM UTC, original submission:

octave-3.2.4.exe:10> bincoeff(100,10)
ans = 17310309456441

The correct result is:

octave-3.2.4.exe:13> 101198971931231394
ans = 17310309456440

Matlab gives the correct result with nchoosek(100,10).
Maybe bincoeff can be improved for large input?

Best regards and thanks
Marcus

Anonymous

 

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

Attach File(s):
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by jordigh (Posted a comment)
  • -unavailable- added by None (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 4 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Fri 24 Feb 2012 04:31:33 PM UTCrik5StatusConfirmed=>Wont Fix
      Open/ClosedOpen=>Closed
    Thu 23 Feb 2012 03:39:08 PM UTCjordighStatusNone=>Confirmed
      Release3.2.4=>dev

    Back to the top


    Powered by Savane 3.1-cleanup1