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