Pythonで乱数を含む無限連根式の期待値を求めるプログラム

問題の概要

整数 n が与えられているとします。ここで、x = rand() mod n と定義します。rand() は 0 以上 10^100 以下の整数を一様な確率で生成する関数です。さらに、次のような無限に入れ子になった平方根(連根式)を考えます。

$$Y = \sqrt{x+\sqrt{x+\sqrt{x+\sqrt{x+...}}}}$$

このとき、Y の期待値を求めるのが課題です。n は 1 以上 5×10^6 以下の範囲にあるものとします。

たとえば、入力が n = 5 の場合、出力は 1.696 になります。

解法のアイデア:閉じた形への変換

一見すると果てしなく続く式ですが、Y の定義自体が再帰的な構造を持っているため、二次方程式に帰着できます。Y = √(x + Y) が成り立つので、両辺を二乗すると

$$Y^2 = x + Y$$

となり、Y は非負であることから、その解は

$$Y = \frac{1+\sqrt{4x+1}}{2}$$

と表せます。また、x は 0 から n−1 までの各値をそれぞれ確率 1/n で等しく取るため、期待値は次の総和として書けます。

$$E[Y]=\frac{1}{2n}\sum_{k=0}^{n-1}\left(1+\sqrt{4k+1}\right)$$

計算手順

この総和を毎回素朴に計算すると時間がかかるため、次の方針で高速化します。

• err := 2235.023971557617(事前に求めた誤差補正定数)を用意する
• max_n := 5 × 10^6 とする
• 累積和リスト pref を、初期値 0 のみを含むリストとして作成する
• i を 1 から 5 × 10^6 まで動かしながら、pref の末尾に「直前の要素 +(1 + √(4i+1)) × 0.5」を追加していく
• n < max_n の場合:前計算済みの累積和を利用し、pref[n − 1] / n を返す
• それ以外の場合:積分近似に基づく漸近式を利用する
  • total := (4 × (n − 1) + 5)^1.5 ÷ 6 − 5^1.5 ÷ 6 − err
  • ans := 0.5 + total ÷ (2 × n)
  • ans を返す

この手法の利点は、小さな n に対しては累積和による正確な値が得られ、より大きな n に対しては O(1) の漸近式で即座に答えを求められる点です。

実装例(Python)

理解を深めるために、以下の実装を見てみましょう。

def solve(n):
    err = 2235.023971557617
    max_n = 5 * 10**6

    pref = [0]
    for i in range(1, 5 * 10**6):
        pref.append(pref[-1] + (1 + (4 * i + 1)**0.5) * 0.5)

    if n < max_n:
        return pref[n - 1] / n
    else:
        total = (4 * (n - 1) + 5)**1.5 / 6 - 5**1.5 / 6 - err
        ans = 0.5 + total / (2 * n)
        return ans

n = 5
print(solve(n))

入力

5

出力

1.69647248786

まとめ

無限連根式は「式全体が自分自身を含む」という性質を利用すると二次方程式へ変換できます。本記事の例では、Y = (1 + √(4x+1)) / 2 への帰着によって、複雑に見える期待値計算を明快な総和の問題に落とし込めました。範囲内の n には累積和による前計算を、それ以上の n には誤差補正付きの漸近近似を組み合わせることで、幅広い入力に対して高速かつ高精度に期待値を求めることが可能です。