#!/usr/bin/env python3
"""N17 — 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 main():
    print('=== row 17a ===')
    print('Scorecard (verbatim): m μ tree (MeV) (φ/√5) 8/3 ·√2/10·m τ theory 106.05 experiment 105.6583755 deviation +0.37% Structural / loaded fit gap')
    tau, ratio, tree, delta, physical = lepton_chain()
    print("tau_anchor_MeV = " + shown(tau,2))
    print("mass_ratio = " + shown(ratio,12))
    print("exponent = 8/3; normalization = sqrt(2)/10")
    print("normalization_value = " + shown(math.sqrt(2.0)/10.0,9))
    print("theory_at_cell_precision_MeV = " + shown(tree,2))
    compare("full_tree_MeV", tree, "105.6583755")
    compare("printed_tree_MeV", "106.05", "105.6583755")
    print()
    print('=== row 17b ===')
    print('Scorecard (verbatim): m μ phys (MeV) tree / (1+δ QED ) theory 105.72 experiment 105.6583755 deviation +0.06% same —')
    tau, ratio, tree, delta, physical = lepton_chain()
    print("supplied_delta_QED = " + shown(delta,6) + "; percent = " + shown(100*delta,4)+"%")
    print("full_tree_input_MeV = " + shown(tree,12))
    print("theory_at_cell_precision_MeV = " + shown(physical,2))
    print("alternative_from_rounded_tree_MeV = " + shown(106.05/(1+delta),12))
    compare("full_matching_chain_MeV", physical, "105.6583755")
    compare("printed_physical_cell_MeV", "105.72", "105.6583755")
    print("source_finer_percent_display = +0.058%; rounded-cell rather than full-chain convention")
    print()


if __name__ == "__main__":
    main()