#!/usr/bin/env python3 # -*- coding: utf-8 -*- """第834コマ — Lean の loI を Python で完全再現し、窓ごとの深さ配分を設計する。 目的: Lean 内の逆数和の下界を d=2（3.0085286163）から 3.00853 台へ上げる。 方針: 一律 d=3（葉 40 万）ではなく、窓ごとに深さを変える（利得/葉 の比で配る）。 Lean 側の定義（Erdos774.T / Erdos814.mxv / BlkF / loI / Erdos815.cntTab）と一対一。 """ import re, sys, time from functools import lru_cache import os # [配布版の変更点 / changed for distribution] 元は書き手の機体の絶対パスを指していた箇所を、このディレクトリからの相対に変えた。計算の中身は無改変。 LEAN = os.path.join(os.path.dirname(os.path.abspath(__file__)), "Shioriproofs") + "/" # main() でだけ使う（wins50/wins100b を Lean ファイルから読む）。erdos1110f4 からの import では使わない SC = 10 ** 100 # ---------- Lean の定義の再現 ---------- def bb(p): return (p - 1) // 2 @lru_cache(maxsize=None) def T(p, k): b = bb(p) if b < k: return (k - b - 1) * (k - b) // 2 else: return (b - k) * (b + 1 - k) // 2 @lru_cache(maxsize=None) def mxv(p, n): if n == 0: return 0 return (p - 1) * (2 * p - 1) ** (n - 1) + mxv(p, n - 1) @lru_cache(maxsize=None) def cntTab(p, n): """n 桁・各桁 1 and out[-1] == 0: out.pop() return tuple(out) def cfT(p, n, r): row = cntTab(p, n) return row[r] if r < len(row) else 0 def loI(p, d, a, n, r): """Lean の loI p SC (cfT p) d a n r。整数。""" if d == 0 or n == 0: return cfT(p, n, r) * (SC // (a + mxv(p, n))) s = 0 base = (2 * p - 1) ** (n - 1) for k in range(p): t = T(p, k) if t <= r: s += loI(p, d - 1, a + k * base, n - 1, r - t) return s @lru_cache(maxsize=None) def leaves(p, d, n, r): """loI の葉の数（枝刈り T<=r を反映）。""" if d == 0 or n == 0: return 1 c = 0 for k in range(p): t = T(p, k) if t <= r: c += leaves(p, d - 1, n - 1, r - t) return c # ---------- 窓の一覧を Lean のファイルから抜く ---------- def parse_wins(fname, name): src = open(LEAN + fname, encoding="utf-8").read() i = src.index("def " + name) j = src.index("[", i) k = src.index("]\n", j) body = src[j + 1:k] nums = re.findall(r"\d+", body) assert len(nums) % 4 == 0, (name, len(nums)) ws = [] for m in range(0, len(nums), 4): ws.append(tuple(int(x) for x in nums[m:m + 4])) # (p, n, r, a) return ws def main(): t0 = time.time() w50 = parse_wins("Erdos815b.lean", "wins50") w100b = parse_wins("Erdos815c.lean", "wins100b") wins = w50 + w100b print("[窓] wins50 = %d, wins100b = %d, 合計 = %d" % (len(w50), len(w100b), len(wins))) # --- (A) d=2 の既知の値を再現できるか（壁8''：道具の一致確認を先に） --- KNOWN50 = 43186594426047909037485327285886450590907364648618651131920393748013061546682393613679261480657280 KNOWN100b = 16484531567385067650538697947409403393169352146330921496288812445275830441250329638192807936 s50 = sum(loI(p, 2, a, n, r) for (p, n, r, a) in w50) s100b = sum(loI(p, 2, a, n, r) for (p, n, r, a) in w100b) print("[再現] wins50 d=2 :", "一致" if s50 == KNOWN50 else "不一致", s50 == KNOWN50) print("[再現] wins100b d=2:", "一致" if s100b == KNOWN100b else "不一致", s100b == KNOWN100b) if s50 != KNOWN50 or s100b != KNOWN100b: print(" Python=", s50, s100b); sys.exit(1) HEAD = 3004209955304941 * 10 ** 85 # 頭 Z0 の下界（SC 倍） RECORD = 300849 * 10 ** 95 tot2 = HEAD + s50 + s100b print("[d=2 一律] 総和 = %.13f" % (tot2 / SC)) # --- (B) 窓ごとに深さ 0..8 の値と葉の数 --- DMAX = 6 tbl = [] # (idx, p, n, r, a, [val_d], [leaf_d]) for idx, (p, n, r, a) in enumerate(wins): vals, lvs = [], [] for d in range(0, DMAX + 1): if d > 0 and d > n: vals.append(vals[-1]); lvs.append(lvs[-1]); continue lv = leaves(p, d, n, r) if lv > 25000: # これ以上は測らない（設計上使わない） vals.append(None); lvs.append(lv); continue vals.append(loI(p, d, a, n, r)); lvs.append(lv) tbl.append((idx, p, n, r, a, vals, lvs)) if idx % 40 == 0: print(" ...窓 %d / %d %.1f 秒" % (idx, len(wins), time.time()-t0)); sys.stdout.flush() print("[表] 完了 %.1f 秒" % (time.time() - t0)) # --- (C) 一律 d の総和 --- print() print("| 一律 d | 葉の総数 | 総和（頭込み） | 記録 3.00849 との差 |") print("|---|---|---|---|") for d in range(0, 5): if any(v[5][d] is None for v in tbl): print("| %d | (未測) | | |" % d); continue s = HEAD + sum(v[5][d] for v in tbl) L = sum(v[6][d] for v in tbl) print("| %d | %d | %.13f | +%.3e |" % (d, L, s / SC, (s - RECORD) / SC)) # --- (D) 混合深さ：利得/葉 の比で貪欲に配る --- print() print("=== 混合深さの設計（葉の予算ごと） ===") # 各窓について「d を 1 段深くしたときの利得 Δ と 追加の葉 ΔL」の階段を作る steps = [] # (gain_per_leaf, idx, d_from, d_to, gain, dleaf) base_d = 2 for (idx, p, n, r, a, vals, lvs) in tbl: cur = base_d while cur < DMAX and vals[cur + 1] is not None: g = vals[cur + 1] - vals[cur] dl = lvs[cur + 1] - lvs[cur] if dl <= 0: cur += 1; continue steps.append((g / dl, idx, cur, cur + 1, g, dl)) cur += 1 steps.sort(key=lambda x: -x[0]) for budget in (30000, 60000, 100000, 150000, 200000, 300000, 500000): depth = {i: base_d for i in range(len(wins))} L = sum(v[6][base_d] for v in tbl) gain = 0 for (rate, idx, df, dt, g, dl) in steps: if depth[idx] != df: continue # 階段は順に踏む if L + dl > budget: continue depth[idx] = dt; L += dl; gain += g s = tot2 + gain nd = {} for i in depth.values(): nd[i] = nd.get(i, 0) + 1 print("葉 ≤ %7d : 実際 %7d 葉 総和 %.13f 深さ分布 %s" % (budget, L, s / SC, sorted(nd.items()))) # --- (E) 目標 3.00853 に必要な最小の葉 --- print() TARGET = 300853 * 10 ** 95 need = TARGET - tot2 print("目標 3.00853 に必要な上げ幅 = %.4e" % (need / SC)) depth = {i: base_d for i in range(len(wins))} L = sum(v[6][base_d] for v in tbl); gain = 0 order = [] for (rate, idx, df, dt, g, dl) in steps: if depth[idx] != df: continue depth[idx] = dt; L += dl; gain += g order.append((L, gain)) if gain > need * 1.15: break print("目標到達に要る葉 = %d（総和 %.13f）" % (L, (tot2 + gain) / SC)) # --- (F) 利得の分布（上位 20 窓が占める割合） --- print() per = sorted(((tbl[i][5][3] - tbl[i][5][2]) if tbl[i][5][3] is not None else 0, i) for i in range(len(wins)))[::-1] tot_gain23 = sum(g for g, _ in per) top20 = sum(g for g, _ in per[:20]) print("d=2→3 の総利得 = %.4e、上位 20 窓 = %.4e（%.1f%%）" % (tot_gain23 / SC, top20 / SC, 100.0 * top20 / tot_gain23)) print("上位 10 窓 (p,n,r, a の桁数, 利得, d=2→3 の追加葉):") for g, i in per[:10]: (idx, p, n, r, a, vals, lvs) = tbl[i] print(" #%3d p=%2d n=%2d r=%3d a~1e%d 利得 %.3e 葉 %d→%d" % (i, p, n, r, len(str(a)) - 1, g / SC, lvs[2], lvs[3])) print() print("経過 %.1f 秒" % (time.time() - t0)) if __name__ == "__main__": main()