#!/usr/bin/env python3 """Independent audit of the coefficient of (I4') (Appendix B of the paper), from Johnston-Starichkova's Theorem 4.2; run with python -B. Arb at 256 bits. The original chain is executed from a local, unmodified copy, with its sole exponent correction made in memory. The direct route below is a separate transcription, written after reading third/js_cert/js_cert.py, not a blind implementation. Adaptive box endpoints are exact rationals; every accepted box is retained. All output files are relative to this directory. No remote imports or writes. """ import json from pathlib import Path from flint import arb as A, fmpq, ctx from independent_check import constants, dump ctx.prec = 256 ctx.threads = 1 gamma = A.const_euler() delta = A('1.3') eps = A('0.0005543') rho = A('0.000105') eps1 = A(1) / 11 b = 2 * A('2.0452') * 7 L0 = fmpq(136843, 50) L1 = fmpq(3000) def old_chain(correct): p = Path('original_js_uniform.py').resolve() src = p.read_text() src = src.replace("save('js_uniform.json',out)", 'pass') src = src.replace("print(Path('js_uniform.json').read_text())", 'pass') if correct: assert src.count('L**(1+delta)') == 1 src = src.replace('L**(1+delta)', 'L**delta') scope = {'__file__': str(p)} exec(compile(src, str(p), 'exec'), scope) return scope['out'] def wfun(L): return L - 5 * L.log() def data(L): w = wfun(L) lw = w.log() nu = 1 / (b * lw) es = 100 * w**A('0.35') / (delta * lw)**2 ep = w * nu j = A('1.3841') * 7 * lw / (7 * lw).log() mult = 1 + 1 / (L**2 * w**3) + 1 / ((1 - 4 / w) * L) ec = (4 * w**A('4.5') / (w - 5 * lw)**7 + 4 / w**A('2.5') + 18 * w * (-w / 12).exp() + A('2.5') * w**A('5.5') * (-w / 6).exp()) small = (-w / 2).exp() / (2 * A(2).log() * w**3) + A('0.4') * w**5 * (-w).exp() add = (A('0.67') + 2 * (-w / 6).exp()) / w**5 deletion = A('0.9') * (w / 2 - L).exp() * L**4 / (w**7 * L.log()) c1 = (A('3.98') + w**2 * ((-ep).exp() / (1 - nu) + w**7 * (A('1.02') * (-w / 2).exp() + 3 * (-2 * w / 3).exp()) + A('9.4') * w**A('1.515') * (-A('0.8274') * w.sqrt()).exp())) c2 = c1 * mult + 1 / L**2 + A('1.3841') * L**4 * (-L).exp() / L.log() return w, lw, nu, es, ep, j, mult, ec, small, add, deletion, c2 def a_majorant(left, c2, tp): L = A(left) w = wfun(L) ell = L.log().log() qlog = (7 * w.log()).log() M3 = ((A('1.5') / L + (L + 7 * L.log()) / (L - 7 * L.log())) * (L / w)**delta * (gamma.exp() * (1 + delta.log() / ell) + A('2.5') / ((delta * w.log()).log() * ell)) + A('1.3841') * L**(2 + delta) * (-L).exp() / (L.log() * ell)) G = A('1.3841') * 7 * L.log() / (L**A('0.7') * ell * qlog) af = ell * (gamma.exp() * (A(7).log() + ell) + A('2.5') / qlog) / (tp * L**delta) return (c2 * G + M3) * af def direct(left, right, tp, C1=113): # Union contains the exact rational endpoints and therefore their interval. L = A(left).union(A(right)) w, lw, nu, es, ep, j, mult, ec, small, add, deletion, c2 = data(L) a = a_majorant(left, c2, tp) c = 40 * (A('0.25') - 14 * L.log() / L) assert c > 7 h = 3 * (-c).exp() / c E1, E2 = eps * C1 * A(2).exp() * h, eps * 114 * A(2).exp() * h brace = 1 - rho - eps1 * rho - (1 + eps1) * E2 - (3 * eps1 + a) * (rho + E1) - a - 1602 / L**2 assert brace > 0 T = 2 * (-gamma).exp() * tp def remainder(exponent, nu_used): p1 = A('1.1') * 7 * lw * (A('3.98') + w**2 * (-exponent).exp() / (1 - nu_used)) + 27 * w**2 * ec + small + add return (p1 * mult + A('2.2') / L**2 + deletion) / L # Case (a) uses min(nuP, nuS). Freeze 1/beta at its safe upper # bound 1/(1-nuP), then take the larger remainder from the two exponents. ra = remainder(ep, nu).max(remainder(es, nu)) rb = remainder(ep, nu) * j ca = 40 * (1 - 801 / L**2) * (1 - rho - eps * 114 * A(2).exp() * A('0.3') * A(-10).exp()) - ra / T cb = 40 * brace - rb / T return ca, cb, a, rb, brace def tail(left, tp, C1=113): L = A(left) w, lw, nu, es, ep, j, mult, ec, small, add, deletion, c2 = data(L) def normalized(exponent): return ((A('1.1') * 7 * lw * (A('3.98') / w + w * (-exponent).exp() / (1 - nu)) + 27 * w * ec + (small + add) / w) * mult + A('2.2') / L**3 + deletion / L) T = 2 * (-gamma).exp() * tp a = a_majorant(left, c2, tp) c = 40 * (A('0.25') - 14 * L.log() / L) h = 3 * (-c).exp() / c E1, E2 = eps * C1 * A(2).exp() * h, eps * 114 * A(2).exp() * h brace = 1 - rho - eps1 * rho - (1 + eps1) * E2 - (3 * eps1 + a) * (rho + E1) - a - 1602 / L**2 assert brace > 0 cb = 40 * brace - normalized(ep) * j / T ca = 40 * (1 - 801 / L**2) * (1 - rho - eps * 114 * A(2).exp() * A('0.3') * A(-10).exp()) - normalized(ep).max(normalized(es)) / T return {'case_a': ca, 'case_b': cb, 'a': a, 'brace': brace} def cover(tp, C1=113): pending, accepted = [(L0, L1)], [] while pending: left, right = pending.pop() ca, cb, a, rem, brace = direct(left, right, tp, C1) if ca > A('38.4') and cb > A('9.8036'): accepted.append({'left': str(left), 'right': str(right), 'case_a': ca, 'case_b': cb, 'case_a_lower': ca.lower(), 'case_b_lower': cb.lower(), 'case_b_upper': cb.upper(), 'a': a, 'remainder': rem, 'brace': brace}) else: assert right - left > fmpq(1, 10**8) mid = (left + right) / 2 pending.extend([(mid, right), (left, mid)]) assert fmpq(accepted[0]['left']) == L0 assert fmpq(accepted[-1]['right']) == L1 assert all(a['right'] == b['left'] for a, b in zip(accepted, accepted[1:])) worst = min(accepted, key=lambda row: row['case_b'].lower()) return {'boxes': accepted, 'count': len(accepted), 'worst_case_b_box': worst, 'min_lower_a': min(row['case_a_lower'] for row in accepted), 'min_lower_b': min(row['case_b_lower'] for row in accepted)} def rounded_chain(): L = A(L0) a, h = A('0.00379644'), A('0.000082075421') E = eps * 114 * A(2).exp() * h brace = 1 - rho - eps1 * rho - (1 + eps1) * E - (3 * eps1 + a) * (rho + E) - a - 1602 / L**2 cb = 40 * brace - A('22.582845') / A('0.7413082243') assert cb > A('9.3682') return {'brace': brace, 'case_b': cb} if __name__ == '__main__': const = constants() assert const['twin_constant'] > A('0.6601618158') assert const['U_low'] > A('0.7413082243') printed, corrected = old_chain(False), old_chain(True) assert printed['case_b_uniform_coefficient'] > A('9.52') assert corrected['case_b_uniform_coefficient'] > A('9.3682') assert corrected['case_b_uniform_coefficient'] < A('9.52') # Conservative decimal below the finite-product lower bound of third/js_cert/tp_bound.py. tp = A('0.6601617536') assert const['twin_constant'] > tp result = {'precision_bits': ctx.prec, 'constants': const, 'original_printed_chain': printed, 'corrected_original_chain': corrected, 'corrected_rounded_chain': rounded_chain(), 'direct_tp_lower_bound': tp, 'direct_endpoint_C1_113': direct(L0, L0, tp), 'direct_cover_C1_113': cover(tp), 'direct_cover_C1_114': cover(tp, 114), 'tail_at_3000_C1_113': tail(L1, tp), 'tail_at_3000_C1_114': tail(L1, tp, 114), 'route_A_at_L0_C1_113': tail(L0, tp)} for C1 in (113, 114): assert result[f'tail_at_3000_C1_{C1}']['case_a'] > A('39.28') assert result[f'tail_at_3000_C1_{C1}']['case_b'] > A('26.21') assert result[f'direct_cover_C1_{C1}']['min_lower_a'] > A('38.417') Path('audit_js.json').write_text(json.dumps(dump(result), indent=2) + '\n') summary = {k: v for k, v in result.items() if not k.startswith('direct_cover')} for C1 in (113, 114): cov = result[f'direct_cover_C1_{C1}'] summary[f'cover_C1_{C1}_count'] = cov['count'] summary[f'cover_C1_{C1}_minimum_a'] = cov['min_lower_a'] summary[f'cover_C1_{C1}_minimum_b'] = cov['min_lower_b'] summary[f'cover_C1_{C1}_worst'] = cov['worst_case_b_box'] print(json.dumps(dump(summary), indent=2))