# On the numerical evaluation of factorials IX

Number nine, who would have thought?!

As it is well known and laid down elsewhere, even by myself in a different post, it is quite easy to multiply and divide integers when their prime factorization is known. It is too complicated for almost all integers but it can be done very fast for factorials.
We have implemented algorithms to compute the binomial coefficient and the Catalan numbers but we did all of it inline and not as a generally applicable function. This generally applicable function needs to be able to handle negative exponents for example and we are going to implement one now.
Why now? Because of the Rule “One, two, many” and we used it three times already. The usage in the factorial function is excused for reasons of speed, so we are back down to two if we ignore the double-factorial. The next function I want to implement is the so called Narayana number:

${\displaystyle{ N(n,k) = \frac{1}{n}{ \binom nk}{\binom n{k-1}} }}$

Which can be written with factorials as

${\displaystyle{ \frac{1}{n} \; {{n!}\over{k!\,\left(n-k\right)!}} \;{{n!}\over{\left(k-1\right)!\,\left(n-k+1\right)!}} }}$

Shoving all together

${\displaystyle{ {{n!^2}\over{\left(k-1\right)!\,k!\,n\,\left(n-k\right)!\,\left(n-k +1\right)!}} }}$

It is most probably faster to use the first formula with the binomials but it is a good example with many factorials to multiply, divide and even use exponentiation (with an integer).

The first function is the one we need to decompose the factorial. It is not much different from the one in the factorial or binomial function but I took the chance to clean it up a little bit.

#include <tommath.h>
#ifndef LN_113
#define LN_113 1.25505871293247979696870747618124469168920275806274
#endif
#include <math.h>
int mp_factor_factorial(unsigned long n, unsigned long start,
long **prime_list,unsigned long *length){

unsigned long pix=0,prime,K;
mp_bitset_t *bst;

/* To be able to handle negative values in subtraction and to be sure that
n is positive*/
if(n > LONG_MAX){
return MP_VAL;
}

bst = malloc(sizeof(mp_bitset_t));
if (bst == NULL) {
return MP_MEM;
}
mp_bitset_alloc(bst, n+1);
mp_eratosthenes(bst);
pix = (unsigned long) (LN_113 *  n/log(n) );

*prime_list = malloc(sizeof(long) * (pix) * 2);
if (prime_list == NULL) {
return MP_MEM;
}
if(start <= 2){
prime = 2;
}
else{
prime = start;
}
K=0;
do {
(*prime_list)[K]   = (long)prime;
(*prime_list)[K+1] = (long)mp_prime_divisors(n,prime);
prime              = mp_bitset_nextset(bst,prime+1);
K+=2;
}while(prime <= n);
*prime_list = realloc(*prime_list,sizeof(long) * K);
if (prime_list == NULL) {
return MP_MEM;
}
free(bst);
*length = K;
return MP_OKAY;
}


The library we use is libtommath with my extensions available at my Github account and spread over this blog.
The function itself has n’t changed: make list of primes, walk over the primes and calculate for every prime the factorization. The result is put in a list in order and with the exponent following the prime, starting with the prime 2.
The mess of different types is due to a desperate try to keep all of the stuff compatible.

What next? Subtraction, OK.

These lists of primes and their exponents are nothing but simple polynomials and subtraction can be done likewise. One little exception: the polynomials are multi-variabled. That means that we need all primes in a consecutive row and we need as much as the bigger one of subtrahend and minuend respectively have.
It is a good idea to save some memory and skip any primes that have an exponent of zero as the result of the subtraction but only as long as the result doesn’t get squeezed another time through the function. The rule of consecutiveness is broken then, so we have to take care to check for such blank spots and fill them in if necessary.

The version without this checks first:

int mp_subtract_factored_factorials( long *subtrahend,
unsigned long l_subtrahend,
long *minuend,
unsigned long l_minuend,
long **difference,
unsigned long /* Vive */ *l_difference /*!*/
){
unsigned long k, counter=0;
unsigned long max_length = MAX(l_subtrahend,l_minuend );
unsigned long min_length = MIN(l_subtrahend,l_minuend );
long d,p,d_1,d_2;
/* TODO: check for sizes > 0 here */

/* Things work a bit different from ordinary arithmetic from now on */
*difference = malloc(sizeof(unsigned long)*(max_length+1));
if(*difference == NULL){
return MP_MEM;
}

/* Loop over smaller chunk first */
for(k=0;k<min_length;k+=2){
/* both have same length, we can take the value of the prime from any */
p = subtrahend[k];
/*d = subtrahend[k+1] - minuend[k+1];*/
d_1 = subtrahend[k+1];
d_2 = minuend[k+1];
/* handle over/underflow */
if((d_2 > 0) && (d_1 < LONG_MIN + d_2)){
return MP_VAL;
}
if((d_2 < 0) && (d_1 > LONG_MAX + d_2)){
return MP_VAL;
}
d = d_1 - d_2;
/* We only keep results that are different from zero */
if(d!=0){
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}
/* We need to find out which one is the smaller array and we have basically
two ways to approach the problem:
a) complicated and error-prone pointer juggling
b) just two loops
Here I have chosen b.
*/

/* If both arrays have the same length we can just stop here */
if(max_length == min_length){
/* We mad nothing dirty, so there's nothing to clean up here, let's just
grab our stuff and run */
return MP_OKAY;
}

/* If the subtrahend is bigger we subtract zeros, so simply copy */
if(l_subtrahend >= l_minuend ){
for(k=min_length;k<max_length;k+=2){
p = subtrahend[k];
d = subtrahend[k+1];
/* We still check, one never knows where things came from */
if(d!=0){
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}
}
/* If the minuend is bigger we subtract from zeros, so negate before copy*/
else{
for(k=min_length;k<max_length;k+=2){
p = minuend[k];
/* Yes, even negation can overflow */
d_1 = minuend[k+1];
if(d_1 == LONG_MIN){
return MP_VAL;
}
d = -d_1;
if(d!=0){
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}
}
/* Take care of the little difference the length might have */
*l_difference = counter;
return MP_OKAY;
}


Yes, I even checked for under/overflow this time 😉
To implement the test for consecutiveness we have to check every step. Let’s pick out the first loop and put a comment where it has to change

  for(k=0;k<min_length;k+=2){
/* Both have same length, we can take the value of the prime from any?
No, that is not true, we have to check and insert if none is there.
We don't have to insert every prime, just enough that both have the
same amount
*/
p_1 = subtrahend[k];
p_2 = minuend[k];
if(p_1 != p_2){
/* where is the gap? */
if(p_1 < p_2){
/* the gap is in the minuend, fill up subtrahend  */
}
if(p_1 > p_2){
/* the gap is in the subtrahend, fill up minuend */
}
}
/* handle difference in indices */
/*d = subtrahend[k+1] - minuend[k+1];*/
d_1 = subtrahend[k+1];
d_2 = minuend[k+1];
/* handle over/underflow */
if((d_2 > 0) && (d_1 < LONG_MIN + d_2)){
return MP_VAL;
}
if((d_2 < 0) && (d_1 > LONG_MAX + d_2)){
return MP_VAL;
}
d = d_1 - d_2;
/* We only keep results that are different from zero */
if(d!=0){
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}


As you might have guessed already, it is not very straight forward. We cannot simply fill-up the minuend because we have to subtract the exponents, so we have to fill-up with negated exponents. To keep up with the difference in the indices we have to know which one, subtrahend or minuend differs and how much. Which results in this (untested! I tend to be one off in these situations the first time, so please check with a debugger before you try it) loop:

  unsigned long gap_in_subtrahend = 0;
unsigned long gap_in_minuend    = 0;
for(k=0;k<min_length;k+=2){
p_1 = subtrahend[k+gap_in_minuend];
p_2 = minuend[k+gap_in_subtrahend];
if(p_1 != p_2){

if(p_1 < p_2){
for(i=0;p_1 < p_2;i+=2){
p_1 = subtrahend[k+i+gap_in_minuend];
(*difference)[counter]   = 1;
(*difference)[counter+1] = 0;
}
gap_in_minuend += i;
}
else{
for(i=0;p_1 > p_2;i+=2){
p_2 = minuend[k+i+gap_in_subtrahend];
d_2 = minuend[k+i+1+gap_in_subtrahend];
if(d_1 == LONG_MIN){
return MP_VAL;
}
d = -d_1;
(*difference)[counter]   = p_2;
(*difference)[counter+1] = d;
}
gap_in_subtrahend += 0;
}
}
/*d = subtrahend[k+1] - minuend[k+1];*/
d_1 = subtrahend[k+1+gap_in_minuend];
d_2 = minuend[k+1+gap_in_subtrahend];
/* handle over/underflow */
if((d_2 > 0) && (d_1 < LONG_MIN + d_2)){
return MP_VAL;
}
if((d_2 < 0) && (d_1 > LONG_MAX + d_2)){
return MP_VAL;
}
d = d_1 - d_2;
/* We only keep results that are different from zero */
if(d!=0){
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}


There is still amiss a check for the bounds of the participating arrays. And all we would gain is some saved memory. And how often does it happen? It does happen only if you manipulate the arrays by hand e.g.: divide by a small integer which means subtracting the prime factorization of a small integer or after a couple of divisions with factorials.
Two small examples:

${\displaystyle{ \frac{\left(\frac{20!}{10!}\right)}{10!} = 2^2 \cdot 3^0 \cdot 5^0 \cdot 7^0 \cdot 11^1 \cdot 13^1 \cdot 19^1 = 184\,756 }}$

And

${\displaystyle{ \frac{10!}{25} = 2^8 \cdot 3^4 \cdot 5^0 \cdot 7^1 = 145\,152 }}$

So it does happen and doesn’t happen rarely, we cannot ignore it.
Another one of the basic rules for programming is:”Make it work first, optimize later!” where “works unusable slow” is to be considered as non-working, too. (One example might be the too optimistic starting value for the nth-root function in libtommath which makes it unusable slow).
Skip the array compression? At what cost? Two examples

${\displaystyle{ \frac{20!}{19!} = 2^2 \cdot 3^0 \cdot 5^1 \cdot 7^0 \cdot 11^0 \cdot 13^0 \cdot 19^0 = 20 }}$

${\displaystyle{ \frac{20!}{18!} = 2^2 \cdot 3^0 \cdot 5^1 \cdot 7^0 \cdot 11^0 \cdot 13^0 \cdot 19^1 = 380 }}$

If we divide two factorials $\frac{m!}{n!}$ the gaps start to build if $n>\frac{m}{2}$ (proof omitted) and if we divide by a small integer $\frac{m!}{k}$ with $k \frac{}{m/2}$, it just needs a place to keep the original, un-evaluated factorials. Maybe at the end of the array with the prime-part set to one which does not change the behavior of the final multiplication ($1^x = 1$). It should be treated, because it would loop several times for nothing (max. 32 or 64 times depending on the bit size of the unsigned long  type but it does loop over the whole array every time).
So it is the check for $n>\frac{m}{2}$ we can fall back to if the result without it is too slow (unlikely, the heaviest burden computationally is the multiplication and squaring in the final multiplication, but cache misses can hurt. It is to be tested) or uses too much memory (more likely, because we don’t do this for small factorials but for large ones, too large to be handled otherwise, something in the range of $n > 10^6$ for $n!$).
So (I am such a so-so man 😉 ) without the memory optimizations fro know. That makes the code look a bit different. Not much but, well, a bit. We just strip the checks if the difference is zero and leave the rest untouched.

int mp_subtract_factored_factorials( long *subtrahend,
unsigned long l_subtrahend,
long *minuend,
unsigned long l_minuend,
long **difference,
unsigned long /* Vive */ *l_difference /*!*/
){
unsigned long k, counter=0;
unsigned long max_length = MAX(l_subtrahend,l_minuend );
unsigned long min_length = MIN(l_subtrahend,l_minuend );
long d,p,d_1,d_2;
/* TODO: check for sizes > 0 here */

/* Things work a bit different from ordinary arithmetic from now on */
*difference = malloc(sizeof(unsigned long)*(max_length+1));
if(*difference == NULL){
return MP_MEM;
}

/* Loop over smaller chunk first */
for(k=0;k<min_length;k+=2){
/* both have same length, we can take the value of the prime from any */
p = subtrahend[k];
/*d = subtrahend[k+1] - minuend[k+1];*/
d_1 = subtrahend[k+1];
d_2 = minuend[k+1];
/* handle over/underflow */
if((d_2 > 0) && (d_1 < LONG_MIN + d_2)){
return MP_VAL;
}
if((d_2 < 0) && (d_1 > LONG_MAX + d_2)){
return MP_VAL;
}
d = d_1 - d_2;
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
/* We need to find out which one is the smaller array and we have basically
two ways to approach the problem:
a) complicated and error-prone pointer juggling
b) just two loops
Here I have chosen b.
*/

/* If both arrays have the same length we can just stop here */
if(max_length == min_length){
/* We made nothing dirty, so there's nothing to clean up here, let's just
grab our stuff and run */
return MP_OKAY;
}

/* If the subtrahend is bigger we subtract zeros, so simply copy */
if(l_subtrahend >= l_minuend ){
for(k=min_length;k<max_length;k+=2){
p = subtrahend[k];
d = subtrahend[k+1];
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}
/* If the minuend is bigger we subtract from zeros, so negate before copy*/
else{
for(k=min_length;k<max_length;k+=2){
p = minuend[k];
/* Yes, even negation can overflow */
d_1 = minuend[k+1];
if(d_1 == LONG_MIN){
return MP_VAL;
}
d = -d_1;
(*difference)[counter]   = p;
(*difference)[counter+1] = d;
counter += 2;
}
}
/* this is now the same as max_length */
*l_difference = counter;
return MP_OKAY;
}


Addition follows the example of subtraction

int mp_add_factored_factorials(  long *summand_1,
unsigned long l_summand_1,
long *summand_2,
unsigned long l_summand_2,
long **sum,
unsigned long *l_sum){
unsigned long k, counter=0;
unsigned long max_length = MAX(l_summand_1,l_summand_2);
unsigned long min_length = MIN(l_summand_1,l_summand_2);
long s,p, s_1,s_2;
/* For more detailed comments see mp_subtract_factored_factorials() */
*sum = malloc(sizeof(unsigned long)*(max_length+1));
if(*sum == NULL){
return MP_MEM;
}
for(k=0;k<min_length;k+=2){
p = summand_1[k];
/* Over/underflow possible! */
/*s = summand_1[k+1] + summand_2[k+1];*/
s_1 = summand_1[k+1];
s_2 = summand_2[k+1];

if((s_2 > 0) && (s_1 > LONG_MAX - s_2)){
/* overflow */
return MP_VAL;
}
if((s_2 < 0) && (s_1 < LONG_MIN - s_2)){
/* underflow */
return MP_VAL;
}
s = s_1 + s_2;
(*sum)[counter]   = p;
(*sum)[counter+1] = s;
counter += 2;
}
if(l_summand_1 >= l_summand_2){
for(k=0;k<max_length;k+=2){
p = summand_1[k];
s = summand_1[k+1];
(*sum)[counter]   = p;
(*sum)[counter+1] = s;
counter += 2;
}
}
else{
for(k=0;k<max_length;k+=2){
p = summand_2[k];
s = summand_2[k+1];
(*sum)[counter]   = p;
(*sum)[counter+1] = s;
counter += 2;
}
}
*l_sum = counter;
return MP_OKAY;
}


Multiplication and division of the polynomials make not much sense, because of their multi-variability but multiplication with a small integer makes. Example

{\displaystyle{ \begin{aligned}20! &= 2^{18} \cdot 3^8 \cdot 5^4 \cdot 7^2 \cdot 11^1 \cdot 13^1 \cdot 19^1 \\ (20!)^2 &= 2^{36} \cdot 3^{16} \cdot 5^8 \cdot 7^4 \cdot 11^2 \cdot 13^2 \cdot 19^3 \end{aligned} }}

We need to check for over/underflow at multiplications here where we have to alternatives: use a data type that has twice the size of the data types of the multiplicands or do it the hard way.
The hard, but safe and portable way is to check every possible combination and test by trial division.

if (a > 0) {
if (b > 0) {
if (a > (SIGNED_TYPE_MAX / b)) {
return ERROR;
}
}
else {
if (b < (SIGNED_TYPE_MIN / a)) {
return ERROR;
}
}
}
else {
if (b > 0) {
if (a < (SIGNED_TYPE_MIN / b)) {
return ERROR;
}
}
else {
if ( (a != 0) && (b < (SIGNED_TYPE_MAX / a))) {
return ERROR;
}
}
}


The other alternative is to use long long and hope for the best or use the explicit data types in stdint.h throughout.
We use libtommath here where a data type mp_word is defined to be able to hold the multiplication of two mp_digit numbers but mp_digit can be as small as a char.

int mp_power_factored_factorials(  long *input,
unsigned long l_input,
long multiplicator,
long **product,
unsigned long *l_product){

unsigned long k, counter=0;

long prod,p, p_1,p_2;
#ifdef MP_USE_LONG_LONG_AND_HOPE_FOR_THE_BEST
long long temp;
#endif
*product = malloc(sizeof(unsigned long)*(l_input+1));
if(*product == NULL){
return MP_MEM;
}
p_2 = multiplicator;
for(k=0;k<l_input;k+=2){
p = input[k];
/* Over/underflow possible! */
/*prod = input[k+1] * multiplicator;*/
p_1 = input[k+1];
/* Two alternatives: use "long long" and hope for the best or do it the
hard way */
#ifdef MP_USE_LONG_LONG_AND_HOPE_FOR_THE_BEST
temp = p * (long long)p_1;
if ((temp > LONG_MAX) || (tmp < LONG_MIN)) {
return MP_VAL;
}
prod = (long) temp;
#else
if (p > 0) {
if (p_1 > 0) {
if (p > (LONG_MAX / p_1)) {
return MP_VAL;
}
}
else {
if (p_1 < (LONG_MIN / p)) {
return MP_VAL;
}
}
}
else {
if (p_1 > 0) {
if (p < (LONG_MIN / p_1)) {
return MP_VAL;
}
}
else {
if ( (p != 0) && (p_1 < (LONG_MAX / p))) {
return MP_VAL;
}
}
}
prod = p * p_1;
#endif
(*product)[counter]   = p;
(*product)[counter+1] = prod;
counter += 2;
}

*l_product = counter;
return MP_OKAY;
}


We can produce the reciprocal of a number if we negate all exponents of its prime factorization.

int mp_negate_factored_factorials(  long *input,
unsigned long l_input,
long **output,
unsigned long *l_output){
unsigned long k, counter=0;
long neg;

*output = malloc(sizeof(unsigned long)*(l_input+1));
if(*output == NULL){
return MP_MEM;
}
for(k=0;k<l_input;k+=2){
neg = input[k+1];
if(neg == LONG_MIN){
return MP_VAL;
}
(*output)[counter]   = input[k];
(*output)[counter+1] = -neg;
counter += 2;
}
*l_output = counter;
return MP_OKAY;
}


One can do the above inline, of course and save a copy but doing things inline is always a rich source of unexpected side-effects and an abundantly gushing spring of feet with gunshot injuries, so pleeease be careful.

The final step, computing the polynomials, is different from the function mp_compute_factored_factorials() in that we have to be able to handle negative exponents, too.
One simple thing to handle negative exponents is to keep them in a second array with absolute values of the exponents and divide at the end. Another alternative might be to return both arrays and let the user decide. No matter which way, we have to check and if there are at least one negative exponent decide what to do.

It is also the place where we have to decide what to do with the fact that the useful function mp_mul_d() accepts only mp_digit data types as the small integers. If we approximate the maximum exponent of a factorial as the exponent of two with a magnitude of roughly $n$ for $n!$ with the next exponents roughly half of the magnitude of their predecessor than the biggest possible exponent for the smallest possible mp_digit data type is 255 which makes the largest processable factorial $256!$ if I calculated it correctly (the biggest prime smaller than 255 is 251). That is not very much.

From here on I assume DIGIT_BIT == 28 and check for it. We will handle the exponent of two separately which will allow us to compute $n!$ for $n<2^{29} = 536\,870\,912$. The factorial of half a billion is about $1.6631e4\,453\,653\,131$, a number with nearly four and a half billion digits. Close to the age of the Earth in years (4.54 ± 0.05 billion according to Wikipedia), not so close to the age in minutes (3.15 billion, according to one James Ushher).
What, me and digressing? Ow, c'mon!

static long local_highbit(unsigned  long n){
long r=0;
while (n >>= 1) {
r++;
}
return r;
}
/*
To compute the resulting polynomials we need an approach different from
mp_compute_factored_factorials() because we have to be able to handle negative
exponents, too.
*/
int mp_multiply_factored_factorials(long *f, unsigned long f_length,
mp_int *c, mp_int *r){
unsigned long start=0,i,count=0,neg_count = 0,k=0;
long shift = 0;
long bit=0,neg_bit=0,hb=0,neg_hb=0,*p=NULL;
mp_int temp;
int e;

if(f[0] == 2){
shift = f[1];
if(f_length == 2){
if(shift > 0){
if( (e = mp_set_int(c,1LU<<(unsigned long)f[1]) ) != MP_OKAY){
return e;
}
}
return MP_OKAY;
}
}
if(shift){
start+=2;
k=2;
}
/*
This loop is expensive and should be done inline instead but for now...
*/
for(k=start;k<f_length;k+=2){
if(f[k+1] < 0){
/* count how many to save some memory */
count++;
}
/* Only primes up to 2^DIGIT_BIT because of mp_mul_d() */
if(f[k] >= 1<<DIGIT_BIT){
return MP_VAL;
}
}
/* all are negative */
//if(count && count == f_length/2){
/* if the alternative with two outputs has been chosen just skip the
computation of the polynomial with all positive exponents and set
*/
/* goto: negative_exponents; */

/* The other alternative would compute 1/x that gives always 0 in integer
divisions if x>1. But it would need exactly one exponent with zero which
has been filtered out already. It is debatable now if that was a good
decision.*/
//return MP_OKAY;
//}

if(count){
p = malloc( sizeof(long)* (count * 2) +1);
if(p == NULL){
return MP_MEM;
}
}

for(k=start;k<f_length;k+=2){
/* won't branch if count == 0 */
if(f[k+1] < 0){
p[neg_count]   = f[k];
p[neg_count+1] = abs(f[k+1]);
neg_count +=2 ;
/* And while we are already here ... */
neg_hb = local_highbit(abs(f[k+1]));
if(neg_bit < neg_hb)neg_bit = neg_hb;
/* This allows for some optimization, mainly w.r.t memory. Not done yet */
f[k]   = 1;
/* Set to zero to let the main loop skip it */
f[k+1] = 0;
}
else{
hb = local_highbit(f[k+1]);
if(bit < hb)bit = hb;
}
}

/* DIGIT_BIT can be as small as 8 */
if(bit >(long) DIGIT_BIT || neg_bit >(long) DIGIT_BIT){
return MP_VAL;
}
/* Anything times zero is zero, so set the result to one in advance */
if( (e = mp_set_int(c,1) ) != MP_OKAY){
return e;
}

/* Everything is in place, so let's start by computing with the positive
exponents first */

/* The other function mp_compute_factored_factorials() has a branch
but the cases where it makes sense are quite rare. Feel free to add it
yourself */
for( ; bit>=0; bit--) {
if( (e = mp_sqr(c, c) ) != MP_OKAY){
return e;
}
for(i=start; i<f_length; i+=2) {
if ((f[i+1] & (1LU << (unsigned long)bit)) != 0) {
/* This is problematic if DIGIT_BIT is too small. checked above */
if( (e = mp_mul_d(c,(mp_digit)f[i], c) ) != MP_OKAY){
return e;
}
}
}
}

/* Compute c * 2^n if n!=0 */
if(shift && shift > 0){
if(shift < 1<<DIGIT_BIT){
/* The choice of types is a bit questionable. Sometimes. */
if( (e = mp_mul_2d (c, (mp_digit) shift, c) ) != MP_OKAY){
return e;
}
}
else{
long multiplicator = 0;
another_round:
while(shift > 1<<DIGIT_BIT){
shift >>= 1;
multiplicator++;
}
if( (e = mp_mul_2d (c, (mp_digit) shift, c) ) != MP_OKAY){
return e;
}
if(multiplicator< DIGIT_BIT){
if( (e = mp_mul_2d (c, (mp_digit)(1<<multiplicator), c) ) != MP_OKAY){
return e;
}
}
else{
shift = 1<<multiplicator;
multiplicator = 0;
goto another_round;
}
}
}

/* None are negative*/
if(count == 0){
/* Clean up */
/* Nothing to clean up */
return MP_OKAY;
}

/* Compute with the negative eponents */
if( (e = mp_init(&temp) ) != MP_OKAY){
return e;
}
if( (e = mp_set_int(&temp,1) ) != MP_OKAY){
return e;
}
for( ; neg_bit>=0; neg_bit--) {
if( (e = mp_sqr(&temp, &temp) ) != MP_OKAY){
return e;
}
/* The exponents of 2 have been stripped already so we start at zero.
"count" counts only the occurrences, the size itself is twice as large.
*/
for(i=0; i<count*2; i+=2) {
if ((p[i+1] & (1LU << (unsigned long)neg_bit)) != 0) {
if( (e = mp_mul_d(&temp,(mp_digit)p[i], &temp) ) != MP_OKAY){
return e;
}
}
}
}

/* Compute c * 2^n if n!=0  for negative n*/
if(shift && shift < 0){
/* check for overflow */
if(shift == LONG_MIN){
return MP_VAL;
}
shift = -shift;
if(shift < 1<<DIGIT_BIT){
/* The choice of types is a bit questionable. Sometimes. */
if( (e = mp_mul_2d (&temp, (mp_digit) shift, &temp) ) != MP_OKAY){
return e;
}
}
else{
long multiplicator = 0;
another_round:
while(shift > 1<<DIGIT_BIT){
shift >>= 1;
multiplicator++;
}
if( (e = mp_mul_2d (&temp, (mp_digit) shift, &temp) ) != MP_OKAY){
return e;
}
if(multiplicator< DIGIT_BIT){
if( (e = mp_mul_2d (&temp, (mp_digit)(1<<multiplicator), &temp) ) != MP_OKAY){
return e;
}
}
else{
shift = 1<<multiplicator;
multiplicator = 0:
goto another_round;
}
}
}
#ifdef BN_MP_DO_LAST_DIVISION
if( (e = mp_div(c,&temp,c,r) ) != MP_OKAY){
return e;
}
#else
if(r != NULL){
if( (e = mp_copy(&temp,r) ) != MP_OKAY){
return e;
}
}
#endif
free(p);
mp_clear(&temp);
return MP_OKAY;
}


It compiles and works for small inputs but it cannot be called well tested. Use with car and drop me a note if you found a bug.

How to compute Narayana numbers? As I said in the beginning: computing them with the formula

${\displaystyle{ N(n,k) = \frac{1}{n}{ \binom nk}{\binom n{k-1}} }}$

would most probably be the fastest way and if not, it is the easiest to implement at least.

int mp_narayana(unsigned long n, unsigned long k, mp_int *c){
mp_int temp;
int e;
if( (e = mp_init(&temp) ) != MP_OKAY){
return e;
}
if( (e = mp_binomial(n,k,c) ) != MP_OKAY){
return e;
}
if( (e = mp_binomial(n,k-1,&temp) ) != MP_OKAY){
return e;
}
if( (e = mp_mul(c,&temp,c) ) != MP_OKAY){
return e;
}
if( (e = mp_set_int(&temp,n) ) != MP_OKAY){
return e;
}
if( (e = mp_div(c,&temp,c,NULL) ) != MP_OKAY){
return e;
}
mp_clear(&temp);
return MP_OKAY;
}


The above code will be merged as mp_narayana() into my fork of libtommath in the next days. The rest only after some thorough tests.