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