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.IPrimeSieve;
   7:  import de.luschny.math.primes.PrimeSieve;
   8:   
   9:  public class FactorialPrimeSchoenhage implements IFactorialFunction
  10:  {
  11:   
  12:      private int[] primeList;
  13:      private int[] multiList;
  14:   
  15:      public FactorialPrimeSchoenhage()
  16:      {
  17:      }
  18:   
  19:      public String getName()
  20:      {
  21:          return "PrimeSchoenhage   ";
  22:      }
  23:   
  24:      public Xint factorial(int n)
  25:      {
  26:          // For very small n the 'NaiveFactorial' is ok.
  27:          if (n < 20)    { return Xmath.Factorial(n); }
  28:   
  29:          int log2n = 31 - Integer.numberOfLeadingZeros(n);
  30:          int piN = 2 + (15 * n) / (8 * (log2n - 1));
  31:   
  32:          primeList = new int[piN];
  33:          multiList = new int[piN];
  34:   
  35:          int len = primeFactors(n);
  36:          return nestedSquare(len).shiftLeft(n - Integer.bitCount(n));
  37:      }
  38:   
  39:      private Xint nestedSquare(int len)
  40:      {
  41:          if (len == 0)
  42:          {
  43:              return Xint.ONE;
  44:          }
  45:   
  46:          int i = 0, mult = multiList[0];
  47:   
  48:          while (mult > 1)
  49:          {
  50:              if ((mult & 1) == 1) // is mult odd ?
  51:              {
  52:                  primeList[len++] = primeList[i];
  53:              }
  54:   
  55:              multiList[i++] = mult / 2;
  56:              mult = multiList[i];
  57:          }
  58:          if (len <= i)
  59:          {
  60:              return nestedSquare(i).square();
  61:          }
  62:   
  63:          return Xint.product(primeList, i, len - i).multiply(nestedSquare(i).square());
  64:      }
  65:   
  66:      private int primeFactors(int n)
  67:      {
  68:          IPrimeSieve sieve = new PrimeSieve(n);
  69:          IPrimeIteration pIter = sieve.getIteration(3, n);
  70:   
  71:          int maxBound = n / 2, count = 0;
  72:   
  73:          for (int prime : pIter)
  74:          {
  75:              int m = prime > maxBound ? 1 : 0;
  76:   
  77:              if (prime <= maxBound)
  78:              {
  79:                  int q = n;
  80:                  while (q >= prime)
  81:                  {
  82:                      m += q /= prime;
  83:                  }
  84:              }
  85:   
  86:              primeList[count] = prime;
  87:              multiList[count++] = m;
  88:          }
  89:          return count;
  90:      }
  91:  } // endOfFactorialPrimeSchoenhage