Opened 3 years ago

Closed 21 months ago

#12448 closed defect (fixed)

The binomial implementation does a quotient of gamma values, which is wrong

Reported by: Snark Owned by: AlexGhitza
Priority: minor Milestone: sage-5.7
Component: basic arithmetic Keywords:
Cc: Merged in: sage-5.7.beta4
Authors: Punarbasu Purkayastha Reviewers: Julien Puydt
Report Upstream: N/A Work issues:
Branch: Commit:
Dependencies: Stopgaps:

Description (last modified by ppurka)

There are special formulas to apply for quotients of gamma values, since those quotients can be pretty reasonable even though the numerator and denominator are too big, for example:

sage: binomial(float(1000),float(1001))
nan
sage: binomial(1000,1001)
0

and:

sage: binomial(float(1001),float(1000))
nan
sage: binomial(1001,1000)
1001

Apply to sage/devel

  1. trac_12448-fix_binomial.patch

Attachments (2)

bug_12448.patch (1.4 KB) - added by Snark 22 months ago.
Patch to make the situation better
trac_12448-fix_binomial.patch (3.9 KB) - added by ppurka 21 months ago.
Apply to sage/devel

Download all attachments as: .zip

Change History (21)

Changed 22 months ago by Snark

Patch to make the situation better

comment:1 Changed 22 months ago by Snark

  • Status changed from new to needs_review

comment:2 Changed 22 months ago by Snark

  • Status changed from needs_review to needs_work

I'm not too happy with my patch after all: binomial(float(1001),float(1)) is still a nan ; I should make it a little smarter.

comment:3 Changed 22 months ago by ppurka

The current code (not your patch) looks overly complicated for a binomial function. What is the definition being used? Can't it simply be defined as (assuming m is an integer):

x * (x-1)/2 * (x-2)/3 * ... * (x-m+1)/m

Doesn't this hold irrespective of whether x is real or integer? The above formula also gives a simple way of implementing the binomial. For x integer, the computations from left to right remains an integer all the time. Secondly the computation remains bounded at all times (if the binomial is not large) even if the values of x and m are large. There is no need to invoke the gamma function. One could also start from the other end to keep the values during the computations even less (note that at each step integral values of x and m implies that the computations remain integers; so we are not breaking anything):

(x-m+1) * (x-m+1)/2 * (x-m+3)/3 * ... * x/m

Am I missing some definition of the binomial here, which contains the gamma function necessarily?

comment:4 Changed 22 months ago by Snark

There are several things to say about the matter:

  1. for m an integer, indeed the formula you give is correct.
  2. but it's possible to define the binomial with m non-integer using gamma-functions, and I think that's what the one who coded the current floating code had in mind.
  3. my patch is about partially unwinding the gamma factors like you mention (with step by step simplification which avoids overflows) before finally calling the gamma function with small parameters.
  4. but my patch should only start looping when the parameters are indeed big (or it will be slow) and use symmetry so binomial(float(1001),float(1)) works.

comment:5 Changed 22 months ago by ppurka

The case of m noninteger is already properly handled in the code. The code checks whether either m or x-m is an integer and only then proceeds. That's why I think calling the gamma function is not necessary.

comment:6 Changed 21 months ago by ppurka

  • Description modified (diff)
  • Milestone set to sage-5.7
  • Status changed from needs_work to needs_review

Added a patch which uses pari now, as was indicated as a future fix in a comment. This speeds up all the floating point binomials by about a factor of 5, and fixes the nan that was being returned from large floats. Needs review now :)

Patchbot apply trac_12448-fix_binomial.patch

comment:7 Changed 21 months ago by Snark

Patching and compiling are ok on amd64 and armv7l.

On amd64, doctesting is ok, but on arm, the line 3143:

a = binomial(RR(1140000.78), 42000000)

makes the line 3213 (added by the patch):

return P(pari(x).binomial(m))

raise an error "ParriError?: (15)".

comment:8 Changed 21 months ago by ppurka

Can you check which part raises the pari error? Is it pari(RR(1140000.78)) or the .binomial() part, or the coercion part? I don't have an arm box, nor do I know pari, so I can only guess.

comment:9 Changed 21 months ago by ppurka

  • Status changed from needs_review to needs_work
  • Work issues set to pari error on arm

comment:10 follow-up: Changed 21 months ago by Snark

It's the .binomial part which raises an error : binomial(RR(1140000.78), 23310031) works, but binomial(RR(1140000.78), 23310032) raises the error (dichotomy).

comment:11 in reply to: ↑ 10 Changed 21 months ago by ppurka

Replying to Snark:

It's the .binomial part which raises an error : binomial(RR(1140000.78), 23310031) works, but binomial(RR(1140000.78), 23310032) raises the error (dichotomy).

So, it seems to be some overflow problems in arm? Should we just change the example

a = binomial(RR(1140000.78), 42000000)

to the one below which works on arm? After all, what this is testing is that the binomial is fast, which is well tested even with something that is half the size of the previous one

a = binomial(RR(1140000.78), 23310031)

Or, does this point to some problem in pari?

comment:12 Changed 21 months ago by Snark

Well, if some piece of code works on a platform but doesn't on another, then there's obviously something to investigate. Is it easy to write a C pari-using test case? I'd gladly give it a try.

comment:13 Changed 21 months ago by ppurka

I would put this problem on sage-devel, if you agree. I don't know pari and neither can I fix architecture specific idiosyncrasies.

comment:14 Changed 21 months ago by Snark

According to sage-devel, indeed the doctest should be modified so it works on both 64 bits and 32 bits architectures.

comment:15 Changed 21 months ago by ppurka

  • Authors set to Punarbasu Purkayastha
  • Reviewers set to Julien Puydt
  • Status changed from needs_work to needs_review
  • Work issues pari error on arm deleted

Thanks for following up in sage-devel. I have updated the patch with a smaller number.

Patchbot apply trac_12448-fix_binomial.patch

comment:16 Changed 21 months ago by ppurka

  • Status changed from needs_review to needs_work
  • Work issues set to fix binomial(1000,1001)

Oops. I just found a failing case - the first example.

Changed 21 months ago by ppurka

Apply to sage/devel

comment:17 Changed 21 months ago by ppurka

  • Status changed from needs_work to needs_review
  • Work issues fix binomial(1000,1001) deleted

Should be fixed now.

Patchbot apply trac_12448-fix_binomial.patch

comment:18 Changed 21 months ago by Snark

  • Status changed from needs_review to positive_review

The new patch applies, compiles and passes all doctests in arith.py on armv7l and amd64.

comment:19 Changed 21 months ago by jdemeyer

  • Merged in set to sage-5.7.beta4
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.