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