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 FactorialPrimeVardi : IFactorialFunction
8: {
9: private PrimeSieve sieve;
10:
11: public FactorialPrimeVardi() { }
12:
13: public string Name
14: {
15: get { return "PrimeVardi "; }
16: }
17:
18: public Xint Factorial(int n)
19: {
20: if (n < 20) { return Xmath.Factorial(n); }
21:
22: sieve = new PrimeSieve(n);
23:
24: return RecFactorial(n);
25: }
26:
27: private Xint RecFactorial(int n)
28: {
29: if (n < 2) return Xint.One;
30:
31: if ((n & 1) == 1)
32: {
33: return RecFactorial(n - 1) * n;
34: }
35:
36: return MiddleBinomial(n) * Xint.Pow(RecFactorial(n / 2), 2);
37: }
38:
39: private Xint MiddleBinomial(int n) // assuming n = 2k
40: {
41: if (n < 50) return (Xint) binomial[n / 2];
42:
43: int k = n / 2, pc = 0, pp = 0, exp;
44: int rootN = Xmath.FloorSqrt(n);
45:
46: Xint bigPrimes = sieve.GetPrimorial(k + 1, n);
47: Xint smallPrimes = sieve.GetPrimorial(k / 2 + 1, n / 3);
48:
49: IPrimeCollection primes = sieve.GetPrimeCollection(rootN + 1, n / 5);
50: int[] primeList = new int[primes.NumberOfPrimes];
51:
52: foreach (int prime in primes)
53: {
54: if ((n / prime & 1) == 1) // if n/prime is odd...
55: {
56: primeList[pc++] = prime;
57: }
58: }
59: Xint prodPrimes = Xmath.Product(primeList, 0, pc);
60:
61: primes = sieve.GetPrimeCollection(1, rootN);
62: Xint[] primePowers = new Xint[primes.NumberOfPrimes];
63:
64: foreach (int prime in primes)
65: {
66: if ((exp = ExpSum(prime, n)) > 0)
67: {
68: primePowers[pp++] = Xint.Pow((Xint)prime, exp);
69: }
70: }
71: Xint powerPrimes = Xmath.Product(primePowers, 0, pp);
72:
73: return bigPrimes * smallPrimes * prodPrimes * powerPrimes;
74: }
75:
76: private static int ExpSum(int p, int n)
77: {
78: int exp = 0, q = n / p;
79:
80: while (0 < q)
81: {
82: exp += q & 1;
83: q /= p;
84: }
85:
86: return exp;
87: }
88:
89: private static long[] binomial = {
90: 1,2,6,20,70,252,924,3432,12870,48620,184756,705432,2704156,
91: 10400600,40116600,155117520,601080390,2333606220L,9075135300L,
92: 35345263800L,137846528820L,538257874440L,2104098963720L,
93: 8233430727600L,32247603683100L };
94: }
95: }