1: package de.luschny.math.factorial;
2:
3: import de.luschny.math.arithmetic.Xint;
4:
5: public class FactorialSwingRationalDouble implements IFactorialFunction
6: {
7:
8: public FactorialSwingRationalDouble()
9: {
10: }
11:
12: private long D, N, g, h;
13: private int i;
14:
15: public String getName()
16: {
17: return "SwingRationalDbl ";
18: }
19:
20: public Xint factorial(int n)
21: {
22: if (n < 0)
23: {
24: throw new ArithmeticException("Factorial: n has to be >= 0, but was " + n);
25: }
26:
27: return recFactorial(n);
28: }
29:
30: private Xint recFactorial(int n)
31: {
32: if (n < 2)
33: {
34: return Xint.ONE;
35: }
36:
37: return recFactorial(n / 2).square().multiply(swing(n));
38: }
39:
40: private Xint swing(int n)
41: {
42: boolean oddN = (n & 1) == 1;
43: boolean div = false;
44: h = n / 2;
45:
46: switch ((n / 2) % 4)
47: {
48: case 0:
49: h = oddN ? h + 1 : 1;
50: break;
51: case 1:
52: h = oddN ? 2 * (h + 2) : 2;
53: break;
54: case 2:
55: h = oddN ? 2 * (h + 1) * (h + 3) : 2 * (h + 1);
56: div = n > 7;
57: break;
58: case 3:
59: h = oddN ? 4 * (h + 2) * (h + 4) : 4 * (h + 2);
60: div = n > 7;
61: break;
62: }
63:
64: g = div ? n / 4 : 1;
65: N = 2 * (n + 3 + (n & 1));
66: D = -1;
67: i = n / 8;
68:
69: return product(i + 1).getNumerator();
70: }
71:
72: private Rational product(int l)
73: {
74: if (l > 1)
75: {
76: int m = l / 2;
77: return product(m).multiply(product(l - m));
78: }
79:
80: if (i-- > 0)
81: {
82: N -= 8;
83: D += 2;
84: return new Rational(N * (N - 4), D * (D + 1));
85: }
86:
87: return new Rational(h, g);
88: }
89:
90: // --------------------------------------------------------
91: // A minimalistic rational arithmetic *only* for the use of
92: // FactorialSwingRational. The sole purpose for inclusion
93:
94: // here is to make the description of the algorithm more
95: // independent and more easy to port.
96: // ---------------------------------------------------------
97: private class Rational
98: {
99:
100: private Xint num; // Numerator
101: private Xint den; // Denominator
102:
103: public Rational(long _num, long _den)
104: {
105: long c = gcd(_num, _den);
106: num = Xint.valueOf(_num / c);
107: den = Xint.valueOf(_den / c);
108: }
109:
110: public Rational(Xint _num, Xint _den)
111: {
112: // If the arithmetic supports a *real* fast gcd
113: // this would lead to a speed up:
114: // Xint cd = _num.gcd(_den);
115: // num = _num.divide(cd);
116: // den = _den.divide(cd);
117: num = _num;
118: den = _den;
119: }
120:
121: public Xint getNumerator()
122: {
123: Xint cd = num.gcd(den);
124: return num.divide(cd);
125: }
126:
127: public Rational multiply(Rational r)
128: {
129: return new Rational(num.multiply(r.num), den.multiply(r.den));
130: }
131:
132: private long gcd(long a, long b)
133: {
134: long x, y;
135:
136: if (a >= b)
137: {
138: x = a;
139: y = b;
140: }
141: else
142: {
143: x = b;
144: y = a;
145: }
146:
147: while (y != 0)
148: {
149: long t = x % y;
150: x = y;
151: y = t;
152: }
153: return x;
154: }
155: } // endOfRational
156: } // endOfSwingRational