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