1: // Author & Copyright : Peter Luschny
2: // Created: 2010-01-15
3: // License: LGPL version 2.1 or (at your option)
4: // Creative Commons Attribution-ShareAlike 3.0
5:
6: #include "schoenhage.h"
7: #include "xmath.h"
8:
9: void Schoenhage::Factorial(Xint result, ulong n)
10: {
11: lmp::SetUi(result, 1);
12:
13: if (n < 2) { return; }
14:
15: slong iterLen = 0, m = n;
16: ulong lim = n / 2, i;
17: while (m > 0) { m >>= 1; iterLen++; }
18:
19: ulong* primes = 0;
20: ulong piN = Xmath::PrimeSieve(&primes, n);
21: ulong* exponents = lmp::MallocUi(piN);
22: ulong* factors = lmp::MallocUi(piN);
23:
24: exponents[0] = n - Xmath::BitCount(n);
25:
26: for(i = 1; i < piN; i++)
27: {
28: ulong p = primes[i];
29: ulong exp = 1;
30: if (p <= lim)
31: {
32: slong N = n / p;
33: exp = N;
34: while (N > 0)
35: {
36: N /= p;
37: exp += N;
38: }
39: }
40: exponents[i] = exp;
41: }
42:
43: Xint prod; lmp::Init(prod);
44:
45: for (; iterLen >= 0; iterLen--)
46: {
47: ulong count = 0;
48: for (ulong m = 1; m < piN ; m++)
49: {
50: if (((exponents[m] >> iterLen) & 1) == 1)
51: {
52: factors[count++] = primes[m];
53: }
54: }
55:
56: lmp::Pow2(result, result);
57: if (count > 0)
58: {
59: Xmath::Product(prod, factors, 0, count);
60: lmp::Mul(result, result, prod);
61: }
62: }
63:
64: lmp::Mul2Exp(result, result, exponents[0]);
65:
66: lmp::FreeUi(factors, piN);
67: lmp::FreeUi(primes, piN);
68: lmp::FreeUi(exponents, piN);
69: lmp::Clear(prod);
70: }
71: