1: using Xint = System.Numerics.BigInteger;
2: using Xmath = Luschny.Math.MathUtils.Xmath;
3: using Luschny.Math.Primes;
4:
5: namespace Luschny.Math.Factorial
6: {
7: public class FactorialPrimeBorwein : IFactorialFunction
8: {
9: public FactorialPrimeBorwein() { }
10:
11: public string Name
12: {
13: get { return "PrimeBorwein "; }
14: }
15:
16: private int[] primeList;
17: private int[] multiList;
18:
19: public Xint Factorial(int n)
20: {
21: if (n < 0)
22: {
23: throw new System.ArgumentOutOfRangeException("n",
24: Name + ": n >= 0 required, but was " + n);
25: }
26:
27: if (n < 20)
28: {
29: return Xmath.Factorial(n);
30: }
31:
32: int lgN = Xmath.FloorLog2(n);
33: int piN = 2 + (15 * n) / (8 * (lgN - 1));
34:
35: primeList = new int[piN];
36: multiList = new int[piN];
37:
38: int len = PrimeFactors(n);
39: int exp2 = n - Xmath.BitCount(n);
40:
41: return RepeatedSquare(len, 1) << exp2;
42: }
43:
44: private Xint RepeatedSquare(int len, int k)
45: {
46: if (len == 0) return Xint.One;
47:
48: int i = 0, mult = multiList[0];
49:
50: while (mult > 1)
51: {
52: if ((mult & 1) == 1) // is mult odd ?
53: {
54: primeList[len++] = primeList[i];
55: }
56:
57: multiList[i++] = mult >> 1;
58: mult = multiList[i];
59: }
60:
61: Xint p = Xmath.Product(primeList, i, len - i);
62: return Xint.Pow(p, k) * RepeatedSquare(i, 2 * k);
63: }
64:
65: private int PrimeFactors(int n)
66: {
67: PrimeSieve sieve = new PrimeSieve(n);
68: IPrimeCollection primeCollection = sieve.GetPrimeCollection(3, n);
69:
70: int maxBound = n / 2, count = 0;
71:
72: foreach (int prime in primeCollection)
73: {
74: int m = prime > maxBound ? 1 : 0;
75:
76: if (prime <= maxBound)
77: {
78: int q = n;
79: while (q >= prime)
80: {
81: m += q /= prime;
82: }
83: }
84: primeList[count] = prime;
85: multiList[count++] = m;
86: }
87: return count;
88: }
89: }
90:
91: } // endOfFactorialPrimeBorwein