Back to the Homepage of Factorial Algorithms. 1: /// -------- ToujoursEnBeta
2: /// Author & Copyright : Peter Luschny
3: /// Created: 2010-03-01
4: /// License: LGPL version 2.1 or (at your option)
5: /// Creative Commons Attribution-ShareAlike 3.0
6:
7: using System;
8: using System.Collections.Generic;
9: using Xint = System.Numerics.BigInteger;
10:
11: namespace Combinatoricum {
12:
13: #region Swing
14:
15: public class Combinatorics
16: {
17: public Combinatorics() { }
18:
19: private Primes prime;
20: private Factors factors;
21:
22: // OEIS: A056040 Swinging factorial
23: public Xint Swing(uint n)
24: {
25: if (n >= smallOddSwing.Length)
26: {
27: prime = new Primes(n);
28: factors = new Factors(n);
29: }
30: return OddSwing(n) << (int)MathFun.BitCount(n);
31: }
32:
33: private Xint OddSwing(uint n)
34: {
35: if (n < smallOddSwing.Length) return smallOddSwing[n];
36:
37: uint rootN = MathFun.FloorSqrt(n);
38: factors.Init();
39:
40: factors.SetMax(rootN);
41: prime.Factorizer(3, rootN, p =>
42: {
43: uint q = n;
44: while ((q /= p) > 0)
45: if ((q & 1) == 1) { factors.Add(p); }
46: });
47:
48: factors.SetMax(n / 3);
49: prime.Factorizer(rootN + 1, n / 3, p =>
50: {
51: if (((n / p) & 1) == 1) { factors.Add(p); }
52: });
53:
54: factors.SetMax(n);
55: prime.Factorizer(n / 2 + 1, n, p =>
56: {
57: factors.Add(p);
58: });
59:
60: return factors.Product();
61: }
62:
63: #endregion
64: #region factorial
65:
66: // OEIS: A000142 Factorial
67: public Xint Factorial(uint n)
68: {
69: int exp2 = (int)(n - MathFun.BitCount(n));
70:
71: if (n < smallOddFactorial.Length)
72: return new Xint(smallOddFactorial[n]) << exp2;
73:
74: if (n >= smallOddSwing.Length)
75: {
76: prime = new Primes(n);
77: factors = new Factors(n);
78: }
79:
80: return OddFactorial(n) << exp2;
81: }
82:
83: private Xint OddFactorial(uint n)
84: {
85: if (n < smallOddFactorial.Length)
86: {
87: return smallOddFactorial[n];
88: }
89:
90: return Xint.Pow(OddFactorial(n / 2), 2) * OddSwing(n);
91: }
92:
93: #endregion
94: #region binomial
95:
96: // OEIS: A007318 Pascal's triangle
97: public Xint Binomial(uint n, uint k)
98: {
99: if (0 > k || k > n) return Xint.Zero;
100: if (k > n / 2) { k = n - k; }
101: if (k < 3)
102: {
103: if (k == 0) return Xint.One;
104: if (k == 1) return new Xint(n);
105: if (k == 2) return new Xint(((ulong)n * (n - 1)) / 2);
106: }
107: if (n == 2 * k) { return Swing(n); }
108:
109: var prime = new Primes(n);
110: var factors = new Factors(n);
111:
112: uint rootN = MathFun.FloorSqrt(n);
113: factors.Init();
114:
115: factors.SetMax(rootN);
116: prime.Factorizer(2, rootN, p =>
117: {
118: uint r = 0, N = n, K = k;
119: while (N > 0)
120: {
121: r = (N % p) < (K % p + r) ? 1u : 0u;
122: if (r == 1)
123: {
124: factors.Add(p);
125: }
126: N /= p; K /= p;
127: }
128: });
129:
130: factors.SetMax(n / 2);
131: prime.Factorizer(rootN + 1, n / 2, p =>
132: {
133: if (n % p < k % p)
134: {
135: factors.Add(p);
136: }
137: });
138:
139: factors.SetMax(n);
140: prime.Factorizer(n - k + 1, n, p =>
141: {
142: factors.Add(p);
143: });
144:
145: return factors.Product();
146: }
147:
148: #endregion
149: #region doublefactorial
150:
151: // OEIS: A006882 Double factorials n!!: a(n) = n*a(n-2).
152: public Xint DoubleFactorial(uint n)
153: {
154: Xint dblFact;
155: uint N = ((n & 1) == 0) ? n / 2 : n + 1;
156:
157: if (n < smallOddDoubleFactorial.Length)
158: {
159: dblFact = (Xint)smallOddDoubleFactorial[n];
160: }
161: else
162: {
163: if (N >= smallOddSwing.Length)
164: {
165: prime = new Primes(N);
166: factors = new Factors(N);
167: }
168: dblFact = OddFactorial(N, n);
169: }
170:
171: if ((n & 1) == 0)
172: {
173: int exp2 = (int)(n - MathFun.BitCount(n / 2));
174: dblFact = dblFact << exp2;
175: }
176: return dblFact;
177: }
178:
179: private Xint OddFactorial(uint n, uint m)
180: {
181: if (n < smallOddFactorial.Length)
182: {
183: return smallOddFactorial[n];
184: }
185:
186: Xint oddFact = OddFactorial(n / 2, m);
187: if (n < m)
188: {
189: oddFact = Xint.Pow(oddFact, 2);
190: }
191:
192: return oddFact * OddSwing(n);
193: }
194:
195: private static ulong[] smallOddSwing = {
196: 1, 1, 1, 3, 3, 15, 5, 35, 35, 315, 63, 693, 231, 3003, 429, 6435, 6435,
197: 109395, 12155, 230945, 46189, 969969, 88179, 2028117, 676039, 16900975,
198: 1300075, 35102025, 5014575, 145422675, 9694845, 300540195, 300540195,
199: 9917826435, 583401555, 20419054425, 2268783825, 83945001525, 4418157975,
200: 172308161025, 34461632205, 1412926920405, 67282234305, 2893136075115,
201: 263012370465, 11835556670925, 514589420475, 24185702762325, 8061900920775,
202: 395033145117975, 15801325804719, 805867616040669, 61989816618513,
203: 3285460280781189, 121683714103007, 6692604275665385, 956086325095055,
204: 54496920530418135, 1879204156221315, 110873045217057585, 7391536347803839,
205: 450883717216034179, 14544636039226909, 916312070471295267, 916312070471295267 };
206:
207: private static ulong[] smallOddFactorial = {
208: 1, 1, 1, 3, 3, 15, 45, 315, 315, 2835, 14175, 155925, 467775, 6081075,
209: 42567525, 638512875, 638512875, 10854718875, 97692469875, 1856156927625,
210: 9280784638125, 194896477400625, 2143861251406875, 49308808782358125,
211: 147926426347074375, 3698160658676859375 };
212:
213: private static ulong[] smallOddDoubleFactorial = {
214: 1, 1, 1, 3, 1, 15, 3, 105, 3, 945, 15, 10395, 45, 135135, 315, 2027025, 315,
215: 34459425, 2835, 654729075, 14175, 13749310575, 155925, 316234143225, 467775,
216: 7905853580625, 6081075, 213458046676875, 42567525, 6190283353629375, 638512875,
217: 191898783962510625, 638512875, 6332659870762850625, 10854718875 };
218: }
219:
220: #endregion
221: #region primes
222:
223: // OEIS: A000040 prime number
224: public class Primes
225: {
226: const int bitsPerInt = 32;
227: const int mask = bitsPerInt - 1;
228: const int log2Int = 5;
229:
230: private static uint[] PrimesOnBits = {
231: 1762821248u, 848611808u, 3299549660u, 2510511646u };
232:
233: private uint[] isComposite;
234: public delegate void Visitor(uint x); // aka function pointer
235:
236: public void Factorizer(uint min, uint max, Visitor visitor)
237: {
238: // isComposite[0] ... isComposite[n] includes
239: // 5 <= prime numbers <= 96*(n+1) + 1
240:
241: if (min <= 2) visitor(2);
242: if (min <= 3) visitor(3);
243:
244: int absPos = (int)((min + (min + 1) % 2) / 3 - 1);
245: int index = absPos / bitsPerInt;
246: int bitPos = absPos % bitsPerInt;
247: bool toggle = (bitPos & 1) == 1;
248: uint prime = (uint)(5 + 3 * (bitsPerInt * index + bitPos) - (bitPos & 1));
249:
250: while (prime <= max)
251: {
252: uint bitField = isComposite[index++] >> bitPos;
253: for (; bitPos < bitsPerInt; bitPos++)
254: {
255: if ((bitField & 1) == 0)
256: {
257: visitor(prime);
258: }
259: prime += (toggle = !toggle) ? 2u : 4u;
260: if (prime > max) return;
261: bitField >>= 1;
262: }
263: bitPos = 0;
264: }
265: }
266:
267: /// Prime number sieve, Eratosthenes (276-194 b. bitPos.)
268: /// This implementation considers only multiples of primes
269: /// greater than 3, so the smallest value has to be mapped to 5.
270: /// Note: There is no multiplication operation in this function
271: /// and *no call to a sqrt* function.
272:
273: public Primes(uint n)
274: {
275: if (n < 386) { isComposite = PrimesOnBits; return; }
276:
277: isComposite = new uint[(n / (3 * bitsPerInt)) + 1];
278: int d1 = 8, d2 = 8, p1 = 3, p2 = 7, s = 7, s2 = 3;
279: int l = 0, c = 1, max = (int)n / 3, inc;
280: bool toggle = false;
281:
282: while (s < max) // -- scan the sieve
283: {
284: // -- if a prime is found ...
285: if ((isComposite[l >> log2Int] & (1u << (l++ & mask))) == 0)
286: {
287: inc = p1 + p2; // -- ... cancel its multiples
288: for (c = s; c < max; c += inc)
289: { // -- ... set c as composite
290: isComposite[c >> log2Int] |= 1u << (c & mask);
291: }
292:
293: for (c = s + s2; c < max; c += inc)
294: {
295: isComposite[c >> log2Int] |= 1u << (c & mask);
296: }
297: }
298:
299: if (toggle = !toggle) { s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2; }
300: else { s += d1; d2 += 8; p1 += 2; p2 += 6; s2 = p1; }
301: }
302: }
303: }
304:
305: #endregion
306: #region support
307:
308: class Factors
309: {
310: private ulong[] factors;
311: private ulong maxProd, prod;
312: private int count;
313:
314: public Factors(uint n)
315: {
316: int mem = (int)(0.63 * n / Math.Log(n));
317: factors = new ulong[mem];
318: }
319:
320: public void Init()
321: {
322: maxProd = 1;
323: prod = 1;
324: count = 0;
325: }
326:
327: public void SetMax(ulong max)
328: {
329: maxProd = ulong.MaxValue / max;
330:
331: if (prod >= maxProd)
332: {
333: factors[count++] = prod;
334: prod = 1;
335: }
336: }
337:
338: public void Add(ulong prime)
339: {
340: if (prod < maxProd) // I learned this clever idea
341: { // from Marco Bodrato.
342: prod *= prime;
343: }
344: else
345: {
346: factors[count++] = prod;
347: prod = prime;
348: }
349: }
350:
351: public Xint Product()
352: {
353: factors[count++] = prod;
354: return MathFun.Product(factors, 0, count);
355: }
356: }
357:
358: public class MathFun
359: {
360: // Calibrate to a Karatsuba treshhold
361: private const int THRESHOLD_PRODUCT_SERIAL = 24;
362:
363: static public Xint ProductSerial(ulong[] seq, int start, int len)
364: {
365: var prod = new Xint(seq[start]);
366:
367: for (int i = start + 1; i < start + len; i++)
368: {
369: prod *= seq[i];
370: }
371:
372: return prod;
373: }
374:
375: static public Xint Product(ulong[] seq, int start, int len)
376: {
377: if (len <= THRESHOLD_PRODUCT_SERIAL)
378: {
379: return ProductSerial(seq, start, len);
380: }
381:
382: int halfLen = len / 2;
383: return Product(seq, start, halfLen) *
384: Product(seq, start + halfLen, len - halfLen);
385: }
386:
387: /// Bit count,
388: /// sometimes referred to as the population count.
389: public static uint BitCount(uint w)
390: {
391: w = w - ((0xaaaaaaaa & w) >> 1);
392: w = (w & 0x33333333) + ((w >> 2) & 0x33333333);
393: w = w + (w >> 4) & 0x0f0f0f0f;
394: w += w >> 8;
395: w += w >> 16;
396:
397: return w & 0xff;
398: }
399:
400: public static uint FloorSqrt(uint n)
401: {
402: uint a, b;
403: a = b = n;
404: do
405: {
406: a = b;
407: b = (n / a + a) / 2;
408: } while (b < a);
409: return a;
410: }
411:
412: private MathFun() { }
413: }
414:
415: #endregion
416: #region test
417:
418: class Test
419: {
420: static void Main()
421: {
422: var combi = new Combinatorics();
423:
424: var swing = combi.Swing(1000);
425: System.Console.WriteLine("Swing: " + swing);
426:
427: var fact = combi.Factorial(1000);
428: System.Console.WriteLine("Factorial: " + fact);
429:
430: var dblfact = combi.DoubleFactorial(1000);
431: System.Console.WriteLine("DoubleFactorial: " + dblfact);
432:
433: var pascal = combi.Binomial(1000, 333);
434: System.Console.WriteLine("Binomial: " + pascal);
435:
436: System.Console.ReadLine();
437: }
438: }
439:
440: #endregion
441: }
Back to the Homepage of Factorial Algorithms.