On the numerical evaluation of factorials XII

Falling Factorial

The algorithm for the falling factorial (x)_n makes shamelessly use of \frac{(x)_n}{n!}  = \binom xn which can be written in terms of factorials as (x)_n  =  \frac{x!}{n! (x-n)!}n! = \frac{x!}{(x-n)!} which is our algorithm for the binomial coefficients without the first loop.
Oh, that was simple!

#ifndef LN_113
#   define LN_113 1.25505871293247979696870747618124469168920275806274
#endif
#include <math.h>

int mp_falling_factorial(unsigned long n, unsigned long k, mp_int * c)
{
  unsigned long *prime_list;
  unsigned long pix = 0, prime, K, diff;
  mp_bitset_t *bst;

  int e;

  if (n < k) {
    mp_set(c, 0);
    return MP_OKAY;
  }
  if (k == 0) {
    mp_set(c, 1);
    return MP_OKAY;
  }
  if (k == 1) {
    if ((e = mp_set_int(c, n)) != MP_OKAY) {
      return e;
    }
    return MP_OKAY;
  }
  if (k == n) {
    if ((e = mp_factorial(n, c)) != MP_OKAY) {
      return e;
    }
    return MP_OKAY;
  }

  bst = malloc(sizeof(mp_bitset_t));
  if (bst == NULL) {
    return MP_MEM;
  }
  mp_bitset_alloc(bst, n + 1);
  mp_eratosthenes(bst);
  /* One could also count the number of primes in the already filled sieve */
  pix = (unsigned long) (LN_113 * n / log(n))+2;

  prime_list = malloc(sizeof(unsigned long) * (pix) * 2);
  if (prime_list == NULL) {
    return MP_MEM;
  }
  prime = 2;
  K = 0;
  do {
    diff = mp_prime_divisors(n, prime) - mp_prime_divisors(n - k, prime);
    if (diff != 0) {
      prime_list[K] = prime;
      prime_list[K + 1] = diff;
      K += 2;
    }
    prime = mp_bitset_nextset(bst, prime + 1);
  } while (prime <= n - k);
  do {
    prime_list[K] = prime;
    prime_list[K + 1] = mp_prime_divisors(n, prime);
    prime = mp_bitset_nextset(bst, prime + 1);
    K += 2;
  } while (prime <= n);
  prime_list = realloc(prime_list, sizeof(unsigned long) * K);
  if (prime_list == NULL) {
    return MP_MEM;
  }
  if ((e = mp_compute_factored_factorial(prime_list, K - 1, c, 0)) != MP_OKAY) {
    return e;
  }
  free(bst);
  free(prime_list);

  return MP_OKAY;
}

But if you think that was simple, wait for the next post đŸ˜‰

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s