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