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

[SnarkBoojum]

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. attachment: trac_12448-fix_binomial.patch

Component: basic arithmetic

Author: Punarbasu Purkayastha

Reviewer: Julien Puydt

Merged: sage-5.7.beta4

Issue created by migration from https://trac.sagemath.org/ticket/12448

[SnarkBoojum]

Attachment: bug_12448.patch.gz

Patch to make the situation better

[SnarkBoojum]

comment:2

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.

[ppurka]

comment:3

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?

[SnarkBoojum]

comment:4

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.

[ppurka]

comment:5

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.

[ppurka]

comment:6

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

[SnarkBoojum]

comment:7

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)".

[ppurka]

comment:8

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.

[ppurka]

Work Issues: pari error on arm

[SnarkBoojum]

comment:10

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).

[ppurka]

comment:11

Replying to @SnarkBoojum:

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?

[SnarkBoojum]

comment:12

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.

[ppurka]

comment:13

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

[SnarkBoojum]

comment:14

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

[ppurka]

Reviewer: Julien Puydt

[ppurka]

comment:15

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

Patchbot apply trac_12448-fix_binomial.patch

[ppurka]

Author: Punarbasu Purkayastha

[ppurka]

Changed work issues from pari error on arm to none

[ppurka]

Work Issues: fix binomial(1000,1001)

[ppurka]

comment:16

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

[ppurka]

Attachment: trac_12448-fix_binomial.patch.gz

Apply to sage/devel

[ppurka]

comment:17

Should be fixed now.

Patchbot apply trac_12448-fix_binomial.patch

[ppurka]

Changed work issues from fix binomial(1000,1001) to none

[SnarkBoojum]

comment:18

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

[jdemeyer]

Merged: sage-5.7.beta4