Pythonで円周率に最も近い分数を求める方法：分母の範囲制約付き近似アルゴリズム

問題概要

2つの長整数値 maximum（最大値）と minimum（最小値）が与えられているとします。このとき、min <= d <= max を満たす分数 n/d のうち、|n/d − π| が最小になるものを見つける必要があります。ここで π = 3.14159265... です。条件を満たす分数が複数存在する場合は、分母が最も小さい分数を返します。

たとえば、minimum = 1、maximum = 10 が入力された場合、出力は 22/7 になります。

解法のアプローチ

この問題はファレイ数列（Farey sequence）の性質を利用すると効率的に解けます。ファレイ数列では、隣接する2つの分数 a/b と c/d に対して、(a+c)/(b+d) という新しい分数が必ず両者の間に現れるという性質があります。これを利用して、π の小数部分に近い分数を段階的に絞り込んでいきます。

具体的な手順は以下の通りです。

• P := 分数 (5706674932067741 / 1816491048114374) − 3（π の高精度近似値から整数部分 3 を引いた値）
• a := 0, b := 1, c := 1, d := 1 で初期化する
• farey := ペアの配列。初期状態として (a, b) と (c, d) の2つのペアを格納する
• 以下を無条件に繰り返します。
  • f := b + d
  • f > maximum − minimum の場合、ループを抜ける
  • e := a + c
  • farey の末尾にペア (e, f) を追加する
  • P < e/f の場合、c := e、d := f と更新する
  • それ以外の場合、a := e、b := f と更新する
• p_min := (P × minimum) の切り捨て値
• minimum <= maximum の間、以下を繰り返します。
  • c := 0, d := 0 でリセットする
  • farey 内の各ペア (a, b) について以下を確認します。
    • minimum + b > maximum の場合、内側のループを抜ける
    • |（p_min + a）/（minimum + b）− P| < |p_min / minimum − P| を満たす場合、c := a、d := b として内側のループを抜ける
  • d が 0 のままの場合、外側のループを抜ける
  • p_min := p_min + c
  • minimum := minimum + d
• 分数 (p_min + 3 × minimum) / minimum を返す

実装例

理解を深めるために、以下のPythonコードを見てみましょう。

from fractions import Fraction

def solve(minimum, maximum):
    # 円周率の高精度近似値から整数部 3 を引いた分数
    P = Fraction(5706674932067741, 1816491048114374) - 3

    a, b, c, d = 0, 1, 1, 1
    farey = [(a, b), (c, d)]

    # ファレイ数列の生成
    while True:
        f = b + d
        if f > maximum - minimum:
            break

        e = a + c
        farey.append((e, f))
        if P < Fraction(e, f):
            c, d = e, f
        else:
            a, b = e, f

    p_min = int(P * minimum)

    # 最良近似分数の探索
    while minimum <= maximum:
        c, d = 0, 0
        for a, b in farey:
            if minimum + b > maximum:
                break
            if abs(Fraction(p_min + a, minimum + b).real - P) < abs(Fraction(p_min, minimum).real - P):
                c, d = a, b
                break
        if d == 0:
            break
        p_min += c
        minimum += d

    return "{}/{}".format(p_min + 3 * minimum, minimum)

minimum = 1
maximum = 10
print(solve(minimum, maximum))

入力

1, 10

出力

22/7

コードのポイント

このコードでは、まず標準ライブラリの fractions.Fraction を使って円周率の小数部分を厳密な有理数として表現しています。これにより、浮動小数点演算による誤差を完全に回避できます。

最初の while ループでは、分母の増分が maximum − minimum を超えない範囲でファレイ数列を構築します。続くループでは、実際の分母の下限 minimum から順に候補分数との誤差を比較しながら、条件を満たす最良の近似を探索します。

その結果、分母が 1〜10 の範囲で π に最も近い分数である 22/7（約 3.142857...）が得られます。これは古くから知られる有名な円周率の近似分数であり、真の値との誤差はわずか約 0.00126 です。