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: #ifndef XMATH_H_
7: #define XMATH_H_
8:
9: #include <assert.h>
10: #include <boost/thread.hpp>
11: #include "lmp.h"
12: #include "primeswing.h"
13:
14: class Xmath {
15: public:
16:
17: static void Product(Xint result, ulong* a, slong start, slong len)
18: {
19: if (len < 4)
20: {
21: lmp::SetUi(result, a[start]);
22: if (len == 1) return;
23: lmp::MulUi(result, result, a[start + 1]);
24: if (len == 2) return;
25: lmp::MulUi(result, result, a[start + 2]);
26: return;
27: }
28:
29: Xint temp1; lmp::Init(temp1);
30: slong halfLen = len / 2;
31:
32: Product(temp1, a, start, halfLen);
33: Product(result, a, start + halfLen, len - halfLen);
34: lmp::Mul(result, result, temp1);
35:
36: lmp::Clear(temp1);
37: }
38:
39: static const slong PITHRESHOLD = 10000;
40:
41: static void ParallelProduct(Xint result, ulong* a, slong start, slong len)
42: {
43: if(len < PITHRESHOLD)
44: {
45: Xmath::Product(result, a, start, len);
46: return;
47: }
48:
49: slong threadCnt = (slong) boost::thread::hardware_concurrency();
50:
51: if(threadCnt < 2)
52: {
53: Xmath::Product(result, a, start, len);
54: return;
55: }
56:
57: slong halfLen = len / 2;
58: Xint temp1; lmp::Init(temp1);
59:
60: boost::thread prod(Xmath::Product, temp1, a, start, halfLen);
61: Xmath::Product(result, a, start + halfLen, len - halfLen);
62: prod.join();
63: lmp::Mul(result, result, temp1);
64:
65: lmp::Clear(temp1);
66: }
67:
68: /* Prime number sieve, Eratosthenes (276-194 b. t.)
69: * This implementation considers only multiples of primes
70: * greater than 3, so the smallest value has to be mapped to 5.
71: * Note: There is no multiplication operation in this function
72: * and no call to a sqrt function.
73: * Returns the prime numbers not exceeding n, i.e. p_1,p_2,...,p_s where
74: * p_1 = 2 and p_s ≤ n. Additionally the number of primes ≤ n,
75: * is returned. The third argument, plist, is boolean. plist = 1 returns
76: * a list of the primes whereas plist = 0 returns only the number of primes
77: * found without allocating a prime list.
78: */
79:
80: #define IS_PRIME(P) ((isComposite[(P) >> 5] & pow2[(P) & (bitsPerInt - 1)]) == 0)
81: #define SET_COMPOSITE(C) (isComposite[(C) >> 5] |= pow2[(C) & (bitsPerInt - 1)])
82:
83: static slong xPrimeSieve(ulong** Primes, ulong n, int plist)
84: {
85: const int bitsPerInt = 32;
86: int toggle = 0, i, t, bitFieldLength;
87: ulong *isComposite; ulong *primes;
88: ulong d1 = 8, d2 = 8, p1 = 3, p2 = 7, s = 7, s2 = 3;
89: ulong p = 0, k = 1, max = n / 3, piN, inc;
90: ulong pow2[bitsPerInt];
91:
92: /* precondition n > 2 */
93: assert(n > 2);
94:
95: bitFieldLength = (n / (3 * bitsPerInt)) + 1;
96: isComposite = lmp::MallocUi(bitFieldLength);
97: for (i = 0; i < bitFieldLength; i++) isComposite[i] = 0;
98:
99: /* bitfield powers of 2 */
100: pow2[0] = 1; for (i = 1; i < bitsPerInt; i++) pow2[i] = k <<= 1;
101:
102: while (s < max) /* -- scan the sieve */
103: {
104: if (IS_PRIME(p)) /* -- if a prime is found */
105: {
106: inc = p1 + p2; /* -- then cancel its multiples */
107: for (k = s; k < max; k += inc) { SET_COMPOSITE(k); }
108: for (k = (s + s2); k < max; k += inc) { SET_COMPOSITE(k); }
109: }
110:
111: if (toggle = ~toggle) { s += d2; d1 += 16; p1 += 2; p2 += 2; s2 = p2; }
112: else { s += d1; d2 += 8; p1 += 2; p2 += 6; s2 = p1; }
113: p++;
114: }
115:
116: toggle = 0; p = 5; s = 0; piN = 2;
117:
118: if(plist == 1) /* Return list of primes. */
119: {
120: /* Create the array of prime numbers. */
121: k = 2;
122: for (t = 0; t < bitFieldLength; t++)
123: k += BitCount(~isComposite[t]);
124:
125: primes = lmp::MallocUi(k);
126: primes[0] = 2; primes[1] = 3;
127:
128: /* Transform the sieve into the sequence of prime numbers. */
129: while (p <= n)
130: {
131: ulong isc = isComposite[s++];
132: for (t = 0; t < bitsPerInt; t++)
133: {
134: if ((isc & 1) == 0) { primes[piN++] = p; }
135: p += (toggle = ~toggle) ? 2u : 4u;
136: if (p > n) break;
137: isc >>= 1;
138: }
139: }
140: *Primes = primes;
141: }
142: else /* Return only the number of primes. */
143: {
144: while (p <= n)
145: {
146: ulong isc = isComposite[s++];
147: for (t = 0; t < bitsPerInt; t++)
148: {
149: if ((isc & 1) == 0) { piN++; }
150: p += (toggle = ~toggle) ? 2u : 4u;
151: if (p > n) break;
152: isc >>= 1;
153: }
154: }
155: *Primes = NULL;
156: }
157:
158: free(isComposite);
159: return piN;
160: }
161:
162: static slong PrimeSieve(ulong** Primes, slong n)
163: {
164: if (n <= 2)
165: {
166: ulong* primes = (ulong*) malloc(sizeof(ulong));
167: if (n == 2) { primes[0] = 2; *Primes = primes; return 1; }
168: else { primes[0] = 0; *Primes = primes; return 0; }
169: }
170:
171: /* precondition of xPrimeSieve is n > 2 */
172: return xPrimeSieve(Primes, n, 1);
173: }
174:
175: static slong NumberOfPrimes(ulong n)
176: {
177: ulong* Primes = NULL;
178: if (n <= 2) {
179: if (n == 2) { return 1; } else {return 0; }
180: }
181:
182: /* precondition of xPrimeSieve is n > 2 */
183: return xPrimeSieve(&Primes, n, 0);
184: }
185:
186: static void NaiveFactorial(Xint result, slong n)
187: {
188: assert(n <= 20); // TODO: Error handling
189: lmp::SetUi(result, 1);
190: if (n < 2) { return; }
191: ulong nat[] = {
192: 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20};
193: Product(result, nat, 0, n);
194: }
195:
196: static void NaiveDoubleFactorial(Xint result, ulong n)
197: {
198: lmp::SetUi(result, 1);
199:
200: while (n > 1)
201: {
202: lmp::MulUi (result, result, n);
203: n -= 2;
204: }
205: }
206:
207: static void NaiveBinomial(Xint result, ulong n, ulong k)
208: {
209: Xint temp1; lmp::Init(temp1);
210: Xint temp2; lmp::Init(temp2);
211:
212: PrimeSwing::Factorial(result, n);
213: PrimeSwing::Factorial(temp1, k);
214: PrimeSwing::Factorial(temp2, n - k);
215: lmp::Mul (temp1, temp1, temp2);
216: lmp::Div (result, result, temp1);
217:
218: lmp::Clear(temp1);
219: lmp::Clear(temp2);
220: }
221:
222: static slong Sqrt(slong n)
223: {
224: slong a, b;
225: a = b = n;
226: do {
227: a = b;
228: b = (n / a + a) / 2;
229: } while (b < a);
230: return a;
231: }
232:
233: static slong BitCount(slong n)
234: {
235: slong bc = n - ((n >> 1) & 0x55555555);
236: bc = (bc & 0x33333333) + ((bc >> 2) & 0x33333333);
237: bc = (bc + (bc >> 4)) & 0x0f0f0f0f;
238: bc += bc >> 8;
239: bc += bc >> 16;
240: bc = bc & 0x3f;
241: return bc;
242: }
243: };
244:
245: #endif // XMATH_H_