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