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:  // Differs only in a single function call from Schoenhage::Factorial.

   7:   

   8:  #include "schoenhage.h"

   9:  #include "xmath.h"

  10:   

  11:  void Schoenhage::ParallelFactorial(Xint result, ulong n)

  12:  {

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

  14:   

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

  16:   

  17:      slong iterLen = 0, m = n;

  18:      ulong lim = n / 2, i;

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

  20:   

  21:      ulong* primes = 0; 

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

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

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

  25:   

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

  27:   

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

  29:      {

  30:          ulong p = primes[i];

  31:          ulong exp = 1;

  32:          if (p <= lim)

  33:          {

  34:              slong N = n / p;

  35:              exp = N;

  36:              while (N > 0)

  37:              {

  38:                  N /= p;

  39:                  exp += N;

  40:              }

  41:          }

  42:          exponents[i] = exp;

  43:      }

  44:   

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

  46:   

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

  48:      {

  49:          ulong count = 0;

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

  51:          {

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

  53:              {

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

  55:              }

  56:          }

  57:   

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

  59:          if (count > 0)

  60:          {

  61:              Xmath::ParallelProduct(prod, factors, 0, count);

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

  63:          }

  64:      }

  65:   

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

  67:   

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

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

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

  71:      lmp::Clear(prod);

  72:  }