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: // Same algorithm as PrimeSwing
9: // but computes swing(n) asynchronous.
10:
11: namespace Sharith.Math.Factorial
12: {
13: using System;
14: using System.Threading.Tasks;
15: using Sharith.Math.Primes;
16: using XMath = Sharith.Math.MathUtils.XMath;
17: using XInt = Sharith.Arithmetic.XInt;
18:
19: public class ParallelPrimeSwing : IFactorialFunction
20: {
21: public ParallelPrimeSwing() { }
22:
23: public string Name
24: {
25: get { return "ParallelPrimeSwing "; }
26: }
27:
28: private const int SMALLSWING = 65;
29: private IAsyncResult[] results;
30: private delegate XInt SwingDelegate(PrimeSieve sieve, int n);
31: private SwingDelegate swingDelegate;
32: private int taskCounter;
33:
34: public XInt Factorial(int n)
35: {
36: if (n < 20) { return XMath.Factorial(n); }
37:
38: var sieve = new PrimeSieve(n);
39: results = new IAsyncResult[XMath.FloorLog2(n)];
40: swingDelegate = Swing; taskCounter = 0;
41: int N = n;
42:
43: // -- It is more efficient to add the big swings
44: // -- first and the small ones later!
45: while (N >= SMALLSWING)
46: {
47: results[taskCounter++] = swingDelegate.BeginInvoke(sieve, N, null, null);
48: N >>= 1;
49: }
50:
51: return RecFactorial(n) << (n - XMath.BitCount(n));
52: }
53:
54: private XInt RecFactorial(int n)
55: {
56: if (n < 2) return XInt.One;
57:
58: XInt recFact = RecFactorial(n / 2);
59: XInt sqrFact = XInt.Pow(recFact, 2);
60: XInt swing;
61:
62: if (n < SMALLSWING)
63: {
64: swing = smallOddSwing[n];
65: }
66: else
67: {
68: swing = swingDelegate.EndInvoke(results[--taskCounter]);
69: }
70: return sqrFact * swing;
71: }
72:
73: private static XInt Swing(PrimeSieve sieve, int n)
74: {
75: var primorial = Task.Factory.StartNew<XInt>(() =>
76: {
77: return sieve.GetPrimorial(n / 2 + 1, n);
78: });
79:
80: int count = 0, rootN = XMath.FloorSqrt(n);
81:
82: var aPrimes = sieve.GetPrimeCollection(3, rootN);
83: var bPrimes = sieve.GetPrimeCollection(rootN + 1, n / 3);
84:
85: int[] primeList = new int[aPrimes.NumberOfPrimes + bPrimes.NumberOfPrimes];
86:
87: foreach (int prime in aPrimes)
88: {
89: int q = n, p = 1;
90:
91: while ((q /= prime) > 0)
92: {
93: if ((q & 1) == 1)
94: {
95: p *= prime;
96: }
97: }
98:
99: if (p > 1)
100: {
101: primeList[count++] = p;
102: }
103: }
104:
105: foreach (int prime in bPrimes)
106: {
107: if (((n / prime) & 1) == 1)
108: {
109: primeList[count++] = prime;
110: }
111: }
112:
113: XInt primeProduct = XMath.Product(primeList, 0, count);
114: return primeProduct * primorial.Result;
115: }
116:
117: private static XInt[] smallOddSwing = {
118: 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003, 429, 6435, 6435,
119: 109395, 12155, 230945, 46189, 969969, 88179, 2028117, 676039, 16900975,
120: 1300075, 35102025, 5014575, 145422675, 9694845, 300540195, 300540195,
121: 9917826435, 583401555, 20419054425, 2268783825, 83945001525, 4418157975,
122: 172308161025, 34461632205, 1412926920405, 67282234305, 2893136075115,
123: 263012370465, 11835556670925, 514589420475, 24185702762325, 8061900920775,
124: 395033145117975, 15801325804719, 805867616040669, 61989816618513,
125: 3285460280781189, 121683714103007, 6692604275665385, 956086325095055,
126: 54496920530418135, 1879204156221315, 110873045217057585, 7391536347803839,
127: 450883717216034179, 14544636039226909, 916312070471295267, 916312070471295267 };
128: }
129: }