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

12.30.2011

How to get current date in Perl

Perl: 現在日付の取得

localtime 関数の戻り値をスライシングすれば、日付だけを取得することができる。
例えば、今日の午前 0:00 を Unix Time で表現する方法は以下のとおり。

use POSIX qw( mktime strftime );
my $time = mktime(0, 0, 0, (localtime)[3..5]);
print strftime('%Y-%m-%d %H:%M:%S', localtime($time));

参考:
http://ebitem.net/program-blog/archives/96

12.27.2011

Etymology of the 'glob'

グロブの語源について

ファイルパスのワイルドカードなどを展開することをグロブとか、グロビングとか言うが
これはどうやらUNIX黎明期の「glob」コマンドが由来らしい。

global を略した言葉のようだ。

http://en.wikipedia.org/wiki/Glob_%28programming%29
http://cm.bell-labs.com/cm/cs/who/dmr/man71.pdf

12.23.2011

Expanding all UNIX environment variables in Perl

Perl: UNIX環境変数を全て展開

正規表現とENVハッシュの組み合わせで実現できる。

sub expand {
    my ($text) = @_;
    $text =~ s{ \$ {? (\w+) }? }
              { exists $ENV{$1} ? $ENV{$1} : $& }gexms;
    return $text;
}

参考:
http://stackoverflow.com/questions/7697216/perl-one-liner-to-expand-unix-environment-variables

12.21.2011

Find out the absolute path from a relative path in Bourne Shell

Bシェル: 相対パスから絶対パスを得る方法

dirname コマンドで得られるディレクトリに cd し、pwd コマンドを行うことで絶対パスを得ることができる。
このときサブシェル内で実行することで、実際のシェルプロセスには影響を与えないようにする。

・例1 自分自身の居場所を特定する

ABSPATH=`(cd \`/usr/bin/dirname $0\` && pwd)`/`/usr/bin/basename $0`

・例2 自分自身の親ディレクトリを特定する

PARENTDIR=`(cd \`/usr/bin/dirname $0\`/.. && pwd)`

プロジェクトごとディレクトリを移動する場合などに便利。

参考:
http://dokonoumanohone.blog47.fc2.com/blog-entry-2.html
http://blog.hansode.org/archives/51481467.html

12.15.2011

XML::TreePP – 100% pure Perl XML parser module

Perl: Pure Perl 実装のXML解析モジュール

Perl モジュールでXML解析を行う場合、XML::Simple では XML::Parsar を要するため、
環境ごとにコンパイルが必要となってしまう。

Pure Perl 実装の XML::TreePP を使えば、.pm ファイルを置くだけで利用することが可能だ。

参考:
http://www.kawa.net/works/perl/treepp/treepp.html
http://search.cpan.org/~kawasaki/XML-TreePP-0.41/lib/XML/TreePP.pm

12.12.2011

Oracle: How to find foreign key constraints

Oracle:外部参照されているテーブルの確認方法

テーブルのトランケートを試み、
エラー「ORA-02266 表には使用可能な外部キーによって参照される一意キー/主キーが含まれています。」
が発生したときに、そのテーブルがどのテーブルから参照されているかを確認するSQL。

トランケートを行うには、外部キー制約を DISABLED にする必要がある。

SELECT
    constraint_name,
    table_name,
    status
FROM
    user_constraints
WHERE
    r_constraint_name IN
        (
        SELECT
            constraint_name
        FROM
            user_constraints
        WHERE
            table_name = ’テーブル名‘
        ); 

参考:
http://yp.xenophy.com/?p=74
http://www.databasejournal.com/features/oracle/article.php/3665591/Finding-Foreign-Key-Constraints-in-Oracle.htm

12.07.2011

Unix time conversion in AWK (with Fliegel-Van Flandern algorithm)

AWK: フリーゲルの公式に基づく Unix time 変換処理

gawk が使える場合を除き、AWKスクリプトにおけるグレゴリオ暦とUnix time(POSIX time)の相互変換処理は自前で準備する必要がある。
そこで今回はフリーゲルの公式を使って実装してみる。

・unix_time.awk ※テスト用コードも実装済み

#!/usr/bin/nawk -f

# Based on Fliegel-Van Flandern algorithm.
# Inputs must be Epoch (1970/01/01 00:00:00) and over.
# Time zone will not be considered.
function to_unix_time(yyyy, mm, dd, hh, mi, ss) {
  if (mm - 2 <= 0) { --yyyy; mm += 12}
  days = int(365.25 * yyyy) + int(yyyy / 400) - int(yyyy / 100) \
      + int(30.59 * (mm - 2)) + dd - 719499  # 1721089 - 2440588(Epoch)
  return ((days * 24 + hh) * 60 + mi) * 60 + ss
}

function to_string(unix_time) {
  ss = unix_time % 60; unix_time = (unix_time - ss) / 60
  mi = unix_time % 60; unix_time = (unix_time - mi) / 60
  hh = unix_time % 24; unix_time = (unix_time - hh) / 24
  j = unix_time + 2472632  # 2440588(Epoch) + 32044
  g = int(j / 146097); dg = j % 146097
  c = int((int(dg / 36524) + 1) * 3 / 4); dc = dg - c * 36524
  b = int(dc / 1461); db = dc % 1461
  a = int((int(db / 365) + 1) * 3 / 4); da = db - a * 365
  y = g * 400 + c * 100 + b * 4 + a
  m = int((da * 5 + 308) / 153) - 2
  d = da - int((m + 4) * 153 / 5) + 122
  yyyy = y - 4800 + int((m + 2) / 12)
  mm = (m + 2) % 12 + 1
  dd = d + 1
  return sprintf("%04d-%02d-%02d %02d:%02d:%02d", yyyy, mm, dd, hh, mi, ss)
}

# test code
{
  x = to_unix_time(substr($1, 1, 4), substr($1, 6, 2), substr($1, 9, 2))
  expected = (NR - 1) * 60 * 60 * 24
  if (x != expected) {
    print "Assertion failed. line " NR " - Expected: " expected ", Received: " x
    exit
  }
  print to_string(x)
}

Python の datetime モジュールの処理結果と比較することで、1970/1/1~9999/12/31 の期間の全日付について連続性および妥当性を検証する。

・unix_time_test.py

#!/usr/bin/env python

import datetime

s = datetime.date(1970, 1, 1)
t = datetime.date(9999, 12, 31)

while True:
  print '%s 00:00:00' % s
  if s >= t: break
  s += datetime.timedelta(days=1)

datetime.date では西暦10000年は表現できないので注意。(Exception が発生)
テスト用のシェルを作るならこのような感じである。

・unix_time_test.sh

#!/usr/bin/ksh

./unix_time_test.py > ./test_in.txt
./unix_time.awk ./test_in.txt > ./test_out.txt
diff ./test_in.txt ./test_out.txt > /dev/null
if [ $? -eq 0 ]; then print 'OK'; else print 'NG'; fi

改行コードがLFの場合、58657940 bytes(≒56MB) のファイルが2個生成される。(test_in.txt, test_out.txt)
Windows 環境なら diff のところを comp で検査。

・出力結果

OK

ただ実際のところ、もう少しわかりやすいコーディングのほうが周囲の理解を得やすいかもしれない。
Perl や Python、Ruby が使える環境であれば尚更である。

参考:
http://www5d.biglobe.ne.jp/~noocyte/Programming/GregorianAndJulianCalendars.html#DateToDayNumber
http://en.wikipedia.org/wiki/Julian_day

12.05.2011

How to show a run-time error message in Eclipse CDT 7.0.2

Eclipse CDT における実行時エラーメッセージの出力

Linux(UNIXも?)環境におけるEclipse Helios CDT 7.0.2のバグ。
例えば以下のようなプログラムを実行(Run)するとSegmentation faultで終了する。

int *x; x[0] = 0;

ところが、Eclipseのコンソール上には何もメッセージが表示されない。

Run configuration を以下のように設定すると実行時のstderrメッセージが表示されるようになる。

Main タブ – Application: /bin/sh
Arguments タブ – Program arguments: –c "Debug/program_name"

Eclipse 自体のバージョンアップによって解消するかは未確認。

参考:
http://stackoverflow.com/questions/3760734/how-to-let-eclipse-cdt-show-runtime-error-e-g-segmentation-fault
http://stackoverflow.com/questions/8099445/eclipse-cdt-running-c-program-show-nothing-in-the-console-why

12.03.2011

Failure probability calculator with Python

Python: 部品故障予測に関する計算

1. 目的

ある期間においてpの確率で故障する部品がn個稼働している。
このとき、期間内の故障発生回数が累計m回よりも多い確率を求めたい。

2. 前提条件

・故障が発生したら、その部品はすぐに新しい部品と交換される。(交換に要する時間は考慮しない)
・交換した部品が故障することもある。その場合も故障回数が加算される。
・故障の発生はポアソン分布とする。

3. 考察

1個の部品の故障率をλとする場合、この部品が期間内にx回故障する確率 P(x) は
ポアソン分布の定義から以下の式で表される。

image … 式(1)

 ※e は数学定数の一つであるネイピア数(Napier's constant)を表す。e = 2.71828…
 ※n! は n の階乗を表す。n! = n × (n-1) × ・・・ × 3 × 2 × 1

部品数がn個になれば、故障の確率はn倍となる。
したがって、故障率pのn個の部品の故障回数の累計がちょうどx回になる確率は以下ように求められる。

image … 式(2)

さて、ここで累計m回よりも多く故障が発生する確率(Q)は

 1-(故障発生回数がm回以下である確率)

と計算できる。つまり、

 故障が一度も発生しない確率(=故障発生回数が0)
 1回だけ故障が発生する確率
 累計2回故障が発生する確率
 累計3回故障が発生する確率
  ・・・
 累計m回故障が発生する確率

の合計を1から差し引いた値が求める確率である。

式で表すとこのようになる。

image … 式(3)

上記の計算を素直に実装すると、階乗の計算が m+1 回発生することとなる。
階乗自体の計算量はO(m)なので、全体の計算量はO(m2)。

mの数値が小さければ問題ないが、10000 程度なら計算不能となってしまうだろう。

 

式(1)をよく見ると、以下の漸化式が成り立つことがわかる。

image … 式(4)

これを利用してDP(Dynamic Programming/動的計画法)を行えば、計算量はO(m)となる。

さらに、このDPは2つの変数を交互に使用することで再利用可能である。
そうすることでメモリ使用量を大幅に節約できる。
これを擬似コードで書くと以下のようになる。

image

4. コード

・関数定義

import math

def FailureProbability(rate, num, times):
  rate *= num
  dp = [math.exp(-rate), 0]
  ret = 1 - dp[0]
  for i in range(1, times + 1):
    dp[i & 1] = dp[i & 1 ^ 1] * rate / i
    ret -= dp[i & 1]
  return ret

・実行例

print FailureProbability(0.03, 18, 0)
print FailureProbability(0.03, 18, 1)
print FailureProbability(0.03, 18, 2)
print FailureProbability(0.03, 18, 3)

 (実行結果)

0.417251747626
0.102567691344
0.0176029961479
0.00230935101263

この例では、故障率3%の部品が18個稼動している場合を想定している。

この結果から、
 41.7%の確率で少なくとも1個の部品が期間内に故障してしまう
 予備部品を2個用意したとしても、1.76%の確率で予備部品不足となる
  (3回以上の故障が発生する)
というようなことが見て取れる。