Advanced Search

Calculate binomial coefficient

/ Published in: C
Save to your folder(s)
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.

Expand | Embed | Plain Text

Copy this code and paste it in your HTML 
  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. }

Report this snippet


Comments

RSS Icon Subscribe to comments

You need to login to post a comment.