1: using System;
2: using Xint = System.Numerics.BigInteger;
3: using Xmath = Luschny.Math.MathUtils.Xmath;
4:
5: namespace Luschny.Math.Factorial
6: {
7: public class FactorialSwing : IFactorialFunction
8: {
9: public FactorialSwing() { }
10:
11: public string Name
12: {
13: get { return "Swing "; }
14: }
15:
16: private Xint oddFactNdiv4, oddFactNdiv2;
17: private const int SMALLSWING = 33;
18: private const int SMALLFACT = 17;
19:
20: public Xint Factorial(int n)
21: {
22: if (n < 0)
23: {
24: throw new ArithmeticException(
25: "Factorial: n has to be >= 0, but was " + n);
26: }
27:
28: oddFactNdiv4 = oddFactNdiv2 = Xint.One;
29:
30: return OddFactorial(n) << (n - Xmath.BitCount(n));
31: }
32:
33: private Xint OddFactorial(int n)
34: {
35: Xint oddFact;
36:
37: if (n < SMALLFACT)
38: {
39: oddFact = smallOddFactorial[n];
40: }
41: else
42: {
43: Xint sqrOddFact = OddFactorial(n / 2);
44:
45: int ndiv4 = n / 4;
46: Xint oddFactNd4 = ndiv4 < SMALLFACT ? smallOddFactorial[ndiv4] : oddFactNdiv4;
47:
48: oddFact = Xint.Pow(sqrOddFact, 2) * OddSwing(n, oddFactNd4);
49: }
50:
51: oddFactNdiv4 = oddFactNdiv2;
52: oddFactNdiv2 = oddFact;
53: return oddFact;
54: }
55:
56: private Xint OddSwing(int n, Xint oddFactNdiv4)
57: {
58: if (n < SMALLSWING) return smallOddSwing[n];
59:
60: int len = (n - 1) / 4;
61: if ((n % 4) != 2) len++;
62: int high = n - ((n + 1) & 1);
63:
64: return Product(high, len) / oddFactNdiv4;
65: }
66:
67: private static Xint Product(int m, int len)
68: {
69: if (len == 1) return new Xint(m);
70: if (len == 2) return new Xint((long) m * (m - 2));
71:
72: int hlen = len >> 1;
73: return Product(m - hlen * 2, len - hlen) * Product(m, hlen);
74: }
75:
76: private static Xint[] smallOddSwing = {
77: 1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,109395,
78: 12155,230945,46189,969969,88179,2028117,676039,16900975,1300075,
79: 35102025,5014575,145422675,9694845,300540195,300540195 };
80:
81: private static Xint[] smallOddFactorial = {
82: 1,1,1,3,3,15,45,315,315,2835,14175,155925,467775,6081075,
83: 42567525,638512875,638512875 };
84:
85: } // endOfFactorialSwingRecursive
86: }