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