範囲for文でcombinationの列挙

これは何?

n個の要素からk個の要素を選び出す組み合わせは{}_n\mathrm{C}_r=\frac{n!}{(n-k)!k!}通りあります。nkを与えられたとき、そのような組み合わせをすべて列挙したい時があるので、C++で実装をしました。

Combinationの表現には、整数のbitを使います。下位n bitを要素に対応させて、選ばれた要素に対応するbitに1を立てることで組み合わせを表現します。例えば、n=5,k=2のときの列挙は次のようになってほしいです。

0b00011
0b00101
0b00110
0b01001
0b01010
0b01100
0b10001
0b10010
0b10100
0b11000

この記事ではさらに、列挙されたcombinationをコンテナとみなし、範囲for文を使って簡潔に列挙を記述できるような構造体を書きました。つまり、次のようなコードが書けます。

for (auto i : combinations(n, k)) {
    // iを使った処理を書く
}

範囲for文の仕様

cpprefjpによれば、C++17以降の範囲for文は、次のように通常のfor文に展開されます。

for ( for-range-declaration : for-range-initializer ) statement

は以下のように展開される:

{
  auto && __range = for-range-initializer;
  auto __begin = begin-expr;
  auto __end = end-expr;     // __begin と __end は異なる型でもよい
  for ( ; __begin != __end; ++__begin ) {
    for-range-declaration = *__begin;
    statement
  }
}

細かい所は置いておいて、どうやらコンテナのイテレータbegin()で初期化して、end()と一致するまでインクリメントで通常のfor文を回しているようです。範囲for文で使いたいだけであれば、都合のいいイテレータを作ってしまえば、実際にコンテナが存在しなくてもよさそうです。

next_combination

辞書順で次に来るcombinationを返す関数です。整数のbitで表現した場合、combinationの辞書順は数値の昇順と同値です。
ネットで調べると、次のような記事が見つかります。
qiita.com
nemutage.hatenablog.jp
一つ目の記事では、次のような実装が紹介されています。

int next_combination(int sub) {
    int x = sub & -sub, y = sub + x;
    return (((sub & ~y) / x) >> 1) | y;
}

二つ目の記事では、一つ目の記事の実装が分かりやすく解説されています。私が書いたnext_combinationも基本的には同じ仕組みで動くので、詳しい解説はこちらを読んでください。

上の実装では、取得したbitをどうにか使いまわして並べ替えるために、除算が使われています。除算は速度が遅いので、除算を使わない実装を考えました。立っているbitの数を数えるpopcountを使うことで、天下り的な実装で除算を回避できます。

int next_combination(int sub) {
    int k = __builtin_popcount(sub);
    sub += sub & -sub;
    return sub + (1 << (k - __builtin_popcount(sub))) - 1;
}

ここでのkはわざわざ求めなくても定数なので、前計算をしておけば演算の数をさらに減らすことができます。

int next_combination(int sub, int k_bits) {
    sub += sub & -sub;
    return sub + (k_bits >> __builtin_popcount(sub));
}

ここで、k_bits = (1 << k) - 1です。

構造体にする

#include <assert.h>

struct combinations {
    using u64 = std::uint64_t;

    uint n;
    uint k;

    combinations(uint n_, uint k_) : n(n_), k(k_) { assert(n_ < 64 && n_ >= k_ && k_ > 0); }

    struct itr {
        u64 s;
        u64 t;

        bool operator!=(itr &right) { return s != right.s; }
        void operator++() {
            s += (s & -s);
            s += t >> __builtin_popcountll(s);
        }
        u64 operator*() { return s; }
    };

    itr begin() { return {(1ull << k) - 1, (1ull << k) - 1}; }
    itr end() { return {(1ull << n) | ((1ull << (k - 1)) - 1), (1ull << k) - 1}; }
};

冒頭で述べたような挙動を実現するには、コンテナと、そのイテレータを自作すればよいです。イテレータはcombinationそのものであるsと、1kbit並んだ演算用の変数であるtしか保持していません。イテレータのインクリメントにnext_combinationを入れれば、combinationがたくさん入ったコンテナのイテレータを一つ進めたような挙動をしてくれます。他にはイテレータへのアクセスと比較、コンテナのbegin()end()を直感通りに実装すればよいです。
私の実装したnext_combinationはk=0の場合をうまく処理できないので、n < 64 && n >= k && k > 0を満たさない引数は初めに弾いています。

速度について

AtCoderのコードテスト機能を使って簡易的な速度の測定を行いました。
n=30,k=15として、以下の三つのコードの速度比較を行いました。
forループは{}_{30}C_{15}=155117520回回ります。

1. こちらのnext_combinationを使ったもの。
ビット演算 (bit 演算) の使い方を総特集! 〜 マスクビットから bit DP まで 〜 #競技プログラミング - Qiita
実行時間:783ms

ソースコード

#include <iostream>
using namespace std;

int next_combination(int sub) {
    int x = sub & -sub, y = sub + x;
    return (((sub & ~y) / x) >> 1) | y;
}

int main() {
    int ans = 0;
    int n = 30, k = 15;
    for (int bit = (1 << k) - 1; bit < (1 << n); bit = next_combination(bit)) {
        ans++;
    }
    cout << ans << "\n";
}

2. 筆者のnext_combinationを使ったもの
実行時間:314ms

ソースコード

#include <iostream>
using namespace std;

int next_combination(int sub, int k_bits) {
    sub += sub & -sub;
    return sub + (k_bits >> __builtin_popcount(sub));
}

int main() {
    int ans = 0;
    int n = 30, k = 15;
    int k_bits = (1 << k) - 1;
    for (int bit = k_bits; bit < (1 << n); bit = next_combination(bit, k_bits)) {
        ans++;
    }
    cout << ans << "\n";
}

3. 構造体に成型したもの
実行時間:359ms

ソースコード

#include <assert.h>
#include <iostream>
using namespace std;

struct combinations {
    using u64 = std::uint64_t;

    uint n;
    uint k;

    combinations(uint n_, uint k_) : n(n_), k(k_) { assert(n_ < 64 && n_ >= k_ && k_ > 0); }

    struct itr {
        u64 s;
        u64 t;

        bool operator!=(itr &right) { return s != right.s; }
        void operator++() {
            s += (s & -s);
            s += t >> __builtin_popcountll(s);
        }
        u64 operator*() { return s; }
    };

    itr begin() { return {(1ull << k) - 1, (1ull << k) - 1}; }
    itr end() { return {(1ull << n) | ((1ull << (k - 1)) - 1), (1ull << k) - 1}; }
};

int main() {
    int ans = 0;
    for (auto i : combinations(30, 15)) {
        ans++;
    }
    cout << ans << "\n";
}


やはり除算を取り除いたことは速度に貢献しているようです。next_combinationを入れ替えると速度が2倍以上になりました。
構造体にしたものはそれより少し遅くなりましたが、どうせ使うときはforループの中の処理も走るため、許容範囲内だと思います。

最後に

筆者が競プロのライブラリのようなものを書いたのはこれが初めてです。はじめは遊びでいじってたんですが、思った以上によさそうなものができたので記事にしてしまいました。おかしい所や改善点など在りましたらご連絡ください。

追記(2024年4月27日)

範囲for文を用いた列挙の実装について、noshi91さんによる部分集合列挙のクリーンな実装を見かけたので、こちらを参考に記事内の実装を更新しました。

12gt 2x2 Spruce Tree Farmの効率期待値

 

12gt 2x2 Spruce Tree Farmとは何か

 

Tree Farm: Minecraft内で木を生産する装置のこと。

2x2 Spruce: 幹の太さが2x2のトウヒの木のこと。

12gt: 木が生成されてから、次の木が生成できるまでの最短の時間が12gt(0.6秒)であること。

つまり

約12gtに1回トウヒの巨木を収穫できる装置のこと。

例:

youtu.be

 

苗木の成長について

空間制約

2x2の木が生成されるためには、北西の苗木を中心として、苗木より上方に5x5x14,苗木と同じ高さに3x3x1の空間*1が必要です*2 。この制約により、苗木に接させることのできるディスペンサーの数は最大4個です。2本の苗木には1つ、1本の苗木には2つのディスペンサーが対応します。

成長段階

苗木には2段階の成長段階があります*3。stageという名前のblock stateになっているので、デバッグ画面で確認することができます。

stage: 0 が成長段階を表す。0か1になる。

苗木は、ディスペンサー化プレイヤーに骨粉を使われた時、またはrandomtick*4を受け取った時に成長段階を1増やすことがあります。成長の確率は、骨粉を使用した場合 45\%、randomtickを受け取った場合 1/7です*5。この記事では、これ以降randomtickについては考慮せず、骨粉を使われた場合を考えます。

成長段階が0のとき

成長段階が0の苗木が成長すると、成長段階が1になります。このとき、上で述べた空間制約は参照されません。この仕様により、Tree Farmがリセットを完了せず、空間制約を満たせていない間も、ディスペンサーを起動することには成長段階を0から1に成長させるメリットがあります。このテクニックのことをpre-bonemealと呼びます。

12gt 2x2 Spruce Tree Farmにおいては、4gt目と8gt目は空間制約を満たせていないため、この時のディスペンサーの起動はpre-bonemealと言えます。詳しい説明は計算段階で行います。

成長段階が1のとき

成長段階が1の苗木が成長するとき、苗木が2x2の苗木の一部で、2x2の木の空間制約を満たしている時、2x2の木を生成します。その苗木以外の3本の苗木の成長段階が1である必要はありません。

2x2の苗木の一部ではなく、1x1の木の空間制約を満たしている時、1x1の木を生成します。当然、2x2のツリーファームではそのような状況は起こらないようになっています。

2x2の苗木の一部で、2x2の木の空間制約を満たしていないとき、その2x2の苗木の成長段階をすべて1に変化させます。これはあまり知られていない事実だと思います。私も計算を始めた時には知らず、実測値とのズレを説明しようと実験をしているうちに気づきました。

2x2の苗木の一部ではなく、1x1の木の空間制約を満たしていないとき、変化は起こりません。但し、内部的には2x2の場合と同様に、成長段階が1の苗木への置き換えが起こっているため、オブザーバー等で検知ができます。BU*6は発生しません。このことを利用して、成長できない苗木をオブザーバーで観察することで、randomtickを1/7の確率で検知することができます。

Tree Farmの動き

ディスペンサーやドロッパーは、最速で4gt(0.2秒)周期で動きます。ツリーファームでは最速の4gt周期で常にディスペンサーを動かします。

苗木を植えるプレイヤーは、クライアントのラグの影響を排除するため、carpet modでサーバー側で操作を行うものとします。この記事では、

/carpet プレイヤー名 use interval 1

を実行した状態での効率を計算します。つまり、プレイヤーは正確に1gtごとに1本苗木を植えます。

細かい動作はFloppyさんによる12gt 2x2 spruce V2に準拠します*7

 

0gt目

前のサイクルの木が成長したtickです。木が成長するTiletickフェーズ,ピストンが動作するBlock Eventフェーズはcarpetのプレイヤーが処理されるEntity phase*8よりも前なので、0gt目のうちに1本目の苗木を植えます*9。植える場所はSE方向の苗木です。

 

1gt目

2本目の苗木をNE方向に植えます。

 

2gt目

レイアウトの都合上、このtickに苗木を植えることはできません。

 

3gt目

3本目の苗木をSW方向に植えます。

ここまでで、ディスペンサーに面する3本の苗木がすべて植わっている状態になります。下の図を参照して下さい。

 

4gt目

まず、Tiletickフェーズにすべてのディスペンサーが動作します。このとき、苗木は3本しか植わっておらず、空間制約が満たされていない状態です。空間制約が満たされていないのは、まだ前の木を収穫している途中だからです。

次に、4本目の苗木をNW方向に植えます。この時初めて2x2の範囲に苗木が植えられたことになります。

苗木を植える順番

 

5gt目~7gt目
ピストンがガシャガシャ動く。

 

8gt目

ディスペンサーが2回目の動作を行います。この時、苗木は2x2すべて植わった状態で、空間制約が満たされていない状態です。

 

9gt目~11gt目
ピストンがガシャガシャ動く。

11gt目には2x2の木の空間制約が満たされている状態になります。

 

12gt目

ディスペンサーが3回目の動作を行います。この時、苗木は2x2をすべて植わった状態で、空間制約が満たされています。つまり、成長段階が1の苗木が成長したとき、2x2の木を生成します。木が生成された場合、次のサイクルの0gt目となります。

 

12+4n gt目

木が生成されなかった場合、木が生成されるまでディスペンサーが4gt周期で動作し続けることになります。木が生成され次第、次のサイクルの0gt目となります。

 

期待値の計算(準備)

木1本あたりの原木の量

2x2のトウヒの木は、基本的には2x2x高さ の大きさの直方体で原木が生成されますが、最上段には1つだけ原木が生成されます。

つまり、原木の量は次の式で表されます。

 Logs=4 \cdot TreeHeight-3

木の高さTreeHeightは、次の式によりランダムに決定されます。ここで、rand(n)は、0以上n以下の自然数を出力する乱数です。

TreeHeight=13+rand(14)+rand(2)

つまり、1本あたりの原木の量の期待値は、期待値の線形性から次のように計算できます。

\begin{eqnarray*} \overline{Logs} &=& 4\cdot\overline{TreeHeight}-3 \\\ &=&49+4\cdot\overline{rand(14)}+4\cdot\overline{rand(2)} \\\ &=&49+4\cdot7+4\cdot1 \\\ &=&81 \end{eqnarray*}

木1本あたりの原木の量の期待値は81個です

 

1時間当たりの木の本数

1時間当たりの木の本数の期待値は次の式で表せます。

 \overline{Trees} = \dfrac{72000}{\overline{duration}}

但し、ここで\overline{duration}は木が生成される間隔(gametick)の期待値を表します。次のように変形すると直感的かもしれません。

\begin{eqnarray*} \dfrac{\overline{Trees}}{72000} &=& \lim_{n\to \infty}\dfrac{\sum_{k=1}^n 1}{\sum_{k=1}^n duration_k} \\\ &=&\lim_{n\to \infty} \dfrac{n}{n \cdot \overline{duration}} \\\ &=&\dfrac{1}{\overline{duration}} \end{eqnarray*}

実際の計測を考えると、\sum_{k=1}^n duration_kは計測時間、\sum_{k=1}^n 1はその時の木の本数を表します。nを限りなく大きくすれば計測時間は限りなく大きくなりますから、この割合は期待値に収束するはずです。

 

効率の期待値

以上より、1時間当たりの原木の効率の期待値は次式で表せます。

\begin{eqnarray*} \overline{Rate}&=&\overline{Logs}\cdot\overline{Trees}\\\ &=&81\cdot\dfrac{72000}{\overline{duration}} \end{eqnarray*}

\overline{duration}の値は次以降で求めていきます。

 

木が生成される間隔の期待値

いよいよ本題です。

状態に対して期待値を対応させる

何を求めたいのかを一度確認しておきましょう。今求めたいのは、木が生成される間隔の期待値でした。これは、次のように言い換えることができます。

0gt目の状態において、次に木が生成されるまでにかかる時間の期待値

なぜこの言いかえをしたのかというと、この記述であれば、期待値をもう少し一般化することができるからです。

ある状態において、次に木が生成されるまでにかかる時間の期待値

0gt目というのはこの記述における特殊状態と言えるわけです。

状態に対して期待値を対応させると何が嬉しいのかというと、期待値に影響しない違いを無視して、ある状態に対する期待値を複数回使いまわせることです。

 

ここで一度Tree Farmのサイクルを思い出しましょう。12gt以降、ディスペンサーは木が生成されるまで4gt周期で骨粉を使い続けるのでした。

次の二つの状況を考えてみます。

12gt目で、ディスペンサーの起動前、すべての苗木の成長段階が1である。

16gt目で、ディスペンサーの起動前、すべての苗木の成長段階が1である。

違うのはそれが何gt目かという点だけです。これらの状態において、次に木が生成されるまでにかかる時間の期待値を考えてみます。二つの状況で期待値が変わらないのが分かるでしょうか。何gt目だろうと、骨粉の成長確率が上がったり、成長段階が勝手に変化したりすることはないため、二つの状況は等価です。

つまり、12gt目以降において、次に木が生成されるまでにかかる時間の期待値は、苗木の成長段階のみに依存するといえます。この時、12gt目以降において考える必要のある状態は2^4通りです。

この状態の数は、さらに減らすことができます。

まず、ディスペンサーに面していない苗木の成長段階は考慮する必要がありません。これは、この苗木をきっかけに木が生成されることがない為です。

次に、ディスペンサーが一つだけ面している苗木同士の区別をつける必要がありません。つまり、片方だけ成長段階が1であるような状態は2通りありますが、どちらの成長段階が1であっても期待値は変わりません。

これらを考えると、考慮する必要のある状態は

  • ディスペンサー1つが面する苗木のうち成長段階が1のものは何本か(3通り)
  • ディスペンサー2つが面する苗木のうち成長段階が1のものは何本か(2通り)

の6通りであることが分かります。

 

と、ここまで書きましたが、ディスペンサーが1つだけ面している苗木同士の区別は付けておいた方が計算式が簡潔に済みます。このため、これ以降次の状態を区別します。

  • SE側(右下)の苗木の成長段階
  • NE側(右上)の苗木の成長段階
  • SW側(左下)の苗木の成長段階

全部で8通りです。この8通りの状態に対応する期待値を表すのに、次の表記を使います。

A_{\alpha,\beta,\gamma}:

12gt目以降であって、SE側(右下)の苗木の成長段階が\alpha,NE側(右上)の苗木の成長段階が\beta,SW側(左下)の苗木の成長段階が\gammaである時、次の木が生成されるまでにかかる時間の期待値

 

12gt目より前の状態についてはとりあえず後回しにします。

 

状態の遷移

前節で12gt目以降について、区別したすべての状態に期待値を対応させました。状態の遷移を考えることで、期待値の線形性を用いて関係式を立式します。

 \displaystyle{A_{\alpha,\beta,\gamma} = \sum_{x_0,x_1,x_2,x_3\in\{0,1\}}0.45^{x_0+x_1+x_2+x_3}\cdot0.55^{4-x_0-x_1-x_2-x_3}\cdot (A_{\alpha+x_0+x_1,\beta+x_2,\gamma+x_3}+4)}

但し、\alpha \geq 2または\beta \geq 2または\gamma \geq 2のとき

 A_{\alpha,\beta,\gamma}=-4

この式の意味を述べます。

まず左辺は注目している状態です。同じgtで、ディスペンサーがすべて起動するような状態での期待値を表しています。ディスペンサーに面する三本の苗木の成長段階は、それぞれ\alpha,\beta,\gammaです。

 

右辺には4gt後の未来をすべて列挙してあります。

 

\displaystyle{\sum_{x_0,x_1,x_2,x_3\in\{0,1\}}}

の意味ですが、x_0,x_1,x_2,x_3が0または1であるような16通りすべての組み合わせについて後ろの数式を足し合わせるという意味です。

 

x_0,x_1,x_2,x_3は、4つのディスペンサーがそのgtに、骨粉による成長を成功させたかを表しています。

 

0.45^{x_0+x_1+x_2+x_3}\cdot0.55^{4-x_0-x_1-x_2-x_3}

は、そのような組み合わせの起こる確率を表します。0.45は骨粉による成長が成功する確率、0.55は失敗する確率です。

 

A_{\alpha+x_0+x_1,\beta+x_2,\gamma+x_3}

は、ディスペンサーがそのような組み合わせで成功/失敗した後に遷移する状態での期待値を表しています。成長に成功した分だけ成長段階が上がるわけですから、元の成長段階にx_0,x_1,x_2,x_3をそれぞれ足し合わせています。

 

このままだと、成長段階が2や3といった、存在しない状態が現れてしまいます。この時、何が起きているかというと、成長段階が1の状態で骨粉による成長が成功しているわけですから、このgtのうちに木が生成されています。

 

逆に、遷移後の成長段階がすべて1以下である場合、また4gt後まで待つ必要があります。よって、この時かかる時間は、遷移後の状態での期待値+4です。これが、A_{\alpha+x_0+x_1,\beta+x_2,\gamma+x_3}に4が足されている理由です。

 

木が生成されたときのことを考えると、この時かかった時間は0であるはずです。整合性を取るため、

\alpha \geq 2または\beta \geq 2または\gamma \geq 2のとき

 A_{\alpha,\beta,\gamma}=-4

とすることで、

A_{\alpha+x_0+x_1,\beta+x_2,\gamma+x_3}+4=0

で変わらずかかった時間を表せるようにしています。

 

式の書き下し

前節では、12gt目以降での木が生成されるまでの時間の期待値を、\alpha,\beta,\gammaといった記号で抽象化して立式しました。このままでは計算できないので、\alpha,\beta,\gammaに代入を行って、総和記号も展開してしまいます。

式を簡潔にするため、以下ではp=0.45,q=0.55とします。

A_{1,1,1}=q^4 (A_{1,1,1}+4)

A_{1,1,0}=q^4 (A_{1,1,0}+4)+q^3p(A_{1,1,1}+4)

A_{1,0,1}=q^4 (A_{1,0,1}+4)+q^3p(A_{1,1,1}+4)

A_{1,0,0}=q^4 (A_{1,0,0}+4)+q^3p(A_{1,1,0}+4)+q^3p(A_{1,0,1}+4)+q^2p^2(A_{1,1,1}+4)

A_{0,1,1}=q^4(A_{0,1,1}+4)+2q^3p(A_{1,1,1}+4)

A_{0,1,0}=q^4(A_{0,1,0}+4)+q^3p(A_{0,1,1}+4)+2q^3p(A_{1,1,0}+4)+2q^2p^2(A_{1,1,1}+4)

A_{0,0,1}=q^4(A_{0,0,1}+4)+q^3p(A_{0,1,1}+4)+2q^3p(A_{1,0,1}+4)+2q^2p^2(A_{1,1,1}+4)

A_{0,0,0}=q^4(A_{0,0,0}+4)+q^3p(A_{0,0,1}+4)+q^3p(A_{0,1,0}+4)+q^2p^2(A_{0,1,1}+4)+2q^3p(A_{1,0,0}+4)+2q^3p(A_{1,0,1}+4)+2q^3p(A_{1,1,0}+4)+2qp^3(A_{1,1,1}+4)

 

よく見ると、このままでは右辺に左辺の内容がまだ残っています。ディスペンサーがすべて失敗した場合は元の状態に戻ってきてしまうためです。これは等式なので、右辺に紛れ込んだ左辺を移行して左辺について解いてしまえばよいです。これを行うと次のようになります。

A_{1,1,1}=\dfrac{4q^4 }{1-q^4}

A_{1,1,0}=\dfrac{4q^4+q^3p(A_{1,1,1}+4)}{1-q^4}

A_{1,0,1}=\dfrac{4q^4+q^3p(A_{1,1,1}+4)}{1-q^4}

A_{1,0,0}=\dfrac{4q^4+q^3p(A_{1,1,0}+4)+q^3p(A_{1,0,1}+4)+q^2p^2(A_{1,1,1}+4)}{1-q^4}

A_{0,1,1}=\dfrac{4q^4+2q^3p(A_{1,1,1}+4)}{1-q^4}

A_{0,1,0}=\dfrac{4q^4+q^3p(A_{0,1,1}+4)+2q^3p(A_{1,1,0}+4)+2q^2p^2(A_{1,1,1}+4)}{1-q^4}

A_{0,0,1}=\dfrac{4q^4+q^3p(A_{0,1,1}+4)+2q^3p(A_{1,0,1}+4)+2q^2p^2(A_{1,1,1}+4)}{1-q^4}

A_{0,0,0}=\dfrac{4q^4+q^3p(A_{0,0,1}+4)+q^3p(A_{0,1,0}+4)+q^2p^2(A_{0,1,1}+4)+2q^3p(A_{1,0,0}+4)+2q^3p(A_{1,0,1}+4)+2q^3p(A_{1,1,0}+4)+2qp^3(A_{1,1,1}+4)}{1-q^4}

 

A_{1,1,1}の式は右辺に定数しかないですから、計算できます。A_{1,1,1}が計算できるので、A_{1,1,0}A_{1,0,1}A_{0,1,1}も計算できます…と続けていくと、上の式の左辺はすべて計算できます。つまり、12gt目以降のすべての状態について、次に木を生成するまでにかかる時間の期待値を計算できたことになります。

 

12gt目より前の状態での期待値

元々知りたかったのは、0gt目での次の木が生成されるまでの期待値でした。しかし、12gt目まで木が生成されることがないにもかかわらず、いちいち4gt目や8gt目の状態を列挙するのは少々面倒です。12gt目のディスペンサーが起動する直前のタイミングで、どのような確率でその成長段階の状態になるのかを直接計算することにしましょう。

P_{\alpha,\beta,\gamma}を12gt目のディスペンサーが起動する直前のタイミングで成長段階が\alpha,\beta,\gammaとなる確率を表すことにします。この時、期待値の線形性より、求める期待値は次の式で表されます。

\displaystyle{\overline{duration}=12+\sum_{\alpha,\beta,\gamma\in \{0,1\}} P_{\alpha,\beta,\gamma}A_{\alpha,\beta,\gamma}}

それぞれの確率を求めましょう。

 

P_{0,0,0}=q^8

成長段階が一つも上がっていなければ、ディスペンサーが成長に成功したことはありません。

 

P_{0,0,1}=2q^7p

P_{0,1,0}=2q^7p

ディスペンサーが1つ面している苗木が1つだけ成長している時、対応するディスペンサーが1度だけ成長に成功したといえます。対応するディスペンサーは4gt目と8gt目の2通り起動できるため、確率は2q^7pです。

ディスペンサーが4gt目と8gt目の両方で成功してしまった場合、8gt目では苗木がすべて植わっている上、成長段階が1の状態で成長しようとするため、すべての苗木の成長段階を1にしてしまいます。したがって、12gt後にこの成長段階の組み合わせにはなりません。

 

P_{0,1,1}=4q^6p^2

1つ上と同じことが、二つの苗木で独立に起こっています。

 

P_{1,0,0}=q^4((1-q^2)q^2+2q^3p)

ディスペンサーが2つ面している苗木の成長段階が1になるとき、2通りのルートが考えられます。

-4gt目に成長段階が1になる。

このとき、8gt目の2回のディスペンサーの起動は両方失敗します。もし成功したら、他の苗木の成長段階が1になってしまうためです。

4gt目の2回のディスペンサーの起動は少なくとも1方が成功してます。こちらは、2回とも成功していても大丈夫です。これは、4gt目の段階では4本目の苗木が植えられておらず、成長段階1の苗木が成長しようとしたとしても、他の苗木の成長段階を1にすることがないからです。

このような確率は、q^6(1-q^2)で表せます。

-8gt目に成長段階が1になる。

このとき、8gt目の2回のディスペンサーの起動のうち、1回のみが成功します。2回とも成功した場合、他の苗木の成長段階が1になってしまうためです。

このような確率は、2q^7pで表せます。起動するディスペンサーが2通りあるためです。

式を簡潔に書き直すと、次のようになります。

P_{1,0,0}=q^6p^2+4q^7p

 

P_{1,0,1}=2q^3p((1-q^2)q^2+2q^3p)

P_{1,1,0}=2q^3p((1-q^2)q^2+2q^3p)

今まで出てきた式の複合です。

それぞれ簡潔に書き直すと、次のようになります。

P_{1,0,1}=2q^5p^3+8q^6p^2

P_{1,1,0}=2q^5p^3+8q^6p^2

 

P_{1,1,1}=1-P_{0,0,0}-P_{0,0,1}-P_{0,1,0}-P_{0,1,1}-P_{1,0,0}-P_{1,0,1}-P_{1,1,0}

最後の項なので余事象を取ります。8gt目にすべての苗木の成長段階を1にされるようなパターンが多いので、愚直に計算するのは大変だと思います。

代入すると、次のようになります。

P_{1,1,1}=1-q^8-8q^7p-21q^6p^2-4q^5p^3

 

計算のまとめ

ここまで計算してきた式で、期待値を計算できます。必要な式を下に列挙します。

p=0.45,q=0.55

 

A_{1,1,1}=\dfrac{4q^4 }{1-q^4}

A_{1,1,0}=\dfrac{4q^4+q^3p(A_{1,1,1}+4)}{1-q^4}

A_{1,0,1}=\dfrac{4q^4+q^3p(A_{1,1,1}+4)}{1-q^4}

A_{1,0,0}=\dfrac{4q^4+q^3p(A_{1,1,0}+4)+q^3p(A_{1,0,1}+4)+q^2p^2(A_{1,1,1}+4)}{1-q^4}

A_{0,1,1}=\dfrac{4q^4+2q^3p(A_{1,1,1}+4)}{1-q^4}

A_{0,1,0}=\dfrac{4q^4+q^3p(A_{0,1,1}+4)+2q^3p(A_{1,1,0}+4)+2q^2p^2(A_{1,1,1}+4)}{1-q^4}

A_{0,0,1}=\dfrac{4q^4+q^3p(A_{0,1,1}+4)+2q^3p(A_{1,0,1}+4)+2q^2p^2(A_{1,1,1}+4)}{1-q^4}

A_{0,0,0}=\dfrac{4q^4+q^3p(A_{0,0,1}+4)+q^3p(A_{0,1,0}+4)+q^2p^2(A_{0,1,1}+4)+2q^3p(A_{1,0,0}+4)+2q^3p(A_{1,0,1}+4)+2q^3p(A_{1,1,0}+4)+2qp^3(A_{1,1,1}+4)}{1-q^4}

 

P_{1,1,1}=1-q^8-8q^7p-21q^6p^2-4q^5p^3

P_{1,1,0}=2q^5p^3+8q^6p^2

P_{1,0,1}=2q^5p^3+8q^6p^2

P_{0,1,1}=4q^6p^2

P_{0,1,0}=2q^7p

P_{0,0,1}=2q^7p

P_{0,0,0}=q^8

 

\displaystyle{\overline{duration}=12+\sum_{\alpha,\beta,\gamma\in \{0,1\}} P_{\alpha,\beta,\gamma}A_{\alpha,\beta,\gamma}}

 

\overline{Rate}=81\cdot\dfrac{72000}{\overline{duration}}

上から順に計算することで、効率の期待値を求められます。

 

結論

desmosに計算式を打ち込んだものがこちらになります。

www.desmos.com

12gt 2x2 spruceの効率の期待値は 463691.576977/hです。これは実測値と100のオーダー程度まで一致しています。

 

おまけ

8gt

世の中には8gt 2x2 Spruce Tree Farmというものもあります。ラグ当たりの効率の観点から、一般的には12gtの方が良いとされていますが、当然最も速いです。

youtu.be

8gt 2x2spruce | Desmos

8gtのバージョンも計算しました。最後の確率分布を少し弄るだけでできます。

 

失敗ver

実は、先ほどの12gtの計算式は第一版ではありません。最初に計算したときは、効率の期待値は464161.120876/hと計算されていました。このずれは、4gt目に2x2の苗木がすべて植え終わっていると勘違いしていたことによるものです。

2x2spruce | Desmos

この微妙な誤差の指摘と、4gt目に2x2の苗木が植え終わっていないせいで成長段階がすべて1に置き換えられる現象が起こらないことの指摘をScorpioさんからいただきました。実は、そもそもこの計算を始めたのも、彼が理論値の計算を探していたことに由来しています。ありがとうございます。

彼は12gt 2x2 spruce tree farmを作り終わったところで、そのうち彼のYoutubeチャンネルから動画が上がるかもしれません。投稿されたら視聴してみてください。

www.youtube.com

 

 

 

分かりづらい所などありましたら、お気軽にご連絡ください。

https://twitter.com/Teitoku_toaru

*1:原木などいくつかの例外があります。

*2:こちらを参照して下さい。Spruce – Minecraft Wiki

*3:余談ですが、ツツジとネザーの木には成長段階がありません。これはオークの原木を集めるうえでツツジのほうがオークより効率が良い一因になっています。

*4:詳しくはこちらを参照して下さい。Tick – Minecraft Wiki
プレイヤーの近くでは、 1/16^3の確率の抽選が1tickあたりに3回行われます。

*5:randomtickについては、成長するために追加で明るさレベルが9以上である必要があります。

*6:wiki内のNC updateと同じものです。Block update – Minecraft Wiki

*7:冒頭で紹介した動画のものの改良版です。Tree HuggersというTree Farmについてのdiscordサーバーで配布されています。

*8:carpet TIS additionのルールfakePlayerTicksLikeRealPlayerを有効化している場合、Player Inputフェーズで動作しますが、いずれにせよTTフェーズとBEフェーズよりも後で処理されます。

*9:Tickフェーズとは何かについては省略します。飽くまで理由なので、分からなければスルーしても大丈夫です。

参考になる資料を置いておきますが、Tickフェーズとは何かを説明してくれるようなものではないです。GitHub - Fallen-Breath/MinecraftGamePhase: Minecraft game phase list and its history