Issue 1:
scipy.stats.binom.sf and scipy.stats.binom_test returns incorrect values for large inputs. The bounds on the inputs is not documented, and similar functions in other libraries (e.g. R's binom.test) do not have this problem.

Reproducing code example:

>>> np.seterr(all='warn')
{'divide': 'warn', 'over': 'warn', 'under': 'ignore', 'invalid': 'warn'}
>>> scipy.stats.binom_test(1e7, 2e7, 0.5, 'greater')
0.5094024848633918
>>> scipy.stats.binom_test(1e8, 2e8, 0.5, 'greater')
0.6860874231168742
>>> scipy.stats.binom_test(1e9, 2e9, 0.5, 'greater')
0.8854076875623922
>>> scipy.stats.binom_test(1e10, 2e10, 0.5, 'greater')
nan

Error message:

None.

Expected behaviour:

All calls in the repro example should return a value near 0.5. For example, in R:

> binom.test(1e7, 2e7, 0.5, 'greater')$p.value
[1] 0.5000892
> binom.test(1e8, 2e8, 0.5, 'greater')$p.value
[1] 0.5000282
> binom.test(1e9, 2e9, 0.5, 'greater')$p.value
[1] 0.5000089
> binom.test(1e10, 2e10, 0.5, 'greater')$p.value
[1] 0.5000028

Notes:

Internally, for the "greater" alternative, scipy calls binom.sf(x - 1, n, p), which then calls scipy.special.bdtrc(floor(x), n, p). This forwards to Cephes' bdtrc, according to scipy's docs.

We can verify that the bug lies within bdtrc:

>>> scipy.special.bdtrc(1e10, 2e10 - 1e10 + 1, 0.5)
nan

Non-broken implementations of the regularized incomplete beta function (such as TensorFlow's tf.math.betainc) will return the expected value (0.5).

This also affects anything else that calls bdtrc, including binom.sf. There may be a case that bdtrc is Working as Intended (because the documentation is explicit that it's just a wrapper for Cephes' buggy implementation), but I think for binom.sf and binom_test it's clear that there is at least a documentation bug if not an opportunity to make things better.

I think there are 4 viable solutions, in what I think is most-preferred to least-preferred:

File a bug with Cephes and wait for them to fix it, documenting the breakage in the meantime.

Use a normal approximation in binom.sf when "k" and "n" are large (perhaps checking that "p" isn't too biased). Fortunately the circumstances when bdtrc breaks usually coincide with when this approximation is good. :)

Implement a regularized incomplete beta function within scipy that isn't broken, instead of using Cephes.

Claim this is Working as Intended, but document the limitation.

I'm happy to send Pull Requests given some direction about which of the choices above to go with.

Scipy/Numpy/Python version information:

>>> import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
1.5.4 1.19.4 sys.version_info(major=3, minor=6, micro=8, releaselevel='final', serial=0)

Issue 2:
The following command correctly returns around 0.5.

scipy.stats.binom.cdf(3000000, 10000000, .3)

The probability that a binomial RV with p=0.3 is less than or equal to 3,000,000 in 10,000,000 draws is around 0.5.

The same should be true for P(X<=30,000,000) in n=100,000,000 draws. (these numbers are 10 times higher than the preceding inputs)

scipy.stats.binom.cdf(30000000, 100000000, .3)

should return around 0.5, but it returns 0.57036.

I also tested scipy.special.betainc, since that's what is used for calculating Binomial's CDF, and it exhibits the same problem.