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