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