1:  package de.luschny.math.factorial;
   2:   
   3:  import de.luschny.math.Xmath;
   4:  import de.luschny.math.arithmetic.Xint;
   5:  import de.luschny.math.primes.IPrimeIteration;
   6:  import de.luschny.math.primes.PrimeSieve;
   7:   
   8:  public class FactorialPrimeLeenstra implements IFactorialFunction
   9:  {
  10:   
  11:      public FactorialPrimeLeenstra()
  12:      {
  13:      }
  14:   
  15:      public String getName()
  16:      {
  17:          return "PrimeLeenstra     ";
  18:      }
  19:   
  20:      public Xint factorial(int n)
  21:      {
  22:          // For very small n the 'NaiveFactorial' is OK.
  23:          if (n < 20)    { return Xmath.Factorial(n); }
  24:   
  25:          int rootN = (int) Math.floor(Math.sqrt(n));
  26:          int log2N = 31 - Integer.numberOfLeadingZeros(n);
  27:          Xint[] expBit = new Xint[log2N + 1];
  28:   
  29:          for (int j = 0; j < expBit.length; j++)
  30:          {
  31:              expBit[j] = Xint.ONE;
  32:          }
  33:   
  34:          PrimeSieve sieve = new PrimeSieve(n);
  35:          IPrimeIteration pIter = sieve.getIteration(3, rootN);
  36:   
  37:          for (int prime : pIter)
  38:          {
  39:              int k = 0, m = 0, q = n;
  40:   
  41:              do
  42:              {
  43:                  m += q /= prime;
  44:   
  45:              }
  46:              while (q >= 1);
  47:   
  48:              while (m > 0)
  49:              {
  50:                  if ((m & 1) == 1)
  51:                  {
  52:                      expBit[k] = expBit[k].multiply(prime);
  53:                  }
  54:                  m = m / 2;
  55:                  k++;
  56:              }
  57:          }
  58:   
  59:          int j = 2, low = n, high;
  60:   
  61:          while (low != rootN)
  62:          {
  63:              high = low;
  64:              low = n / j++;
  65:   
  66:              if (low < rootN)
  67:              {
  68:                  low = rootN;
  69:              }
  70:   
  71:              Xint primorial = sieve.getPrimorial(low + 1, high);
  72:   
  73:              if (!primorial.isONE())
  74:              {
  75:                  int k = 0, m = j - 2;
  76:   
  77:                  while (m > 0)
  78:                  {
  79:                      if ((m & 1) == 1)
  80:                      {
  81:                          expBit[k] = expBit[k].multiply(primorial);
  82:                      }
  83:                      m = m / 2;
  84:                      k++;
  85:                  }
  86:              }
  87:          }
  88:   
  89:          Xint fact = expBit[log2N];
  90:          for (int i = log2N - 1; i >= 0; --i)
  91:          {
  92:              fact = fact.square().multiply(expBit[i]);
  93:          }
  94:   
  95:          return fact.shiftLeft(n - Integer.bitCount(n));
  96:      }
  97:      
  98:  } // endOfFactorialPrimeLeenstra