# 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.