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