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.IPrimeSieve;
7: import de.luschny.math.primes.PrimeSieve;
8:
9: public class FactorialPrimeSchoenhage implements IFactorialFunction
10: {
11:
12: private int[] primeList;
13: private int[] multiList;
14:
15: public FactorialPrimeSchoenhage()
16: {
17: }
18:
19: public String getName()
20: {
21: return "PrimeSchoenhage ";
22: }
23:
24: public Xint factorial(int n)
25: {
26: // For very small n the 'NaiveFactorial' is ok.
27: if (n < 20) { return Xmath.Factorial(n); }
28:
29: int log2n = 31 - Integer.numberOfLeadingZeros(n);
30: int piN = 2 + (15 * n) / (8 * (log2n - 1));
31:
32: primeList = new int[piN];
33: multiList = new int[piN];
34:
35: int len = primeFactors(n);
36: return nestedSquare(len).shiftLeft(n - Integer.bitCount(n));
37: }
38:
39: private Xint nestedSquare(int len)
40: {
41: if (len == 0)
42: {
43: return Xint.ONE;
44: }
45:
46: int i = 0, mult = multiList[0];
47:
48: while (mult > 1)
49: {
50: if ((mult & 1) == 1) // is mult odd ?
51: {
52: primeList[len++] = primeList[i];
53: }
54:
55: multiList[i++] = mult / 2;
56: mult = multiList[i];
57: }
58: if (len <= i)
59: {
60: return nestedSquare(i).square();
61: }
62:
63: return Xint.product(primeList, i, len - i).multiply(nestedSquare(i).square());
64: }
65:
66: private int primeFactors(int n)
67: {
68: IPrimeSieve sieve = new PrimeSieve(n);
69: IPrimeIteration pIter = sieve.getIteration(3, n);
70:
71: int maxBound = n / 2, count = 0;
72:
73: for (int prime : pIter)
74: {
75: int m = prime > maxBound ? 1 : 0;
76:
77: if (prime <= maxBound)
78: {
79: int q = n;
80: while (q >= prime)
81: {
82: m += q /= prime;
83: }
84: }
85:
86: primeList[count] = prime;
87: multiList[count++] = m;
88: }
89: return count;
90: }
91: } // endOfFactorialPrimeSchoenhage