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