Faster Factorial

[FirstName LastName]

The following code implements a binary splitting method to finding a
pochhammer (and the factorial as a specific case).  It balances the
multiplications more effectively than mpz/fac_ui.c


/* pochhammer.c */

#include <stdio.h>
#include <time.h>

#include <gmp.h>


#define MAX_THRESHOLD    14

void pochhammer_basecase(mpz_ptr d, const size_t n0,  const size_t n1,
const
size_t s)
{
 size_t i;
size_t t;


t=1<<s;
 mpz_set_ui(d,n0);
 for(i=n0+t;i<=n1;i+=t)
 {
  mpz_mul_ui(d,d,i);
 }

 return;
}

void pochhammer_recursive(mpz_ptr d, const size_t n0,  const size_t n1,
const size_t s)
{
 mpz_t d2;


 if(((n1-n0)>>s) < MAX_THRESHOLD)
 {
 pochhammer_basecase( d, n0, n1, s);

 return;
 }

 pochhammer_recursive( d, n0, n1, s+1);

 mpz_init(d2);
 pochhammer_recursive( d2, n0+(1<<s), n1, s+1);

 mpz_mul(d,d,d2);
 mpz_clear(d2);

return ;
}

void pochhammer(mpz_ptr d, const size_t n0, const size_t n1)
{
 pochhammer_recursive(d, n0, n1, 0);

return;
}