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 FactorialPrimeBorwein implements IFactorialFunction
10: {
11:
12: private int[] primeList;
13: private int[] multiList;
14:
15: public FactorialPrimeBorwein()
16: {
17: }
18:
19: public String getName()
20: {
21: return "PrimeBorwein ";
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 lgN = 31 - Integer.numberOfLeadingZeros(n);
30: int piN = 2 + (15 * n) / (8 * (lgN - 1));
31:
32: primeList = new int[piN];
33: multiList = new int[piN];
34:
35: int len = primeFactors(n);
36: return repeatedSquare(len, 1).shiftLeft(n - Integer.bitCount(n));
37: }
38:
39: private Xint repeatedSquare(int len, int k)
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: return Xint.product(primeList, i, len - i).toPowerOf(k).multiply(repeatedSquare(i, 2 * k));
59: }
60:
61: private int primeFactors(int n)
62: {
63: IPrimeSieve sieve = new PrimeSieve(n);
64: IPrimeIteration pIter = sieve.getIteration(3, n);
65:
66: int maxBound = n / 2, count = 0;
67:
68: for (int prime : pIter)
69: {
70: int m = prime > maxBound ? 1 : 0;
71:
72: if (prime <= maxBound)
73: {
74: int q = n;
75: while (q >= prime)
76: {
77: m += q /= prime;
78: }
79: }
80:
81: primeList[count] = prime;
82: multiList[count++] = m;
83: }
84: return count;
85: }
86: } // endOfFactorialPrimeBorwein