Calculate binominal coefficient

bico.c

/*
 * calculate binominal coefficient
 * 
 * gcc -std=c99 -g0 -O3 -s bico.c -o bico -lm
 * 
 * (c) 2017 by Heiko Schaefer <[email protected]>
 *
 */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double gammln(double xx) {
  
  double x = xx, y = xx, tmp = x + 5.5, ser = 1.000000000190015;
  
  static const double cof[6] = {
      76.18009172947146,
      -86.50532032941677,
      24.01409824083091,
      -1.231739572450155,
      0.1208650973866179e-2,
      -0.5395239384953e-5
  };
  
  tmp -= (x + 0.5) * log(tmp);
  
  for(int j = 0; j <= 5; ++j) {
    ser += cof[j]/++y;
  }
  
  return -tmp + log(2.5066282746310005 * ser/x);
}

static double factln(int n) {
  
  static double a[101];

  if (n < 0) {
    return NAN;
  }
  
  if (n <= 1) {
    return 0.0;
  }
  
  if (n <= 100) {
    return a[n] ? a[n] : (a[n] = gammln(n + 1.0));
  } else {
    return gammln(n + 1.0);
  }
}

inline static double nr_bico(int n, int k) {
  return floor(0.5 + exp(factln(n) - factln(k) - factln(n - k)));
}

inline static unsigned long ef_bico(unsigned long n, unsigned long k) {
  
  if(k == 0) {
    return 1;
  }
  
  if(2 * k > n) {
    k = n - k;
  }
  
  unsigned long bico = 1;
  
  for(unsigned long i = 1; i <= k; ++i) {
    bico = bico * (n - k + i)/i;
  }
  
  return bico;
}

int main(int argc, char *argv[]) {
  
  int n = atoi(argv[1]), k = atoi(argv[2]);;
  
  printf("nr_bico(%d, %d): %lu\n", n, k, (unsigned long)nr_bico(n, k));
  printf("ef_bico(%d, %d): %lu\n", n, k, (unsigned long)ef_bico(n, k));
  
  return EXIT_SUCCESS;
}