Adventures of a Programmer: Parser Writing Peril IX

Now, that we hae a function to compute the binom…(what? The grammar? Well)…ial coefficient we can use it to compute the Catalan numbers.

{\displaystyle{ \mathrm{C}_n = \frac{\binom {2n}{n}}{n+1}   }}

The term 2n in the binomial coefficient needs some care and even the n+1 causes problems if it is larger than 1<<DIGIT_BIT because the function mp_div_d() only accepts mp_digit as the divisor. These thoughts result in the following short function:

int mp_catalan(unsigned long n, mp_int *c){
  int e;
  if(n == 0 || n == 1){
    if( (e = mp_set_int(c,1LU) ) != MP_OKAY){
      return e;
    }
    return MP_OKAY;
  }
  /*  C(2^31) is already a very large number
      C(2^31) ~ 1.7593e1,292,913,972
      The author's box would need approximately 1 year and seven month
      computing time.
      The number of digits follow loosely the size of n here, e.g.
      C(2^32)  ~ 1.9303e2,585,827,958
      C(2^64)  ~ 2.5891e11,106,046,577,046,714,235
      C(2^128) ~ 1.2728e204,870,398,877,478,727,500,024,218,500,206,465,342
   */
  if(n >= ULONG_MAX/2){
    return MP_VAL;
  }
  if( (e = mp_binomial(2*n,n,c) ) != MP_OKAY){
    return e;
  }
  if((n+1) >= 1<<DIGIT_BIT){
    mp_int temp;
    if( (e = mp_init(&temp) ) != MP_OKAY){
      return e;
    }
    if( (e = mp_set_int(&temp,n+1) ) != MP_OKAY){
      return e;
    }
    if( (e = mp_div(c,&temp,c,NULL) ) != MP_OKAY){
      return e;
    }
    mp_clear(&temp);
  }
  else{
    if( (e = mp_div_d(c,n+1,c,NULL) ) != MP_OKAY){
      return e;
    }
  }
  return MP_OKAY;
}

It has been merged into my fork of libtommath.

Further optimizations possible? Yes.
The binomial coefficient is the central binomial coefficient with the interesting equality:

{\displaystyle{ \binom {2n}{n} = \frac{\left(2n\right)!}{\left(n!\right)^2}  }}

which means that we can skip one part of the calculation of the final prime factorization.

  1. {\displaystyle{  f_{2n} - \left( f_n + f_n\right) } }
  2. {\displaystyle{  f_{2n} } }

The following snippet might be of use for a better understanding.

  prime = 2;K=0;
  do{
    diff = mp_prime_divisors(2*n,prime)-(2*mp_prime_divisors(n,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);
  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 <= 2*n);

The time spend in calculating the prime factorization is negligible in contrast to the final computation but all of the exponents for the primes larger than n are one for 2n, so it might be a good idea to use binary splitting for this primorial. Former experiments showed that the cutoff is above the factorial of two hundred thousand (2,000,000!) for my machine.

Advertisements

One thought on “Adventures of a Programmer: Parser Writing Peril IX

  1. Pingback: On the numerical evaluation of factorials IX | De Amentiae Mundi

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