e1071 に permutations という関数がある。
実装はいろいろあるけど,簡単そうな奥村先生に典拠するアルゴリズムを使って書いてみた。
ただし,n! なんて,すぐにメモリ制限を超えてしまうので,実際の利用はたかだか n <= 10 しかないので,どんなに速いアルゴリズムでも,生涯期待計算時間(ともいうようなもの)を見ると,ゼロに近いだろう。(いろいろなプログラムの実装法を経験するという意味はあるだろうとやってみているだけなのだが...それにしても,本当に利用すべきと言う場合は未だ見つからず)。
src <- '
int N = as(n);
if (N < 1 || N > 10) {
return wrap(NA_REAL);
}
int i, k, t, c[N + 1], p[N];
int size = 1;
for (i = 2; i <= N; i++) {
size *= i;
}
Rcpp::IntegerMatrix res(size, N);
int count = 0;
for (i = 0; i < N; i++) {
p[i] = i + 1;
}
for (i = 1; i <= N; i++) {
c[i] = i;
}
k = 1;
while (k < N) {
if (k & 1) {
i = c[k];
}
else {
i = 0;
}
t = p[k];
p[k] = p[i];
p[i] = t;
for (i = 0; i < N; i++) {
res(count, i) = p[i];
}
count++;
k = 1;
while (c[k] == 0) {
c[k] = k;
k++;
}
c[k]--;
}
return wrap(res);
'
library(inline)
system.time(genperm <- cxxfunction(signature(n="integer"), src, plugin="Rcpp"))
system.time(genperm(10))
library(e1071)
system.time(permutations(10))
実行例
> system.time(genperm <- cxxfunction(signature(n="integer"), src, plugin="Rcpp"))
ユーザ システム 経過
3.146 0.708 8.395 # コンパイルするのに必要な時間
> system.time(genperm(10)) # 10! のリストを計算する時間
ユーザ システム 経過
0.253 0.123 0.373
> library(e1071)
要求されたパッケージ class をロード中です
> system.time(permutations(10)) # perumutations による 10! のリストを計算する時間
ユーザ システム 経過
1.537 0.589 2.109
ほとんど無駄な努力と言える。
※コメント投稿者のブログIDはブログ作成者のみに通知されます