mirror of
				https://github.com/saitohirga/WSJT-X.git
				synced 2025-11-04 05:50:31 -05:00 
			
		
		
		
	
		
			
	
	
		
			510 lines
		
	
	
		
			14 KiB
		
	
	
	
		
			C++
		
	
	
	
	
	
		
		
			
		
	
	
			510 lines
		
	
	
		
			14 KiB
		
	
	
	
		
			C++
		
	
	
	
	
	
| 
								 | 
							
								//  Copyright Jeremy Murphy 2016.
							 | 
						|||
| 
								 | 
							
								//  Use, modification and distribution are subject to the
							 | 
						|||
| 
								 | 
							
								//  Boost Software License, Version 1.0. (See accompanying file
							 | 
						|||
| 
								 | 
							
								//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								#ifdef _MSC_VER
							 | 
						|||
| 
								 | 
							
								#  pragma warning (disable : 4224)
							 | 
						|||
| 
								 | 
							
								#endif
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								#include <boost/math/common_factor_rt.hpp>
							 | 
						|||
| 
								 | 
							
								#include <boost/math/special_functions/prime.hpp>
							 | 
						|||
| 
								 | 
							
								#include <boost/multiprecision/cpp_int.hpp>
							 | 
						|||
| 
								 | 
							
								#include <boost/multiprecision/integer.hpp>
							 | 
						|||
| 
								 | 
							
								#include <boost/random.hpp>
							 | 
						|||
| 
								 | 
							
								#include <boost/array.hpp>
							 | 
						|||
| 
								 | 
							
								#include <iostream>
							 | 
						|||
| 
								 | 
							
								#include <algorithm>
							 | 
						|||
| 
								 | 
							
								#include <numeric>
							 | 
						|||
| 
								 | 
							
								#include <string>
							 | 
						|||
| 
								 | 
							
								#include <tuple>
							 | 
						|||
| 
								 | 
							
								#include <type_traits>
							 | 
						|||
| 
								 | 
							
								#include <vector>
							 | 
						|||
| 
								 | 
							
								#include <functional>
							 | 
						|||
| 
								 | 
							
								#include "fibonacci.hpp"
							 | 
						|||
| 
								 | 
							
								#include "../../test/table_type.hpp"
							 | 
						|||
| 
								 | 
							
								#include "table_helper.hpp"
							 | 
						|||
| 
								 | 
							
								#include "performance.hpp"
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								using namespace std;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								boost::multiprecision::cpp_int total_sum(0);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <typename Func, class Table>
							 | 
						|||
| 
								 | 
							
								double exec_timed_test_foo(Func f, const Table& data, double min_elapsed = 0.5)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								    double t = 0;
							 | 
						|||
| 
								 | 
							
								    unsigned repeats = 1;
							 | 
						|||
| 
								 | 
							
								    typename Table::value_type::first_type sum{0};
							 | 
						|||
| 
								 | 
							
								    stopwatch<boost::chrono::high_resolution_clock> w;
							 | 
						|||
| 
								 | 
							
								    do
							 | 
						|||
| 
								 | 
							
								    {
							 | 
						|||
| 
								 | 
							
								       for(unsigned count = 0; count < repeats; ++count)
							 | 
						|||
| 
								 | 
							
								       {
							 | 
						|||
| 
								 | 
							
								          for(typename Table::size_type n = 0; n < data.size(); ++n)
							 | 
						|||
| 
								 | 
							
								            sum += f(data[n].first, data[n].second);
							 | 
						|||
| 
								 | 
							
								       }
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								        t = boost::chrono::duration_cast<boost::chrono::duration<double>>(w.elapsed()).count();
							 | 
						|||
| 
								 | 
							
								        if(t < min_elapsed)
							 | 
						|||
| 
								 | 
							
								            repeats *= 2;
							 | 
						|||
| 
								 | 
							
								    }
							 | 
						|||
| 
								 | 
							
								    while(t < min_elapsed);
							 | 
						|||
| 
								 | 
							
								    total_sum += sum;
							 | 
						|||
| 
								 | 
							
								    return t / repeats;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <typename T>
							 | 
						|||
| 
								 | 
							
								struct test_function_template
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								    vector<pair<T, T> > const & data;
							 | 
						|||
| 
								 | 
							
								    const char* data_name;
							 | 
						|||
| 
								 | 
							
								    
							 | 
						|||
| 
								 | 
							
								    test_function_template(vector<pair<T, T> > const &data, const char* name) : data(data), data_name(name) {}
							 | 
						|||
| 
								 | 
							
								    
							 | 
						|||
| 
								 | 
							
								    template <typename Function>
							 | 
						|||
| 
								 | 
							
								    void operator()(pair<Function, string> const &f) const
							 | 
						|||
| 
								 | 
							
								    {
							 | 
						|||
| 
								 | 
							
								        auto result = exec_timed_test_foo(f.first, data);
							 | 
						|||
| 
								 | 
							
								        auto table_name = string("gcd method comparison with ") + compiler_name() + string(" on ") + platform_name();
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								        report_execution_time(result, 
							 | 
						|||
| 
								 | 
							
								                            table_name,
							 | 
						|||
| 
								 | 
							
								                            string(data_name), 
							 | 
						|||
| 
								 | 
							
								                            string(f.second) + "\n" + boost_name());
							 | 
						|||
| 
								 | 
							
								    }
							 | 
						|||
| 
								 | 
							
								};
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								boost::random::mt19937 rng;
							 | 
						|||
| 
								 | 
							
								boost::random::uniform_int_distribution<> d_0_6(0, 6);
							 | 
						|||
| 
								 | 
							
								boost::random::uniform_int_distribution<> d_1_20(1, 20);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								T get_prime_products()
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   int n_primes = d_0_6(rng);
							 | 
						|||
| 
								 | 
							
								   switch(n_primes)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								   case 0:
							 | 
						|||
| 
								 | 
							
								      // Generate a power of 2:
							 | 
						|||
| 
								 | 
							
								      return static_cast<T>(1u) << d_1_20(rng);
							 | 
						|||
| 
								 | 
							
								   case 1:
							 | 
						|||
| 
								 | 
							
								      // prime number:
							 | 
						|||
| 
								 | 
							
								      return boost::math::prime(d_1_20(rng) + 3);
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   T result = 1;
							 | 
						|||
| 
								 | 
							
								   for(int i = 0; i < n_primes; ++i)
							 | 
						|||
| 
								 | 
							
								      result *= boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3);
							 | 
						|||
| 
								 | 
							
								   return result;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								T get_uniform_random()
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   static boost::random::uniform_int_distribution<T> minmax(std::numeric_limits<T>::min(), std::numeric_limits<T>::max());
							 | 
						|||
| 
								 | 
							
								   return minmax(rng);
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								inline bool even(T const& val)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   return !(val & 1u);
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates>
							 | 
						|||
| 
								 | 
							
								inline bool even(boost::multiprecision::number<Backend, ExpressionTemplates> const& val)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   return !bit_test(val, 0);
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								T euclid_textbook(T a, T b)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   using std::swap;
							 | 
						|||
| 
								 | 
							
								   if(a < b)
							 | 
						|||
| 
								 | 
							
								      swap(a, b);
							 | 
						|||
| 
								 | 
							
								   while(b)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      T t = b;
							 | 
						|||
| 
								 | 
							
								      b = a % b;
							 | 
						|||
| 
								 | 
							
								      a = t;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   return a;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								T binary_textbook(T u, T v)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   if(u && v)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      unsigned shifts = std::min(boost::multiprecision::lsb(u), boost::multiprecision::lsb(v));
							 | 
						|||
| 
								 | 
							
								      if(shifts)
							 | 
						|||
| 
								 | 
							
								      {
							 | 
						|||
| 
								 | 
							
								         u >>= shifts;
							 | 
						|||
| 
								 | 
							
								         v >>= shifts;
							 | 
						|||
| 
								 | 
							
								      }
							 | 
						|||
| 
								 | 
							
								      while(u)
							 | 
						|||
| 
								 | 
							
								      {
							 | 
						|||
| 
								 | 
							
								         unsigned bit_index = boost::multiprecision::lsb(u);
							 | 
						|||
| 
								 | 
							
								         if(bit_index)
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            u >>= bit_index;
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								         else if(bit_index = boost::multiprecision::lsb(v))
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            v >>= bit_index;
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								         else
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            if(u < v)
							 | 
						|||
| 
								 | 
							
								               v = (v - u) >> 1u;
							 | 
						|||
| 
								 | 
							
								            else
							 | 
						|||
| 
								 | 
							
								               u = (u - v) >> 1u;
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								      }
							 | 
						|||
| 
								 | 
							
								      return v << shifts;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   return u + v;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								//
							 | 
						|||
| 
								 | 
							
								// The Mixed Binary Euclid Algorithm
							 | 
						|||
| 
								 | 
							
								// Sidi Mohamed Sedjelmaci
							 | 
						|||
| 
								 | 
							
								// Electronic Notes in Discrete Mathematics 35 (2009) 169<36>176
							 | 
						|||
| 
								 | 
							
								//
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								T mixed_binary_gcd(T u, T v)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   using std::swap;
							 | 
						|||
| 
								 | 
							
								   if(u < v)
							 | 
						|||
| 
								 | 
							
								      swap(u, v);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   unsigned shifts = 0;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   if(!u)
							 | 
						|||
| 
								 | 
							
								      return v;
							 | 
						|||
| 
								 | 
							
								   if(!v)
							 | 
						|||
| 
								 | 
							
								      return u;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   while(even(u) && even(v))
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      u >>= 1u;
							 | 
						|||
| 
								 | 
							
								      v >>= 1u;
							 | 
						|||
| 
								 | 
							
								      ++shifts;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   while(v > 1)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      u %= v;
							 | 
						|||
| 
								 | 
							
								      v -= u;
							 | 
						|||
| 
								 | 
							
								      if(!u)
							 | 
						|||
| 
								 | 
							
								         return v << shifts;
							 | 
						|||
| 
								 | 
							
								      if(!v)
							 | 
						|||
| 
								 | 
							
								         return u << shifts;
							 | 
						|||
| 
								 | 
							
								      while(even(u)) u >>= 1u;
							 | 
						|||
| 
								 | 
							
								      while(even(v)) v >>= 1u;
							 | 
						|||
| 
								 | 
							
								      if(u < v)
							 | 
						|||
| 
								 | 
							
								         swap(u, v);
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   return (v == 1 ? v : u) << shifts;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								void test_type(const char* name)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   using namespace boost::math::detail;
							 | 
						|||
| 
								 | 
							
								   typedef T int_type;
							 | 
						|||
| 
								 | 
							
								   std::vector<pair<int_type, int_type> > data;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   for(unsigned i = 0; i < 1000; ++i)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      data.push_back(std::make_pair(get_prime_products<T>(), get_prime_products<T>()));
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   std::string row_name("gcd<");
							 | 
						|||
| 
								 | 
							
								   row_name += name;
							 | 
						|||
| 
								 | 
							
								   row_name += "> (random prime number products)";
							 | 
						|||
| 
								 | 
							
								   
							 | 
						|||
| 
								 | 
							
								   typedef pair< function<int_type(int_type, int_type)>, string> f_test;
							 | 
						|||
| 
								 | 
							
								   array<f_test, 5> test_functions{ { 
							 | 
						|||
| 
								 | 
							
								      { Stein_gcd<int_type>, "Stein_gcd" } ,
							 | 
						|||
| 
								 | 
							
								      { Euclid_gcd<int_type>, "Euclid_gcd" },
							 | 
						|||
| 
								 | 
							
								      { binary_textbook<int_type>, "Stein_gcd_textbook" },
							 | 
						|||
| 
								 | 
							
								      { euclid_textbook<int_type>, "gcd_euclid_textbook" },
							 | 
						|||
| 
								 | 
							
								      { mixed_binary_gcd<int_type>, "mixed_binary_gcd" } 
							 | 
						|||
| 
								 | 
							
								   } };
							 | 
						|||
| 
								 | 
							
								   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   data.clear();
							 | 
						|||
| 
								 | 
							
								   for(unsigned i = 0; i < 1000; ++i)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      data.push_back(std::make_pair(get_uniform_random<T>(), get_uniform_random<T>()));
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   row_name.erase();
							 | 
						|||
| 
								 | 
							
								   row_name += "gcd<";
							 | 
						|||
| 
								 | 
							
								   row_name += name;
							 | 
						|||
| 
								 | 
							
								   row_name += "> (uniform random numbers)";
							 | 
						|||
| 
								 | 
							
								   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(data, row_name.c_str()));
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   // Fibonacci number tests:
							 | 
						|||
| 
								 | 
							
								   row_name.erase();
							 | 
						|||
| 
								 | 
							
								   row_name += "gcd<";
							 | 
						|||
| 
								 | 
							
								   row_name += name;
							 | 
						|||
| 
								 | 
							
								   row_name += "> (adjacent Fibonacci numbers)";
							 | 
						|||
| 
								 | 
							
								   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_1<T>(), row_name.c_str()));
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   row_name.erase();
							 | 
						|||
| 
								 | 
							
								   row_name += "gcd<";
							 | 
						|||
| 
								 | 
							
								   row_name += name;
							 | 
						|||
| 
								 | 
							
								   row_name += "> (permutations of Fibonacci numbers)";
							 | 
						|||
| 
								 | 
							
								   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(fibonacci_numbers_permution_2<T>(), row_name.c_str()));
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   row_name.erase();
							 | 
						|||
| 
								 | 
							
								   row_name += "gcd<";
							 | 
						|||
| 
								 | 
							
								   row_name += name;
							 | 
						|||
| 
								 | 
							
								   row_name += "> (Trivial cases)";
							 | 
						|||
| 
								 | 
							
								   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(trivial_gcd_test_cases<T>(), row_name.c_str()));
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								/*******************************************************************************************************************/
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <class T>
							 | 
						|||
| 
								 | 
							
								T generate_random(unsigned bits_wanted)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   static boost::random::mt19937 gen;
							 | 
						|||
| 
								 | 
							
								   typedef boost::random::mt19937::result_type random_type;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   T max_val;
							 | 
						|||
| 
								 | 
							
								   unsigned digits;
							 | 
						|||
| 
								 | 
							
								   if(std::numeric_limits<T>::is_bounded && (bits_wanted == (unsigned)std::numeric_limits<T>::digits))
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      max_val = (std::numeric_limits<T>::max)();
							 | 
						|||
| 
								 | 
							
								      digits = std::numeric_limits<T>::digits;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   else
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      max_val = T(1) << bits_wanted;
							 | 
						|||
| 
								 | 
							
								      digits = bits_wanted;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   unsigned bits_per_r_val = std::numeric_limits<random_type>::digits - 1;
							 | 
						|||
| 
								 | 
							
								   while((random_type(1) << bits_per_r_val) > (gen.max)()) --bits_per_r_val;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   unsigned terms_needed = digits / bits_per_r_val + 1;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   T val = 0;
							 | 
						|||
| 
								 | 
							
								   for(unsigned i = 0; i < terms_needed; ++i)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      val *= (gen.max)();
							 | 
						|||
| 
								 | 
							
								      val += gen();
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   val %= max_val;
							 | 
						|||
| 
								 | 
							
								   return val;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <typename N>
							 | 
						|||
| 
								 | 
							
								N gcd_stein(N m, N n)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   BOOST_ASSERT(m >= static_cast<N>(0));
							 | 
						|||
| 
								 | 
							
								   BOOST_ASSERT(n >= static_cast<N>(0));
							 | 
						|||
| 
								 | 
							
								   if(m == N(0)) return n;
							 | 
						|||
| 
								 | 
							
								   if(n == N(0)) return m;
							 | 
						|||
| 
								 | 
							
								           // m > 0 && n > 0
							 | 
						|||
| 
								 | 
							
								   unsigned d_m = 0;
							 | 
						|||
| 
								 | 
							
								   while(even(m)) { m >>= 1; d_m++; }
							 | 
						|||
| 
								 | 
							
								   unsigned d_n = 0;
							 | 
						|||
| 
								 | 
							
								   while(even(n)) { n >>= 1; d_n++; }
							 | 
						|||
| 
								 | 
							
								           // odd(m) && odd(n)
							 | 
						|||
| 
								 | 
							
								      while(m != n) {
							 | 
						|||
| 
								 | 
							
								      if(n > m) swap(n, m);
							 | 
						|||
| 
								 | 
							
								      m -= n;
							 | 
						|||
| 
								 | 
							
								      do m >>= 1; while(even(m));
							 | 
						|||
| 
								 | 
							
								                  // m == n
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   return m << std::min(d_m, d_n);
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								boost::multiprecision::cpp_int big_gcd(const boost::multiprecision::cpp_int& a, const boost::multiprecision::cpp_int& b)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   return boost::multiprecision::gcd(a, b);
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								namespace boost { namespace multiprecision { namespace backends {
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1>
							 | 
						|||
| 
								 | 
							
								inline typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value>::type
							 | 
						|||
| 
								 | 
							
								   eval_gcd_new(
							 | 
						|||
| 
								 | 
							
								      cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& result, 
							 | 
						|||
| 
								 | 
							
								      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& a, 
							 | 
						|||
| 
								 | 
							
								      const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>& b)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   using default_ops::eval_lsb;
							 | 
						|||
| 
								 | 
							
								   using default_ops::eval_is_zero;
							 | 
						|||
| 
								 | 
							
								   using default_ops::eval_get_sign;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   if(a.size() == 1)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      eval_gcd(result, b, *a.limbs());
							 | 
						|||
| 
								 | 
							
								      return;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   if(b.size() == 1)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      eval_gcd(result, a, *b.limbs());
							 | 
						|||
| 
								 | 
							
								      return;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   int shift;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> u(a), v(b), mod;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   int s = eval_get_sign(u);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   /* GCD(0,x) := x */
							 | 
						|||
| 
								 | 
							
								   if(s < 0)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      u.negate();
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   else if(s == 0)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      result = v;
							 | 
						|||
| 
								 | 
							
								      return;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   s = eval_get_sign(v);
							 | 
						|||
| 
								 | 
							
								   if(s < 0)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      v.negate();
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								   else if(s == 0)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      result = u;
							 | 
						|||
| 
								 | 
							
								      return;
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   /* Let shift := lg K, where K is the greatest power of 2
							 | 
						|||
| 
								 | 
							
								   dividing both u and v. */
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   unsigned us = eval_lsb(u);
							 | 
						|||
| 
								 | 
							
								   unsigned vs = eval_lsb(v);
							 | 
						|||
| 
								 | 
							
								   shift = (std::min)(us, vs);
							 | 
						|||
| 
								 | 
							
								   eval_right_shift(u, us);
							 | 
						|||
| 
								 | 
							
								   eval_right_shift(v, vs);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   // From now on access u and v via pointers, that way we have a trivial swap:
							 | 
						|||
| 
								 | 
							
								   cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>* up(&u), *vp(&v), *mp(&mod);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   do 
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      /* Now u and v are both odd, so diff(u, v) is even.
							 | 
						|||
| 
								 | 
							
								      Let u = min(u, v), v = diff(u, v)/2. */
							 | 
						|||
| 
								 | 
							
								      s = up->compare(*vp);
							 | 
						|||
| 
								 | 
							
								      if(s > 0)
							 | 
						|||
| 
								 | 
							
								         std::swap(up, vp);
							 | 
						|||
| 
								 | 
							
								      if(s == 0)
							 | 
						|||
| 
								 | 
							
								         break;
							 | 
						|||
| 
								 | 
							
								      if(vp->size() <= 2)
							 | 
						|||
| 
								 | 
							
								      {
							 | 
						|||
| 
								 | 
							
								         if(vp->size() == 1)
							 | 
						|||
| 
								 | 
							
								            *up = mixed_binary_gcd(*vp->limbs(), *up->limbs());
							 | 
						|||
| 
								 | 
							
								         else
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            double_limb_type i, j;
							 | 
						|||
| 
								 | 
							
								            i = vp->limbs()[0] | (static_cast<double_limb_type>(vp->limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
							 | 
						|||
| 
								 | 
							
								            j = (up->size() == 1) ? *up->limbs() : up->limbs()[0] | (static_cast<double_limb_type>(up->limbs()[1]) << sizeof(limb_type) * CHAR_BIT);
							 | 
						|||
| 
								 | 
							
								            u = mixed_binary_gcd(i, j);
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								         break;
							 | 
						|||
| 
								 | 
							
								      }
							 | 
						|||
| 
								 | 
							
								      if(vp->size() > up->size() /*eval_msb(*vp) > eval_msb(*up) + 32*/)
							 | 
						|||
| 
								 | 
							
								      {
							 | 
						|||
| 
								 | 
							
								         eval_modulus(*mp, *vp, *up);
							 | 
						|||
| 
								 | 
							
								         std::swap(vp, mp);
							 | 
						|||
| 
								 | 
							
								         eval_subtract(*up, *vp);
							 | 
						|||
| 
								 | 
							
								         if(eval_is_zero(*vp) == 0)
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            vs = eval_lsb(*vp);
							 | 
						|||
| 
								 | 
							
								            eval_right_shift(*vp, vs);
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								         else
							 | 
						|||
| 
								 | 
							
								            break;
							 | 
						|||
| 
								 | 
							
								         if(eval_is_zero(*up) == 0)
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            vs = eval_lsb(*up);
							 | 
						|||
| 
								 | 
							
								            eval_right_shift(*up, vs);
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								         else
							 | 
						|||
| 
								 | 
							
								         {
							 | 
						|||
| 
								 | 
							
								            std::swap(up, vp);
							 | 
						|||
| 
								 | 
							
								            break;
							 | 
						|||
| 
								 | 
							
								         }
							 | 
						|||
| 
								 | 
							
								      }
							 | 
						|||
| 
								 | 
							
								      else
							 | 
						|||
| 
								 | 
							
								      {
							 | 
						|||
| 
								 | 
							
								         eval_subtract(*vp, *up);
							 | 
						|||
| 
								 | 
							
								         vs = eval_lsb(*vp);
							 | 
						|||
| 
								 | 
							
								         eval_right_shift(*vp, vs);
							 | 
						|||
| 
								 | 
							
								      }
							 | 
						|||
| 
								 | 
							
								   } 
							 | 
						|||
| 
								 | 
							
								   while(true);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   result = *up;
							 | 
						|||
| 
								 | 
							
								   eval_left_shift(result, shift);
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								}}}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								boost::multiprecision::cpp_int big_gcd_new(const boost::multiprecision::cpp_int& a, const boost::multiprecision::cpp_int& b)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   boost::multiprecision::cpp_int result;
							 | 
						|||
| 
								 | 
							
								   boost::multiprecision::backends::eval_gcd_new(result.backend(), a.backend(), b.backend());
							 | 
						|||
| 
								 | 
							
								   return result;
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								#if 0
							 | 
						|||
| 
								 | 
							
								void test_n_bits(unsigned n, std::string data_name, const std::vector<pair<boost::multiprecision::cpp_int, boost::multiprecision::cpp_int> >* p_data = 0)
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								   using namespace boost::math::detail;
							 | 
						|||
| 
								 | 
							
								   typedef boost::multiprecision::cpp_int int_type;
							 | 
						|||
| 
								 | 
							
								   std::vector<pair<int_type, int_type> > data, data2;
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   for(unsigned i = 0; i < 1000; ++i)
							 | 
						|||
| 
								 | 
							
								   {
							 | 
						|||
| 
								 | 
							
								      data.push_back(std::make_pair(generate_random<int_type>(n), generate_random<int_type>(n)));
							 | 
						|||
| 
								 | 
							
								   }
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								   typedef pair< function<int_type(int_type, int_type)>, string> f_test;
							 | 
						|||
| 
								 | 
							
								   array<f_test, 2> test_functions{ { /*{ Stein_gcd<int_type>, "Stein_gcd" } ,{ Euclid_gcd<int_type>, "Euclid_gcd" },{ binary_textbook<int_type>, "Stein_gcd_textbook" },{ euclid_textbook<int_type>, "gcd_euclid_textbook" },{ mixed_binary_gcd<int_type>, "mixed_binary_gcd" },{ gcd_stein<int_type>, "gcd_stein" },*/{ big_gcd, "boost::multiprecision::gcd" },{ big_gcd_new, "big_gcd_new" } } };
							 | 
						|||
| 
								 | 
							
								   for_each(begin(test_functions), end(test_functions), test_function_template<int_type>(p_data ? *p_data : data, data_name.c_str(), true));
							 | 
						|||
| 
								 | 
							
								}
							 | 
						|||
| 
								 | 
							
								#endif
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								int main()
							 | 
						|||
| 
								 | 
							
								{
							 | 
						|||
| 
								 | 
							
								    test_type<unsigned short>("unsigned short");
							 | 
						|||
| 
								 | 
							
								    test_type<unsigned>("unsigned");
							 | 
						|||
| 
								 | 
							
								    test_type<unsigned long>("unsigned long");
							 | 
						|||
| 
								 | 
							
								    test_type<unsigned long long>("unsigned long long");
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								    test_type<boost::multiprecision::uint256_t>("boost::multiprecision::uint256_t");
							 | 
						|||
| 
								 | 
							
								    test_type<boost::multiprecision::uint512_t>("boost::multiprecision::uint512_t");
							 | 
						|||
| 
								 | 
							
								    test_type<boost::multiprecision::uint1024_t>("boost::multiprecision::uint1024_t");
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								    /*
							 | 
						|||
| 
								 | 
							
								    test_n_bits(16, "   16 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(32, "   32 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(64, "   64 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(125, "  125 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(250, "  250 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(500, "  500 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(1000, " 1000 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(5000, " 5000 bit random values");
							 | 
						|||
| 
								 | 
							
								    test_n_bits(10000, "10000 bit random values");
							 | 
						|||
| 
								 | 
							
								    //test_n_bits(100000);
							 | 
						|||
| 
								 | 
							
								    //test_n_bits(1000000);
							 | 
						|||
| 
								 | 
							
								
							 | 
						|||
| 
								 | 
							
								    test_n_bits(0, "consecutive first 1000 fibonacci numbers", &fibonacci_numbers_cpp_int_permution_1());
							 | 
						|||
| 
								 | 
							
								    test_n_bits(0, "permutations of first 1000 fibonacci numbers", &fibonacci_numbers_cpp_int_permution_2());
							 | 
						|||
| 
								 | 
							
								    */
							 | 
						|||
| 
								 | 
							
								}
							 |