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 XInt = Sharith.Arithmetic.XInt;
  11:   
  12:      public class SwingDouble : IFactorialFunction  
  13:      {
  14:          public SwingDouble() { }
  15:   
  16:          public string Name
  17:          {
  18:              get { return "SwingDouble         "; }
  19:          }
  20:   
  21:          private XInt f;
  22:          private long gN;
  23:   
  24:          public XInt Factorial(int n)
  25:          {
  26:              if (n < 0)
  27:              {
  28:                  throw new System.ArgumentOutOfRangeException("n",
  29:                  Name + ": n >= 0 required, but was " + n);
  30:              }
  31:   
  32:              gN = 1;
  33:              f = XInt.One;
  34:              return RecFactorial(n);
  35:          }
  36:   
  37:          private XInt RecFactorial(int n)
  38:          {
  39:              if (n < 2) return XInt.One;
  40:   
  41:              return XInt.Pow(RecFactorial(n / 2),2) * Swing(n);
  42:          }
  43:   
  44:          private XInt Swing(long n)
  45:          {
  46:              long s = gN - 1 + ((n - gN + 1) % 4);
  47:              bool oddN = (gN & 1) != 1;
  48:   
  49:              for (; gN <= s; gN++)
  50:              {
  51:                  if (oddN = !oddN) f *= gN;
  52:                  else f = (f * 4) / gN;
  53:              }
  54:   
  55:              if (oddN) for (; gN <= n; gN += 4)
  56:              {
  57:                      long m = ((gN + 1) * (gN + 3)) << 1;
  58:                      long d = (gN * (gN + 2)) >> 3;
  59:   
  60:                      f = (f * m) / d;
  61:              }
  62:              else for (; gN <= n; gN += 4)
  63:              {
  64:                      long m = (gN * (gN + 2)) << 1;
  65:                      long d = ((gN + 1) * (gN + 3)) >> 3;
  66:   
  67:                      f = (f * m) / d;
  68:               }
  69:   
  70:              return f;
  71:          }
  72:      }
  73:  } // endOfSwingDouble