# Ticket #12448(closed defect: fixed)

Opened 16 months ago

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

Reported by: Owned by: Snark AlexGhitza minor sage-5.7 basic arithmetic N/A Julien Puydt Punarbasu Purkayastha sage-5.7.beta4

### Description (last modified by ppurka) (diff)

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

## Change History

### Changed 5 months ago by Snark

Patch to make the situation better

### comment:1 Changed 5 months ago by Snark

• Status changed from new to needs_review

### comment:2 Changed 5 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 5 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 5 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 5 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 4 months ago by ppurka

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

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 4 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 4 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 4 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 4 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 4 months ago by ppurka

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 4 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 4 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 4 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 4 months ago by ppurka

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

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 4 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 4 months ago by ppurka

Apply to sage/devel

### comment:17 Changed 4 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 4 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 3 months ago by jdemeyer

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