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