1: /// -------- ToujoursEnBeta
2: /// Author & Copyright : Peter Luschny
3: /// License: LGPL version 3.0 or (at your option)
4: /// Creative Commons Attribution-ShareAlike 3.0
5: /// Comments mail to: peter(at)luschny.de
6: /// Created: 2010-03-01
7:
8: namespace Sharith.Math.Factorial
9: {
10: using Sharith.Math.Primes;
11: using XMath = Sharith.Math.MathUtils.XMath;
12: using XInt = Sharith.Arithmetic.XInt;
13:
14: public class PrimeSchoenhage : IFactorialFunction
15: {
16: public PrimeSchoenhage() { }
17:
18: public string Name
19: {
20: get { return "PrimeSchoenhage "; }
21: }
22:
23: private int[] primeList;
24: private int[] multiList;
25:
26: public XInt Factorial(int n)
27: {
28: if (n < 20) { return XMath.Factorial(n); }
29:
30: int lgN = XMath.FloorLog2(n);
31: int piN = 2 + (15 * n) / (8 * (lgN - 1));
32:
33: primeList = new int[piN];
34: multiList = new int[piN];
35:
36: int len = PrimeFactors(n);
37: int exp2 = n - XMath.BitCount(n);
38:
39: return NestedSquare(len) << exp2;
40: }
41:
42: private XInt NestedSquare(int len)
43: {
44: if (len == 0) return XInt.One;
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 >> 1;
56: mult = multiList[i];
57: }
58:
59: return XMath.Product(primeList, i, len - i)
60: * XInt.Pow(NestedSquare(i), 2);
61: }
62:
63: private int PrimeFactors(int n)
64: {
65: var sieve = new PrimeSieve(n);
66: var primes = sieve.GetPrimeCollection(3, n);
67:
68: int maxBound = n / 2, count = 0;
69:
70: foreach (int prime in primes)
71: {
72: int m = prime > maxBound ? 1 : 0;
73:
74: if (prime <= maxBound)
75: {
76: int q = n;
77: while (q >= prime)
78: {
79: m += q /= prime;
80: }
81: }
82:
83: primeList[count] = prime;
84: multiList[count++] = m;
85: }
86: return count;
87: }
88: }
89: } // endOfFactorialPrimeSchoenhage