【Qiskit 2.1.1対応】量子モンテカルロ法でπを計算しよう！

はじめにこのノートブックでは、最も有名な数学定数の一つである円周率πを、量子コンピュータ（のシミュレータ）を使って計算します。これは「量子コンピュータは計算が速い」ということを示すためではなく、「量子状態の重ね合わせと測定」という基本的な性質が、どのようにして確率的なシミュレーションに応用できるのかを理解するための、教育的な題材です。
使用するバージョン:
Python 3.11
Qiskit 2.1.1
qiskit-aer (Qiskit 1.0から分離された公式シミュレータ)
matplotlib, numpy

Step 1: 準備 (ライブラリのインポート)

まずは、計算と可視化に必要なライブラリをインポートします。qBraid環境ではこれらはプリインストールされています。

# 必要なライブラリをインポート
import numpy as np
import matplotlib.pyplot as plt
import time

# Qiskit関連のインポート
from qiskit import QuantumCircuit, transpile
from qiskit_aer import AerSimulator

Step 2: 古典的なモンテカルロ法によるπの計算

量子計算と比較するため、まずは古典的なコンピュータでπを計算してみましょう。
方法はシンプルです。一辺が1の正方形を考えます。その中にランダムな点をたくさん打ちます。中心を原点としたとき、半径1の円（実際には1/4円）の内側に入った点の数を数えます。
(円の内側の点の数) / (全ての点の数) ≈ (1/4円の面積) / (正方形の面積) = (π/4) / 1 この関係から、π ≈ 4 * (円の内側の点の数) / (全ての点の数)
となります。

def calculate_pi_classical(n_samples: int):
    """古典的なモンテカルロ法でπを計算する関数"""
    start_time = time.time()
    
    # 0から1の範囲でランダムな(x, y)座標をN個生成
    x = np.random.rand(n_samples)
    y = np.random.rand(n_samples)
    
    # 原点からの距離の2乗が1未満（円の内側）の点の数を数える
    inside_circle_count = np.sum(x**2 + y**2 < 1.0)
    
    # πの近似値を計算
    pi_approx = 4 * inside_circle_count / n_samples
    
    execution_time = time.time() - start_time
    
    print(f"--- 古典計算 ---")
    print(f"サンプル数: {n_samples}")
    print(f"πの近似値: {pi_approx:.6f}")
    print(f"実行時間: {execution_time:.6f} 秒")
    return pi_approx, execution_time

# 実行してみる
classical_samples = 1_000_000
c_pi, c_time = calculate_pi_classical(classical_samples)

Step 3: 量子モンテカルロ法によるπの計算

次に、量子コンピュータで同じことをしてみます
どうやってランダムな点を生成するのでしょうか？
ここで**「重ね合わせ」**を使います。座標の各軸（x, y）を表現するために、複数個の量子ビットを用意します。例えば、x軸に5ビット、y軸に5ビット使えば、2^5=32段階の解像度で座標を表現できます。
全ての量子ビットにアダマール(H)ゲートをかけ、|0⟩と|1⟩が均等に重なった状態を作ります。
この状態で測定すると、各量子ビットは完全にランダムに0か1のどちらかに確定します。
得られたビット列（例: 01101）を座標の値に変換します。これが量子的に生成された「ランダムな点」です。あとは古典的な方法と同じく、円の内側に入った点の割合を数えます。

def calculate_pi_quantum(n_qubits_per_axis: int, shots: int):
    """量子モンテカルロ法でπを計算する関数"""
    start_time = time.time()
    
    total_qubits = 2 * n_qubits_per_axis
    
    # 1. 量子回路の作成
    qc = QuantumCircuit(total_qubits)
    
    # 2. 全ての量子ビットにHゲートをかけて重ね合わせ状態を生成
    qc.h(range(total_qubits))
    
    # 3. 全ての量子ビットを測定
    qc.measure_all() # Qiskit 1.0以降の便利な機能

    # 回路図を描画して確認
    print("使用する量子回路:")
    display(qc.draw('mpl')) # Jupyter Notebook/Lab用の表示方法
    
    # 4. シミュレータで実行
    simulator = AerSimulator()
    # 回路をシミュレータ用に最適化（トランスパイル）
    compiled_circuit = transpile(qc, simulator)
    # 実行して結果を取得
    result = simulator.run(compiled_circuit, shots=shots).result()
    counts = result.get_counts()
    
    # 5. 結果の集計
    inside_circle_count = 0
    max_val = 2**n_qubits_per_axis # 座標の最大値 (例: 5qubitなら32)

    for bitstring, count in counts.items():
        # bitstringをxとyに分割
        x_bits = bitstring[:n_qubits_per_axis]
        y_bits = bitstring[n_qubits_per_axis:]
        
        # 2進数を10進数に変換し、[0, 1)の範囲に正規化
        x = int(x_bits, 2) / max_val
        y = int(y_bits, 2) / max_val
        
        # 円の内側にあるか判定
        if x**2 + y**2 < 1.0:
            inside_circle_count += count
            
    pi_approx = 4 * inside_circle_count / shots
    execution_time = time.time() - start_time
    
    print(f"\n--- 量子計算 ---")
    print(f"座標ごとの量子ビット数: {n_qubits_per_axis} (合計 {total_qubits} q)")
    print(f"ショット数 (＝サンプル数): {shots}")
    print(f"πの近似値: {pi_approx:.6f}")
    print(f"実行時間: {execution_time:.6f} 秒")
    return pi_approx, execution_time

# 実行してみる
# 量子計算はシミュレーションに時間がかかるため、サンプル数は少なめにします
qubits_per_coord = 8 # 座標の解像度 (x,yそれぞれ8bit = 256段階)
quantum_shots = 8192 # 測定回数

q_pi, q_time = calculate_pi_quantum(qubits_per_coord, quantum_shots)

Step 4: 考察

結果を比較してみましょう。
精度: おそらく、遥かに多くのサンプルを使った古典計算の方が、真のπの値に近い結果を出したはずです。量子計算の精度を上げるには、ショット数を増やすか、量子ビット数を増やして座標の解像度を上げる必要がありますが、どちらも計算コストが増大します。
実行時間: 量子シミュレーションは、古典コンピュータ上で量子ビットの状態ベクトル（巨大な複素数配列）を計算しているため、非常に時間がかかります。古典計算の方が圧倒的に速いことがわかります。
結論: このπ計算の例では、量子コンピュータは古典コンピュータに全く敵いません。
この演習の重要な点は、**「量子コンピュータの重ね合わせと測定が、質の高い乱数生成器として機能し、モンテカルロ法のような確率的アルゴリズムの基盤となりうる」**という原理を理解することです。