On the numerical evaluation of factorials VI

Another addendum to the last post.
After playing around with the implementation of the Borwein algorithm I found that the small multiplication eats a lot of time, most of the all-over running time to be a bit more specific. So I added the splitting of the evaluation between the primes with an exponent larger than one and those with one that is one again. Now the upper half of the primes is now treated as a primorial and evaluated with binary splitting to a good effect.

define factorial(n){
  local prime result shift prime_list k k1 k2 expo_list pix cut;
  local s1 e1 r1 s2 e2 r2 primorial; 
    
  result = 1;
  prime      = 2;

  if(!isint(n))  return "false"; ## or gamma(n)?
  if(n < 0)      return "false";
  if(n < 9000)  return n!;

  shift = __CZ__prime_divisors(n,prime);
  prime = 3;
  cut = n//2;
  pix   = pix(cut);
  prime_list = mat[pix];
  expo_list  = mat[pix];
  k = 0;
  do {
    prime_list[k]  = prime;
    expo_list[k++] = __CZ__prime_divisors(n,prime);
    prime          = nextprime(prime);
  }while(prime <= cut);

  /* size of the largest exponent in bits  */
  k1 = highbit(expo_list[0]);
  k2 = size(prime_list)-1;

  for(;k1>=0;k1--){
    /*
       the cut-off for T-C-4 is still to low, using T-C-3 here
       TODO: check cutoff of T-C-4
     */
    result = toomCook3square(result);
    for(k=0; k<=k2; k++) {
      if((expo_list[k] & (1 << k1)) != 0) {
        result *= prime_list[k];
      }
    }

  }
  primorial = primorial( cut, n);
  ##print digits(result,2), digits(primorial,2);
  ##result = toomCook3(result,primorial);
  result *= primorial;
  result <<= shift;
  return result;
}

As you can see, the point where the builtin factorial is faster got down from 21,000 to 9,000. I also used fixed size arrays instead of the lists.

A little bug was in the construction of the prime list, too and gets repaired with the follwing replacement:

/* A simple list to keep things...uhm...simple?*/
global __CZ__primelist;
/* Helper for primorial: fill list with primes in range a,b */
define __CZ__fill_prime_list(a,b)
{
  local k;
  k=a;
  /* reset list every time or cache it */
  __CZ__primelist = list();
  if(isprime(k))k--;
  while(1){
    k = nextprime(k);
    if(k > b) break;
    append(__CZ__primelist,k );
  }

Next: the double factorial.