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 FactorialPrimeLeenstra implements IFactorialFunction
9: {
10:
11: public FactorialPrimeLeenstra()
12: {
13: }
14:
15: public String getName()
16: {
17: return "PrimeLeenstra ";
18: }
19:
20: public Xint factorial(int n)
21: {
22: // For very small n the 'NaiveFactorial' is OK.
23: if (n < 20) { return Xmath.Factorial(n); }
24:
25: int rootN = (int) Math.floor(Math.sqrt(n));
26: int log2N = 31 - Integer.numberOfLeadingZeros(n);
27: Xint[] expBit = new Xint[log2N + 1];
28:
29: for (int j = 0; j < expBit.length; j++)
30: {
31: expBit[j] = Xint.ONE;
32: }
33:
34: PrimeSieve sieve = new PrimeSieve(n);
35: IPrimeIteration pIter = sieve.getIteration(3, rootN);
36:
37: for (int prime : pIter)
38: {
39: int k = 0, m = 0, q = n;
40:
41: do
42: {
43: m += q /= prime;
44:
45: }
46: while (q >= 1);
47:
48: while (m > 0)
49: {
50: if ((m & 1) == 1)
51: {
52: expBit[k] = expBit[k].multiply(prime);
53: }
54: m = m / 2;
55: k++;
56: }
57: }
58:
59: int j = 2, low = n, high;
60:
61: while (low != rootN)
62: {
63: high = low;
64: low = n / j++;
65:
66: if (low < rootN)
67: {
68: low = rootN;
69: }
70:
71: Xint primorial = sieve.getPrimorial(low + 1, high);
72:
73: if (!primorial.isONE())
74: {
75: int k = 0, m = j - 2;
76:
77: while (m > 0)
78: {
79: if ((m & 1) == 1)
80: {
81: expBit[k] = expBit[k].multiply(primorial);
82: }
83: m = m / 2;
84: k++;
85: }
86: }
87: }
88:
89: Xint fact = expBit[log2N];
90: for (int i = log2N - 1; i >= 0; --i)
91: {
92: fact = fact.square().multiply(expBit[i]);
93: }
94:
95: return fact.shiftLeft(n - Integer.bitCount(n));
96: }
97:
98: } // endOfFactorialPrimeLeenstra