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 PrimeSwing : IFactorialFunction
15: {
16: public PrimeSwing() { }
17:
18: private PrimeSieve sieve;
19: private int[] primeList;
20:
21: public string Name
22: {
23: get { return "PrimeSwing "; }
24: }
25:
26: public XInt Factorial(int n)
27: {
28: if (n < 20) { return XMath.Factorial(n); }
29:
30: sieve = new PrimeSieve(n);
31: int pLen = (int)(2.0 * (XMath.FloorSqrt(n)
32: + (double)n / (XMath.Log2(n) - 1)));
33: primeList = new int[pLen];
34:
35: int exp2 = n - XMath.BitCount(n);
36: return RecFactorial(n) << exp2;
37: }
38:
39: private XInt RecFactorial(int n)
40: {
41: if (n < 2) return XInt.One;
42:
43: return XInt.Pow(RecFactorial(n / 2), 2) * Swing(n);
44: }
45:
46: private XInt Swing(int n)
47: {
48: if (n < 33) return smallOddSwing[n];
49:
50: int count = 0, rootN = XMath.FloorSqrt(n);
51:
52: var aPrimes = sieve.GetPrimeCollection(3, rootN);
53: var bPrimes = sieve.GetPrimeCollection(rootN + 1, n / 3);
54:
55: foreach (int prime in aPrimes)
56: {
57: int q = n, p = 1;
58:
59: while ((q /= prime) > 0)
60: {
61: if ((q & 1) == 1)
62: {
63: p *= prime;
64: }
65: }
66:
67: if (p > 1)
68: {
69: primeList[count++] = p;
70: }
71: }
72:
73: foreach (int prime in bPrimes)
74: {
75: if (((n / prime) & 1) == 1)
76: {
77: primeList[count++] = prime;
78: }
79: }
80:
81: XInt primorial = sieve.GetPrimorial(n / 2 + 1, n);
82: return primorial * XMath.Product(primeList, 0, count);
83: }
84:
85: private static XInt[] smallOddSwing = {
86: 1,1,1,3,3,15,5,35,35,315,63,693,231,3003,429,6435,6435,109395,
87: 12155,230945,46189,969969,88179,2028117,676039,16900975,1300075,
88: 35102025,5014575,145422675,9694845,300540195,300540195 };
89: }
90: } // endOfFactorialPrimeSwingLuschny