#!/usr/bin/env python3
"""N18 — arithmetic from D1287 supplied sources.
No network, external packages, fitting operations or shared runtime files.
Printed source cells and unrounded arithmetic are distinct outputs.
"""

import math
from decimal import Decimal, ROUND_HALF_UP, getcontext
from fractions import Fraction

getcontext().prec = 40
PHI = (1.0 + math.sqrt(5.0)) / 2.0


def shown(value, places=9, signed=False):
    """Decimal half-up display; never use display-rounded inputs implicitly."""
    value = Decimal(str(value))
    rounded = value.quantize(Decimal(1).scaleb(-places), rounding=ROUND_HALF_UP)
    return format(rounded, ("+" if signed else "") + "." + str(places) + "f")


def compare(label, theory, reference, lower=None, upper=None):
    """Naive central-value arithmetic, not a likelihood or theory-error model.

    The asymmetric denominator points from the comparator toward the theory:
    lower error below the central value; upper error above it. Any different
    printed convention is audited separately, rather than silently substituted.
    """
    theory, reference = Decimal(str(theory)), Decimal(str(reference))
    if reference == 0:
        raise ValueError("A relative deviation needs a nonzero reference.")
    offset = theory - reference
    percent = 100 * offset / reference
    print(label + ".theory = " + shown(theory))
    print(label + ".reference = " + shown(reference))
    print(label + ".signed_percent = " + shown(percent, signed=True)
          + "%; rounded = " + shown(percent, 2, True) + "%")
    if lower is None or upper is None:
        print(label + ".d = NOT AVAILABLE (no uncertainty supplied for this comparison variable)")
        return
    side = "lower" if offset < 0 else "upper"
    uncertainty = Decimal(str(lower if offset < 0 else upper))
    if uncertainty <= 0:
        raise ValueError("The selected comparator uncertainty must be positive.")
    distance = offset / uncertainty
    print(label + ".uncertainty_used = " + str(uncertainty) + " (" + side + ")")
    print(label + ".d = " + shown(distance, signed=True)
          + "; rounded = " + shown(distance, 2, True))


def lepton_chain():
    # The printed matching correction is an input, not re-derived here.
    tau = 1776.93
    ratio = (PHI / math.sqrt(5.0))**(8.0/3.0) * math.sqrt(2.0)/10.0
    tree = tau * ratio
    delta = 0.003125
    return tau, ratio, tree, delta, tree / (1.0 + delta)


def koide_roots(muon, tau):
    # Solve the supplied K=2/3 equation in x=sqrt(m_e), without a mass fit.
    s = math.sqrt(muon) + math.sqrt(tau)
    c = muon + tau - 4.0 * math.sqrt(muon*tau)
    disc = 4.0*s*s - c
    if disc < 0:
        raise ValueError("The supplied Koide inputs have no real roots.")
    lo, hi = 2.0*s - math.sqrt(disc), 2.0*s + math.sqrt(disc)
    if lo < 0:
        raise ValueError("The light square-root solution is not physical.")
    return lo*lo, hi*hi

def main():
    print('=== row 18 ===')
    print('Scorecard (verbatim): m e (MeV) QED-Koide, K = 2/3 theory 0.5076 experiment 0.51099895069 deviation −0.67% Koide-consistent ext. selector')
    tau, ratio, tree, delta, physical = lepton_chain()
    light, heavy = koide_roots(physical,tau)
    rounded_mu_light, unused = koide_roots(105.72,tau)
    rounded_tree_light, unused = koide_roots(106.05/(1+delta),tau)
    k = (light+physical+tau)/(math.sqrt(light)+math.sqrt(physical)+math.sqrt(tau))**2
    assert abs(k-2.0/3.0) < 1e-13
    print("tau_anchor_MeV = " + shown(tau,2))
    print("unrounded_mu_tree_MeV = " + shown(tree,12))
    print("supplied_delta_QED = " + shown(delta,6))
    print("unrounded_mu_physical_MeV = " + shown(physical,12))
    print("selected_light_root_MeV = " + shown(light,12))
    print("other_positive_root_MeV = " + shown(heavy,9))
    print("Koide_K_check = " + shown(k,12))
    print("theory_at_cell_precision_MeV = " + shown(light,4))
    print("light_root_from_rounded_mu_MeV = " + shown(rounded_mu_light,12)
          + "; rounded = " + shown(rounded_mu_light,4))
    print("light_root_from_rounded_tree_MeV = " + shown(rounded_tree_light,12)
          + "; rounded = " + shown(rounded_tree_light,4))
    compare("full_chain_MeV", light, "0.51099895069")
    compare("printed_cell_MeV", "0.5076", "0.51099895069")
    print()


if __name__ == "__main__":
    main()