1:  // Author & Copyright : Peter Luschny

   2:  // Created: 2010-01-15

   3:  // License: LGPL version 2.1 or (at your option)

   4:  // Creative Commons Attribution-ShareAlike 3.0

   5:   

   6:  #include "schoenhage.h"

   7:  #include "xmath.h"

   8:   

   9:  void Schoenhage::Factorial(Xint result, ulong n)

  10:  {

  11:      lmp::SetUi(result, 1);

  12:   

  13:      if (n < 2) { return; }

  14:   

  15:      slong iterLen = 0, m = n;

  16:      ulong lim = n / 2, i;

  17:      while (m > 0) { m >>= 1; iterLen++; }

  18:   

  19:      ulong* primes = 0; 

  20:      ulong piN = Xmath::PrimeSieve(&primes, n);

  21:      ulong* exponents = lmp::MallocUi(piN);

  22:      ulong* factors = lmp::MallocUi(piN);

  23:   

  24:      exponents[0] = n - Xmath::BitCount(n);

  25:   

  26:      for(i = 1; i < piN; i++)

  27:      {

  28:          ulong p = primes[i];

  29:          ulong exp = 1;

  30:          if (p <= lim)

  31:          {

  32:              slong N = n / p;

  33:              exp = N;

  34:              while (N > 0)

  35:              {

  36:                  N /= p;

  37:                  exp += N;

  38:              }

  39:          }

  40:          exponents[i] = exp;

  41:      }

  42:   

  43:      Xint prod; lmp::Init(prod);

  44:   

  45:      for (; iterLen >= 0; iterLen--)

  46:      {

  47:          ulong count = 0;

  48:          for (ulong m = 1; m < piN ; m++)

  49:          {

  50:              if (((exponents[m] >> iterLen) & 1) == 1)

  51:              {

  52:                  factors[count++] = primes[m];

  53:              }

  54:          }

  55:   

  56:          lmp::Pow2(result, result); 

  57:          if (count > 0)

  58:          {

  59:              Xmath::Product(prod, factors, 0, count);

  60:              lmp::Mul(result, result, prod);

  61:          }

  62:      }

  63:   

  64:      lmp::Mul2Exp(result, result, exponents[0]);

  65:   

  66:      lmp::FreeUi(factors, piN);

  67:      lmp::FreeUi(primes, piN);

  68:      lmp::FreeUi(exponents, piN);

  69:      lmp::Clear(prod);

  70:  }

  71: