1: package de.luschny.math.factorial;
2:
3: import de.luschny.math.Xmath;
4: import de.luschny.math.arithmetic.Xint;
5: import de.luschny.math.primes.IPrimeIteration;
6: import de.luschny.math.primes.PrimeSieve;
7:
8: public class FactorialPrimeVardi implements IFactorialFunction
9: {
10:
11: private PrimeSieve sieve;
12:
13: public FactorialPrimeVardi()
14: {
15: }
16:
17: public final String getName()
18: {
19: return "PrimeVardi ";
20: }
21:
22: public Xint factorial(int n)
23: {
24: // For very small n the 'NaiveFactorial' is ok.
25: if (n < 20) { return Xmath.Factorial(n); }
26:
27: sieve = new PrimeSieve(n);
28:
29: return recFactorial(n);
30: }
31:
32: private Xint recFactorial(int n)
33: {
34: if (n < 2)
35: {
36: return Xint.ONE;
37: }
38:
39: if ((n & 1) == 1)
40: {
41: return recFactorial(n - 1).multiply(n);
42: }
43:
44: return middleBinomial(n).multiply(recFactorial(n / 2).square());
45: }
46:
47: private Xint middleBinomial(int n) // assuming n = 2k
48: {
49: if (n < 50)
50: {
51: return Xint.valueOf(binom[n / 2]);
52: }
53:
54: int k = n / 2, pc = 0, pp = 0, e;
55: int rootN = (int) Math.floor(Math.sqrt(n));
56:
57: Xint bigPrimes = sieve.getPrimorial(k + 1, n);
58: Xint smallPrimes = sieve.getPrimorial(k / 2 + 1, n / 3);
59:
60: IPrimeIteration pIter = sieve.getIteration(rootN + 1, n / 5);
61: int[] primeList = new int[pIter.getNumberOfPrimes()];
62:
63: for (int prime : pIter)
64: {
65: if ((n / prime & 1) == 1) // if n/prime is odd...
66: {
67: primeList[pc++] = prime;
68: }
69: }
70: Xint prodPrimes = Xint.product(primeList, 0, pc);
71:
72: pIter = sieve.getIteration(1, rootN);
73: Xint[] primePowerList = new Xint[pIter.getNumberOfPrimes()];
74:
75: for (int prime : pIter)
76: {
77: if ((e = expSum(prime, n)) > 0)
78: {
79: primePowerList[pp++] = Xint.valueOf(prime).toPowerOf(e);
80: }
81: }
82: Xint powerPrimes = Xint.product(primePowerList, 0, pp);
83:
84: return bigPrimes.multiply(smallPrimes).multiply(prodPrimes).multiply(powerPrimes);
85: }
86:
87: private int expSum(int p, int n)
88: {
89: int exp = 0, q = n / p;
90:
91: while (0 < q)
92: {
93: exp += q & 1;
94: q /= p;
95: }
96:
97: return exp;
98: }
99:
100: private static long[] binom =
101: { 1, 2, 6, 20, 70, 252, 924, 3432, 12870, 48620, 184756, 705432, 2704156,
102: 10400600, 40116600, 155117520, 601080390, 2333606220L, 9075135300L,
103: 35345263800L, 137846528820L, 538257874440L, 2104098963720L, 8233430727600L,
104: 32247603683100L };
105:
106: } // endOfFactorialPrimeVardi