1: // Author & Copyright : Peter Luschny
2: // Created: 2010-01-15
3: // License: LGPL version 2.1 or (at your option)
4: // Creative Commons Attribution-ShareAlike 3.0
5:
6: #include "primeswing.h"
7: #include "xmath.h"
8:
9: void PrimeSwing::ParallelFactorial(Xint fact, ulong n)
10: {
11: if (n < THRESHOLD) { Xmath::NaiveFactorial(fact, n); return; }
12: lmp::SetUi(fact, 1);
13:
14: ulong* primes;
15: ulong piN = Xmath::PrimeSieve(&primes, n);
16: ulong* factors = lmp::MallocUi(piN);
17: ulong iterLen = 0; ulong i = n;
18: slong m = n; while (m > 0) { m >>= 1; iterLen++; }
19: ulong* iter = lmp::MallocUi(iterLen);
20: ulong* lim = lmp::MallocUi(iterLen);
21: m = n; i = iterLen;
22: while (m > 0) { iter[--i] = m; m >>= 1; }
23:
24: Xint primorial; lmp::Init(primorial);
25: Xint swing; lmp::Init(swing);
26:
27: lim[0] = 0;
28: for (i = 1; i < iterLen; i++)
29: lim[i] = iter[i] < SOSLEN / 2 ? 0 :
30: GetIndexOf(primes, iter[i], lim[i-1], piN);
31:
32: for (i = 1; i < iterLen; i++)
33: {
34: ulong N = iter[i];
35:
36: if (N < SOSLEN)
37: {
38: lmp::Pow2(fact, fact);
39: lmp::SetUi(swing, smallOddSwing[N]);
40: }
41: else
42: {
43: boost::thread pow(lmp::Pow2, fact, fact);
44:
45: ulong prime = 3;
46: slong pi = 2, fi = 0;
47: ulong max = Xmath::Sqrt(N);
48:
49: while (prime <= max)
50: {
51: ulong q = N, p = 1;
52: while ((q /= prime) > 0)
53: {
54: if ((q & 1) == 1) { p *= prime; }
55: }
56:
57: if (p > 1) { factors[fi++] = p; }
58: prime = primes[pi++];
59: }
60:
61: max = N / 3;
62: while (prime <= max)
63: {
64: if (((N / prime) & 1) == 1)
65: {
66: factors[fi++] = prime;
67: }
68: prime = primes[pi++];
69: }
70:
71: pi = lim[i] - lim[i-1];
72: memcpy(factors + fi, primes + lim[i-1], pi * sizeof(ulong));
73:
74: Xmath::Product(swing, factors, 0, fi + pi);
75: pow.join();
76: }
77:
78: lmp::Mul(fact, fact, swing);
79: }
80:
81: lmp::Mul2Exp(fact, fact, n - Xmath::BitCount(n));
82:
83: lmp::Clear(swing);
84: lmp::Clear(primorial);
85: lmp::FreeUi(primes, piN);
86: lmp::FreeUi(factors, piN);
87: lmp::FreeUi(iter, iterLen);
88: lmp::FreeUi(lim, iterLen);
89: }
90: