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