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 System;
  11:      using System.Threading.Tasks;
  12:      using XInt = Sharith.Arithmetic.XInt;
  13:      using XMath = Sharith.Math.MathUtils.XMath;
  14:   
  15:      public class ParallelSwing : IFactorialFunction 
  16:      {
  17:          public ParallelSwing() { }
  18:   
  19:          public string Name
  20:          {
  21:              get { return "ParallelSwing         "; }
  22:          }
  23:   
  24:          private XInt oddFactNdiv4, oddFactNdiv2;
  25:          private const int SMALLSWING = 33;
  26:          private const int SMALLFACT = 17;
  27:          private Task<XInt> oddSwingTask;
  28:   
  29:          public XInt Factorial(int n)
  30:          {
  31:              if (n < 0)
  32:              {
  33:                  throw new ArithmeticException(
  34:                  "Factorial: n has to be >= 0, but was " + n);
  35:              }
  36:   
  37:              oddFactNdiv4 = oddFactNdiv2 = XInt.One;
  38:              return OddFactorial(n) << (n - XMath.BitCount(n));
  39:          }
  40:   
  41:          private XInt OddFactorial(int n)
  42:          {
  43:              XInt oddFact;
  44:   
  45:              if (n < SMALLFACT)
  46:              {
  47:                  return smallOddFactorial[n];
  48:              }
  49:              else
  50:              {
  51:                  XInt sqrOddFact = OddFactorial(n / 2);
  52:   
  53:                  if (n < SMALLSWING)
  54:                  {
  55:                      oddFact = XInt.Pow(sqrOddFact, 2) * smallOddSwing[n];
  56:                  }
  57:                  else
  58:                  {
  59:                      int ndiv4 = n / 4;
  60:                      var oddFactNd4 = ndiv4 < SMALLFACT ? smallOddFactorial[ndiv4] : oddFactNdiv4;
  61:                      oddSwingTask = Task.Factory.StartNew<XInt>( () => 
  62:                      {
  63:                          return OddSwing(n, oddFactNd4); 
  64:                      });
  65:   
  66:                      sqrOddFact = XInt.Pow(sqrOddFact, 2);
  67:                      oddFact = sqrOddFact * oddSwingTask.Result;
  68:                  }
  69:   
  70:                  oddFactNdiv4 = oddFactNdiv2;
  71:                  oddFactNdiv2 = oddFact;
  72:                  return oddFact;
  73:              }
  74:          }
  75:   
  76:          private XInt OddSwing(int n, XInt oddFactNdiv4)
  77:          {
  78:              int len = (n - 1) / 4;
  79:              if ((n % 4) != 2) len++;
  80:   
  81:              //-- if type(n, odd) then high = n else high = n-1.
  82:              int high = n - ((n + 1) & 1);
  83:   
  84:              return Product(high, len) / oddFactNdiv4;
  85:          }
  86:   
  87:          private static XInt Product(int m, int len)
  88:          {
  89:              if (len == 1) return new XInt(m);
  90:              if (len == 2) return new XInt((long)m * (m - 2));
  91:   
  92:              int hlen = len >> 1;
  93:              return Product(m - hlen * 2, len - hlen) * Product(m, hlen);
  94:          }
  95:   
  96:          private static XInt[] smallOddSwing = {
  97:              1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,109395,
  98:              12155,230945,46189,969969,88179,2028117,676039,16900975,1300075,
  99:              35102025,5014575,145422675,9694845,300540195,300540195 };
 100:   
 101:          private static XInt[] smallOddFactorial = {
 102:              1,1,1,3,3,15,45,315,315,2835,14175,155925,467775,6081075,
 103:              42567525,638512875,638512875 };
 104:   
 105:      } // endOfFactorialParallelSwing
 106:  }