CoxIter 1.3
CoxIter - Computing invariants of hyperbolic Coxeter groups
Loading...
Searching...
No Matches
math_tools.h
Go to the documentation of this file.
1/*
2Copyright (C) 2013-2017
3Rafael Guglielmetti, rafael.guglielmetti@unifr.ch
4*/
5
6/*
7This file is part of CoxIter and AlVin.
8
9CoxIter is free software: you can redistribute it and/or modify
10it under the terms of the GNU General Public License as
11published by the Free Software Foundation, either version 3 of the
12License, or (at your option) any later version.
13
14CoxIter is distributed in the hope that it will be useful,
15but WITHOUT ANY WARRANTY; without even the implied warranty of
16MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17GNU General Public License for more details.
18
19You should have received a copy of the GNU General Public License
20along with CoxIter. If not, see <http://www.gnu.org/licenses/>.
21*/
22
30#ifndef __MATH_TOOLS_H__
31#define __MATH_TOOLS_H__
32
33#include <map>
34#include <string>
35#include <vector>
36#ifdef _USE_LOCAL_GMP_
37#include "gmpxx.h"
38#else
39#include <gmpxx.h>
40#endif
41
42using namespace std;
43
44namespace MathTools {
45extern unsigned int iSmallPrimeNumbers[1229];
46
53bool isPrime(unsigned int n);
54
62template <typename Type>
63typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
64integerSqrt(Type n) {
65 Type place =
66 (Type)1
67 << (sizeof(Type) * 8 -
68 2); // calculated by precompiler = same runtime as: place = 0x40000000
69 while (place > n)
70 place /= 4; // optimized by complier as place >>= 2
71
72 Type root = 0;
73
74 while (place) {
75 if (n >= root + place) {
76 n -= root + place;
77 root += place * 2;
78 }
79
80 root /= 2;
81 place /= 4;
82 }
83
84 return root;
85}
86
93template <typename Type>
94typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
95sqrtSup(Type n) {
96 if (n < 2) // 0 or 1
97 return n;
98
99 Type place =
100 (Type)1
101 << (sizeof(Type) * 8 -
102 2); // calculated by precompiler = same runtime as: place = 0x40000000
103 n--;
104 while (place > n)
105 place /= 4; // optimized by complier as place >>= 2
106
107 Type root = 0;
108
109 while (place) {
110 if (n >= root + place) {
111 n -= root + place;
112 root += place * 2;
113 }
114
115 root /= 2;
116 place /= 4;
117 }
118
119 return (root + 1);
120}
121
129template <typename Type>
130vector<typename std::enable_if<std::is_unsigned<Type>::value, Type>::type>
131listDivisors(const Type &n, const bool &nonTrivialOnly = false) {
132#ifdef _MSC_VER
133 static vector<vector<Type>> divisors_(
134 vector<vector<Type>>(60, vector<Type>(0)));
135 static vector<vector<Type>> divisors_nonTrivialOnly_(
136 vector<vector<Type>>(60, vector<Type>(0)));
137
138 if (n <= 60) {
139 if (nonTrivialOnly && divisors_nonTrivialOnly_[n - 1].size())
140 return divisors_nonTrivialOnly_[n - 1];
141
142 if (!nonTrivialOnly && divisors_[n - 1].size())
143 return divisors_[n - 1];
144 }
145#else // Remark: the following two initialisers don't work on Visual C++ 2013
146 // and 2015
147 static vector<vector<Type>> divisors_ = vector<vector<Type>>(
148 {vector<Type>({1}),
149 vector<Type>({1, 2}),
150 vector<Type>({1, 3}),
151 vector<Type>({1, 2, 4}),
152 vector<Type>({1, 5}),
153 vector<Type>({1, 2, 3, 6}),
154 vector<Type>({1, 7}),
155 vector<Type>({1, 2, 4, 8}),
156 vector<Type>({1, 3, 9}),
157 vector<Type>({1, 2, 5, 10}),
158 vector<Type>({1, 11}),
159 vector<Type>({1, 2, 3, 4, 6, 12}),
160 vector<Type>({1, 13}),
161 vector<Type>({1, 2, 7, 14}),
162 vector<Type>({1, 3, 5, 15}),
163 vector<Type>({1, 2, 4, 8, 16}),
164 vector<Type>({1, 17}),
165 vector<Type>({1, 2, 3, 6, 9, 18}),
166 vector<Type>({1, 19}),
167 vector<Type>({1, 2, 4, 5, 10, 20}),
168 vector<Type>({1, 3, 7, 21}),
169 vector<Type>({1, 2, 11, 22}),
170 vector<Type>({1, 23}),
171 vector<Type>({1, 2, 3, 4, 6, 8, 12, 24}),
172 vector<Type>({1, 5, 25}),
173 vector<Type>({1, 2, 13, 26}),
174 vector<Type>({1, 3, 9, 27}),
175 vector<Type>({1, 2, 4, 7, 14, 28}),
176 vector<Type>({1, 29}),
177 vector<Type>({1, 2, 3, 5, 6, 10, 15, 30}),
178 vector<Type>({1, 31}),
179 vector<Type>({1, 2, 4, 8, 16, 32}),
180 vector<Type>({1, 3, 11, 33}),
181 vector<Type>({1, 2, 17, 34}),
182 vector<Type>({1, 5, 7, 35}),
183 vector<Type>({1, 2, 3, 4, 6, 9, 12, 18, 36}),
184 vector<Type>({1, 37}),
185 vector<Type>({1, 2, 19, 38}),
186 vector<Type>({1, 3, 13, 39}),
187 vector<Type>({1, 2, 4, 5, 8, 10, 20, 40}),
188 vector<Type>({1, 41}),
189 vector<Type>({1, 2, 3, 6, 7, 14, 21, 42}),
190 vector<Type>({1, 43}),
191 vector<Type>({1, 2, 4, 11, 22, 44}),
192 vector<Type>({1, 3, 5, 9, 15, 45}),
193 vector<Type>({1, 2, 23, 46}),
194 vector<Type>({1, 47}),
195 vector<Type>({1, 2, 3, 4, 6, 8, 12, 16, 24, 48}),
196 vector<Type>({1, 7, 49}),
197 vector<Type>({1, 2, 5, 10, 25, 50}),
198 vector<Type>({1, 3, 17, 51}),
199 vector<Type>({1, 2, 4, 13, 26, 52}),
200 vector<Type>({1, 53}),
201 vector<Type>({1, 2, 3, 6, 9, 18, 27, 54}),
202 vector<Type>({1, 5, 11, 55}),
203 vector<Type>({1, 2, 4, 7, 8, 14, 28, 56}),
204 vector<Type>({1, 3, 19, 57}),
205 vector<Type>({1, 2, 29, 58}),
206 vector<Type>({1, 59}),
207 vector<Type>({1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60})});
208 static vector<vector<Type>> divisors_nonTrivialOnly_ =
209 vector<vector<Type>>({vector<Type>(0),
210 vector<Type>(0),
211 vector<Type>(0),
212 vector<Type>({2}),
213 vector<Type>(0),
214 vector<Type>({2, 3}),
215 vector<Type>(0),
216 vector<Type>({2, 4}),
217 vector<Type>({3}),
218 vector<Type>({2, 5}),
219 vector<Type>(0),
220 vector<Type>({2, 3, 4, 6}),
221 vector<Type>(0),
222 vector<Type>({2, 7}),
223 vector<Type>({3, 5}),
224 vector<Type>({2, 4, 8}),
225 vector<Type>(0),
226 vector<Type>({2, 3, 6, 9}),
227 vector<Type>(0),
228 vector<Type>({2, 4, 5, 10}),
229 vector<Type>({3, 7}),
230 vector<Type>({2, 11}),
231 vector<Type>(0),
232 vector<Type>({2, 3, 4, 6, 8, 12}),
233 vector<Type>({5}),
234 vector<Type>({2, 13}),
235 vector<Type>({3, 9}),
236 vector<Type>({2, 4, 7, 14}),
237 vector<Type>(0),
238 vector<Type>({2, 3, 5, 6, 10, 15}),
239 vector<Type>(0),
240 vector<Type>({2, 4, 8, 16}),
241 vector<Type>({3, 11}),
242 vector<Type>({2, 17}),
243 vector<Type>({5, 7}),
244 vector<Type>({2, 3, 4, 6, 9, 12, 18}),
245 vector<Type>(0),
246 vector<Type>({2, 19}),
247 vector<Type>({3, 13}),
248 vector<Type>({2, 4, 5, 8, 10, 20}),
249 vector<Type>(0),
250 vector<Type>({2, 3, 6, 7, 14, 21}),
251 vector<Type>(0),
252 vector<Type>({2, 4, 11, 22}),
253 vector<Type>({3, 5, 9, 15}),
254 vector<Type>({2, 23}),
255 vector<Type>(0),
256 vector<Type>({2, 3, 4, 6, 8, 12, 16, 24}),
257 vector<Type>({7}),
258 vector<Type>({2, 5, 10, 25}),
259 vector<Type>({3, 17}),
260 vector<Type>({2, 4, 13, 26}),
261 vector<Type>(0),
262 vector<Type>({2, 3, 6, 9, 18, 27}),
263 vector<Type>({5, 11}),
264 vector<Type>({2, 4, 7, 8, 14, 28}),
265 vector<Type>({3, 19}),
266 vector<Type>({2, 29}),
267 vector<Type>(0),
268 vector<Type>({2, 3, 4, 5, 6, 10, 12, 15, 20, 30})});
269
270 if (n <= 60) {
271 if (nonTrivialOnly)
272 return divisors_nonTrivialOnly_[n - 1];
273 else
274 return divisors_[n - 1];
275 }
276#endif
277
278 vector<Type> divisors;
279
280 if (!nonTrivialOnly) {
281 divisors.push_back(1);
282 divisors.push_back(n);
283 }
284
285 Type max(integerSqrt(n)), temp;
286 for (Type i(2); i <= max; i++) {
287 if (!(n % i)) {
288 temp = n / i;
289 divisors.push_back(i);
290 if (i != temp)
291 divisors.push_back(temp);
292 }
293 }
294
295#ifdef _MSC_VER
296 if (nonTrivialOnly)
297 divisors_nonTrivialOnly_[n - 1] = iDivisors;
298 else
299 divisors_[n - 1] = iDivisors;
300#endif
301
302 return divisors;
303}
304
312template <typename Type>
313typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
314ugcd(Type u, Type v) {
315 while (v != 0) {
316 unsigned int r = u % v;
317 u = v;
318 v = r;
319 }
320
321 return u;
322}
323
331int jacobiSymbol(int a, unsigned int b);
332
339vector<unsigned int> primeFactorsWithoutSquares(unsigned int n);
340
347template <typename Type> vector<Type> primeFactors(const Type &n) {
348 Type nWork(n < 0 ? -n : n);
349 vector<Type> primes;
350
351 if (nWork == 1)
352 return vector<Type>();
353
354 if (!(nWork % 2)) {
355 primes.push_back(2);
356 nWork /= 2;
357
358 while (!(nWork % 2))
359 nWork /= 2;
360 }
361
362 for (unsigned int i(1); i < 1229 && nWork > 1; i++) {
363 if (!(nWork % iSmallPrimeNumbers[i])) {
364 primes.push_back(iSmallPrimeNumbers[i]);
365 nWork /= iSmallPrimeNumbers[i];
366
367 while (!(nWork % iSmallPrimeNumbers[i]))
368 nWork /= iSmallPrimeNumbers[i];
369 }
370 }
371
372 Type iDivisor(10007);
373 while (nWork > 1) {
374 if (!(nWork % iDivisor)) {
375 primes.push_back(iDivisor);
376 nWork /= iDivisor;
377
378 while (!(nWork % iDivisor))
379 nWork /= iDivisor;
380 }
381
382 iDivisor += 2;
383 }
384
385 return primes;
386}
387
394template <typename Type>
395typename std::enable_if<std::is_unsigned<Type>::value,
396 map<Type, unsigned int>>::type
397primeDecomposition(Type n) {
398 map<Type, unsigned int> primes;
399
400 if (n == 1)
401 return map<Type, unsigned int>();
402
403 if (!(n % 2)) {
404 primes[2] = 1;
405 n /= 2;
406
407 while (!(n % 2)) {
408 n /= 2;
409 primes[2]++;
410 }
411 }
412
413 for (unsigned int i(1); i < 1229 && n > 1; i++) {
414 if (!(n % iSmallPrimeNumbers[i])) {
415 primes[iSmallPrimeNumbers[i]] = 1;
416 n /= iSmallPrimeNumbers[i];
417
418 while (!(n % iSmallPrimeNumbers[i])) {
419 n /= iSmallPrimeNumbers[i];
420 primes[iSmallPrimeNumbers[i]]++;
421 }
422 }
423 }
424
425 Type iDivisor(10007);
426 while (n > 1) {
427 if (!(n % iDivisor)) {
428 primes[iDivisor] = 1;
429 n /= iDivisor;
430
431 while (!(n % iDivisor)) {
432 n /= iDivisor;
433 primes[iDivisor]++;
434 }
435 }
436
437 iDivisor += 2;
438 }
439
440 return primes;
441}
442
449template <typename Type> Type iRemoveSquareFactors(Type n) {
450 Type nWork(n > 0 ? n : -n), res(1), square;
451
452 if (nWork < 3)
453 return n;
454
455 for (unsigned int i(0); i < 1229 && nWork > 1; i++) {
456 square = iSmallPrimeNumbers[i] * iSmallPrimeNumbers[i];
457
458 while (nWork % square == 0)
459 nWork /= square;
460
461 if (nWork % iSmallPrimeNumbers[i] == 0) {
462 nWork /= iSmallPrimeNumbers[i];
463 res *= iSmallPrimeNumbers[i];
464 }
465 }
466
467 Type iDivisor(10007);
468 while (nWork > 1) {
469 square = iDivisor * iDivisor;
470
471 while (nWork % square == 0)
472 nWork /= square;
473
474 if (nWork % iDivisor == 0) {
475 nWork /= iDivisor;
476 res *= iDivisor;
477 }
478
479 iDivisor += 2;
480 }
481
482 return (n < 0 ? -res : res);
483}
484
492template <typename Type>
493inline Type ceilQuotient(const Type &numerator, const Type &denominator) {
494 return ((numerator % denominator) ? numerator / denominator + 1
495 : numerator / denominator);
496}
497
506template <typename Type>
507typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
508sqrtQuotient(const Type &numerator, const Type &denominator) {
509 Type tRes(integerSqrt(numerator / denominator));
510
511 while (denominator * (tRes + 1) * (tRes + 1) <= numerator)
512 tRes++;
513
514 return tRes;
515}
516
517mpz_class sqrtQuotient(const mpz_class &numerator,
518 const mpz_class &denominator);
519mpz_class sqrtSupQuotient(const mpz_class &numerator,
520 const mpz_class &denominator);
521
530template <typename Type>
531typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
532sqrtSupQuotient(const Type &numerator, const Type &denominator) {
533 if (!numerator)
534 return 0;
535
536 Type tRes((numerator % denominator) != 0 ? numerator / denominator + 1
537 : numerator / denominator);
538 tRes = sqrtSup(tRes);
539
540 while (numerator <= denominator * (tRes - 1) * (tRes - 1))
541 tRes--;
542
543 return tRes;
544}
545} // namespace MathTools
546#endif
MathTools
Definition math_tools.cpp:25
MathTools::sqrtQuotient
mpz_class sqrtQuotient(const mpz_class &numerator, const mpz_class &denominator)
Definition math_tools.cpp:137
MathTools::integerSqrt
std::enable_if< std::is_unsigned< Type >::value, Type >::type integerSqrt(Type n)
Definition math_tools.h:64
MathTools::primeFactors
vector< Type > primeFactors(const Type &n)
Definition math_tools.h:347
MathTools::listDivisors
vector< typename std::enable_if< std::is_unsigned< Type >::value, Type >::type > listDivisors(const Type &n, const bool &nonTrivialOnly=false)
Definition math_tools.h:131
MathTools::jacobiSymbol
int jacobiSymbol(int a_, unsigned b)
Definition math_tools.cpp:43
MathTools::ugcd
std::enable_if< std::is_unsigned< Type >::value, Type >::type ugcd(Type u, Type v)
Definition math_tools.h:314
MathTools::iSmallPrimeNumbers
unsigned int iSmallPrimeNumbers[1229]
Definition math_tools.cpp:166
MathTools::sqrtSup
std::enable_if< std::is_unsigned< Type >::value, Type >::type sqrtSup(Type n)
Definition math_tools.h:95
MathTools::ceilQuotient
Type ceilQuotient(const Type &numerator, const Type &denominator)
Definition math_tools.h:493
MathTools::primeDecomposition
std::enable_if< std::is_unsigned< Type >::value, map< Type, unsignedint > >::type primeDecomposition(Type n)
Definition math_tools.h:397
MathTools::primeFactorsWithoutSquares
vector< unsigned int > primeFactorsWithoutSquares(unsigned int n)
Definition math_tools.cpp:100
MathTools::isPrime
bool isPrime(unsigned int n)
Definition math_tools.cpp:27
MathTools::sqrtSupQuotient
mpz_class sqrtSupQuotient(const mpz_class &numerator, const mpz_class &denominator)
Definition math_tools.cpp:148
MathTools::iRemoveSquareFactors
Type iRemoveSquareFactors(Type n)
Definition math_tools.h:449
Generated by doxygen 1.9.8