Efficient Computation of Binomial Coefficients and Applications

Abstract

1  Introduction

bk= k å i=0  Cki ai,    bk= n å i=k  Cik ai,    bk= k å i=0  Cn-ik-iai, ¼    (k=0, 1, ¼, n)

where,

Cni= æ ç è n i ö ÷ ø = n! (n-i)! i!

are called the binomial coefficients. Assume a table of binomial coefficients is given. Then, a straight forward method to compute, for example, b_k=å_i=0^k C_kⁱ p_i (k=1,2,...,n) is to multiply p_i by C_kⁱ and then add them together. However, even with the supply of binomial coefficients, this simple method requires n(n-1)/2 multiplications and the same number of additions. In this paper, we shall derive an efficient algorithm to compute b_k. Using our method, only n(n-1)/2 additions are needed to compute all b_k=å_i=0^k C_kⁱ p_i (k=0, 1, ..., n) in the absence of a table of binomial coefficients.

Our methods are derived analogously to the construction of Pascal triangle, the well-known algorithm to compute binomial coefficients. Some practical applications of our methods are given in this paper.

2  Computing b_k=å_i=0^k C_kⁱ a_i and b_k=å_i=0^k C_kⁱ a_ibⁱ

   c(0,0) = 1
   do k=1, n
      c(k,0) = 1
      c(k,k) = 1
      do i=1, k-1
         c(k,i) = c(k-1,i-1) + c(k-1,i)
      enddo
   enddo

1 1 1 1 2 1 1 3 3 1

We now consider computation of b_k=å_i=0^k C_kⁱ a_i. For n=3, we have

æ ç ç ç ç è b0 b1 b2 b3 ö ÷ ÷ ÷ ÷ ø = æ ç ç ç ç è 1 0 0 0 1 1 0 0 1 2 1 0 1 3 3 1 ö ÷ ÷ ÷ ÷ ø æ ç ç ç ç è a0 a1 a2 a3 ö ÷ ÷ ÷ ÷ ø = æ ç ç ç ç è a0 a0+a1 a0+2a1+a2 a0+3a1+3a2+a3 ö ÷ ÷ ÷ ÷ ø .

It is noted that the transformation matrix resembles the property of Pascal triangle. Therefore, analogous to the construction of Pascal triangle, we may derive the following method to compute b_k

   Copy a(i) into b(i)
   do k=1, n
   do i=n, k, -1
      b(i) = b(i) + b(i-1)
   enddo

   (k = 1)              (k = 2)                   (k = 3)
   b3=b3+b2=a3+a2       b3=b3+b2=a3+2a2+a1        b3=b3+b2=a3+3a2+3a1+a0
   b2=b2+b1=a2+a1       b2=b2+b1=a2+2a1+a0
   b1=b1+b0=a1+a0 

We may extend the above method to compute b_k=å_i=0^k C_kⁱ a_ibⁱ:

   Copy a(i) into b(i)
   do k=1, n
   do i=n, k, -1
      b(i) = b(i) * beta + b(i-1)
   enddo

We now consider the application of the above method. In computational geometry, a Bézier polynomial is commonly used to describe the shape of geometric entity. A Bézier polynomial of degree n may be represented in terms of Bernstein polynomials [2-3]:

P(t)= n å i=0  Cni(1-t)n-iti pi.

For efficient evaluation of Bézier polynomial and its derivatives, we often convert it into a canonical polynomial so that Horner method can be employed. By the binomial expansion theorem, we have

n å i=0  Cni(1-t)n-itipi = n å i=0  n-i å j=0  (-1)jCniCn-ijpiti+j.

Then, we want to find b_k such that

n å k=0  bktk= n å i=0  n-i å j=0  (-1)jCniCn-ijpiti+j.

Since the coefficients of the same degree terms at both sides should be equal, we thus have

bk= k å i=0  k-i å j=k-i  (-1)jCniCn-ijpiti+j = k å i=0  (-1)k-iCniCn-ik-ipi.

From the relation C_nⁱC_n-i^k-i=C_n^kC_kⁱ, we obtain

bk= k å i=0  (-1)k-iCnkCkipi = (-1)kCnk k å i=0  (-1)iCkipi.

To compute b_k efficiently, we first evaluate å_i=0^k(-1)ⁱC_kⁱp_i using the above method, then multiply the results by (-1)^kC_n^k (noting that C_n^k=[(n-k+1)/k]C_n^k-1). The pseudo code is given as follows:

  do i=0, n, 2
     b(i) = a(i)
  enddo
  do i=1, n, 2
     b(i) = -a(i)
  enddo

  do k=1, n
  do i=n, k, -1
     b(i) = b(i) + b(i-1)
  enddo

  temp = 1.0
  do i=1, n
     temp = -(n + 1 - i) * temp / i
     b(i) = b(i) * temp
  enddo

3  Compute b_k=å_i=kⁿ C_i^ka_i, b_k=å_i=kⁿ C_i^ka_ib^i-k

æ ç ç ç ç è b0 b1 b2 b3 ö ÷ ÷ ÷ ÷ ø = æ ç ç ç ç è 1 1 1 1 0 1 2 3 0 0 1 3 0 0 0 1 ö ÷ ÷ ÷ ÷ ø æ ç ç ç ç è a0 a1 a2 a3 ö ÷ ÷ ÷ ÷ ø = æ ç ç ç ç è a0+a1+a2+a3 a1+2a2+3a3 a2+3a3 a3 ö ÷ ÷ ÷ ÷ ø .

From this special case, we can see that the transpose of a transformation matrix resembles the property of Pascal triangle. Therefore, analogous to the construction of Pascal triangle we can derive the following algorithm to compute b_k=å_i=kⁿ C_i^ka_i:

   Copy a(i) into b(i)
   do k=0, n-1
   do i=n-1, k, -1
      b(i) = b(i) + b(i+1)
   enddo

   Copy a(i) into b(i)
   do k=0, n-1
   do i=n-1, k, -1
      b(i) = b(i) + b(i+1) * beta
   enddo

To demonstrate an application of the above method, we consider reparametrising a canonical polynomial

P(t)= n å i=0  ai ti,       t Î [a, b]

by t=[(b+a)+(b-a)u]/2 such that P(u) is defined over [-1,1]. Let a = (b-a)/2 and b = (b+a)/(b-a). Then, t=a(b+u). Therefore,

P(u)= n å i=0  ai ai(b+u)i.

By the theorem of binomial expansion, we have

P(u)= n å i=0  i å j=0  aiai Cij bj ui-j.

We want to find b_k such that

n å k=0  bk uk= n å i=0  i å j=0  aiai Cij bj ui-j.

Since the coefficients of the same degree terms on both side should be equal, we have

bk= n å i=k  i-k å j=i-k  aiai Cij bj = n å i=k  aiai Cii-k bi-k = n å i=k  Cik(aiai)bi-k.

Therefore, an efficient method to compute b_k would be

    b(0) = a(0)
    temp = alpha
    do i=1, n
       b(i) = a(i) * temp
       temp = temp * alpha
    enddo

    do k=0, n-1
    do i=n-1, k, step=-1
       b(i) = b(i) + b(i+1) * beta
    enddo

4  Compute b_k=å_i=0^k C_n-i^k-ia_i, b_k=å_i=0^k C_n-i^k-ia_ib^k-i

æ ç ç ç ç è b0 b1 b2 b3 ö ÷ ÷ ÷ ÷ ø = æ ç ç ç ç è 1 0 0 0 3 1 0 0 3 2 1 0 1 1 1 1 ö ÷ ÷ ÷ ÷ ø æ ç ç ç ç è a0 a1 a2 a3 ö ÷ ÷ ÷ ÷ ø = æ ç ç ç ç è a0 3a0+a1 3a0+2a1+a2 a0+a1+a2+a3 ö ÷ ÷ ÷ ÷ ø .

From this special case, we can see that the transpose of a transformation matrix resembles the property of Pascal triangle. Therefore, analogous to the construction of Pascal triangle, we give the following method to compute b_k:

   Copy a(k) into b(k)
   do k=0, n
   do i=1, n-k
      b(i) = b(i) + b(i-1)
   enddo

   Copy a(k) into b(k)
   do k=0, n
   do i=1, n-k
      b(i) = b(i) + b(i-1) * beta
   enddo

r(u)=r(t),    r¢(u)|u=0 = ar¢(t)|t=0.

Such reparametrisation is required whenever we want change the derivatives at the end points without affecting the geometric shape of the curve. In addition to the requirement of modifying the derivative at t=0, we also want to keep the degree of curve and the parameter interval unchanged. To meet all these requirements, the following bilinear transformation is a choice

t= au (a-1)u+1 = au bu+1        (u Î [0,1], t Î [0,1]).

By the chain rule of differentiation we have r^¢(u)=r^¢(t)t^¢(u). Then, it is readily checked that r^¢(u)|_u=0=ar^¢(t)|_t=0.

We denote r(t)=(x(t),y(t),z(t)) in homogeneous coordinate system by R(t)=(X(t),Y(t),Z(t),W(t)) and write R(t) in the canonical form R(t)=å_k=0ⁿa_k t^k, where a_k=(x_k,y_k,z_k,w_k) are vectors in homogeneous coordinates. Replacing t in R(t) by t=t(u), we have

R(u)= n å i=0  ai æ ç è au bu+1 ö ÷ ø i   = 1 (bu+1)n n å i=0  ai(au)i(bu+1)n-i.

Since R(u) is represented in homogeneous coordinates, the common factor 1/(bu+1)ⁿ can be discarded. Therefore, we have

R(u) = n å i=0  ai(au)i(bu+1)n-i.

By the binomial expansion theorem we obtain

R(u)= n å i=0  n-i å j=0  aiaiCn-ijbj ui+j.

We now want to find b_k such that

n å k=0  bkuk= n å i=0  n-i å j=0  aiai Cn-ijbj ui+j.

Since the same degree terms at both side should be equal, we can then derive

bk= k å i=0  k-i å j=k-i  aiaiCn-ijbj = k å i=0  aiaiCn-ik-ibk-i.

Let [`(a<sub>i</sub>)]=a<sub>i</sub>a<sup>i</sup>. Then, b<sub>k</sub>=å<sub>i=0</sub><sup>k</sup> C<sub>n-i</sub><sup>k-i</sup>[`(a_i)]b^k-i (k=0,1,¼,n), which can be computed very efficiently using the method we discussed earlier.

5  Conclusions

References

Davis, P.J. (1975), Interpolation & Approximation, Dover Pub., Inc.

Hoschek, J. and Lasser, D. (1993), Computer Aided Geometric Design, A.K. Peters, Ltd.

Faux, I.D. and Pratt, M.J. (1979), Computational Geometry for Design and Manufacture, Ellis Horwood Ltd., London.