`   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`