30 #ifndef __MATH_TOOLS_H__ 31 #define __MATH_TOOLS_H__ 36 #ifdef _USE_LOCAL_GMP_ 62 template <
typename Type>
63 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
iSQRT( Type iN )
65 Type place = (Type)1 << (
sizeof (Type) * 8 - 2 );
73 if( iN >= root+place )
92 template <
typename Type>
93 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
iSQRTsup( Type iN )
98 Type place = (Type)1 << (
sizeof (Type) * 8 - 2 );
107 if( iN >= root+place )
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 )
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 ) ) );
136 if( bNonTrivialOnly && iDivisors_bNonTrivialOnly_[ iN - 1 ].size() )
137 return iDivisors_bNonTrivialOnly_[ iN - 1 ];
139 if( !bNonTrivialOnly && iDivisors_[ iN - 1 ].size() )
140 return iDivisors_[ iN - 1 ];
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})} );
148 if( bNonTrivialOnly )
149 return iDivisors_bNonTrivialOnly_[ iN - 1 ];
151 return iDivisors_[ iN - 1 ];
155 vector< Type > iDivisors;
157 if( !bNonTrivialOnly )
159 iDivisors.push_back( 1 );
160 iDivisors.push_back( iN );
163 Type iMax(
iSQRT( iN ) ), iTemp;
164 for( Type i(2); i <= iMax; i++ )
169 iDivisors.push_back( i );
171 iDivisors.push_back( iTemp );
176 if( bNonTrivialOnly )
177 iDivisors_bNonTrivialOnly_[ iN - 1 ] = iDivisors;
179 iDivisors_[ iN - 1 ] = iDivisors;
192 template <
typename Type>
193 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
ugcd( Type u, Type v )
197 unsigned int r = u % v;
228 template <
typename Type>
231 Type iWork( iN < 0 ? -iN : iN );
232 vector< Type > iPrimes;
235 return vector< Type >( );
239 iPrimes.push_back(2);
242 while( !( iWork % 2 ) )
246 for(
unsigned int i(1); i < 1229 && iWork > 1; i++ )
248 if( !( iWork % iSmallPrimeNumbers[i] ) )
250 iPrimes.push_back( iSmallPrimeNumbers[i] );
251 iWork/= iSmallPrimeNumbers[i];
253 while( !( iWork % iSmallPrimeNumbers[i] ) )
254 iWork/= iSmallPrimeNumbers[i];
258 Type iDivisor( 10007 );
261 if( !( iWork % iDivisor ) )
263 iPrimes.push_back( iDivisor );
266 while( !( iWork % iDivisor ) )
282 template <
typename Type>
283 typename std::enable_if<std::is_unsigned<Type>::value, map< Type, unsigned int >>::type
iPrimeDecomposition( Type iN )
285 map< Type, unsigned int > iPrimes;
288 return map< Type, unsigned int >( );
302 for(
unsigned int i(1); i < 1229 && iN > 1; i++ )
304 if( !( iN % iSmallPrimeNumbers[i] ) )
306 iPrimes[iSmallPrimeNumbers[i]] = 1;
307 iN/= iSmallPrimeNumbers[i];
309 while( !( iN % iSmallPrimeNumbers[i] ) )
311 iN/= iSmallPrimeNumbers[i];
312 iPrimes[iSmallPrimeNumbers[i]]++;
317 Type iDivisor( 10007 );
320 if( !( iN % iDivisor ) )
322 iPrimes[iDivisor] = 1;
325 while( !( iN % iDivisor ) )
344 template <
typename Type>
347 Type iWork( iN > 0 ? iN : -iN ), iRes( 1 ), iSquare;
352 for(
unsigned int i(0); i < 1229 && iWork > 1; i++ )
354 iSquare = iSmallPrimeNumbers[i] * iSmallPrimeNumbers[i];
356 while( iWork % iSquare == 0 )
359 if( iWork % iSmallPrimeNumbers[i] == 0 )
361 iWork /= iSmallPrimeNumbers[i];
362 iRes *= iSmallPrimeNumbers[i];
366 Type iDivisor( 10007 );
369 iSquare = iDivisor * iDivisor;
371 while( iWork % iSquare == 0 )
374 if( iWork % iDivisor == 0 )
383 return ( iN < 0 ? -iRes : iRes );
393 template <
typename Type>
394 inline Type
iCeilQuotient(
const Type& iNumerator,
const Type& iDenominator )
396 return ( (iNumerator % iDenominator) ? iNumerator / iDenominator + 1 : iNumerator / iDenominator );
406 template <
typename Type>
407 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
iSQRTQuotient(
const Type& iNumerator,
const Type& iDenominator )
409 Type tRes(
iSQRT( iNumerator / iDenominator ) );
411 while( iDenominator * ( tRes + 1 ) * ( tRes + 1 ) <= iNumerator )
417 mpz_class
iSQRTQuotient(
const mpz_class& iNumerator,
const mpz_class& iDenominator );
418 mpz_class
iSQRTsupQuotient(
const mpz_class& iNumerator,
const mpz_class& iDenominator );
427 template <
typename Type>
428 typename std::enable_if<std::is_unsigned<Type>::value, Type>::type
iSQRTsupQuotient(
const Type& iNumerator,
const Type& iDenominator )
433 Type tRes( (iNumerator % iDenominator) != 0 ? iNumerator / iDenominator + 1 : iNumerator / iDenominator );
436 while( iNumerator <= iDenominator * ( tRes - 1 ) * ( tRes - 1 ) )
type ugcd(type u, type v)
Definition: mpz_rational.cpp:25