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