Posted By

s1m55r on 06/25/10

Calculate binomial coefficient

/ Published in: C   Program to calculate large binomial coefficients efficiently. As long as the value fits into an unsigned long long int, the program should find it in a relatively short amount of time.

1. /**
2.  * Program to calculate large binomial coefficients
3.  * efficiently. As long as the value fits into an
4.  * unsigned long long int, the program should find
5.  * it in a relatively short amount of time.
6.  *
7.  * Author: Sameer Vijaykar <[email protected]>
8.  *
9.  * References:
10.  * http://blog.plover.com/math/choose.html
11.  * http://en.wikipedia.org/wiki/Binary_GCD_algorithm
12.  *
13.  */
14.
15. #include<stdio.h>
16.
17. typedef unsigned long long int LL;
18.
19. /**
20.  * Sourced from Wikipedia
21.  * http://en.wikipedia.org/wiki/Binary_GCD_algorithm#Implementation_in_C
22.  */
23. LL gcd(LL u, LL v)
24. {
25. int shift;
26.
27. /* GCD(0,x) := x */
28. if (u == 0 || v == 0)
29. return u | v;
30.
31. /* Let shift := lg K, where K is the greatest power of 2
32.   dividing both u and v. */
33. for (shift = 0; ((u | v) & 1) == 0; ++shift) {
34. u >>= 1;
35. v >>= 1;
36. }
37.
38. while ((u & 1) == 0)
39. u >>= 1;
40.
41. /* From here on, u is always odd. */
42. do {
43. while ((v & 1) == 0) /* Loop X */
44. v >>= 1;
45.
46. /* Now u and v are both odd, so diff(u, v) is even.
47.   Let u = min(u, v), v = diff(u, v)/2. */
48. if (u < v) {
49. v -= u;
50. } else {
51. LL diff = u - v;
52. u = v;
53. v = diff;
54. }
55. v >>= 1;
56. } while (v != 0);
57.
58. return u << shift;
59. }
60.
61. int main()
62. {
63. int n, k, g, i, i1;
64. LL f;
65.
66. printf("Enter integers n and k: ");
67. scanf("%d", &n);
68. scanf("%d", &k);
69.
70. // Check for trivial cases
71. if ((n == k) || (n == 1) || (k == 1))
72. {
73. printf("nCk : 1\n");
74. return 0;
75. }
76.
77. /* Beyond this check, everything going into the
78.   loop will at least be 1 */
79.
80. f = 1;
81.
82. /* Pick the smaller denominator */
83. k = (k < (n-k) ? k : (n-k));
84.
85. for (i = 1; i <= k; i++)
86. {
87. // Compute f = f * (n/i)
88.
89. if (n % i == 0)
90. {
91. f *= (n/i);
92. }
93. else if (f % i == 0)
94. {
95. f = (f/i) * n;
96. }
97. else
98. {
99. g = gcd(f, i);
100. f /= g;
101. i1 = i/g;
102. if (n % i1 == 0)
103. {
104. f *= (n/i1);
105. }
106. else
107. {
108. f = (f*n) / i1;
109. }
110. }
111.
112. n--;
113. }
114.
115. printf("nCk : %lld\n", f);
116. } Subscribe to comments