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 "primeswing.h"

   7:  #include "xmath.h"

   8:   

   9:  void PrimeSwing::ParallelFactorial(Xint fact, ulong n)

  10:  {

  11:      if (n < THRESHOLD) { Xmath::NaiveFactorial(fact, n); return; }

  12:      lmp::SetUi(fact, 1);

  13:   

  14:      ulong* primes;

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

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

  17:      ulong iterLen = 0; ulong i = n;

  18:      slong m = n; while (m > 0) { m >>= 1; iterLen++; }

  19:      ulong* iter = lmp::MallocUi(iterLen);

  20:      ulong* lim  = lmp::MallocUi(iterLen);

  21:      m = n; i = iterLen;

  22:      while (m > 0) { iter[--i] = m; m >>= 1; }

  23:   

  24:      Xint primorial; lmp::Init(primorial);

  25:      Xint swing; lmp::Init(swing);

  26:   

  27:      lim[0] = 0;

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

  29:          lim[i] = iter[i] < SOSLEN / 2 ?  0 :

  30:          GetIndexOf(primes, iter[i], lim[i-1], piN);

  31:   

  32:      for (i = 1; i < iterLen; i++)

  33:      {

  34:          ulong N = iter[i];

  35:   

  36:          if (N < SOSLEN)

  37:          {

  38:              lmp::Pow2(fact, fact);

  39:              lmp::SetUi(swing, smallOddSwing[N]);

  40:          }

  41:          else

  42:          {

  43:              boost::thread pow(lmp::Pow2, fact, fact);

  44:   

  45:              ulong prime = 3;

  46:              slong pi = 2, fi = 0;

  47:              ulong max = Xmath::Sqrt(N);

  48:   

  49:              while (prime <= max)

  50:              {

  51:                  ulong q = N, p = 1;

  52:                  while ((q /= prime) > 0)

  53:                  {

  54:                      if ((q & 1) == 1) { p *= prime; }

  55:                  }

  56:   

  57:                  if (p > 1) { factors[fi++] = p; }

  58:                  prime = primes[pi++];

  59:              }

  60:   

  61:              max = N / 3;

  62:              while (prime <= max)

  63:              {

  64:                  if (((N / prime) & 1) == 1)

  65:                  {

  66:                      factors[fi++] = prime;

  67:                  }

  68:                  prime = primes[pi++];

  69:              }

  70:   

  71:              pi = lim[i] - lim[i-1];

  72:              memcpy(factors + fi, primes + lim[i-1], pi * sizeof(ulong));

  73:   

  74:              Xmath::Product(swing, factors, 0, fi + pi);

  75:              pow.join();

  76:          }

  77:   

  78:          lmp::Mul(fact, fact, swing);

  79:      }

  80:   

  81:      lmp::Mul2Exp(fact, fact, n - Xmath::BitCount(n));

  82:   

  83:      lmp::Clear(swing);

  84:      lmp::Clear(primorial);

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

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

  87:      lmp::FreeUi(iter, iterLen);

  88:      lmp::FreeUi(lim, iterLen);

  89:  }

  90: