This site is supported by donations to The OEIS Foundation.
User:Peter Luschny/ExtensionsOfTheBinomial
Blaise Pascal, Traité du triangle arithmétique (Cambridge University Library, Licensed under (CC BY-NC 3.0))
Extensions of the Binomial
Contents
The basic conditions
The binomial is a function ( is the set of all integers) based on a quotient of factorials
def b(n, k): return factorial(n) / (factorial(k) * factorial(n-k))
and ruled by three conditions which describe the nonzero values of this function:
- the Pascal condition,
- the Riordan condition (a.k.a. reflection),
- the Taylor condition (a.k.a. upper negation).
def binomial(n, k):
if 0 <= k <= n: return b(n, k) # Pascal
if k <= n < 0: return b(-k - 1, n - k) * (-1)^(n - k) # Riordan
if n < 0 <= k: return b(-n + k - 1, k) * (-1)^k # Taylor
return 0
|
|
|
||||||||||||||||||||||||||||||||||||
|
|
|||||||||||||||||||||||||||||||||||||
|
|
|
Alternatively the binomial function can be defined as the limit of the quotient of Gamma values:
def limit_binomial(n, k):
return limit(gamma(n + x) / (gamma(k + x)*gamma(n - k + x)), x = 1)
for n in (-5..5): print [limit_binomial(n, k) for k in (-5..5)]
Both definitions yield the same result. This proves that the three conditions form a coherent unit and are anything but an arbitrary collection of formulas.
There is also a basic representation of the binomial function defined with the Pascal and the Riordan condition but without the Taylor condition in terms of the matrix exponential: it is the matrix exponential of the matrix which has ..., -3, -2, -1, 0, 1, 2, 3, ... as subdiagonal and zero everywhere else.
matrix(SR,11,[ 0,0,0,0,0,0,0,0,0,0,0, -4,0,0,0,0,0,0,0,0,0,0, 0,-3,0,0,0,0,0,0,0,0,0, 0,0,-2,0,0,0,0,0,0,0,0, 0,0,0,-1,0,0,0,0,0,0,0, 0,0,0,0, 0,0,0,0,0,0,0, 0,0,0,0,0, 1,0,0,0,0,0, 0,0,0,0,0,0, 2,0,0,0,0, 0,0,0,0,0,0,0, 3,0,0,0, 0,0,0,0,0,0,0,0, 4,0,0, 0,0,0,0,0,0,0,0,0, 5,0 ]).exp()
[ 1 0 0 0 0 0 0 0 0 0 0] [-4 1 0 0 0 0 0 0 0 0 0] [ 6 -3 1 0 0 0 0 0 0 0 0] [-4 3 -2 1 0 0 0 0 0 0 0] [ 1 -1 1 -1 1 0 0 0 0 0 0] [ 0 0 0 0 0 1 0 0 0 0 0] [ 0 0 0 0 0 1 1 0 0 0 0] [ 0 0 0 0 0 1 2 1 0 0 0] [ 0 0 0 0 0 1 3 3 1 0 0] [ 0 0 0 0 0 1 4 6 4 1 0] [ 0 0 0 0 0 1 5 10 10 5 1]
In some sense this shows that the Riordan condition is more fundamental than the Taylor condition as it reveals a fundamental reflection principle (or duality) of the binomial function similar to the reflections of the Stirling numbers {n,k} = [-k,-n] or of the Lah numbers |n,k| = |-k,-n| whereas the Taylor condition mainly is a convenience relation suggested by Taylor's theorem from calculus.
Maple's and Mathematica's binomial functions behave like the binomial defined here. It is a sad fact that SageMath has only implemented a crippled form of the binomial function without the Riordan condition. (The same is true for some other open source math libraries.) Because of the ubiquity and importance of the binomial function we therefore cannot recommend SageMath for the use of implementing sequences on the OEIS.
A plot of log(abs(binomial(x,y))) for x, y in [-5..5] compared to the unmotivated cropped version of SageMath/GKP.
plot3d((x,y) -> log(abs(binomial(x,y))), -5..5, -5..5, grid=[60,60]); # Maple! plot3d((x,y) -> `if`(x<0 and y<0,0,log(abs(binomial(x,y)))), -5..5, -5..5, grid=[60,60]);
Sage users which want to switch to the extended form as described above can use the following function based on Sage's binomial.
def Binomial(n, k):
if n in ZZ and k in ZZ:
if n >= 0:
return binomial(n, k) # Pascal
if k >= 0:
return (-1)^k*binomial(-n+k-1, k) # Taylor
if k <= n:
return (-1)^(n-k)*binomial(-k-1, n-k) # Riordan
return 0
return binomial(n, k)
However if n or k is not in ZZ this function still has its problems. For example
for n in (-5..5): print [Binomial(n+1/2, k+1/2) for k in (-5..5)]
will compute as it should but refuse its service in the case of Binomial(n+1/2,k+1/3). To make it quite clear what response I expect from a CAS in this case: binomial(n+1/2, k+1/3) = C*r(n,k) where r(n,k) is a rational number and C is the constant sqrt(3)*Gamma(2/3)*Gamma(5/6)/Pi^(3/2).
References
- [BP] Blaise Pascal; Traité du triangle arithmétique; G. Desprez, 1665. Facsimile
- [ADAL] S. Anelli, Ernesto Damiani, Ottavio D'Antona, Daniel E. Loeb; Getting results with negative thinking; arXiv:9502214 (1995), p. 11.
- [LMMS] Ana Luzón, Donatella Merlini, Manuel A. Morón, Renzo Sprugnoli; Identities induced by Riordan arrays; Linear Algebra and its Applications, 436 (2012) 631-647.
- [WP] Pascal's triangle
- [MJK] M. J. Kronenburg; The Binomial Coefficient for Negative Arguments; arXiv:1105.3689v2
- [RS] Renzo Sprugnoli; Negation of binomial coefficients; Discrete Mathematics 308 (2008) 5070–5077.
An application of the reflection formula
The situation is quite confusing: depending on the conditions implemented the math tools will not only give wrong results relative to our definition but also overshadow valid identities. To make it clear as of which variant we talk it is convenient to control the conditions via parameters.
def _binomial(n,k, pascal = true, riordan = true, taylor = true):
def b(n,k): return factorial(n)/(factorial(k)*factorial(n-k))
if pascal and n >= 0 and k >= 0 and k <= n: return b(n,k)
if riordan and n < 0 and k < 0 and k <= n: return (-1)^(n-k)*b(-k-1,n-k)
if taylor and n < 0 and k >= 0: return (-1)^k*b(k-n-1,k)
return 0
For example let us consider in which case the reflection identity is true:
pascal=true, riordan=false, taylor=true [SageMath]
for n,k in Z: False
for n,k in N\{0}: False
pascal=true, riordan=true, taylor=true [Maple, Mathematica]
for n,k in Z: False
for n,k in N\{0}: True
pascal=true, riordan=true, taylor=false [Pure reflection]
for n,k in Z: True
for n,k in N\{0}: True
In my last blog post I introduced two number triangles which give raise to new identities related to the Stirling numbers.
![{\displaystyle \left[{n \atop k}\right]^{\star }={\frac {(2n)!}{(n+k)^{\underline {n}}}}\,\sum _{m=0}^{k}(-1)^{m+k}{\binom {n+k}{n+m}}\,\left[{n+m \atop m}\right]}](https://en.wikipedia.org/api/rest_v1/media/math/render/svg/626c170efee6f1c267463261e0f18bdf98085af8)
The numerical values of these numbers can be found in A268437 and A268438. The identities display a nice formal symmetry:
![{\displaystyle {\binom {n}{k}}\sum _{i=0}^{k}{\binom {k}{i}}\left\{{n-k \atop i}\right\}^{\star }={\binom {-k}{-n}}\sum _{i=0}^{n-k}{\binom {-n}{i}}\left[{n-k \atop i}\right]^{\star }}](https://en.wikipedia.org/api/rest_v1/media/math/render/svg/fd7a25fceff49f89d7b39ac99e40a6aaf4943201)
![{\displaystyle {\binom {n}{k}}\sum _{i=0}^{k}{\binom {k}{i}}\left[{n-k \atop i}\right]^{\star }\ \ ={\binom {-k}{-n}}\sum _{i=0}^{n-k}{\binom {-n}{i}}\left\{{n-k \atop i}\right\}^{\star }}](https://en.wikipedia.org/api/rest_v1/media/math/render/svg/b216708302915fe3f10db6ee85426ab5126c2b59)
Let's check these the formulas for some small values.
C = lambda n,k: Binomial(n,k)
R = lambda n,k: rising_factorial(n,k)
F = lambda n,k: falling_factorial(n,k)
N = lambda n,k: factorial(2*n)/F(n+k,n)
S = lambda n,k: stirling_number2(n,k)
s = lambda n,k: stirling_number1(n,k)
B = lambda n,k: N(n,k)*sum((-1)^(m+k)*C(n+k,n+m)*S(n+m,m) for m in (0..k))
b = lambda n,k: N(n,k)*sum((-1)^(m+k)*C(n+k,n+m)*s(n+m,m) for m in (0..k))
A1 = lambda n,k: R(n-k+1,n-k)*S(n,k)
A2 = lambda n,k: C(n,k)*sum(C(k,i)*B(n-k,i) for i in (0..k))
A3 = lambda n,k: C(-k,-n)*sum(C(-n,i)*b(n-k,i) for i in (0..n-k))
B1 = lambda n,k: R(n-k+1,n-k)*s(n,k)
B2 = lambda n,k: C(n,k)*sum(C(k,i)*b(n-k,i) for i in (0..k))
B3 = lambda n,k: C(-k,-n)*sum(C(-n,i)*B(n-k,i) for i in (0..n-k))
def test(f):
for n in (0..6): print [f(n,k) for k in (0..n)]
test(A1);test(A2);test(A3)
test(B1);test(B2);test(B3)
From our discussion it is clear that these identities fail when checked with SageMath but hold with our definition or checked with Maple or Mathematica. I, for my part, rather abandon SageMath than these beautiful identities.
How to compute the binomial?
A nice implementation of Pascal's triangle can be found on rosettacode (GNU FD-License 1.2) and is based on summing the previous row content.
def pascal(n):
"""Prints out n rows of Pascal's triangle.
It returns False for failure and True for success.
"""
row = [1]
k = [0]
for x in range(max(n,0)):
print row
row = [l + r for l,r in zip(row + k,k + row)]
return n >= 1
Clearly this algorithm is not well suited to compute the binomial coefficient for large n and k. We will show an efficient algorithm implemented with SageMath below.
Amazingly an efficient algorithm of the binomial coefficients was not long ago unavailable in the world's fastest arithmetic library, the GNU Multiple Precision library GMP. Indeed in Feb. 2010 I wrote on my homepage:
Here is a mini-benchmark. It compares our implementation with the library functions of GMP-5.0.1 (The timings for MPIR 1.3.0 are very similar.) The speedup is remarkable. For some input values the binomial function is several hundred or even several thousand times faster than the 5.0.1 equivalents. Testing binomial: (6400000,2133333) ParallelBinomial: 0.422 sec PrimeBinomial: 0.610 sec GMP5Binomial: 1733.750 sec
When the lead developer of GMP Torbjorn Granlund read this he decided to replace the GMP code for the binomial by the algorithm below. Interestingly not only the performance increased but also the code became much simpler. Marco Bodrato finally carried out the (highly optimized) 'C'-implementation.
The algorithm is not based on the quotient of factorials n!/(k!(n-k)!) but on an old idea of Ernst Kummer which says that the largest integer k such that p^k (p prime) divides the binomial coefficient C(n,m) is equal to the number of carries when m is added to n - m in base p. (Ernst Kummer, Über die Ergänzungssätze zu den allgemeinen Reciprocitätsgesetzen, 1852, Journal für die reine und angewandte Mathematik 44: 93–146.)
The implementation (without reflection) can also be found in Sympy's "sympy/functions/combinatorial/factorials.py" without any notes on their origin.
def kummer(n, k):
if k == 0 or k == n: return 1
if k > n // 2: k = n - k
nk = n - k
pfact = []
rootN = isqrt(n)
for prime in primes(2, n):
if prime > nk:
pfact.append(prime)
continue
if prime > n // 2:
continue
if prime > rootN:
if n % prime < k % prime:
pfact.append(prime)
continue
r, N, K, p = 0, n, k, 1
while N > 0:
r = 1 if N % prime < K % prime + r else 0
if r == 1: p *= prime
N = // prime
K = K // prime
if p > 1: pfact.append(p)
return prod(pfact)
def Binomial(n, k):
if 0 <= k <= n: return kummer(n, k)
if k <= n < 0: return kummer(-k - 1, n - k) * (-1)^(n-k)
if n < 0 <= k: return kummer(-n + k - 1, k) * (-1)^k
return 0
Basically all we use here from SageMath is the function primes(2, n). Thus adding a Sieve of Eratosthenes would make our implementation stand-alone.
Extensions of factorial products
There are two other basic functions which have extensions from to : the falling product
fp(n, k) = prod([(n-j) for j in range(k)])
and the rising product
rp(n, k) = prod([(n+j) for j in range(k)])
which are defined for integer n ≥ 0 and k ≥ 0 and extended to the falling factorial respectively the rising factorial which are functions .
|
|
|
||||||||||||||||||||||||||||||||||||
|
|
|||||||||||||||||||||||||||||||||||||
|
|
|
def fp(n, k):
return prod([(n-j) for j in range(k)])
def falling_factorial(n, k):
if k >= 0: return fp(n, k)
if n < k or n >= 0:
return Rational((-1)^k, fp(-n-1, -k))
return 0 # (or ComplexInfinity)
|
|
|
||||||||||||||||||||||||||||||||||||
|
|
|||||||||||||||||||||||||||||||||||||
|
|
|
def rp(n,k):
return prod([(n+j) for j in range(k)])
def rising_factorial(n, k):
if k >= 0: return rp(n, k)
if k > -n or n <= 0:
return Rational((-1)^k, rp(-n+1, -k))
return 0 # (or ComplexInfinity)
In our context it is interesting to see that both functions observe conditions which could be called reflection laws. This is obvious from studying the tables or looking at the above implementations.
The definitions can be extended to the whole complex plane. The only difference in the tables would be that the '0's in the tables at the left hand side (the case when k < 0) have to be replaced by 'complex infinity'. So nothing is gained from our integer point of view but things are certainly more complex now as the codomain of the functions is no longer .
With our implementations we want to go one further step of optimization. First of all we want to make clear that a product which is implemented by the naive approach to multiply one factor after the other is a very poor one and should not be used.
def product(a, b): # inefficient, do not use.
p = 1
for k in range(a, b + 1):
p *= k
return p
Instead we need an implementation which employs a divide-and-conquer strategy to exploit the asymptotical speed of the recursion.
def rec_product(a, b):
"""
Input:
a -- integer, lower bound of interval (inclusive)
b -- integer, upper bound of interval (inclusive)
Output:
If a <= b then return
a * (a+1) * ... * (b-1) * b
else return 1.
"""
n = b - a
if n < 25:
p = 1
for k in range(a, b + 1):
p *= k
return p
m = a + n // 2
return rec_product(a, m) * rec_product(m + 1, b)
Now we can implement our final version of the falling and rising factorial:
def falling_factorial(n, k):
if k >= 0:
return rec_product(n - k + 1, n)
if n < k or n >= 0:
p = rec_product(-n + k, -n - 1)
return Rational(-1 if k % 2 else 1, p)
return 0
def rising_factorial(n, k):
if k >= 0:
return rec_product(n, n + k - 1)
if k > -n or n <= 0:
p = rec_product(-n + 1, -n - k)
return Rational(-1 if k % 2 else 1, p)
return 0
Assorted triangles using reflection
S1 Stirling cycle numbers A132393 and S2 Stirling set numbers A048993. L denotes the unsigned Lah numbers A271703, E1 the Eulerian numbers A173018 and E2 the second-order Eulerian numbers A201637.
Sum of powers
Sum_{j=0..n} n^j*C(-j,-n) = A055897(n).
Sum_{j=0..n} n^j*C(-j,-n)*(-1)^(n-j) = A053506(n+1) (n \ge 1)
Sum_{j=0..n} n^j*C(-j-1,-n-1)*(-1)^(n-j) = A000169(n)
Sum_{j=0..n} j^n*C(-j,-n) = A001710(n+1) (n \ge 1).
Sum_{j=0..n} j^n*C(-j,-n)*(-1)^(n-j) = A195242(n).
Sum_{j=0..n} j^n*C(-j-1,-n-1)*(-1)^(n-j) = A072034(n)
Stirling cycle numbers S1 A132393
Sum_{j=0..n} S1(j,k)*C(-j,-n) = A269954(n,k)
Sum_{j=0..n} S1(j,k)*C(-j,-n)*(-1)^(n-j) = A269951(n,k)
Sum_{j=0..n} S1(j,k)*C(-j-1,-n-1) = A269953(n,k)
Sum_{j=0..n} S1(j,k)*C(-j-1,-n-1)*(-1)^(n-j) = A094816(n,k)
Sum_{j=0..n} S1(k,j)*C(-j-1,-n-1)*(-1)^(n-j) = A271700(n,k)
Stirling set numbers S2 A048993
Sum_{j=0..n} S2(j,k)*C(-j,-n) = S2(n-1,k-1) with S2(-1,-1) = 1
Sum_{j=0..n} S2(j,k)*C(-j,-n)*(-1)^(n-j) = A269952(n,k) = S2(n+1,k+1)-S2(n,k+1)
Sum_{j=0..n} S2(j,k)*C(-j-1,-n-1) = A105794(n,k)
Sum_{j=0..n} S2(j,k)*C(-j-1,-n-1)*(-1)^(n-j) = S2(n+1,k+1)
Sum_{j=0..n} S2(k,j)*C(-j-1,-n-1)*(-1)^(n-j) = A271702(n,k)
Unsigned Lah numbers L A271703
Sum_{j=0..n} L(j,k)*C(-j,-n) = A216154(n,k)
Sum_{j=0..n} L(j,k)*C(-j,-n)*(-1)^(n-j) = A271704(n,k)
Sum_{j=0..n} L(j,k)*C(-j-1,-n-1) = A271706(n,k)
Sum_{j=0..n} L(j,k)*C(-j-1,-n-1)*(-1)^(n-j) = A271705(n,k)
Eulerian numbers E1 A173018
Sum_{j=0..n} E1(j,k)*C(-j,-n) = A271698(n,k)
Sum_{j=0..n} E1(j,k)*C(-j-1,-n-1) = A271697(n,k)
Sum_{j=0..n} E1(j,k)*C(-j-1,-n-1)*(-1)^(n-j) = A272098(n,k)
