CoxIter  1.2
CoxIter - Computing invariants of hyperbolic Coxeter groups
math_tools.h
Go to the documentation of this file.
1 /*
2 Copyright (C) 2013, 2014, 2015, 2016
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 <string>
34 #include <vector>
35 #include <map>
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 {
46  extern unsigned int iSmallPrimeNumbers[1229];
47 
54  bool bIsPrime( unsigned int iN );
55 
62  template <typename Type>
63  typename std::enable_if<std::is_unsigned<Type>::value, Type>::type iSQRT( Type iN )
64  {
65  Type place = (Type)1 << ( sizeof (Type) * 8 - 2 ); // calculated by precompiler = same runtime as: place = 0x40000000
66  while( place > iN )
67  place /= 4; // optimized by complier as place >>= 2
68 
69  Type root = 0;
70 
71  while( place )
72  {
73  if( iN >= root+place )
74  {
75  iN -= root+place;
76  root += place * 2;
77  }
78 
79  root /= 2;
80  place /= 4;
81  }
82 
83  return root;
84  }
85 
92  template <typename Type>
93  typename std::enable_if<std::is_unsigned<Type>::value, Type>::type iSQRTsup( Type iN )
94  {
95  if( iN < 2 ) // 0 or 1
96  return iN;
97 
98  Type place = (Type)1 << ( sizeof (Type) * 8 - 2 ); // calculated by precompiler = same runtime as: place = 0x40000000
99  iN--;
100  while( place > iN )
101  place /= 4; // optimized by complier as place >>= 2
102 
103  Type root = 0;
104 
105  while( place )
106  {
107  if( iN >= root+place )
108  {
109  iN -= root+place;
110  root += place * 2;
111  }
112 
113  root /= 2;
114  place /= 4;
115  }
116 
117  return (root+1);
118  }
119 
127  template <typename Type>
128  vector< typename std::enable_if<std::is_unsigned<Type>::value, Type>::type > iListDivisors( const Type& iN, const bool& bNonTrivialOnly = false )
129  {
130  #ifdef _MSC_VER
131  static vector< vector< Type > > iDivisors_( vector< vector< Type > >( 60, vector< Type >( 0 ) ) );
132  static vector< vector< Type > > iDivisors_bNonTrivialOnly_( vector< vector< Type > >( 60, vector< Type >( 0 ) ) );
133 
134  if( iN <= 60 )
135  {
136  if( bNonTrivialOnly && iDivisors_bNonTrivialOnly_[ iN - 1 ].size() )
137  return iDivisors_bNonTrivialOnly_[ iN - 1 ];
138 
139  if( !bNonTrivialOnly && iDivisors_[ iN - 1 ].size() )
140  return iDivisors_[ iN - 1 ];
141  }
142  #else // Remark: the following two initialisers don't work on Visual C++ 2013 and 2015
143  static vector< vector< Type > > iDivisors_ = vector< vector< Type > >( {vector<Type>({1}),vector<Type>({1,2}),vector<Type>({1,3}),vector<Type>({1,2,4}),vector<Type>({1,5}),vector<Type>({1,2,3,6}),vector<Type>({1,7}),vector<Type>({1,2,4,8}),vector<Type>({1,3,9}),vector<Type>({1,2,5,10}),vector<Type>({1,11}),vector<Type>({1,2,3,4,6,12}),vector<Type>({1,13}),vector<Type>({1,2,7,14}),vector<Type>({1,3,5,15}),vector<Type>({1,2,4,8,16}),vector<Type>({1,17}),vector<Type>({1,2,3,6,9,18}),vector<Type>({1,19}),vector<Type>({1,2,4,5,10,20}),vector<Type>({1,3,7,21}),vector<Type>({1,2,11,22}),vector<Type>({1,23}),vector<Type>({1,2,3,4,6,8,12,24}),vector<Type>({1,5,25}),vector<Type>({1,2,13,26}),vector<Type>({1,3,9,27}),vector<Type>({1,2,4,7,14,28}),vector<Type>({1,29}),vector<Type>({1,2,3,5,6,10,15,30}),vector<Type>({1,31}),vector<Type>({1,2,4,8,16,32}),vector<Type>({1,3,11,33}),vector<Type>({1,2,17,34}),vector<Type>({1,5,7,35}),vector<Type>({1,2,3,4,6,9,12,18,36}),vector<Type>({1,37}),vector<Type>({1,2,19,38}),vector<Type>({1,3,13,39}),vector<Type>({1,2,4,5,8,10,20,40}),vector<Type>({1,41}),vector<Type>({1,2,3,6,7,14,21,42}),vector<Type>({1,43}),vector<Type>({1,2,4,11,22,44}),vector<Type>({1,3,5,9,15,45}),vector<Type>({1,2,23,46}),vector<Type>({1,47}),vector<Type>({1,2,3,4,6,8,12,16,24,48}),vector<Type>({1,7,49}),vector<Type>({1,2,5,10,25,50}),vector<Type>({1,3,17,51}),vector<Type>({1,2,4,13,26,52}),vector<Type>({1,53}),vector<Type>({1,2,3,6,9,18,27,54}),vector<Type>({1,5,11,55}),vector<Type>({1,2,4,7,8,14,28,56}),vector<Type>({1,3,19,57}),vector<Type>({1,2,29,58}),vector<Type>({1,59}),vector<Type>({1,2,3,4,5,6,10,12,15,20,30,60})} );
144  static vector< vector< Type > > iDivisors_bNonTrivialOnly_ = vector< vector< Type > >( {vector<Type>(0),vector<Type>(0),vector<Type>(0),vector<Type>({2}),vector<Type>(0),vector<Type>({2,3}),vector<Type>(0),vector<Type>({2,4}),vector<Type>({3}),vector<Type>({2,5}),vector<Type>(0),vector<Type>({2,3,4,6}),vector<Type>(0),vector<Type>({2,7}),vector<Type>({3,5}),vector<Type>({2,4,8}),vector<Type>(0),vector<Type>({2,3,6,9}),vector<Type>(0),vector<Type>({2,4,5,10}),vector<Type>({3,7}),vector<Type>({2,11}),vector<Type>(0),vector<Type>({2,3,4,6,8,12}),vector<Type>({5}),vector<Type>({2,13}),vector<Type>({3,9}),vector<Type>({2,4,7,14}),vector<Type>(0),vector<Type>({2,3,5,6,10,15}),vector<Type>(0),vector<Type>({2,4,8,16}),vector<Type>({3,11}),vector<Type>({2,17}),vector<Type>({5,7}),vector<Type>({2,3,4,6,9,12,18}),vector<Type>(0),vector<Type>({2,19}),vector<Type>({3,13}),vector<Type>({2,4,5,8,10,20}),vector<Type>(0),vector<Type>({2,3,6,7,14,21}),vector<Type>(0),vector<Type>({2,4,11,22}),vector<Type>({3,5,9,15}),vector<Type>({2,23}),vector<Type>(0),vector<Type>({2,3,4,6,8,12,16,24}),vector<Type>({7}),vector<Type>({2,5,10,25}),vector<Type>({3,17}),vector<Type>({2,4,13,26}),vector<Type>(0),vector<Type>({2,3,6,9,18,27}),vector<Type>({5,11}),vector<Type>({2,4,7,8,14,28}),vector<Type>({3,19}),vector<Type>({2,29}),vector<Type>(0),vector<Type>({2,3,4,5,6,10,12,15,20,30})} );
145 
146  if( iN <= 60 )
147  {
148  if( bNonTrivialOnly )
149  return iDivisors_bNonTrivialOnly_[ iN - 1 ];
150  else
151  return iDivisors_[ iN - 1 ];
152  }
153  #endif
154 
155  vector< Type > iDivisors;
156 
157  if( !bNonTrivialOnly )
158  {
159  iDivisors.push_back( 1 );
160  iDivisors.push_back( iN );
161  }
162 
163  Type iMax( iSQRT( iN ) ), iTemp;
164  for( Type i(2); i <= iMax; i++ )
165  {
166  if( !( iN % i ) )
167  {
168  iTemp = iN / i;
169  iDivisors.push_back( i );
170  if( i != iTemp )
171  iDivisors.push_back( iTemp );
172  }
173  }
174 
175  #ifdef _MSC_VER
176  if( bNonTrivialOnly )
177  iDivisors_bNonTrivialOnly_[ iN - 1 ] = iDivisors;
178  else
179  iDivisors_[ iN - 1 ] = iDivisors;
180  #endif
181 
182  return iDivisors;
183  }
184 
192  template <typename Type>
193  typename std::enable_if<std::is_unsigned<Type>::value, Type>::type ugcd( Type u, Type v )
194  {
195  while ( v != 0)
196  {
197  unsigned int r = u % v;
198  u = v;
199  v = r;
200  }
201 
202  return u;
203  }
204 
212  int iJacobiSymbol( int iA, unsigned int iB );
213 
220  vector< unsigned int > iPrimeFactorsWithoutSquares( unsigned int iN );
221 
228  template <typename Type>
229  vector<Type> iPrimeFactors( const Type& iN )
230  {
231  Type iWork( iN < 0 ? -iN : iN );
232  vector< Type > iPrimes;
233 
234  if( iWork == 1 )
235  return vector< Type >( );
236 
237  if( !( iWork % 2 ) )
238  {
239  iPrimes.push_back(2);
240  iWork /= 2;
241 
242  while( !( iWork % 2 ) )
243  iWork/= 2;
244  }
245 
246  for( unsigned int i(1); i < 1229 && iWork > 1; i++ )
247  {
248  if( !( iWork % iSmallPrimeNumbers[i] ) )
249  {
250  iPrimes.push_back( iSmallPrimeNumbers[i] );
251  iWork/= iSmallPrimeNumbers[i];
252 
253  while( !( iWork % iSmallPrimeNumbers[i] ) )
254  iWork/= iSmallPrimeNumbers[i];
255  }
256  }
257 
258  Type iDivisor( 10007 );
259  while( iWork > 1 )
260  {
261  if( !( iWork % iDivisor ) )
262  {
263  iPrimes.push_back( iDivisor );
264  iWork/= iDivisor;
265 
266  while( !( iWork % iDivisor ) )
267  iWork/= iDivisor;
268  }
269 
270  iDivisor += 2;
271  }
272 
273  return iPrimes;
274  }
275 
282  template <typename Type>
283  typename std::enable_if<std::is_unsigned<Type>::value, map< Type, unsigned int >>::type iPrimeDecomposition( Type iN )
284  {
285  map< Type, unsigned int > iPrimes;
286 
287  if( iN == 1 )
288  return map< Type, unsigned int >( );
289 
290  if( !( iN % 2 ) )
291  {
292  iPrimes[2] = 1;
293  iN /= 2;
294 
295  while( !( iN % 2 ) )
296  {
297  iN/= 2;
298  iPrimes[2]++;
299  }
300  }
301 
302  for( unsigned int i(1); i < 1229 && iN > 1; i++ )
303  {
304  if( !( iN % iSmallPrimeNumbers[i] ) )
305  {
306  iPrimes[iSmallPrimeNumbers[i]] = 1;
307  iN/= iSmallPrimeNumbers[i];
308 
309  while( !( iN % iSmallPrimeNumbers[i] ) )
310  {
311  iN/= iSmallPrimeNumbers[i];
312  iPrimes[iSmallPrimeNumbers[i]]++;
313  }
314  }
315  }
316 
317  Type iDivisor( 10007 );
318  while( iN > 1 )
319  {
320  if( !( iN % iDivisor ) )
321  {
322  iPrimes[iDivisor] = 1;
323  iN/= iDivisor;
324 
325  while( !( iN % iDivisor ) )
326  {
327  iN/= iDivisor;
328  iPrimes[iDivisor]++;
329  }
330  }
331 
332  iDivisor += 2;
333  }
334 
335  return iPrimes;
336  }
337 
344  template <typename Type>
345  Type iRemoveSquareFactors( Type iN )
346  {
347  Type iWork( iN > 0 ? iN : -iN ), iRes( 1 ), iSquare;
348 
349  if( iWork < 3 )
350  return iN;
351 
352  for( unsigned int i(0); i < 1229 && iWork > 1; i++ )
353  {
354  iSquare = iSmallPrimeNumbers[i] * iSmallPrimeNumbers[i];
355 
356  while( iWork % iSquare == 0 )
357  iWork /= iSquare;
358 
359  if( iWork % iSmallPrimeNumbers[i] == 0 )
360  {
361  iWork /= iSmallPrimeNumbers[i];
362  iRes *= iSmallPrimeNumbers[i];
363  }
364  }
365 
366  Type iDivisor( 10007 );
367  while( iWork > 1 )
368  {
369  iSquare = iDivisor * iDivisor;
370 
371  while( iWork % iSquare == 0 )
372  iWork /= iSquare;
373 
374  if( iWork % iDivisor == 0 )
375  {
376  iWork /= iDivisor;
377  iRes *= iDivisor;
378  }
379 
380  iDivisor += 2;
381  }
382 
383  return ( iN < 0 ? -iRes : iRes );
384  }
385 
393  template <typename Type>
394  inline Type iCeilQuotient( const Type& iNumerator, const Type& iDenominator )
395  {
396  return ( (iNumerator % iDenominator) ? iNumerator / iDenominator + 1 : iNumerator / iDenominator );
397  }
398 
406  template <typename Type>
407  typename std::enable_if<std::is_unsigned<Type>::value, Type>::type iSQRTQuotient( const Type& iNumerator, const Type& iDenominator )
408  {
409  Type tRes( iSQRT( iNumerator / iDenominator ) );
410 
411  while( iDenominator * ( tRes + 1 ) * ( tRes + 1 ) <= iNumerator )
412  tRes++;
413 
414  return tRes;
415  }
416 
417  mpz_class iSQRTQuotient( const mpz_class& iNumerator, const mpz_class& iDenominator );
418  mpz_class iSQRTsupQuotient( const mpz_class& iNumerator, const mpz_class& iDenominator );
419 
427  template <typename Type>
428  typename std::enable_if<std::is_unsigned<Type>::value, Type>::type iSQRTsupQuotient( const Type& iNumerator, const Type& iDenominator )
429  {
430  if( !iNumerator )
431  return 0;
432 
433  Type tRes( (iNumerator % iDenominator) != 0 ? iNumerator / iDenominator + 1 : iNumerator / iDenominator );
434  tRes = iSQRTsup( tRes );
435 
436  while( iNumerator <= iDenominator * ( tRes - 1 ) * ( tRes - 1 ) )
437  tRes--;
438 
439  return tRes;
440  }
441 }
442 #endif
type ugcd(type u, type v)
Definition: mpz_rational.cpp:25
Type iCeilQuotient(const Type &iNumerator, const Type &iDenominator)
Definition: math_tools.h:394
Definition: math_tools.cpp:25
vector< Type > iPrimeFactors(const Type &iN)
Definition: math_tools.h:229
unsigned int iSmallPrimeNumbers[1229]
Definition: math_tools.cpp:175
int iJacobiSymbol(int iA, unsigned int iB)
vector< typename std::enable_if< std::is_unsigned< Type >::value, Type >::type > iListDivisors(const Type &iN, const bool &bNonTrivialOnly=false)
Definition: math_tools.h:128
std::enable_if< std::is_unsigned< Type >::value, Type >::type iSQRTsup(Type iN)
Definition: math_tools.h:93
vector< unsigned int > iPrimeFactorsWithoutSquares(unsigned int iN)
Definition: math_tools.cpp:108
Type iRemoveSquareFactors(Type iN)
Definition: math_tools.h:345
bool bIsPrime(unsigned int iN)
Definition: math_tools.cpp:28
std::enable_if< std::is_unsigned< Type >::value, Type >::type iSQRTsupQuotient(const Type &iNumerator, const Type &iDenominator)
Definition: math_tools.h:428
std::enable_if< std::is_unsigned< Type >::value, Type >::type iSQRTQuotient(const Type &iNumerator, const Type &iDenominator)
Definition: math_tools.h:407
std::enable_if< std::is_unsigned< Type >::value, Type >::type iSQRT(Type iN)
Definition: math_tools.h:63
std::enable_if< std::is_unsigned< Type >::value, map< Type, unsigned int > >::type iPrimeDecomposition(Type iN)
Definition: math_tools.h:283