1: using System;
2: using System.Threading.Tasks;
3: using Xint = System.Numerics.BigInteger;
4: using Xmath = Luschny.Math.MathUtils.Xmath;
5: using Luschny.Math.Primes;
6:
7: namespace Luschny.Math.Factorial
8: {
9: public class FactorialParallelPrimeSplit : IFactorialFunction
10: {
11: private PrimeSieve sieve;
12:
13: public FactorialParallelPrimeSplit() { }
14:
15: public string Name
16: {
17: get { return "ParallelPrimeBinSplit "; }
18: }
19:
20: private delegate Xint SwingDelegate(int n);
21:
22: public Xint Factorial(int n)
23: {
24: if (n < 20) { return Xmath.Factorial(n); }
25:
26: sieve = new PrimeSieve(n);
27: int log2n = Xmath.FloorLog2(n);
28:
29: SwingDelegate swingDelegate = Swing;
30: IAsyncResult[] result = new IAsyncResult[log2n];
31:
32: int h = 0, shift = 0, taskCounter = 0;
33:
34: // -- It is more efficient to add the big intervals
35: // -- first and the small ones later!
36: while (h != n)
37: {
38: shift += h;
39: h = n >> log2n--;
40: if (h > 2)
41: {
42: result[taskCounter++] = swingDelegate.BeginInvoke(h, null, null);
43: }
44: }
45:
46: Xint p = Xint.One, r = Xint.One, rl = Xint.One;
47:
48: for (int i = 0; i < taskCounter; i++)
49: {
50: var t = rl * swingDelegate.EndInvoke(result[i]);
51: p = p * t;
52: rl = r;
53: r = r * p;
54: }
55:
56: return r << shift;
57: }
58:
59: private Xint Swing(int n)
60: {
61: if (n < 33) return (Xint)smallOddSwing[n];
62:
63: var primorial = Task.Factory.StartNew<Xint>(() =>
64: { return sieve.GetPrimorial(n / 2 + 1, n); });
65:
66: int count = 0, rootN = Xmath.FloorSqrt(n);
67:
68: IPrimeCollection aPrimes = sieve.GetPrimeCollection(3, rootN);
69: IPrimeCollection bPrimes = sieve.GetPrimeCollection(rootN + 1, n / 3);
70: int piN = aPrimes.NumberOfPrimes + bPrimes.NumberOfPrimes;
71: int[] primeList = new int[piN];
72:
73: foreach (int prime in aPrimes)
74: {
75: int q = n, p = 1;
76:
77: while ((q /= prime) > 0)
78: {
79: if ((q & 1) == 1)
80: {
81: p *= prime;
82: }
83: }
84:
85: if (p > 1)
86: {
87: primeList[count++] = p;
88: }
89: }
90:
91: foreach (int prime in bPrimes)
92: {
93: if (((n / prime) & 1) == 1)
94: {
95: primeList[count++] = prime;
96: }
97: }
98:
99: Xint primeProduct = Xmath.Product(primeList, 0, count);
100:
101: return primeProduct * primorial.Result;
102: }
103:
104: private static int[] smallOddSwing = {
105: 1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,109395,
106: 12155,230945,46189,969969,88179,2028117,676039,16900975,1300075,
107: 35102025,5014575,145422675,9694845,300540195,300540195 };
108: }
109: }