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.

The term in the binomial coefficient needs some care and even the 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:

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

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 are one for , 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.

### Like this:

Like Loading...

*Related*

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