/cln-1.3.2/src/integer/misc/combin/cl_I_binomial.cc

 1// binomial().
 2
 3// General includes.
 4#include "base/cl_sysdep.h"
 5
 6// Specification.
 7#include "cln/integer.h"
 8
 9
10// Implementation.
11
12#include "integer/misc/combin/cl_I_combin.h"
13
14namespace cln {
15
16const cl_I binomial (uintL n, uintL k)
17{
18	// Method:
19	// Ensure 0 <= k <= n: If n<k, return 0.
20	// Ensure 0 <= k <= n/2: If k>n/2, replace k by n-k.
21	// Compute the product (n-k+1)...n and divide by k!.
22	// When computing the product m < i <= n, split into the odd part
23	// and the even part. The even part is 2^(ord2(n!)-ord2(m!)),
24	// and recall that ord2(n!) = n - logcount(n). The odd part is
25	// computed recursively: It is the product of the odd part of
26	// m/2 < i <= n/2, times the product of m < 2*i+1 <= n. The
27	// recursion goes until floor(m/2^j) = floor(n/2^j) or floor(n/2^j) < 2.
28	if (n < k)
29		return 0;
30	// Now 0 <= k <= n.
31	if (2*k > n)
32		k = n-k;
33	// Now 0 <= k <= n/2.
34	var cl_I prod = 1;
35	var uintL m = n-k;
36	var uintL j = 0;
37	{
38		var uintL a = m;
39		var uintL b = n;
40		while (a < b && b > 1) {
41			a = floor(a,2); b = floor(b,2);
42			j++;
43		}
44	}
45	while (j > 0) {
46		j--;
47		var uintL a = m>>j;
48		var uintL b = n>>j;
49		// Compute product(a < i <= b and i odd, i)
50		a = floor(a-1,2);
51		b = floor(b-1,2);
52		// = product(a < i <= b, 2i+1)
53		if (a < b)
54			prod = prod * cl_I_prod_ungerade(a,b);
55	}
56	prod = prod << (k + logcount(m) - logcount(n));
57	return exquopos(prod,factorial(k));
58}
59
60}  // namespace cln