CxxlMan: 排列組合產生器

整理來自 http://www.programmer-club.com.tw/ShowSameTitleN/c/46015.html 的討論，討論結果大致可分為遞迴法、計算法及量子演算法三種。但遞迴法有 stack overflow 的隱憂，所以改寫成模擬遞迴法；而計算法則因須要使用大數才能完整，所以加入了 goost 的 multiprecision；量子演算法要有量子電腦才能真正發揮實力，這裡使用 thread 去模擬，過這些做法象徵意義大於實質意義。

模擬遞迴法

相較於計算法，在單一執行緒上處理速度較快，但也只能用單一執行緒來跑，此例中使用 std::list、std::vector 及 size_t 使得最大處理數量受其限制

排列(Permutation)

#include <list>
#include <vector>
#include <iostream>
#include <iomanip>

template <typename T> 
std::ostream& operator <<(std::ostream& os, std::vector<T>& vt)
{
  os << "{ ";
  for(std::vector<T>::iterator cit = vt.begin(); 
      cit != vt.end();
      ++cit)
  {
    os << *cit << ' ';
  }
  os << "}";
  return os;
}

void Permutation(const size_t n, const size_t r,
                 void (*pFn)(std::vector<size_t>&))
{
  if(n == 0 || n < r) return;

  std::list<size_t> Samples; // 要排列的樣本
  std::vector<size_t> Index(r);
  for(size_t i = 0; i < n; ++i)
    Samples.push_back(i+1);

  std::vector<size_t> result(r);
  size_t Level = 1;
  bool up = true;

  while(true)
  {
    if(Level == r)
    {
      for(size_t i = Samples.size(); i > 0; --i)
      {
        result[Level-1] = Samples.front();
        Samples.pop_front();
        pFn(result);
        Samples.push_back(result[Level-1]);
      }
      --Level;
      up = false;
    }
    else if(Level == 0)
    {
      return;
    }
    else
    {
      if(up)
      {
        Index[Level-1] = Samples.size();
        result[Level-1] = Samples.front();
        Samples.pop_front();
        ++Level;
      }
      else
      {
        --Index[Level-1];
        Samples.push_back(result[Level-1]);
        if(Index[Level-1] > 0)
        {
          result[Level-1] = Samples.front();
          Samples.pop_front();
          ++Level;
          up = true;
        }
        else
          --Level;
      }
    }
  }
}

unsigned long long m = 1;

// 回叫函數
void Fn(std::vector<size_t>& result)
{
  std::cout << std::setw(10) << m++ << ": " << result << std::endl;
}

int main(int argc, char* argv[])
{
  size_t n, r;
  std::cout << "Permutation nPr" << std::endl;
  while (std::cout << "Enter n r? " << std::flush, std::cin >> n >> r)
  {
    m = 1;
    Permutation(n, r, Fn);
  }

  return 0;
}

組合(Combination)

#include <vector>
#include <iostream>
#include <iomanip>

template <typename T> 
std::ostream& operator <<(std::ostream& os, std::vector<T>& vt)
{
  os << "{ ";
  for(std::vector<T>::iterator cit = vt.begin(); 
      cit != vt.end();
      ++cit)
  {
    os << *cit << ' ';
  }
  os << "}";
  return os;
}

void Combination(const size_t n, const size_t r,
                 void (*pFn)(std::vector<size_t>& vt))
{
  if(n == 0 || n < r) return;

  std::vector<size_t> result(r);
  bool up = true;
  size_t Level = 0;

  while(Level < r)
  {
    result[Level] = Level + 1;
    ++Level;
  }

  size_t b = n-r+1;
  while(true)
  {
    if(Level == r)
    {
      pFn(result);
      --Level;
      up = false;
    }
    else if(Level == 0)
    {
      if(up == false)
      {
        if(result[0] < b)
        {
          up = true;
          ++result[0];
          ++Level;
        }
        else
          return;
      }
    }
    else
    {
      if(up == false)
      {
        if(result[Level] < (b+Level))
        {
          up = true;
          ++result[Level];
          ++Level;
        }
        else
          --Level;
      }
      else // if(up == true)
      {
        result[Level] = result[Level-1]+1;
        ++Level;
      }
    }
  }
}

unsigned long long m = 1;

// 回叫函數
void Fn(std::vector<size_t>& result)
{
  std::cout << std::setw(10) << m++ << ": " << result << std::endl;
}

int main(int argc, char* argv[])
{  
  size_t n, r;
  std::cout << "Calculate nCr" << std::endl;
  while(std::cout << "Enter n r? " << std::flush, std::cin >> n >> r)
  {
    m = 1;
    Combination(n, r, Fn);
  }

  return 0;
}

計算法

這方法雖然處理費時，但每組計算可獨立運作，很適合用平行處理來加速，即可用多個 CPU 或多台電腦一起計算。另外因可以指定要取得哪一個組合的這個特點，可作為不重覆取得亂數的應用。

組合(Combination)

這方法先由 sflam(Raymond) 提出，也做了一般版及大數版，但他的計算法我看不懂，所以我也做了一套，不過我沒寫大數版。另外 sflam(Raymond) 也提供 Python 版，可在討論中找到

sflam(Raymond) 的一般版

--------------------------------------

[參考資料: https://msdn.microsoft.com/en-us/library/aa289166.aspx]
組合程式. 程式用了一些 C++11 的語法如 auto, 編譯器如不支援自己轉
成 std::vector<T>::const_iterator. 
因為用的是內建類型, 計算會有上限, unsigned long long 應該可以
到 64C32, 再大就 overflow 了. 

[code: C++11]
#include <vector>
#include <string>
#include <iostream>
#include <iomanip>

template <typename T> 
std::ostream& operator <<(std::ostream& os, const std::vector<T>& vt)
{
  os << "{ ";
  for (auto cit = vt.cbegin(); cit != vt.cend(); ++cit)
  {
    os << *cit << ' ';
  }
  os << "}";
  return os;
}

namespace combination
{
  unsigned long long choose(const size_t n, const size_t r)
  {
    if(n < r)
      return 0;
    else if (n == r)
      return 1;

    // choose(100, 97) ==> choose(100, 3)
    size_t k2 = (r < n - r) ? r : n - r;

    // n * (n-1) * (n-2) ... (n-k2+1)
    // -------------------------------
    // 1 * 2 * ... k2
    unsigned long long ans = n;
    for (size_t i = 2; i <= k2; ++i)
    {
      ans = (ans * (n - i + 1))/i;
    }
    return ans;
  }

  size_t LargestV(size_t a, size_t b, unsigned long long x)
  {
    long v = a - 1;
    while (choose(v, b) > x)
    --v;
    return v;
  }

  void GetMthElement(const size_t n, const size_t r, unsigned long long m,
                    std::vector<size_t>& result)
  {
    result.resize(r);
    size_t a = n;
    size_t b = r;
    unsigned long long x = (choose(n, r) - 1) - m;
    for (size_t i = 0; i < r; ++i)
    {
      result[i] = LargestV(a, b, x);
      x = x - choose(result[i], b);
      a = result[i];
      b = b - 1;
    }
    for (size_t i = 0; i < r; ++i)
    {
      result[i] = n - result[i];
    }
  }

  void C(const size_t n, const size_t r)
  {
    std::vector<size_t> result(r);
    const unsigned long long mMax = choose(n, r);
    for (unsigned long long m = 0; m < mMax; ++m)
    {
      GetMthElement(n, r, m, result);
      std::cout << std::setw(10) << m+1 << ": " << result << std::endl;
    }
  }
}

int main()
{
  size_t n, r;
  std::cout << "Calculate nCr" << std::endl;
  while (std::cout << "Enter n r? " << std::flush, std::cin >> n >> r)
  {
    combination::C(n, r);
  }
} 
[/code: C++11]

[執行:]
Calculate nCr
Enter n r? 5 3
     1: { 1 2 3 }
     2: { 1 2 4 }
     3: { 1 2 5 }
     4: { 1 3 4 }
     5: { 1 3 5 }
     6: { 1 4 5 }
     7: { 2 3 4 }
     8: { 2 3 5 }
     9: { 2 4 5 }
     10: { 3 4 5 }
Enter n r?

sflam(Raymond) 的大數版
--------------------------------------

[code: C++11]
#include <vector>
#include <string>
#include <iostream>
#include <iomanip>
#include <boost/multiprecision/cpp_int.hpp>

using namespace boost::multiprecision;
namespace combination
{
  cpp_int choose(const size_t n, const size_t r)
  {
    if (n < r)
      return 0;
    else if (n == r)
      return 1;

    // choose(100, 97) ==> choose(100, 3)
    size_t k2 = (r < n - r) ? r : n - r;

    // n * (n-1) * (n-2) ... (n-k2+1)
    // -------------------------------
    // 1 * 2 * ... k2
    cpp_int ans = n;
    for (size_t i = 2; i <= k2; ++i)
    {
      ans = (ans * (n - i + 1))/i;
    }
    return ans;
  }

  size_t LargestV(size_t a, size_t b, cpp_int x)
  {
    size_t v = a - 1;
    while (choose(v, b) > x)
      --v;
    return v;
  }

  void GetMthElement(const size_t n, const size_t r, cpp_int m,
                     std::vector<size_t>& result)
  {
    result.resize(r);
    size_t a = n;
    size_t b = r;
    cpp_int x = (choose(n, r) - 1) - m;
    for (size_t i = 0; i < r; ++i)
    {
      result[i] = LargestV(a, b, x);
      x = x - choose(result[i], b);
      a = result[i];
      b = b - 1;
    }
    for (size_t i = 0; i < r; ++i)
    {
      result[i] = n - result[i];
    }
  }

  void C(const size_t n, const size_t r)
  {
    std::vector<size_t> result(r);
    const cpp_int mMax = choose(n, r);
    for (cpp_int m = 0; m < mMax; ++m)
    {
      GetMthElement(n, r, m, result);
      std::cout << std::setw(10) << m+1 << ": " << result << std::endl;
    }
  }
}
[/code: C++11]

我的版本
--------------------------------------

#include <vector>
#include <string>
#include <iostream>
#include <iomanip>

template <typename T> 
std::ostream& operator <<(std::ostream& os, std::vector<T>& vt)
{
  os << "{ ";
  for (std::vector<T>::iterator cit = vt.begin(); 
      cit != vt.end();
      ++cit)
  {
    os << *cit << ' ';
  }
  os << "}";
  return os;
}

unsigned long long choose(const size_t n, const size_t r)
{
  if (n < r)
    return 0;
  else if (n == r)
    return 1;

  // choose(100, 97) ==> choose(100, 3)
  size_t k2 = (r < n - r) ? r : n - r;

  // n * (n-1) * (n-2) ... (n-k2+1)
  // -------------------------------
  // 1 * 2 * ... k2
  unsigned long long ans = n;
  for (size_t i = 2; i <= k2; ++i)
  {
    ans = (ans * (n - i + 1))/i;
  }
  return ans;
}

size_t V(size_t a, size_t b, long long &m)
{
  if(b == 1)
    return (size_t)m+1;

  size_t l = 0;
  while(true)
  {
    long long z = choose(a-l-1,b-1);
    if(z < (m+1))
    {
      m = m - z;
      ++l;
    }
    else 
      return l+1;
  }
}

void GetMthElement(const size_t n, const size_t r, unsigned long long mm, 
                   std::vector<size_t>& result)
{
  long long m = mm;
  result.resize(r);
  size_t a = n;
  size_t b = r;
  size_t c = 0;
  for(size_t i = 0; i < r; ++i)
  {
    result[i] = V(a, b, m) + c;
    c = result[i];
    a = n - c;
    --b;
  }
}

void C(const size_t n, const size_t r)
{
  std::vector<size_t> result(r);
  const unsigned long long mMax = choose(n, r);
  for (unsigned long long m = 0; m < mMax; ++m)
  {
    GetMthElement(n, r, m, result);
    std::cout << std::setw(10) << m+1 << ": " << result << std::endl;
  }
}

int main(int argc, char* argv[])
{
  size_t n, r;
  std::cout << "Calculate nCr" << std::endl;
  while (std::cout << "Enter n r? " << std::flush, std::cin >> n >> r)
  {
    C(n, r);
  }

  return 0;
}

排列(Permutation)

Permutation 的 C++ 程式碼, 配合之前的 combination, 就可以打出 nPr 所有
組合的排列. NextPermutation() 的邏輯參
考 http://en.wikipedia.org/wiki/Permutation#Generation_in_lexicographic_order.
[code: C++]
#include <algorithm>
#include <iterator>
namespace permutation
{
  bool NextPermutation(std::vector<size_t>& result)
  {
    if (result.size() < 2)
      return false;

    size_t k = result.size() - 2;
    while (k != 0 && result[k] >= result[k+1])
    {
      --k;
    }
    if (k == 0 && result[k] >= result[k+1])
      return false;

    size_t l = result.size() - 1;
    while (l > k && result[k] >= result[l])
    {
      --l;
    }
    std::swap(result[k], result[l]);
    std::reverse(result.begin() + k + 1, result.end());
    return true;
  }

  void P(const size_t n, const size_t r)
  {
    std::vector<size_t> result(r);
    const unsigned long long mMax = combination::choose(n, r);
    for (unsigned long long line = 1, m = 0; m < mMax; ++m)
    {
      combination::GetMthElement(n, r, m, result);
      do 
      {
        std::cout << std::setw(10) << line++ << ": " << 
            result << std::endl;
      } while (NextPermutation(result));
    }
  }
}

int main()
{
  size_t n, r;
  std::cout << "Calculate nPr" << std::endl;
  while (std::cout << "Enter n r? " << std::flush, std::cin >> n >> r)
  {
    permutation::P(n, r);
  }
}
[/code: C++]

[執行:]
Enter n r? 4 3
     1: { 1 2 3 }
     2: { 1 3 2 }
     3: { 2 1 3 }
     4: { 2 3 1 }
     5: { 3 1 2 }
     6: { 3 2 1 }
     7: { 1 2 4 }
     8: { 1 4 2 }
     9: { 2 1 4 }
     10: { 2 4 1 }
     11: { 4 1 2 }
     12: { 4 2 1 }
     13: { 1 3 4 }
     14: { 1 4 3 }
     15: { 3 1 4 }
     16: { 3 4 1 }
     17: { 4 1 3 }
     18: { 4 3 1 }
     19: { 2 3 4 }
     20: { 2 4 3 }
     21: { 3 2 4 }
     22: { 3 4 2 }
     23: { 4 2 3 }
     24: { 4 3 2 }
Enter n r?

量子演算法：

這裡的思路是以 thread 代替遞廻呼叫，若 thread 的數量可以沒上限，而且 thread 的處理又得到硬體的支援速度極快(比如量子電腦)，那麼一瞬間就可以得到所有答案。

因沒有量子電腦所以只能像徵性的意思一下，先由遞廻版本了解這個演算法的實作方法。

遞廻版：

#include <iostream>

#include <iomanip>

struct resul_t
{
  resul_t(resul_t &previous_arg, size_t value_arg)
    :previous(previous_arg)
  {
    value = value_arg;
  }
  resul_t &previous;
  size_t value;
};

class Combination
{
  static void(*m_pFn)(const resul_t &);

  static void v(resul_t &pre, const size_t n, const size_t r)
  {
    if (r == 1)
    {
      for (size_t i = n; i > 0; --i)
        m_pFn(resul_t(pre, i - 1));
      return;
    }

    for (size_t i = n - 1; i >= r - 1; --i)
      v(resul_t(pre, i), i, r - 1);
  }
public:
  static void C(const size_t n, const size_t r, void(*pFn)(const resul_t &))
  {
    m_pFn = pFn;
    v(*(resul_t*)NULL, n, r);
  }
};

void(*Combination::m_pFn)(const resul_t &);

unsigned long long m = 1;

void Fn(const resul_t &resul_arg)
{
  std::cout << std::setw(10) << m++ << '=' << '{';
  for (const resul_t *resul = &resul_arg; resul != (resul_t*)NULL; 
    resul = &resul->previous
  )
    std::cout << resul->value << ' ';
  std::cout << '}' << std::endl;
}

int main(int argc, char* argv[])
{
  size_t n, r;
  std::cout << "Calculate nCr" << std::endl;
  while (std::cout << "Enter n r? " << std::flush, std::cin >> n >> r)
  {
    m = 1;
    Combination::C(n, r, Fn);
  }

  return 0;
}

接著用 std::thread 將遞廻版改寫。

std::thread 版：


#include <mutex>
#include <iostream>
#include <iomanip>
#include <memory>
#include <thread>

using namespace std;
struct result_t
{
  result_t(const shared_ptr<result_t> &previous_arg, const size_t value_arg)
    :previous_ptr(previous_arg), value(value_arg)
  {
  }
  const shared_ptr<result_t> previous_ptr;
  const size_t value;
};

class Combination
{
  static void(*m_pFn)(const shared_ptr<result_t> &);
  static void v(const shared_ptr<result_t> &pre,
    const size_t n, const size_t r)
  {
    if (r == 1)
    {
      for (size_t i = n; i > 0; --i)
      {
        thread m(m_pFn, shared_ptr<result_t>(new result_t(pre, i - 1)));
        m.detach();
      }
      return;
    }
    for (size_t i = n - 1; i >= r - 1; --i)
    {
      thread m(v, shared_ptr<result_t>(new result_t(pre, i)), i, r - 1);
      m.detach();
    }
  }
public:
  static void C(const size_t n, const size_t r, 
    void(*pFn)(const shared_ptr<result_t> &))
  {
    m_pFn = pFn;
    v(shared_ptr<result_t>((result_t*)NULL), n, r);
  }
};

void(*Combination::m_pFn)(const shared_ptr<result_t> &);

unsigned long long m = 1;
void Fn(const shared_ptr<result_t> &result_arg)
{
  static mutex io_mutex;
  lock_guard<mutex> lk(io_mutex);
  cout << setw(10) << m++ << '=' << '{';
  for (shared_ptr<result_t> result_ptr = result_arg; 
   result_ptr.get() != (result_t*)NULL; 
   result_ptr = result_ptr->previous_ptr
  )
    cout << result_ptr->value << ' ';

cout << '}' << endl;
}

int main(int argc, char* argv[])
{
  size_t n, r;
  cout << "Calculate nCr" << endl;
  cout << "Enter n r?" << endl;
  while (cin >> n >> r)
  {
    m = 1;
    Combination::C(n, r, Fn);
  }
  return 0;
}

因 std::thread 可同時產生的數量是有上限的，所以若超過上限，程式會被強制結束，因此面對這種大量吃 thread 的東西，應該先經過一層控管機制來產生 thread 才能確保安全。

所以再用 CxxlMan2 程式庫裡的 ThreadLimit 改寫一次

ThreadLimit 版：


#include <THREADMGR.HPP> // CxxlMan2::ThreadLimit

#include <iostream>
#include <iomanip>


using namespace std;
using namespace CxxlMan2;

struct result_t
{
  result_t(const shared_ptr<result_t> &previous_arg,
    const size_t value_arg
  )
    :previous_ptr(previous_arg), value(value_arg)
  {
  }
  const shared_ptr<result_t> previous_ptr;
  const size_t value;
};

class Combination
{
  static void(*m_pFn)(const shared_ptr<result_t> &);
  static ThreadLimit TL; 

  static void v(const shared_ptr<result_t> &pre, const size_t n,
       const size_t r)
  {
    if (r == 1)
    {
      for (size_t i = n; i > 0; --i)
        TL.Thread(m_pFn, shared_ptr<result_t>(new result_t(pre, i - 1)));
      return;
    }

    for (size_t i = n - 1; i >= r - 1; --i)
      TL.Thread(v, shared_ptr<result_t>(new result_t(pre, i)), i, r - 1);
  }
public:
  static void C(const size_t n, const size_t r, 
    void(*pFn)(const shared_ptr<result_t> &))
  {
    m_pFn = pFn;
    v(shared_ptr<result_t>((result_t*)NULL), n, r);
  }
};

void(*Combination::m_pFn)(const shared_ptr<result_t> &);

// 視自己電腦的 CPU 做修改，這裡假設是 4 核心，
// 設定太大不會比較快
ThreadLimit Combination::TL(4);

unsigned long long m = 1;

void Fn(const shared_ptr<result_t> &result_arg)
{
  static mutex io_mutex;
  lock_guard<mutex> lk(io_mutex);

  cout << setw(10) << m++ << '=' << '{';
  for (shared_ptr<result_t> result_ptr = result_arg;
   result_ptr.get() != (result_t*)NULL;
   result_ptr = result_ptr->previous_ptr
  )
    cout << result_ptr->value << ' ';
  cout << '}' << endl;
}



int main()
{
  size_t n, r;
  cout << "Calculate nCr" << endl;
  cout << "Enter n r?" << endl;
  while (cin >> n >> r)
  {
    m = 1;
    Combination::C(n, r, Fn);
  }

  
  return 0;
}

CxxlMan

2015年3月28日星期六

排列組合產生器

模擬遞迴法

排列(Permutation)

組合(Combination)

計算法

組合(Combination)

排列(Permutation)

量子演算法：

遞廻版：

std::thread 版：

ThreadLimit 版：

沒有留言:

張貼留言

2015年3月28日 星期六

排列組合產生器

模擬遞迴法

排列(Permutation)

組合(Combination)

計算法

組合(Combination)

排列(Permutation)

量子演算法：

遞廻版：

std::thread 版：

ThreadLimit 版：

沒有留言:

張貼留言

2015年3月28日星期六