1: using System.Collections.Generic;
2: using Xint = System.Numerics.BigInteger;
3: using Xmath = Luschny.Math.MathUtils.Xmath;
4: using Luschny.Math.Primes;
5:
6: namespace Luschny.Math.Factorial
7: {
8: public class FactorialPrimeSwingCache : IFactorialFunction
9: {
10: private Dictionary<int, CachedPrimorial> cache;
11: private PrimeSieve sieve;
12: private int[] primeList;
13:
14: public FactorialPrimeSwingCache() { }
15:
16: public string Name
17: {
18: get { return "PrimeSwingCache "; }
19: }
20:
21: public Xint Factorial(int n)
22: {
23: if (n < 20)
24: {
25: return Xmath.Factorial(n);
26: }
27:
28: cache = new Dictionary<int, CachedPrimorial>();
29: sieve = new PrimeSieve(n);
30:
31: int pLen = (int)(2.0 * (Xmath.FloorSqrt(n)
32: + (double)n / (Xmath.Log2(n) - 1)));
33: primeList = new int[pLen];
34:
35: int exp2 = n - Xmath.BitCount(n);
36: return RecFactorial(n) << exp2;
37: }
38:
39: private Xint RecFactorial(int n)
40: {
41: if (n < 2) return Xint.One;
42:
43: //-- Not commutative!!
44: return Swing(n) * Xint.Pow(RecFactorial(n / 2),2);
45: }
46:
47: private Xint Swing(int n)
48: {
49: if (n < 33) return smallOddSwing[n];
50:
51: int count = 0, rootN = Xmath.FloorSqrt(n);
52: int j = 1, low, high;
53:
54: Xint prod = Xint.One;
55:
56: while (true)
57: {
58: high = n / j++;
59: low = n / j++;
60:
61: if (low < rootN) { low = rootN; }
62: if (high - low < 32) break;
63:
64: Xint primorial = GetPrimorial(low + 1, high);
65: prod *= primorial;
66: }
67:
68: IPrimeCollection primes = sieve.GetPrimeCollection(3, high);
69:
70: foreach (int prime in primes)
71: {
72: int q = n, p = 1;
73:
74: while ((q /= prime) > 0)
75: {
76: if ((q & 1) == 1)
77: {
78: p *= prime;
79: }
80: }
81:
82: if (p > 1)
83: {
84: primeList[count++] = p;
85: }
86: }
87:
88: return prod * Xmath.Product(primeList, 0, count);
89: }
90:
91: Xint GetPrimorial(int low, int high)
92: {
93: CachedPrimorial lowPrimorial;
94: Xint primorial;
95:
96: if (cache.TryGetValue(low, out lowPrimorial))
97: {
98: //-- This is the catch! The intervals expand.
99: int mid = lowPrimorial.High + 1;
100: Xint highPrimorial = sieve.GetPrimorial(mid, high);
101: primorial = highPrimorial * lowPrimorial.Value;
102: }
103: else
104: {
105: primorial = sieve.GetPrimorial(low, high);
106: }
107:
108: cache[low] = new CachedPrimorial(high, primorial);
109: return primorial;
110: }
111:
112: private static Xint[] smallOddSwing = {1, 1, 1, 3, 3, 15, 5, 35, 35,
113: 315, 63, 693, 231, 3003, 429, 6435, 6435, 109395, 12155, 230945,
114: 46189, 969969, 88179, 2028117, 676039, 16900975, 1300075, 35102025,
115: 5014575, 145422675, 9694845, 300540195, 300540195};
116:
117: private struct CachedPrimorial
118: {
119: public int High; // class { get; set; }
120: public Xint Value; // class { get; set; }
121:
122: public CachedPrimorial(int highBound, Xint val)
123: {
124: High = highBound;
125: Value = val;
126: }
127: }
128: }
129: } // endOfFactorialPrimeSwingCache
130: