The Binomial (n over k)

The binomial function computes the binomial coefficients. If the arguments are both positive integers with 0 <= k <= n, then

binomial(n, k) = n! / (k! * (n-k)!)

which is the number of ways of choosing k objects from n distinct objects. Otherwise an exception is thrown.

For each prime p <= n, we compute the power of pe dividing the binomial coefficient binomial(n, k). We use the fact that e is the number of borrows in the subtraction n-k in base p. 

See P. Goetgheluck, Computing Binomial Coefficients, American Math. Monthly 94 (1987), 360-365.

However, our implementation pays special attention to the memory consumption. In fact it uses only PrimePi(n) longs. As far as I know there is no mention in the literature to a binomial algorithm with this feature.

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.* equivalents.

Testing binomial: (400000,133333)
ParallelBinomial:       0.015 sec
PrimeBinomial:          0.026 sec
GMP5Binomial:           4.938 sec

Testing binomial: (1600000,533333)
ParallelBinomial:       0.078 sec
PrimeBinomial:          0.125 sec
GMP5Binomial:          91.063 sec

Testing binomial: (6400000,2133333)
ParallelBinomial:        0.422 sec
PrimeBinomial:           0.610 sec
GMP5Binomial:         1733.750 sec

C# (similar JAVA)

   1:  /// Author & Copyright : Peter Luschny
   2:  /// Created: 2010-02-01
   3:  /// License: LGPL version 2.1 or (at your option)
   4:  /// Creative Commons Attribution-ShareAlike 3.0
   5:   
   8:  using System;
   9:  using System.Numerics;
  10:  using Math.Primes;
  11:  using FastFactorial;
  12:  using MathUtils;
  13:   
  14:  namespace Math.FastFactorial {
  15:   
  16:  class Program {
  17:   
  18:  static BigInteger NaiveBinomial(int n, int k)
  19:  {
  20:      var f = new FactorialPrimeSwing();
  21:      var fn = f.Factorial(n);
  22:      var fk = f.Factorial(k);
  23:      var fnk = f.Factorial(n - k);
  24:      return fn / (fk * fnk);
  25:  }
  26:   
  27:  static void Main()
  28:  {
  29:      for(int n=0; n < 100;n++)
  30:          for (int k = 0; k <= n; k++)
  31:          {
  32:              BigInteger b = NaiveBinomial(n, k);
  33:              BigInteger c = FastBinomial.Binomial(n, k);
  34:              if (b != c)
  35:              {
  36:                  Console.WriteLine("NB(" + n + "," + k + ") = " + b);
  37:                  Console.WriteLine("FB(" + n + "," + k + ") = " + c);
  38:                  Console.ReadLine();
  39:              }
  40:              Console.WriteLine("C(" + n + "," + k + ") = " + c);
  41:          }
  42:          Console.ReadLine();
  43:  }}
  44:   
  45:   
  46:  public class FastBinomial {
  47:   
  48:  private FastBinomial() { }
  49:   
  50:  public static BigInteger Binomial(int n, int k)
  51:  {
  52:      if (0 > k || k > n)
  53:      {
  54:          throw new System.ArgumentException(
  55:           "Binomial: 0 <= k and k <= n required, but n was "
  56:           + n + " and k was " + k);
  57:      }
  58:   
  59:      if ((k == 0) || (k == n)) return BigInteger.One;
  60:      if (k > n / 2) { k = n - k; }
  61:      int fi = 0, nk = n - k;
  62:   
  63:      int rootN = (int)System.Math.Floor(System.Math.Sqrt(n));
  64:      int[] primes = new PrimeSieve(n).GetPrimeCollection(2, n).ToArray();
  65:   
  66:      foreach (int prime in primes) // Equivalent to a nextPrime() function.
  67:      {
  68:          if (prime > nk)
  69:          {
  70:              primes[fi++] = prime;
  71:              continue;
  72:          }
  73:   
  74:          if (prime > n / 2)
  75:          {
  76:              continue;
  77:          }
  78:   
  79:          if (prime > rootN)
  80:          {
  81:              if (n % prime < k % prime)
  82:              {
  83:                  primes[fi++] = prime;
  84:              }
  85:              continue;
  86:          }
  87:   
  88:          int r = 0, N = n, K = k, p = 1;
  89:   
  90:          while (N > 0)
  91:          {
  92:              r = (N % prime) < (K % prime + r) ? 1 : 0;
  93:              if (r == 1)
  94:              {
  95:                  p *= prime;
  96:              }
  97:              N /= prime;
  98:              K /= prime;
  99:          }
 100:          if (p > 1) primes[fi++] = p;
 101:      }
 102:   
 103:      return Xmath.Product(primes, 0, fi);
 104:  }
 105:  }}
 106:   

C++ with GMP or MPIR

 108:   
 109:  #include <assert.h>
 110:  #include "binomial.h"
 111:  #include "xmath.h"
 112:   
 113:  void Binomial::PrimeBinomial(Xint binom, ulong n, ulong k)
 114:  {
 115:      assert(k <= n);  // TODO: Error handling
 116:   
 117:      lmp::SetUi(binom, 1);
 118:      if ((k == 0) || (k == n)) return;
 119:   
 120:      if (k > n / 2) { k = n - k; }
 121:      ulong fi = 0, nk = n - k;
 122:   
 123:      ulong rootN = Xmath::Sqrt(n);
 124:      ulong* primes;
 125:      ulong piN = Xmath::PrimeSieve(&primes, n);
 126:   
 127:      for(ulong i = 0; i < piN; i++) 
 128:      {
 129:          ulong prime = primes[i];
 130:   
 131:          if (prime > nk)
 132:          {
 133:              primes[fi++] = prime;
 134:              continue;
 135:          }
 136:   
 137:          if (prime > n / 2)
 138:          {
 139:              continue;
 140:          }
 141:   
 142:          if (prime > rootN)
 143:          {
 144:              if (n % prime < k % prime)
 145:              {
 146:                  primes[fi++] = prime;
 147:              }
 148:              continue;
 149:          }
 150:   
 151:          ulong N = n, K = k, p = 1, r = 0;
 152:   
 153:          while (N > 0)
 154:          {
 155:              r = (N % prime) < (K % prime + r) ? 1 : 0;
 156:              if (r == 1)
 157:              {
 158:                  p *= prime;
 159:              }
 160:              N /= prime;
 161:              K /= prime;
 162:          }
 163:          if (p > 1) primes[fi++] = p; 
 164:      }
 165:   
 166:      Xmath::Product(binom, primes, 0, fi);
 167:      free(primes);
 168:  }

       The c++ sources, which can also be browsed here (see the c++ and h tags).
previous Back to the Homepage of Factorial Algorithms.     2010/Feb/08.