CoxIter  1.3
CoxIter - Computing invariants of hyperbolic Coxeter groups
math_tools.h
Go to the documentation of this file.
1 /*
2 Copyright (C) 2013-2017
3 Rafael Guglielmetti, rafael.guglielmetti@unifr.ch
4 */
5 
6 /*
7 This file is part of CoxIter and AlVin.
8 
9 CoxIter is free software: you can redistribute it and/or modify
10 it under the terms of the GNU General Public License as
11 published by the Free Software Foundation, either version 3 of the
12 License, or (at your option) any later version.
13 
14 CoxIter is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 GNU General Public License for more details.
18 
19 You should have received a copy of the GNU General Public License
20 along 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 
42 using namespace std;
43 
44 namespace MathTools {
45 extern unsigned int iSmallPrimeNumbers[1229];
46 
53 bool isPrime(unsigned int n);
54 
62 template <typename Type>
63 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
64 integerSqrt(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 
93 template <typename Type>
94 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
95 sqrtSup(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 
129 template <typename Type>
130 vector<typename std::enable_if<std::is_unsigned<Type>::value, Type>::type>
131 listDivisors(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 
312 template <typename Type>
313 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
314 ugcd(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 
331 int jacobiSymbol(int a, unsigned int b);
332 
339 vector<unsigned int> primeFactorsWithoutSquares(unsigned int n);
340 
347 template <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 
394 template <typename Type>
395 typename std::enable_if<std::is_unsigned<Type>::value,
396  map<Type, unsigned int>>::type
397 primeDecomposition(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 
449 template <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 
492 template <typename Type>
493 inline Type ceilQuotient(const Type &numerator, const Type &denominator) {
494  return ((numerator % denominator) ? numerator / denominator + 1
495  : numerator / denominator);
496 }
497 
506 template <typename Type>
507 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
508 sqrtQuotient(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 
517 mpz_class sqrtQuotient(const mpz_class &numerator,
518  const mpz_class &denominator);
519 mpz_class sqrtSupQuotient(const mpz_class &numerator,
520  const mpz_class &denominator);
521 
530 template <typename Type>
531 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
532 sqrtSupQuotient(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
ugcd
type ugcd(type u, type v)
Definition: mpz_rational.cpp:25
MathTools
Definition: math_tools.cpp:25
MathTools::sqrtSup
std::enable_if< std::is_unsigned< Type >::value, Type >::type sqrtSup(Type n)
Definition: math_tools.h:95
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 int b)
MathTools::iSmallPrimeNumbers
unsigned int iSmallPrimeNumbers[1229]
Definition: math_tools.cpp:166
MathTools::primeDecomposition
std::enable_if< std::is_unsigned< Type >::value, map< Type, unsigned int > >::type primeDecomposition(Type n)
Definition: math_tools.h:397
MathTools::ceilQuotient
Type ceilQuotient(const Type &numerator, const Type &denominator)
Definition: math_tools.h:493
MathTools::sqrtQuotient
std::enable_if< std::is_unsigned< Type >::value, Type >::type sqrtQuotient(const Type &numerator, const Type &denominator)
Definition: math_tools.h:508
MathTools::primeFactorsWithoutSquares
vector< unsigned int > primeFactorsWithoutSquares(unsigned int n)
Definition: math_tools.cpp:100
MathTools::sqrtSupQuotient
std::enable_if< std::is_unsigned< Type >::value, Type >::type sqrtSupQuotient(const Type &numerator, const Type &denominator)
Definition: math_tools.h:532
MathTools::integerSqrt
std::enable_if< std::is_unsigned< Type >::value, Type >::type integerSqrt(Type n)
Definition: math_tools.h:64
MathTools::isPrime
bool isPrime(unsigned int n)
Definition: math_tools.cpp:27
MathTools::iRemoveSquareFactors
Type iRemoveSquareFactors(Type n)
Definition: math_tools.h:449
MathTools::primeFactors
vector< Type > primeFactors(const Type &n)
Definition: math_tools.h:347
Generated on Sat Oct 21 2023 17:38:30 for CoxIter by doxygen 1.9.1