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 Sharith.Math.Primes;
11: using XInt = Sharith.Arithmetic.XInt;
12: using XMath = Sharith.Math.MathUtils.XMath;
13:
14: public class PrimeSwingList : IFactorialFunction
15: {
16: public PrimeSwingList() { }
17:
18: private int[][] primeList;
19: private int[] listLength;
20: private int[] tower;
21: private int[] bound;
22:
23: public string Name
24: {
25: get { return "PrimeSwingList "; }
26: }
27:
28: public XInt Factorial(int n)
29: {
30: if (n < 20) { return XMath.Factorial(n); }
31:
32: int log2n = XMath.FloorLog2(n);
33: int j = log2n;
34: int hN = n;
35:
36: primeList = new int[log2n][];
37: listLength = new int[log2n];
38: bound = new int[log2n];
39: tower = new int[log2n + 1];
40:
41: while (true)
42: {
43: tower[j] = hN;
44: if (hN == 1) break;
45: bound[--j] = hN / 3;
46: int pLen = hN < 4 ? 6 : (int)(2.0 * (XMath.FloorSqrt(hN)
47: + (double) hN / (XMath.Log2(hN) - 1)));
48: primeList[j] = new int[pLen];
49: hN >>= 1;
50: }
51:
52: tower[0] = 2;
53: PrimeFactors(n);
54:
55: int init = listLength[0] == 0 ? 1 : 3;
56: var oddFactorial = new XInt(init);
57:
58: for (int i = 1; i < log2n; i++)
59: {
60: oddFactorial = XInt.Pow(oddFactorial, 2)
61: * XMath.Product(primeList[i], 0, listLength[i]);
62: }
63: return oddFactorial << (n - XMath.BitCount(n));
64: }
65:
66: private void PrimeFactors(int n)
67: {
68: int maxBound = n / 3;
69: int lastList = primeList.Length - 1;
70: int start = tower[1] == 2 ? 1 : 0;
71: var sieve = new PrimeSieve(n);
72:
73: for (int section = start; section < primeList.Length; section++)
74: {
75: var primes = sieve.GetPrimeCollection(tower[section] + 1, tower[section + 1]);
76:
77: foreach (int prime in primes)
78: {
79: primeList[section][listLength[section]++] = prime;
80: if (prime > maxBound) continue;
81:
82: int np = n;
83: do
84: {
85: int k = lastList;
86: int q = np /= prime;
87:
88: do if ((q & 1) == 1) //if q is odd
89: {
90: primeList[k][listLength[k]++] = prime;
91: }
92: while (((q >>= 1) > 0) && (prime <= bound[--k]));
93:
94: } while (prime <= np);
95: }
96: }
97: }
98: }
99: } // endOfFactorialPrimeSwingList