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 SwingRationalDouble : IFactorialFunction 
  13:      {
  14:          public SwingRationalDouble() { }
  15:   
  16:          public string Name
  17:          {
  18:              get { return "SwingRationalDouble "; }
  19:          }                
  20:   
  21:          private long den, num, g, h;
  22:          private int i;
  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:              return RecFactorial(n);
  33:          }
  34:   
  35:          private XInt RecFactorial(int n)
  36:          {
  37:              if (n < 2) return XInt.One;
  38:   
  39:              return XInt.Pow(RecFactorial(n / 2),2) * Swing(n);
  40:          }
  41:   
  42:          private XInt Swing(int n)
  43:          {
  44:              bool oddN = (n & 1) == 1;
  45:              bool div = false;
  46:              h = n / 2;
  47:   
  48:              switch ((n / 2) % 4)
  49:              {
  50:                  case 0: h = oddN ? h + 1 : 1; break;
  51:                  case 1: h = oddN ? 2 * (h + 2) : 2; break;
  52:                  case 2: h = oddN ? 2 * (h + 1) * (h + 3) : 2 * (h + 1);
  53:                      div = n > 7; break;
  54:                  case 3: h = oddN ? 4 * (h + 2) * (h + 4) : 4 * (h + 2);
  55:                      div = n > 7; break;
  56:              }
  57:   
  58:              g = div ? n / 4 : 1;
  59:              num = 2 * (n + 3 + (n & 1));
  60:              den = -1;
  61:              i = n / 8;
  62:   
  63:              return Product(i + 1).Numerator;
  64:          }
  65:   
  66:          private Rational Product(int l)
  67:          {
  68:              if (l > 1)
  69:              {
  70:                  int m = l / 2;
  71:                  return Product(m) * Product(l - m);
  72:              }
  73:   
  74:              if (i-- > 0)
  75:              {
  76:                  num -= 8;
  77:                  den += 2;
  78:                  return new Rational(num * (num - 4), den * (den + 1));
  79:              }
  80:   
  81:              return new Rational(h, g);
  82:          }
  83:   
  84:          //----------------------------------------------------------
  85:          // A minimalistic rational arithmetic *only* for the use of
  86:          // SwingRational. The sole purpose for inclusion
  87:          // here is to make the description of the algorithm more
  88:          // independent and more easy to port.
  89:          //---------------------------------------------------------
  90:          private class Rational
  91:          {
  92:              private XInt num; // Numerator
  93:              private XInt den; // Denominator
  94:   
  95:              public XInt Numerator
  96:              {
  97:                  get
  98:                  {
  99:                      XInt cd = XInt.GreatestCommonDivisor(num, den);
 100:                      return num / cd;
 101:                  }
 102:              }
 103:   
 104:              public Rational(long _num, long _den)
 105:              {
 106:                  long g = Gcd(_num, _den);
 107:                  num = new XInt(_num / g);
 108:                  den = new XInt(_den / g);
 109:              }
 110:   
 111:              public Rational(XInt _num, XInt _den)
 112:              {
 113:                  //  If (and only if) the arithmetic supports a
 114:                  //  *real* fast Gcd this would lead to a speed up:
 115:                  num = _num;
 116:                  den = _den;
 117:              }
 118:   
 119:              public static Rational operator *(Rational a, Rational r)
 120:              {
 121:                  return new Rational(a.num * r.num, a.den * r.den);
 122:              }
 123:   
 124:              private static long Gcd(long a, long b)
 125:              {
 126:                  long x, y;
 127:   
 128:                  if (a >= b) { x = a; y = b; }
 129:                  else { x = b; y = a; }
 130:   
 131:                  while (y != 0)
 132:                  {
 133:                      long t = x % y; x = y; y = t;
 134:                  }
 135:                  return x;
 136:              }
 137:          } // endOfRational
 138:      }
 139:   
 140:  } // endOfSwingRationalDouble