Opened 2 years ago
Closed 18 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
Attachments (2)
Change History (21)
Changed 19 months ago by Snark
comment:1 Changed 19 months ago by Snark
- Status changed from new to needs_review
comment:2 Changed 19 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 19 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 19 months ago by Snark
There are several things to say about the matter:
- for m an integer, indeed the formula you give is correct.
- 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.
- 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.
- 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 19 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 18 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 18 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 18 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 18 months ago by ppurka
- Status changed from needs_review to needs_work
- Work issues set to pari error on arm
comment:10 follow-up: ↓ 11 Changed 18 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 18 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 18 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 18 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 18 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 18 months ago by ppurka
- 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 18 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.
comment:17 Changed 18 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 18 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 18 months ago by jdemeyer
- Merged in set to sage-5.7.beta4
- Resolution set to fixed
- Status changed from positive_review to closed
Patch to make the situation better