定義に沿って計算すると、行列累乗をアレンジした感じの後述の謎演算に帰着できて、二分法で高速化できます。以下ご説明していきます。

事象 \(A _ { i, j } ( s )\)と確率 \(p _ { i, j } ( s )\)、確率変数 \(X _ { i, j } ( s )\) の定義です。

長さ固定の（ある確率分布に従う）ランダム文字列 \(s\) が各文字が 0, 1 のいずれかを取るとしましょう。このとき、\(s\) の先頭、末尾の文字がそれぞれ \(i, \ j\) であるという事象を \(A _ { i, j } ( s )\) と置きます。この事象が起こる確率を \(p _ { i, j } ( s )\) と置きます。また確率変数 \(X _ { i, j } ( s )\) を、事象 \(A _ { i, j } ( s )\) が起きるときには必要な操作回数の期待値、さもなくば \(0\) であると定めます。

期待値 \(E [ X _ { i, j } ( s ) ]\) の計算です。

このとき、長さ固定の\(2\) つのランダム文字列 \(s, t\) について、それを並べてできるランダム文字列を \(s + t\) とおくと、

\[ p _ { i, l } ( s+ t ) = \sum _ { j = 0 } ^ 1 \sum _ { k = 0 } ^ 1 p _ { i, j } ( s ) p _ { k, l } ( t ) \\ E [ X _ { i, l } ( s + t ) ] = \sum _ { j = 0 } ^ 1 \sum _ { k = 0 } ^ 1 E [ X _ { i, j } ( s ) ] p _ { k, l } ( t ) + p _ { i, j } ( s ) E [ X _ { k, l } ( t ) ] + f ( i, j, k, l )\\ \]

（ただし \(f\) は

\[ f ( i, j, k, l ) = \begin{cases} 1 & ( i = j, j \neq k, k = l ),\\ -1 & ( i \neq j, j =k, k \neq l) \\ \end{cases} \]

で定まる関数です。）と書けますから、これで計算すると良いです。

たとえば、このように実装することができます。（ただしあまりを取る操作を考慮から外すためにすべて f64 型で書いてありますので、参考にされる場合は適宜よろしくおねがいです。）（これは Rust です。）

fn mul(a: [[[f64; 2]; 2]; 2], b: [[[f64; 2]; 2]; 2]) -> [[[f64; 2]; 2]; 2] {
    let mut c = [[[0.0; 2]; 2]; 2];
    for (((i, j), k), l) in (0..2)
        .cartesian_product(0..2)
        .cartesian_product(0..2)
        .cartesian_product(0..2)
    {
        let [p0, e0] = a[i][j];
        let [p1, e1] = b[k][l];
        c[i][l][0] += p0 * p1;
        c[i][l][1] += e0 * p1
            + e1 * p0
            + p0 * p1 * ((i == j && j != k && k == l) as f64 - (i != j && j == k && k != l) as f64);
    }
    c
}

二分法による高速化でトドメです。

ところでこの演算は定義より結合的ですから、入力の \(S\) について計算したあと、\(K\) 回の繰り返しは、二分累乗により高速化することができ、\(\Theta ( \lvert S \rvert \lg K )\) が達成できます。