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