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