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