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
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:
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).
Back to the Homepage of Factorial Algorithms. 2010/Feb/08.