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 XMath = Sharith.Math.MathUtils.XMath;
12: using XInt = Sharith.Arithmetic.XInt;
13:
14: public class PrimeLeenstra : IFactorialFunction
15: {
16: public PrimeLeenstra() { }
17:
18: public string Name
19: {
20: get { return "PrimeLeenstra "; }
21: }
22:
23: public XInt Factorial(int n)
24: {
25: if (n < 20) { return XMath.Factorial(n); }
26:
27: int rootN = XMath.FloorSqrt(n);
28: int log2N = XMath.FloorLog2(n);
29: var section = new XInt[log2N + 1];
30:
31: for (int i = 0; i < section.Length; i++)
32: {
33: section[i] = XInt.One;
34: }
35:
36: var sieve = new PrimeSieve(n);
37: var primes = sieve.GetPrimeCollection(3, rootN);
38:
39: foreach (int prime in primes)
40: {
41: int k = 0, m = 0, q = n;
42:
43: do
44: {
45: m += q /= prime;
46:
47: } while (q >= 1);
48:
49: while (m > 0)
50: {
51: if ((m & 1) == 1)
52: {
53: section[k] *= prime;
54: }
55: m = m / 2;
56: k++;
57: }
58: }
59:
60: int j = 2, low = n, high;
61:
62: while (low != rootN)
63: {
64: high = low;
65: low = n / j++;
66:
67: if (low < rootN) { low = rootN; }
68:
69: XInt primorial = sieve.GetPrimorial(low + 1, high);
70:
71: if (primorial != XInt.One)
72: {
73: int k = 0, m = j - 2;
74:
75: while (m > 0)
76: {
77: if ((m & 1) == 1)
78: {
79: section[k] *= primorial;
80: }
81: m = m / 2;
82: k++;
83: }
84: }
85: }
86:
87: XInt factorial = section[log2N];
88: for (int i = log2N - 1; i >= 0; --i)
89: {
90: factorial = XInt.Pow(factorial,2) * section[i];
91: }
92:
93: int exp2N = n - XMath.BitCount(n);
94: return factorial << exp2N;
95: }
96: }
97: } // endOfFactorialPrimeLeenstra