#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""第110便 — f(4) の数。Walker の K(S,55)+1 を尺度 Z で切った頭に、Wróblewski ブロック
B(p,q,r)（3-AP-free）を倍加税 l_j = 2 r_{j-1} + 1 だけで継ぐ（窓長は自由：第107便 補題 2）。
逆数和は整数厳密（固定小数 SC = 10^100、方向付き丸め）で挟む。lo の式は erdos834.loI（Lean）と同一。

  python3 erdos1110f4.py bench            費用の見積り
  python3 erdos1110f4.py head             Walker の H(K+1) の再現と、切り所ごとの頭 [lo,hi]
  python3 erdos1110f4.py k3 [cap]         較正：k=3（Szekeres 頭・p≤14・110 段）で 3.008532672324 を再現
  python3 erdos1110f4.py scan             切り所 Z を浮動小数の見積りで走査
  python3 erdos1110f4.py k4 log10Z [cap]  本番（厳密）
"""
import sys, time, math, re
from math import comb
from fractions import Fraction
import os
# [配布版の変更点 / changed for distribution] 元は書き手の機体の絶対パスを指していた箇所を、このディレクトリからの相対に変えた。計算の中身は無改変。
_HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, _HERE)
from erdos834 import T, mxv, cntTab, loI, SC
import numpy as np

def rlo(n): return SC // n
def rhi(n): return -((-SC) // n)
def show(x, d=12):
    q, rr = divmod(x, SC); return f"{q}.{(rr*10**d)//SC:0{d}d}"

# ---------------- ブロック表（Kronecker 詰めの厳密 DP） ----------------
class Rows:
    """cntTab(p, n), n=0..q を一本の整数に詰めて持ち、要る行だけ展開する。"""
    def __init__(self, p, q):
        self.p, self.q = p, q
        Tk = [T(p, k) for k in range(p)]
        hist = {}
        for t in Tk: hist[t] = hist.get(t, 0) + 1
        self.hist = hist; self.Tmax = max(Tk)
        self.Wb = (p ** q).bit_length() // 8 + 2      # スロット幅（バイト）。p^q < 2^(8Wb-8)
        W = 8 * self.Wb
        self.packed = [1]; self.lens = [1]
        cur = 1
        for n in range(1, q + 1):
            nxt = 0
            for t, m in hist.items(): nxt += (cur * m) << (t * W)
            cur = nxt
            self.packed.append(cur); self.lens.append(self.lens[-1] + self.Tmax)
        Ts = sorted(hist); self.supp = [1]
        for n in range(1, q + 1):
            s = 0
            for t in Ts: s |= self.supp[-1] << t
            self.supp.append(s)
        self._rows = {}; self._lv = {}
    def row(self, n):
        if n not in self._rows:
            L = self.lens[n]; Wb = self.Wb
            b = self.packed[n].to_bytes(L * Wb, 'little')
            r = [int.from_bytes(b[i*Wb:(i+1)*Wb], 'little') for i in range(L)]
            while len(r) > 1 and r[-1] == 0: r.pop()
            self._rows[n] = tuple(r)
        return self._rows[n]
    def cf(self, n, r):
        row = self.row(n); return row[r] if 0 <= r < len(row) else 0
    def ok(self, n, r): return r >= 0 and (self.supp[n] >> r) & 1
    def argmax(self, n=None):
        row = self.row(self.q if n is None else n); return max(range(len(row)), key=lambda i: row[i])
    def maxelt(self, r, n=None):
        p, base = self.p, 2 * self.p - 1; v = 0
        for m in range((self.q if n is None else n), 0, -1):
            for k in range(p - 1, -1, -1):
                t = T(p, k)
                if t <= r and self.ok(m - 1, r - t):
                    v = v * base + k; r -= t; break
            else: raise ValueError("unreachable")
        assert r == 0
        return v
    def leaves(self, d, n, r):
        if d == 0 or n == 0: return 1
        key = (d, n, r)
        if key in self._lv: return self._lv[key]
        c = 0
        for k in range(self.p):
            t = T(self.p, k)
            if t <= r and self.ok(n - 1, r - t): c += self.leaves(d - 1, n - 1, r - t)
        self._lv[key] = c; return c
    def bracket(self, d, a, n, r):
        """Σ_{s∈B(p,n,r)} 1/(a+s) の [lo,hi]（SC 倍整数）。lo は erdos834.loI と同式。"""
        if d == 0 or n == 0:
            c = self.cf(n, r)
            if c == 0: return 0, 0
            return c * rlo(a + mxv(self.p, n)), c * rhi(a)
        lo = hi = 0; base = (2 * self.p - 1) ** (n - 1)
        for k in range(self.p):
            t = T(self.p, k)
            if t <= r and self.ok(n - 1, r - t):
                l, h = self.bracket(d - 1, a + k * base, n - 1, r - t); lo += l; hi += h
        return lo, hi
    def depth_for(self, r, cap, n=None):
        n = self.q if n is None else n; d = 0
        while d < n and self.leaves(d + 1, n, r) <= cap: d += 1
        return d

def float_table(pmin, pmax, qmax, lim, qmin=4):
    """政策用の浮動小数の候補表。(p,q) -> (r0, |B(p,q,r0)| の近似)。"""
    tab = {}
    for p in range(pmin, pmax + 1):
        Tk = [T(p, k) for k in range(p)]
        hist = {}
        for t in Tk: hist[t] = hist.get(t, 0) + 1
        Tmax = max(Tk); row = np.array([1.0])
        for n in range(1, qmax + 1):
            if (2 * p - 1) ** n > lim: break
            nxt = np.zeros(len(row) + Tmax)
            for t, m in hist.items(): nxt[t:t + len(row)] += m * row
            row = nxt
            if n >= qmin:
                r0 = int(row.argmax()); tab[(p, n)] = (r0, float(row[r0]))
    return tab

class Supp:
    """台だけ（maxelt 用）。"""
    def __init__(self, p, qmax):
        self.p = p; Ts = sorted(set(T(p, k) for k in range(p))); self.supp = [1]
        for n in range(1, qmax + 1):
            s = 0
            for t in Ts: s |= self.supp[-1] << t
            self.supp.append(s)
    def ok(self, n, r): return r >= 0 and (self.supp[n] >> r) & 1
    def maxelt(self, q, r):
        p, base = self.p, 2 * self.p - 1; v = 0
        for m in range(q, 0, -1):
            for k in range(p - 1, -1, -1):
                t = T(p, k)
                if t <= r and self.ok(m - 1, r - t):
                    v = v * base + k; r -= t; break
            else: raise ValueError("unreachable")
        assert r == 0
        return v

# ---------------- Walker の Kempner 集合 K(S,b)+1 ----------------
S55 = (0,1,2,4,5,9,10,11,14,16,17,18,21,24,30,37,39,41,42,45,47)
class Kempner:
    def __init__(self, S, b, mmax, I=14):
        assert I % 2 == 0
        self.S, self.b, self.I = tuple(S), b, I; self.Sset = set(S)
        D = [sum(d ** j for d in S) for j in range(I + 1)]
        M = [[1] + [0] * I]
        for m in range(1, mmax + 1):
            prev = M[-1]
            M.append([sum(comb(i, j) * b ** j * prev[j] * D[i - j] for j in range(i + 1)) for i in range(I + 1)])
        self.M = M; dmax = max(S)
        self.R = [dmax * (b ** m - 1) // (b - 1) for m in range(mmax + 1)]
        # V_i(m) = Σ_{s∈S^m} (2s − R_m)^i
        self.V = []
        for m in range(mmax + 1):
            R = self.R[m]
            self.V.append([sum(comb(i, j) * (2 ** j) * M[m][j] * ((-R) ** (i - j)) for j in range(i + 1)) for i in range(I + 1)])
            assert self.V[m][I] >= 0
    def series(self, x0, m):
        """Σ_{s∈S^m} 1/(x0+s) の [lo,hi]（SC 倍整数）。展開点は区間の中央、I 項＋剰余の符号付き評価。"""
        I = self.I; R = self.R[m]; X = 2 * x0 + R; V = self.V[m]
        lo = hi = 0
        for i in range(I):
            N = 2 * (-1) ** i * V[i] * SC; den = X ** (i + 1)
            lo += N // den; hi += -((-N) // den)
        hi += -((-2 * V[I] * SC) // (X ** I * (X - R)))     # 0 ≤ 剰余 ≤ 2V_I/(X^I (X−R))
        return lo, hi
    def free_block(self, P, m):
        """Σ_{s∈S^m} 1/(1 + P·b^m + s)。P ≥ 1。P < b なら次の桁で割ってから級数。"""
        if m == 0: return rlo(1 + P), rhi(1 + P)
        if P < self.b:
            lo = hi = 0
            for e in self.S:
                l, h = self.free_block(P * self.b + e, m - 1); lo += l; hi += h
            return lo, hi
        return self.series(1 + P * self.b ** m, m)
    def head(self, Zk):
        """Σ_{k∈K, k≤Zk} 1/(k+1) の [lo,hi]（SC 倍整数）と rmax = 1 + max{k∈K: k≤Zk}。"""
        b = self.b; digs = []; v = Zk
        while v > 0: digs.append(v % b); v //= b
        digs.reverse(); N = len(digs)
        lo, hi = rlo(1), rhi(1)
        for n in range(1, N):
            for e in self.S:
                if e == 0: continue
                l, h = self.free_block(e, n - 1); lo += l; hi += h
        Pt = 0; tight = True
        for i in range(N):
            m = N - 1 - i
            for e in self.S:
                if e < digs[i] and not (i == 0 and e == 0):
                    l, h = self.free_block(Pt * b + e, m); lo += l; hi += h
            if digs[i] in self.Sset:
                Pt = Pt * b + digs[i]
            else:
                tight = False; break
        if tight: lo += rlo(Zk + 1); hi += rhi(Zk + 1)
        Pt = 0; kmax = None
        for i in range(N):
            m = N - 1 - i; e = max(x for x in self.S if x <= digs[i])
            if e < digs[i]: kmax = (Pt * b + e) * b ** m + self.R[m]; break
            Pt = Pt * b + e
        if kmax is None: kmax = Zk
        return lo, hi, kmax + 1
    def full(self, N):
        """H(K+1) の [lo,hi]：N 桁までを厳密に、N+1 桁以上は 20ρ^N/(1−ρ) で抑える（|S|−1 = 20 は先頭桁）。"""
        lo, hi, _ = self.head(self.b ** N - 1)
        rho = Fraction(len(self.S), self.b); lead = len(self.S) - 1
        tail = Fraction(lead) * rho ** N / (1 - rho)
        hi += -((-tail.numerator * SC) // tail.denominator)
        return lo, hi

def selftest():
    t0 = time.time()
    # (a) Rows vs erdos834.cntTab
    for (p, q) in ((4, 9), (6, 10), (9, 7), (14, 6), (15, 5)):
        R = Rows(p, q)
        for n in range(q + 1):
            assert R.row(n) == cntTab(p, n), (p, q, n)
            assert all((R.ok(n, r) != 0) == (R.row(n)[r] > 0) for r in range(len(R.row(n))))
    # (b) maxelt vs erdos775-maxelt.txt
    src = open(os.path.join(_HERE, 'erdos775-maxelt.txt')).read()
    ents = re.findall(r'\((\d+), (\d+), (\d+), (\d+),', src)
    chk = 0
    for (p, q, r, mx) in ents[:60]:
        p, q, r, mx = map(int, (p, q, r, mx))
        R = Rows(p, q); assert R.maxelt(r) == mx, (p, q, r, mx, R.maxelt(r)); chk += 1
    # (c) maxelt vs 全列挙（p=4,q=9）
    p, q = 4, 9; base = 2 * p - 1; R = Rows(p, q)
    best = {}
    for v in range(base ** q):
        x = v; w = 0; okd = True
        for _ in range(q):
            d = x % base; x //= base
            if d >= p: okd = False; break
            w += T(p, d)
        if okd: best[w] = max(best.get(w, -1), v)
    for w, mx in best.items(): assert R.maxelt(w) == mx
    # (d) bracket.lo == loI
    for (p, q, r, a, d) in ((4, 9, 4, 43046723, 3), (6, 10, 13, 10 ** 9 + 7, 2), (12, 8, 44, 10 ** 12 + 1, 2)):
        R = Rows(p, q); lo, hi = R.bracket(d, a, q, r)
        assert lo == loI(p, d, a, q, r); assert lo <= hi
    # (e) Kempner の級数 vs 全列挙（Fraction）
    K = Kempner(S55, 55, 6)
    for (P, m) in ((56, 3), (1, 3), (1234, 2)):
        lo, hi = K.free_block(P, m)
        ex = Fraction(0)
        for v in range(21 ** m):
            x = v; s = 0
            for i in range(m):
                s += S55[x % 21] * 55 ** i; x //= 21
            ex += Fraction(1, 1 + P * 55 ** m + s)
        assert Fraction(lo, SC) <= ex <= Fraction(hi, SC), (P, m)
        assert hi - lo < 10 ** 75, (P, m, hi - lo)
    # (f) head(Zk) vs 全列挙（小さい Zk）
    Zk = 55 ** 2 + 1234
    lo, hi, rmax = K.head(Zk)
    ex = Fraction(0); mx = 0
    for v in range(21 ** 3):
        x = v; s = 0
        for i in range(3):
            s += S55[x % 21] * 55 ** i; x //= 21
        if s <= Zk: ex += Fraction(1, 1 + s); mx = max(mx, s)
    assert Fraction(lo, SC) <= ex <= Fraction(hi, SC) and rmax == mx + 1
    print(f"[selftest] OK（maxelt 照合 {chk} 件・全列挙 4 種）{time.time()-t0:.1f} 秒")

# ---------------- k=4 の鎖 ----------------
RMIN, RMAX = 0.01, 1000.0
def pick4(z, cand):
    best = None
    for (p, q, r, t, n) in cand:
        if not (RMIN * z <= t <= RMAX * z): continue
        a = 2 * z + 1
        g = (n / t) * math.log1p(t / a); dl = math.log((a + t) / z)
        if best is None or g / dl > best[0]: best = (g / dl, p, q, r, t, n, g)
    return best

def chain4(r0, cand, nst, exact, cap, supp=None, log=None):
    z = r0; lo = hi = 0; est = 0.0; stages = []
    for j in range(nst):
        b = pick4(z, cand)
        if b is None: break
        _, p, q, r, t, n, g = b
        a = 2 * z + 1
        d = None; l = h = 0
        if exact:
            t1 = time.time()
            R = Rows(p, q); r = R.argmax(); t = R.maxelt(r); n = R.cf(q, r)
            d = R.depth_for(r, cap); l, h = R.bracket(d, a, q, r)
            lo += l; hi += h
            g = (n / t) * math.log1p(t / a)
            if log: log(f"  段{j+1:2d} B({p},{q},{r}) |B|={n:.3e} t=1e{math.log10(t):.2f} a=1e{math.log10(a):.2f} d={d} "
                        f"[{l/SC:.4e}, {h/SC:.4e}] 幅 {(h-l)/SC:.1e}  {time.time()-t1:.1f}s")
        assert a > 2 * z and t > 0
        est += g
        stages.append(dict(p=p, q=q, r=r, t=t, a=a, n=n, lo=l, hi=h, d=d, zprev=z, g=g))
        z = a + t
    return lo, hi, est, stages, z

def build_cand4(pmax=200, qmax=26, lim=10 ** 82, log=print):
    t0 = time.time(); tab = float_table(3, pmax, qmax, lim)
    log(f"[候補表] 浮動小数 {len(tab)} 対 (p≤{pmax}, q≤{qmax}, (2p−1)^q≤1e{round(math.log10(lim))}) {time.time()-t0:.1f} 秒")
    t0 = time.time(); cand = []; sp = {}
    for (p, q), (r0, n0) in sorted(tab.items()):
        if p not in sp: sp[p] = Supp(p, qmax)
        cand.append((p, q, r0, sp[p].maxelt(q, r0), n0))
    log(f"[候補表] maxelt（台の bitset で厳密） {time.time()-t0:.1f} 秒")
    return cand

def main():
    mode = sys.argv[1] if len(sys.argv) > 1 else 'bench'
    T0 = time.time()
    if mode == 'bench':
        selftest()
        for (p, q) in ((115, 15), (200, 20), (83, 15)):
            t0 = time.time(); R = Rows(p, q); t1 = time.time(); r0 = R.argmax(); t2 = time.time()
            mx = R.maxelt(r0); d = R.depth_for(r0, 30000); lv = R.leaves(d, q, r0); t3 = time.time()
            lo, hi = R.bracket(d, 10 ** 33, q, r0); t4 = time.time()
            print(f"[bench] Rows({p},{q}) 詰め {t1-t0:.2f}s 展開+argmax {t2-t1:.2f}s maxelt+depth {t3-t2:.2f}s "
                  f"bracket(d={d}, 葉 {lv}) {t4-t3:.2f}s  |B|={R.cf(q,r0):.3e} 幅/値={(hi-lo)/max(lo,1):.1e}")
        t0 = time.time(); tab = float_table(3, 60, 22, 10 ** 75); print(f"[bench] float_table p≤60: {time.time()-t0:.1f}s ({len(tab)} 対)")
        t0 = time.time(); tab = float_table(150, 200, 22, 10 ** 75); print(f"[bench] float_table 150≤p≤200: {time.time()-t0:.1f}s ({len(tab)} 対)")
        t0 = time.time(); K = Kempner(S55, 55, 40); print(f"[bench] Kempner 冪和 m≤40: {time.time()-t0:.1f}s")
        t0 = time.time(); lo, hi, r0 = K.head(10 ** 32); print(f"[bench] head(1e32): {time.time()-t0:.1f}s  [{show(lo,14)}, {show(hi,14)}] rmax=1e{math.log10(r0):.3f}")
    elif mode == 'head':
        K = Kempner(S55, 55, 90)
        lo, hi = K.full(85)
        print(f"[H(K+1)] 85 桁まで厳密＋尾の上界: [{show(lo,15)}, {show(hi,15)}]  幅 {(hi-lo)/SC:.1e}")
        print(f"         Walker の 4.43975 との差: lo−4.43975 = {(lo - 443975*SC//100000)/SC:+.3e}")
        HH = hi
        print()
        print("| 切り所 Zk | rmax | 頭 [lo, hi] | 捨てた尾 H−頭 |")
        for e in [28, 29, 30, 30.5, 31, 31.5, 32, 32.5, 33, 33.5, 34, 35, 36]:
            Zk = int(10 ** e)
            l, h, r0 = K.head(Zk)
            print(f"| 1e{e} | 1e{math.log10(r0):.3f} | [{show(l,14)}, {show(h,14)}] | {(lo-h)/SC:.3e} .. {(HH-l)/SC:.3e} |")
        print(f"経過 {time.time()-T0:.1f} 秒")
    elif mode == 'k3':
        cap = int(sys.argv[2]) if len(sys.argv) > 2 else 250000
        run3(cap)
    elif mode == 'scan':
        K = Kempner(S55, 55, 60); Hlo, Hhi = K.full(60)
        cand = build_cand4()
        print(f"H(K+1) ∈ [{show(Hlo,12)}, {show(Hhi,12)}]")
        print("| log10 Z | rmax | 頭 lo | 尾（見積り） | 総和（見積り） | H との差 | 段数 |")
        for e in [29, 30, 30.5, 31, 31.5, 32, 32.5, 33, 33.5, 34, 35]:
            Zk = int(10 ** e); hl, hh, r0 = K.head(Zk)
            lo, hi, est, st, z = chain4(r0, cand, 60, False, 0)
            tot = hl / SC + est
            print(f"| {e} | 1e{math.log10(r0):.3f} | {show(hl,10)} | {est:.4e} | {tot:.10f} | {tot - Hhi/SC:+.3e} | {len(st)} (z=1e{math.log10(z):.1f}) |")
        print(f"経過 {time.time()-T0:.1f} 秒")
    elif mode == 'k4':
        e = float(sys.argv[2]); cap = int(sys.argv[3]) if len(sys.argv) > 3 else 200000
        nst = int(sys.argv[4]) if len(sys.argv) > 4 else 60
        selftest()
        K = Kempner(S55, 55, 90); Hlo, Hhi = K.full(85)
        print(f"[H(K+1)] [{show(Hlo,15)}, {show(Hhi,15)}]")
        cand = build_cand4()
        Zk = int(10 ** e); hl, hh, r0 = K.head(Zk)
        print(f"[頭] Zk=1e{e}  rmax = {r0} (1e{math.log10(r0):.4f})  頭 ∈ [{show(hl,15)}, {show(hh,15)}]  捨てた尾 ≤ {(Hhi-hl)/SC:.4e}")
        lo, hi, est, st, z = chain4(r0, cand, nst, True, cap, log=print)
        TL, TH = hl + lo, hh + hi
        print()
        print(f"[段] {len(st)} 段、最終 z = 1e{math.log10(z):.2f}、尾の厳密和 ∈ [{lo/SC:.6e}, {hi/SC:.6e}]（見積り {est:.4e}）")
        print(f"[倍加税] 全段 a_j = 2 r_(j−1) + 1 > 2 r_(j−1): {all(s['a'] == 2*s['zprev']+1 for s in st)}")
        print(f"[総和] 厳密な下界 = {show(TL,15)}")
        print(f"[総和] 厳密な上界 = {show(TH,15)}")
        print(f"[総和] 幅 = {(TH-TL)/SC:.3e}")
        print(f"[比較] 下界 − H(K+1) の上界 = {(TL-Hhi)/SC:+.4e}   下界 − 4.43975 = {(TL-443975*SC//100000)/SC:+.4e}")
        print(f"[比較] 超えた: {TL > Hhi}")
        with open('erdos1110f4-stages-%s.txt' % sys.argv[2], 'w') as f:
            f.write(f"# Zk=1e{e} rmax={r0} head_lo={hl} head_hi={hh} H_lo={Hlo} H_hi={Hhi}\n")
            f.write("# j p q r t a n d lo hi\n")
            for j, s in enumerate(st, 1):
                f.write(f"{j} {s['p']} {s['q']} {s['r']} {s['t']} {s['a']} {s['n']} {s['d']} {s['lo']} {s['hi']}\n")
            f.write(f"# TOTAL_LO={TL}\n# TOTAL_HI={TH}\n")
        print(f"経過 {time.time()-T0:.1f} 秒")

# ---------------- k=3 較正（erdos732b/733b の政策をそのまま） ----------------
def sep_ok(iv):
    n = len(iv)
    for i in range(n):
        for j in range(n):
            for k in range(n):
                if i == j and j == k: continue
                L = iv[i][0] + iv[k][0]; Rr = iv[i][1] + iv[k][1]
                a, b = 2 * iv[j][0], 2 * iv[j][1]
                if not (Rr < a or b < L): return False
    return True
def offsets_M1(z, t): return [max(2 * z, z + t) + 1]
def offsets_M2(z, t):
    a = offsets_M1(z, t)[0]; lowb = max(2 * z, a + 2 * t) + 1; c = []
    if lowb + t + z < 2 * a: c.append([a, lowb])
    c.append([a, 2 * a + 2 * t + 1])
    def sc(o):
        num = den = 1
        for x in o: num *= (x + t); den *= x
        return Fraction(num, den)
    return max(c, key=sc)
def offsets_M3(z, t):
    m = max(z, t); o = [m + z + 1, 3 * m + 2 * t + z + 3, 3 * m + 4 * t + z + 4]
    return o if sep_ok([(0, z)] + [(a, a + t) for a in o]) else None
def offsets(M, z, t): return {1: offsets_M1, 2: offsets_M2}[M](z, t) if M in (1, 2) else offsets_M3(z, t)

def run3(jmax):
    t0 = time.time(); Kd = 16; head = []
    for mask in range(1 << Kd):
        v = 1
        for i in range(Kd):
            if (mask >> i) & 1: v += 3 ** i
        head.append(v)
    z = max(head); lo = sum(rlo(n) for n in head); hi = sum(rhi(n) for n in head)
    print(f"[k3 頭] |Z0|={len(head)} max={z} Σ ∈ [{show(lo,15)}, {show(hi,15)}]")
    cand = []; RW = {}
    for p in range(3, 15):
        qmax = 4
        while (2 * p - 1) ** (qmax + 1) <= 10 ** 100 and qmax < 49: qmax += 1
        R = Rows(p, qmax); RW[p] = R
        for q in range(4, qmax + 1):
            row = R.row(q); r0 = max(range(len(row)), key=lambda i: row[i])
            for r in (r0 - 1, r0, r0 + 1):
                if 0 <= r < len(row) and row[r] > 0: cand.append((p, q, r, R.maxelt(r, q), row[r]))
    print(f"[k3 候補] {len(cand)} 個 {time.time()-t0:.1f} 秒")
    Ms = {}
    for step in range(110):
        best = None
        for (p, q, r, t, n) in cand:
            if not (0.05 * z <= t <= 40 * z): continue
            for M in (1, 2, 3):
                o = offsets(M, z, t)
                if o is None or not sep_ok([(0, z)] + [(a, a + t) for a in o]): continue
                g = (n / t) * sum(math.log(1 + t / a) for a in o); dl = math.log((o[-1] + t) / z)
                if best is None or g / dl > best[0]: best = (g / dl, (p, q, r, M), o, t)
        (_, (p, q, r, M), o, t) = best
        assert sep_ok([(0, z)] + [(a, a + t) for a in o])
        R = RW[p]; j = 1
        while j < q and p ** (j + 1) <= jmax: j += 1
        for a in o:
            l, h = R.bracket(j, a, q, r); lo += l; hi += h
        z = o[-1] + t; Ms[M] = Ms.get(M, 0) + 1
        if step < 6 or step % 20 == 19:
            print(f"   段{step+1:3d} B({p},{q},{r}) M={M} 累計 [{show(lo,13)}, {show(hi,13)}] {time.time()-t0:.0f}s")
    print(f"[k3 結果] 110 段 M の内訳 {sorted(Ms.items())}  最終 z = 1e{math.log10(z):.1f}")
    print(f"[k3 結果] 厳密な下界 {show(lo,13)}  上界 {show(hi,13)}  幅 {(hi-lo)/SC:.2e}")
    ref = 3008532672324 * (SC // 10 ** 12)
    print(f"[k3 較正] ノートの 110 段の値 3.008532672324 との差: lo−ref = {(lo-ref)/SC:+.3e}, hi−ref = {(hi-ref)/SC:+.3e}")
    print(f"経過 {time.time()-t0:.1f} 秒")

if __name__ == '__main__':
    main()