12.31.2011

Generating integer partitions with C++

C++: 自然数の分割を列挙する

正の整数 n の分割(integer partition)とは、与えられた正整数 n を正整数の和として表す方法をいう。

例えば、「6」 はこのように11通りの異なる方法で分割できる。

　6　　5+1　　4+2　　4+1+1　　3+3　　3+2+1　　3+1+1+1　　2+2+2　　2+2+1+1　　2+1+1+1+1　　1+1+1+1+1+1

この成分(part)の組み合わせを std::vector<int> として順次生成するクラスを作りたい。
以下、この関数を P(n) と呼ぶ。

1. 再帰を利用したアルゴリズム(Python)

Python のジェネレータ機能を使ったコードが、こちらに投稿されている。
http://code.activestate.com/recipes/218332-generator-for-integer-partitions/

処理時間計測用のコードと併せて記載する。

def P(n):
  # base case of recursion: zero is the sum of the empty list
  if n == 0:
    yield []
    return

  # modify partitions of n-1 to form partitions of n
  for p in P(n - 1):
    p.append(1)
    yield p
    p.pop()
    if p and (len(p) < 2 or p[-2] > p[-1]):
      p[-1] += 1
      yield p

def main():
  import timeit
  for n in (15, 30, 45, 60, 75, 90):
    code = 'for _ in P(%d): pass' % n
    print 'n=%d  elapsed=%f' % (
        cnt,
        timeit.Timer(code, 'from __main__ import P').timeit(1))

if __name__ == '__main__':
  main()

手元のPC（Intel Core2Duo E4600 @ 2.40GHz）での実行結果は以下のとおり。（処理時間の単位は秒）

n=15  elapsed=0.000423
n=30  elapsed=0.015092
n=45  elapsed=0.290839
n=60  elapsed=3.473330
n=75  elapsed=31.071480
n=90  elapsed=236.823416

このコードのアルゴリズムは以下のようになる。

①P(n) の計算にあたり、P(n-1) を列挙する。　※P(0) は空のリストを返す。

②P(n-1) の各出力に対して、末尾に新規の成分「1」を追加したリストを出力する。

③P(n-1) の各出力を k個の成分を持つリスト x1,x2, … ,xk とした場合、
　成分の個数が1個(k=1)であれば、その唯一の成分を インクリメントする。（x1 ← x1 + 1）
　末尾の2成分(xk-1, xk)について、xk-1　> xk が成立する場合、末尾の成分をインクリメントする。（xk ← xk + 1）
　インクリメントしたリストを出力する。
　（言い換えれば、末尾の成分をインクリメントしても要素の降順が守られる場合にのみインクリメント実行）

この結果出力される各成分は降順にソートされている。また、出力結果は辞書順となる。
P(0)～P(5) の処理の様子はこのようになる。

image
緑色の箇所はP(n-1)に対して新規の成分「1」を追加したこと、ピンク色の箇所は既存の成分をインクリメントしたことを表している。

2. 再帰を利用しないアルゴリズム(C++)

以下の文献を参考に、C++で処理を実装してみる。

Antoine Zoghbiu and Ivan Stojmenovic. Fast Algorithms for Generating Integer Partitions. Intern. J. Computer Math., 70, 319-332 (1998) (PDF)
http://www.site.uottawa.ca/~ivan/F49-int-part.pdf

出力は辞書の逆順とするほうがコストが小さい。そのアルゴリズムはこのような感じだ。
尚、この場合も1つの出力における各成分は降順にソートされるようにする。

①n = 0 の場合は空のリストを出力して終了。

②成分が1個だけのリスト {n} を作成し、それを出力する。

③2以上の成分の数を h とおき、1 を代入する。（h ← 1）

＃　以下、④～⑦を繰り返す

④最初の成分が1 （x1 = 1） ならば終了。

⑤2以上の成分のうち最も右側にあるもの（xh）をデクリメントし、その値を r とする。（r ← xh-1; xh ← r）

⑥xh より右側にあった成分を全て結合し、{r, r, … , d}  (0 ＜ d ≦ r) という形に再構成する。
　xh より右側の成分は、定義より全て 1 であるから、現在のリストの個数を|{x}|とすると
　|{x}| h + 1 = c × r + d  (0 < d ≦ r)
　という変換を行えばよいことが分かる。

⑦r = 1 の場合は h をデクリメントする。（h ← h - 1）
　r > 1 かつ d = 1 の場合は h に c を加算。（h ← h + c）
　r > 1 かつ d ≧ 1 の場合は h に c+1 を加算する。（h ← h + c + 1）

＃　④へ戻る

P(6) の処理の様子を以下に示す。

  ・リスト {n} を作成
・h ← 1
  ・xh = x1 をデクリメント
・1 を 5 単位で再構成：　1 = 0 × 5 + 1
・2 以上の成分の個数(h)は変わらず
  ・xh = x1 をデクリメント
・2 を 4 単位で再構成：　2 = 0 × 4 + 2
・2 以上の成分ができたので h を増加
 

・xh = x2 をデクリメント
・1 を 1以下に再構成
・2 以上の成分が減ったので h を減少

 

・xh = x1 をデクリメント
・3 を 3 単位で再構成：　3 = 1 × 3 + 0
・2 以上の成分ができたので h を増加

 

・xh = x2 をデクリメント
・1 を 2 単位で再構成：　1 = 0 × 2 + 1

 

・xh = x2 をデクリメント
・2 を 1以下に再構成
・2 以上の成分が減ったので h を減少

 

・xh = x1 をデクリメント
・4 を 2 単位で再構成：　4 = 2 × 2 + 0
・2 以上の成分が2つできたので h を2つ増加

 

・xh = x3 をデクリメント
・1 を 1以下に再構成
・2 以上の成分が減ったので h を減少

 

・xh = x2 をデクリメント
・3 を 1以下に再構成
・2 以上の成分が減ったので h を減少

 

・xh = x1 をデクリメント
・5 を 1以下に再構成
・2 以上の成分が減ったので h を減少

・これで、最初の成分(x1)が 1 になったので終了

黄色の箇所は 2以上の数値（＝h としてカウントされるもの）を、青色の箇所はデクリメントした結果1となったこと（＝hのデクリメントを伴ったこと）を表している。

コーディングにあたっては、簡単のため h を 1-indexed から 0-indexed へ変換する。
xh = 2 のときの再配分処理はリストの末尾に 1 を追加するだけでよい。
また、出力したリストの個数を count 変数として保持することにした。

今回もやはり時間計測用のコードと一緒に掲載する。

#include <vector>
#include <ctime>
#include <cstdio>

// Generator for integer partition.
class Partition {
 public:
  int count;

  Partition(int n) {
    x_ = std::vector<int>(std::min(1, n), n);
    h_ = n > 1 ? 0 : -1;
    count = 0;
  }

  std::vector<int> &next() {
    if (x_.empty()) return x_;
    if (count > 0) {
      if (h_ < 0) {
        x_.clear();
        return x_;
      }
      if (x_[h_] == 2) {
        --x_[h_--];
        x_.push_back(1);
      }
      else {
        int r = --x_[h_];
        int t = x_.size() - h_;
        while (t >= r) {
          x_[++h_] = r;
          t -= r;
        }
        if (t == 0) {
          x_.resize(h_ + 1);
        }
        else {
          x_.resize(h_ + 2, 1);
          if (t > 1) x_[++h_] = t;
        }
      }
    }
    ++count;
    return x_;
  }

 private:
  std::vector<int> x_;
  int h_;
};

int main() {
  // Print each P(6).
  Partition p(6);
  std::vector<int> v;
  while ((v = p.next()).size()) {
    printf("%2d: [", p.count);
    for (int i = 0; i < (int)v.size(); ++i) {
      printf("%s%d", i == 0 ? "" : ", ", v[i]);
    }
    printf("]\n");
  }

  // Measuring time elapsed.
  int n[] = {15, 30, 45, 60, 75, 90};
  for (int i = 0; i < static_cast<int>(sizeof(n) / sizeof(n[0])); ++i) {
    clock_t t = clock();
    Partition p(n[i]);
    while (p.next().size()) {}
    clock_t s = clock();

    printf("n=%d  elapsed=%.6f  count=%d\n", n[i],
        static_cast<double>(s - t) / CLOCKS_PER_SEC, p.count);
  }
  return 0;
}

実行結果は以下のとおり。（処理時間の単位は秒）

 1: [6]
 2: [5, 1]
 3: [4, 2]
 4: [4, 1, 1]
 5: [3, 3]
 6: [3, 2, 1]
 7: [3, 1, 1, 1]
 8: [2, 2, 2]
 9: [2, 2, 1, 1]
10: [2, 1, 1, 1, 1]
11: [1, 1, 1, 1, 1, 1]
n=15  elapsed=0.000000  count=176
n=30  elapsed=0.000000  count=5604
n=45  elapsed=0.010000  count=89134
n=60  elapsed=0.105000  count=966467
n=75  elapsed=0.818000  count=8118264
n=90  elapsed=5.612000  count=56634173

参考：
http://handasse.blogspot.com/2010/01/blog-post.html

Wikipedia
http://en.wikipedia.org/wiki/Partition_%28number_theory%29
http://ja.wikipedia.org/wiki/%E8%87%AA%E7%84%B6%E6%95%B0%E3%81%AE%E5%88%86%E5%89%B2

投稿者 mog project 時刻: 18:41
ラベル: ,

0 件のコメント:

コメントを投稿

登録: コメントの投稿 (Atom)