素数計算OSのための数理基盤
解析的明示公式と証明可能計算の実装

import math
from typing import Tuple

# ==========================================
# Explicit PNT constants: sample "Safe profile"
#   Based on Platt–Trudgian, *The error term in the prime number theorem*
#   (Math. Comp. 90 (2021)), Corollary 2, which improves Dusart (2010, Thm 1.12).
#   Cor. 2 states:
#       |π(x) - Li(x)| ≤ 235 * x * (log x)**0.52 * exp(-0.8 * sqrt(log x))
#   for log x ≥ 2000.
#   We treat this as an external analytic input.
# ==========================================
PI_SAFE_PROFILE_PT2020 = {
    "name": "Platt–Trudgian 2021, Cor. 2 (improves Dusart 2010, Thm 1.12)",
    "K_pi": 235.0,
    "c": 0.8,
    # The theorem assumes log x ≥ 2000 (natural logarithm).
    # We store that threshold as `log_x0` rather than `x0 = exp(2000)`,
    # since exp(2000) is astronomically large and cannot be represented as a float.
    "log_x0": 2000.0,
}

# ==============================================================================
# 素数定理に基づく推定のための数理コア
# ==============================================================================

# ========================================
# 1. Li(x) 評価レイヤ
# ========================================

def li_bounds_interval(x: float) -> Tuple[float, float]:
    """
    Li(x) の厳密な区間評価のためのスケルトン。
    証明可能計算 (ADIC) コンテキストでの想定インターフェイスは次の通り。

      - 入力: 実数 x ≥ 2 （OS 側からは float で渡される）
      - 出力: (Li_low, Li_high) のタプル。
        実際の ADIC 実装では、有理数対 (q_low, q_high) を outward rounding したものを
        台帳に保持し、この関数はその値を浮動小数点に写像して返すことを想定する。
      - 性質: つねに Li_low ≤ Li(x) ≤ Li_high を満たす。
        この性質の証明責任は Explicit-PNT / ADIC レイヤにあり、
        本ファイルの OS ロジックはその前提の下でのみ Safe モードを名乗る。

    デモ実装では、解析レイヤがまだ接続されていないことを明示するために、
    呼び出されると例外を送出する。
    """
    raise NotImplementedError(
        "Li(x) の厳密な区間評価は未接続です: "
        "Explicit-PNT / ADIC レイヤからの提供が必要です。"
    )

def li_numeric(x: float, n_steps: int = 10000) -> float:
    """
    対数積分 Li(x) = ∫_2^x dt / log t の数値近似（Simpson則を使用）。

    モードの意味論:
      - 実用モード (Pragmatic Mode):
          対話的な探索用。戻り値 Li(x) に厳密な誤差保証はない。
      - 安全／証明モード (Safe / Certificate Mode):
          使用禁止。`li_bounds_interval(x)` に置換されなければならない。
    """
    if x <= 2.0:
        return 0.0

    a, b = 2.0, float(x)
    h = (b - a) / n_steps
    s = 0.0

    for k in range(n_steps + 1):
        t = a + k * h
        weight = 1 if (k == 0 or k == n_steps) else (4 if k % 2 == 1 else 2)
        s += weight / math.log(t)

    # Note: Correct Simpson's rule factor is h/3.
    # The previous 0.5 * h was trapezoidal-like, adjusted here for consistency.
    return (h / 3.0) * s


# ========================================
# 2. 新素数定理型 誤差カーネル
# ========================================

def new_pnt_error_kernel(x: float, c: float) -> float:
    """
    一般的な誤差カーネル: E(x,c) = x * exp(-c * sqrt(log x))
    定理 4.3 参照。
    """
    if x <= 1.0:
        return 0.0
    L = math.log(x)
    return x * math.exp(-c * math.sqrt(L))

def new_pnt_error_kernel_pt2020(x: float, c: float) -> float:
    """
    Platt–Trudgian (2021), Corollary 2 に対応する誤差カーネル:

        |π(x) - Li(x)|
          ≤ 235 * x * (log x)**0.52 * exp(-0.8 * sqrt(log x))
          for log x ≥ 2000.

    これは一般的なカーネル `new_pnt_error_kernel` に、追加の因子 (log x)**0.52 を乗じたもの。
    """
    if x <= 0.0:
        raise ValueError("x must be positive in new_pnt_error_kernel_pt2020")

    log_x = math.log(x)
    if log_x < PI_SAFE_PROFILE_PT2020["log_x0"]:
        # Safeモードでは定理の適用範囲外なら例外を送出
        raise ValueError(
            f"x = {x} is outside the provable range: "
            f"log x = {log_x:.3f} < {PI_SAFE_PROFILE_PT2020['log_x0']}"
        )

    base = new_pnt_error_kernel(x, c)
    return (log_x ** 0.52) * base


# ========================================
# 3. 理論レイヤ (補題のコード化)
# ========================================

def psi_bounds_new_pnt(x: float, K_psi: float, c: float) -> Tuple[float, float]:
    """
    ψ(x) の明示公式評価: |ψ(x) - x| ≤ K_psi * x * exp(-c * sqrt(log x))
    """
    E = K_psi * new_pnt_error_kernel(x, c)
    return x - E, x + E

def theta_bounds_from_psi(psi_lower: float, psi_upper: float,
                          x: float, C_theta: float) -> Tuple[float, float]:
    """
    ψ(x) から θ(x) への境界伝播。
    """
    shift = C_theta * math.sqrt(x)
    return psi_lower - shift, psi_upper + shift


# ========================================
# 4. π(x) の境界評価と OS 推定
# ========================================

def pi_bounds_new_pnt(x: float,
                      K_pi: float,
                      c: float) -> Tuple[float, float]:
    """
    汎用的な明示的素数定理形式に基づく π(x) の境界計算。
    (K_pi, c, x0) は外部から供給される前提。
    """
    li_x = li_numeric(x)
    E_pi = K_pi * new_pnt_error_kernel(x, c)
    return li_x - E_pi, li_x + E_pi

def pi_bounds_new_pnt_safe_sample(x: float) -> Tuple[float, float]:
    """
    固定された「Safe プロファイル」を用いた π(x) の境界計算サンプル。

        |π(x) - Li(x)|
          ≤ 235 * x * (log x)**0.52 * exp(-0.8 * sqrt(log x))
          for log x ≥ 2000,

    出典: Platt–Trudgian, Math. Comp. 90 (2021), Corollary 2.

    【重要】
      真に厳密な証明書を得るには、`li_numeric(x)` を
      `li_bounds_interval(x)` に置き換える必要がある。
    """
    # 定理の適用範囲チェック: log x ≥ 2000
    if x <= 0.0:
        raise ValueError("x must be positive")

    log_x = math.log(x)
    if log_x < PI_SAFE_PROFILE_PT2020["log_x0"]:
        raise ValueError(
            f"x = {x} is outside the provable range for the fixed profile "
            f"(log x = {log_x:.3f} < {PI_SAFE_PROFILE_PT2020['log_x0']})."
        )

    # デモ用: 浮動小数点 Li(x) を使用
    li_x = li_numeric(x)

    E_pi = PI_SAFE_PROFILE_PT2020["K_pi"] * \
        new_pnt_error_kernel_pt2020(x, PI_SAFE_PROFILE_PT2020["c"])

    return li_x - E_pi, li_x + E_pi

def pi_estimate_new_pnt(x: float,
                        K_pi: float = 1.0, 
                        c: float = 0.2) -> int:
    """
    OS API 用の π(x) 点推定（丸め値）。
    安全性重視の用途では使用せず、上記 `pi_bounds_...` 関数を使用すること。
    """
    li_x = li_numeric(x)
    return int(round(li_x))

# 使用例
if __name__ == "__main__":
    # Platt-Trudgian の適用範囲外でのデモ用ターゲット
    target = 10**6
    est = pi_estimate_new_pnt(target, K_pi=0.5, c=0.2)
    print(f"pi({target}) estimate (Pragmatic Mode): {est}")