1:  /// -------- ToujoursEnBeta
   2:  /// Author & Copyright : Peter Luschny
   3:  /// License: LGPL version 3.0 or (at your option)
   4:  /// Creative Commons Attribution-ShareAlike 3.0
   5:  /// Comments mail to: peter(at)luschny.de
   6:  /// Created: 2010-03-01
   7:   
   8:  namespace Sharith.Math.Factorial 
   9:  {    
  10:      using Sharith.Math.Primes;
  11:      using XMath = Sharith.Math.MathUtils.XMath;
  12:      using XInt = Sharith.Arithmetic.XInt;
  13:   
  14:      public class PrimeSwing : IFactorialFunction 
  15:      {
  16:          public PrimeSwing() { }
  17:   
  18:          private PrimeSieve sieve;
  19:          private int[] primeList;
  20:   
  21:          public string Name
  22:          {
  23:              get { return "PrimeSwing          "; }
  24:          }   
  25:   
  26:          public XInt Factorial(int n)
  27:          {
  28:              if (n < 20) { return XMath.Factorial(n); }
  29:   
  30:              sieve = new PrimeSieve(n);
  31:              int pLen = (int)(2.0 * (XMath.FloorSqrt(n)
  32:                       + (double)n / (XMath.Log2(n) - 1)));
  33:              primeList = new int[pLen];
  34:              
  35:              int exp2 = n - XMath.BitCount(n);
  36:              return RecFactorial(n) << exp2;
  37:          }
  38:   
  39:          private XInt RecFactorial(int n)
  40:          {
  41:              if (n < 2) return XInt.One;
  42:   
  43:              return XInt.Pow(RecFactorial(n / 2), 2) * Swing(n);
  44:          }
  45:   
  46:          private XInt Swing(int n)
  47:          {
  48:              if (n < 33) return smallOddSwing[n];
  49:   
  50:              int count = 0, rootN = XMath.FloorSqrt(n);
  51:   
  52:              var aPrimes = sieve.GetPrimeCollection(3, rootN);
  53:              var bPrimes = sieve.GetPrimeCollection(rootN + 1, n / 3);
  54:   
  55:              foreach (int prime in aPrimes)
  56:              {
  57:                  int q = n, p = 1;
  58:   
  59:                  while ((q /= prime) > 0)
  60:                  {
  61:                      if ((q & 1) == 1)
  62:                      {
  63:                          p *= prime;
  64:                      }
  65:                  }
  66:   
  67:                  if (p > 1)
  68:                  {
  69:                      primeList[count++] = p;
  70:                  }
  71:              }
  72:   
  73:              foreach (int prime in bPrimes)
  74:              {
  75:                  if (((n / prime) & 1) == 1)
  76:                  {
  77:                      primeList[count++] = prime;
  78:                  }
  79:              }
  80:   
  81:              XInt primorial = sieve.GetPrimorial(n / 2 + 1, n);
  82:              return primorial * XMath.Product(primeList, 0, count);
  83:          }
  84:   
  85:          private static XInt[] smallOddSwing = {
  86:              1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,109395,
  87:              12155,230945,46189,969969,88179,2028117,676039,16900975,1300075,
  88:              35102025,5014575,145422675,9694845,300540195,300540195 };
  89:      }
  90:  } // endOfFactorialPrimeSwingLuschny