# SPDX-License-Identifier: MIT

# physics — check quantities, units and dimensions, and calculate classical
# motion.
#
# Two halves that share nothing but the subject.
#
# CHECKING. You supply quantities as facts (`physics.Quantity`), each a value
# in a unit, and the units they are written in (`physics.Unit`), each with its
# dimension as powers of mass, length, time, electric current and temperature,
# and what one of it is in SI. `si_units()` writes a table of common units,
# each with the source of its value:
#
#     import physics
#     physics.si_units()
#     physics.quantity("impulse", 100.0, "lbf_s")
#     physics.expect("impulse", "N_s")
#     print(physics.checks())
#
# The standing rules then conclude each quantity's value in SI, whether an
# equation written as a product of quantities is dimensionally consistent and
# what is left over when it is not, a quantity handed over in a unit other
# than the one its reader expects, and the resistance of resistors joined in
# series or in parallel. Each conclusion answers `why`, and withdraws itself
# when a fact it rests on is corrected. Where two sources disagree about a
# fact, both values are reasoned about, and `unsettled()` names what they
# disagree about.
#
# An equation's dimensions and a group's resistance are sums, and a sum adds
# every value it reads: a resistor whose ohms a second source writes, even the
# same ohms, would be counted twice, 400 ohms where every source says 300. So
# while a member is bound to its sum twice, as a second value binds it, that
# sum is disputed and not concluded, and neither is what rests on it: a factor
# whose power two sources give leaves each dimension it carries undecided, and
# `dimension()` says "disputed"; a resistor of two values leaves its group's
# resistance undecided, and `resistance()` refuses it. `checks()` names each
# disputed equation and group, and `why_not` on the sum names the dispute. A
# sum both sources agree on stands. Forgetting the claim that is not to stand
# settles it. A factor whose equation, or a resistor whose group, two sources
# disagree on counts in each a value names, so each of those sums is one
# source's. No other query refuses.
#
# CALCULATION. Value scripts that return a value and write nothing: a
# pendulum's period by integration and exactly, an orbit under three
# integrators, the heat equation, a linear system by Gaussian elimination and
# the bridge circuit it solves, and the Ising model of a magnet.

import math

# ── Quantities as facts ──────────────────────────────────────────────────────

class Unit:
    # A unit, by a name of letters, digits and underscores: "N", "lbf_s" for
    # a pound-force second, "m_per_s2" for a metre per second squared. Its
    # dimension is the powers of the five base dimensions it is made of, so a
    # newton, a kilogram metre per second squared, is mass 1, length 1 and
    # time -2. `factor` is what one of it is in the SI unit of the same
    # dimension: 4.4482216152605 for a pound-force, in newtons. A unit whose
    # zero is not SI's, such as a degree Celsius, is not a factor and is not
    # written here.
    mass: int
    length: int
    time: int
    current: int
    temperature: int
    factor: real

class Quantity:
    # A value in a unit, by a name you choose. `expects` is the unit its
    # reader takes it in, where it is handed from one system to another. The
    # rules conclude `factor` and the five powers from the unit, and
    # `efactor` and the expected unit's powers, `emass` to `etemperature`,
    # from `expects`. `converted` takes the ratio of the two factors first, so
    # a value expected in its own unit comes back as it was, and it is a real
    # wherever the value in the expected unit is one, for any factors from
    # 10^-150 to 10^150.
    value: real
    unit: string
    expects: string
    factor: real
    mass: int
    length: int
    time: int
    current: int
    temperature: int
    efactor: real
    emass: int
    elength: int
    etime: int
    ecurrent: int
    etemperature: int
    si: real = value * factor
    converted: real = value * (factor / efactor)
    # What the rules conclude about the quantity: a unit the table does not
    # hold, and how a hand-off stands, "convert" when the expected unit has
    # the same dimension and "incompatible" when it has another.
    unknown: bool
    unexpected: bool
    handoff: string

class Factor:
    # One quantity in an equation written as a product of powers: the
    # equation, the quantity and its power. Write every power on the left side
    # positive and every power on the right side negative, so `F = m a` is
    # F to the 1, m to the -1 and a to the -1. A pure number such as the 1/2
    # of `E = m v^2 / 2` has no dimension and needs no factor. A square root
    # is written by doubling every power in the equation. The rules conclude
    # the quantity's five powers, and the computed fields are this factor's
    # share of each.
    equation: string
    quantity: string
    power: int
    mass: int
    length: int
    time: int
    current: int
    temperature: int
    unknown: bool
    smass: int = power * mass
    slength: int = power * length
    stime: int = power * time
    scurrent: int = power * current
    stemperature: int = power * temperature

class Equation:
    # An equation, by a name you choose. The rules sum its factors' shares of
    # each base dimension, which are all 0 when the equation is consistent,
    # and conclude `consistent`, or `unbalanced` when any sum is not 0. They
    # also conclude, for each dimension, how many factors hold a share of it
    # (`mass_factors` and so on) and how many claims of that share they hold
    # (`mass_claims`); where the claims are more than the factors, a factor is
    # bound to the equation twice by that share, the dimension is disputed
    # (`mass_disputed`), and its sum is not concluded.
    mass_factors: int
    length_factors: int
    time_factors: int
    current_factors: int
    temperature_factors: int
    mass_claims: int
    length_claims: int
    time_claims: int
    current_claims: int
    temperature_claims: int
    mass_disputed: bool
    length_disputed: bool
    time_disputed: bool
    current_disputed: bool
    temperature_disputed: bool
    mass: int
    length: int
    time: int
    current: int
    temperature: int
    consistent: bool
    unbalanced: bool

# ── Circuits as facts ────────────────────────────────────────────────────────

class Resistor:
    # A resistor, by name: its resistance in ohms, and the group it is joined
    # in, if any.
    ohms: real
    group: string
    siemens: real = 1.0 / ohms

class Group:
    # Resistors joined "series" or "parallel". The rules sum a series group's
    # resistances into `ohms`, and a parallel group's conductances into
    # `siemens`, from which `parallel` is its resistance. A group's members are
    # resistors; a circuit that is not built of series and parallel groups,
    # such as a bridge, is solved with `solve()`. The rules also conclude how
    # many resistors carry a resistance into a series group's sum
    # (`ohms_resistors`) and how many claims of it they hold (`ohms_claims`),
    # and the same of conductance for a parallel group (`siemens_resistors`,
    # `siemens_claims`); a resistor bound to the group twice disputes that sum
    # (`ohms_disputed`, `siemens_disputed`), and it is not concluded.
    joined: string
    ohms_resistors: int
    siemens_resistors: int
    ohms_claims: int
    siemens_claims: int
    ohms_disputed: bool
    siemens_disputed: bool
    ohms: real
    siemens: real
    parallel: real = 1.0 / siemens

# ── Standing rules: a quantity's unit ────────────────────────────────────────

@rule
def factor_of():
    for q in Quantity, u in Unit:
        if q.unit == u:
            q.factor = u.factor

@rule
def mass_of():
    for q in Quantity, u in Unit:
        if q.unit == u:
            q.mass = u.mass

@rule
def length_of():
    for q in Quantity, u in Unit:
        if q.unit == u:
            q.length = u.length

@rule
def time_of():
    for q in Quantity, u in Unit:
        if q.unit == u:
            q.time = u.time

@rule
def current_of():
    for q in Quantity, u in Unit:
        if q.unit == u:
            q.current = u.current

@rule
def temperature_of():
    for q in Quantity, u in Unit:
        if q.unit == u:
            q.temperature = u.temperature

@rule
def unknown_unit():
    # A unit the table does not hold, so no factor reaches the quantity.
    for q in Quantity:
        if q.unit != "" and absent(q.factor):
            q.unknown = True

# ── Standing rules: a quantity handed over ───────────────────────────────────
#
# The unit the reader expects, and its dimension, beside the unit the
# quantity is written in. Same dimension and a different unit: the value must
# be converted, and `converted` is what it reads in the expected unit. A
# different dimension: no conversion makes it the quantity the reader wants.

@rule
def efactor_of():
    for q in Quantity, u in Unit:
        if q.expects == u:
            q.efactor = u.factor

@rule
def emass_of():
    for q in Quantity, u in Unit:
        if q.expects == u:
            q.emass = u.mass

@rule
def elength_of():
    for q in Quantity, u in Unit:
        if q.expects == u:
            q.elength = u.length

@rule
def etime_of():
    for q in Quantity, u in Unit:
        if q.expects == u:
            q.etime = u.time

@rule
def ecurrent_of():
    for q in Quantity, u in Unit:
        if q.expects == u:
            q.ecurrent = u.current

@rule
def etemperature_of():
    for q in Quantity, u in Unit:
        if q.expects == u:
            q.etemperature = u.temperature

@rule
def unexpected_unit():
    # An expected unit the table does not hold.
    for q in Quantity:
        if q.expects != "" and absent(q.efactor):
            q.unexpected = True

@rule
def handoff_convert():
    for q in Quantity:
        if q.unit != q.expects and q.mass == q.emass and q.length == q.elength and q.time == q.etime and q.current == q.ecurrent and q.temperature == q.etemperature:
            q.handoff = "convert"

@rule
def handoff_mass():
    for q in Quantity:
        if q.mass != q.emass:
            q.handoff = "incompatible"

@rule
def handoff_length():
    for q in Quantity:
        if q.length != q.elength:
            q.handoff = "incompatible"

@rule
def handoff_time():
    for q in Quantity:
        if q.time != q.etime:
            q.handoff = "incompatible"

@rule
def handoff_current():
    for q in Quantity:
        if q.current != q.ecurrent:
            q.handoff = "incompatible"

@rule
def handoff_temperature():
    for q in Quantity:
        if q.temperature != q.etemperature:
            q.handoff = "incompatible"

# ── Standing rules: dimensions of an equation ────────────────────────────────
#
# Each factor takes its quantity's dimension from the quantity's unit, its
# share is that dimension times its power, and the equation sums the shares.

@rule
def factor_mass():
    for f in Factor, q in Quantity, u in Unit:
        if f.quantity == q and q.unit == u:
            f.mass = u.mass

@rule
def factor_length():
    for f in Factor, q in Quantity, u in Unit:
        if f.quantity == q and q.unit == u:
            f.length = u.length

@rule
def factor_time():
    for f in Factor, q in Quantity, u in Unit:
        if f.quantity == q and q.unit == u:
            f.time = u.time

@rule
def factor_current():
    for f in Factor, q in Quantity, u in Unit:
        if f.quantity == q and q.unit == u:
            f.current = u.current

@rule
def factor_temperature():
    for f in Factor, q in Quantity, u in Unit:
        if f.quantity == q and q.unit == u:
            f.temperature = u.temperature

@rule
def factor_unknown():
    # A factor whose quantity or unit is not recorded has no dimension to
    # share, so the equation it is in is not judged consistent.
    for f in Factor:
        if f.quantity != "" and absent(f.mass):
            f.unknown = True

@rule
def equation_mass_factors():
    for e in Equation, f in Factor:
        if f.equation == e and f.smass > -1000000000:
            e.mass_factors = distinct(f)

@rule
def equation_length_factors():
    for e in Equation, f in Factor:
        if f.equation == e and f.slength > -1000000000:
            e.length_factors = distinct(f)

@rule
def equation_time_factors():
    for e in Equation, f in Factor:
        if f.equation == e and f.stime > -1000000000:
            e.time_factors = distinct(f)

@rule
def equation_current_factors():
    for e in Equation, f in Factor:
        if f.equation == e and f.scurrent > -1000000000:
            e.current_factors = distinct(f)

@rule
def equation_temperature_factors():
    for e in Equation, f in Factor:
        if f.equation == e and f.stemperature > -1000000000:
            e.temperature_factors = distinct(f)

@rule
def equation_mass_claims():
    for e in Equation, f in Factor:
        if f.equation == e and f.smass > -1000000000:
            e.mass_claims = count(f)

@rule
def equation_length_claims():
    for e in Equation, f in Factor:
        if f.equation == e and f.slength > -1000000000:
            e.length_claims = count(f)

@rule
def equation_time_claims():
    for e in Equation, f in Factor:
        if f.equation == e and f.stime > -1000000000:
            e.time_claims = count(f)

@rule
def equation_current_claims():
    for e in Equation, f in Factor:
        if f.equation == e and f.scurrent > -1000000000:
            e.current_claims = count(f)

@rule
def equation_temperature_claims():
    for e in Equation, f in Factor:
        if f.equation == e and f.stemperature > -1000000000:
            e.temperature_claims = count(f)

@rule
def equation_mass_disputed():
    # A sum adds every value it reads, so a factor bound to the equation twice
    # would be counted twice, which no source gives.
    for e in Equation:
        if e.mass_claims > e.mass_factors:
            e.mass_disputed = True

@rule
def equation_length_disputed():
    for e in Equation:
        if e.length_claims > e.length_factors:
            e.length_disputed = True

@rule
def equation_time_disputed():
    for e in Equation:
        if e.time_claims > e.time_factors:
            e.time_disputed = True

@rule
def equation_current_disputed():
    for e in Equation:
        if e.current_claims > e.current_factors:
            e.current_disputed = True

@rule
def equation_temperature_disputed():
    for e in Equation:
        if e.temperature_claims > e.temperature_factors:
            e.temperature_disputed = True

@rule
def equation_mass():
    for e in Equation, f in Factor:
        if f.equation == e and absent(e.mass_disputed):
            e.mass = sum(f.smass)

@rule
def equation_length():
    for e in Equation, f in Factor:
        if f.equation == e and absent(e.length_disputed):
            e.length = sum(f.slength)

@rule
def equation_time():
    for e in Equation, f in Factor:
        if f.equation == e and absent(e.time_disputed):
            e.time = sum(f.stime)

@rule
def equation_current():
    for e in Equation, f in Factor:
        if f.equation == e and absent(e.current_disputed):
            e.current = sum(f.scurrent)

@rule
def equation_temperature():
    for e in Equation, f in Factor:
        if f.equation == e and absent(e.temperature_disputed):
            e.temperature = sum(f.stemperature)

@rule
def consistent():
    for e in Equation:
        if e.mass == 0 and e.length == 0 and e.time == 0 and e.current == 0 and e.temperature == 0 and none(f in Factor: f.equation == e and f.unknown == True):
            e.consistent = True

@rule
def unbalanced_mass():
    for e in Equation:
        if e.mass != 0 and none(f in Factor: f.equation == e and f.unknown == True):
            e.unbalanced = True

@rule
def unbalanced_length():
    for e in Equation:
        if e.length != 0 and none(f in Factor: f.equation == e and f.unknown == True):
            e.unbalanced = True

@rule
def unbalanced_time():
    for e in Equation:
        if e.time != 0 and none(f in Factor: f.equation == e and f.unknown == True):
            e.unbalanced = True

@rule
def unbalanced_current():
    for e in Equation:
        if e.current != 0 and none(f in Factor: f.equation == e and f.unknown == True):
            e.unbalanced = True

@rule
def unbalanced_temperature():
    for e in Equation:
        if e.temperature != 0 and none(f in Factor: f.equation == e and f.unknown == True):
            e.unbalanced = True

# ── Standing rules: series and parallel ──────────────────────────────────────

@rule
def group_ohms_resistors():
    for g in Group, r in Resistor:
        if g.joined == "series" and r.group == g and r.ohms > -1.0e300:
            g.ohms_resistors = distinct(r)

@rule
def group_siemens_resistors():
    for g in Group, r in Resistor:
        if g.joined == "parallel" and r.group == g and r.siemens > -1.0e300:
            g.siemens_resistors = distinct(r)

@rule
def group_ohms_claims():
    for g in Group, r in Resistor:
        if g.joined == "series" and r.group == g and r.ohms > -1.0e300:
            g.ohms_claims = count(r)

@rule
def group_siemens_claims():
    for g in Group, r in Resistor:
        if g.joined == "parallel" and r.group == g and r.siemens > -1.0e300:
            g.siemens_claims = count(r)

@rule
def group_ohms_disputed():
    # A sum adds every value it reads, so a resistor bound to the group twice
    # would be counted twice, which no source gives.
    for g in Group:
        if g.ohms_claims > g.ohms_resistors:
            g.ohms_disputed = True

@rule
def group_siemens_disputed():
    for g in Group:
        if g.siemens_claims > g.siemens_resistors:
            g.siemens_disputed = True

@rule
def series():
    # Resistors in series carry one current, so their resistances add.
    for g in Group, r in Resistor:
        if g.joined == "series" and r.group == g and absent(g.ohms_disputed):
            g.ohms = sum(r.ohms)

@rule
def parallel():
    # Resistors in parallel share one voltage, so their conductances add.
    for g in Group, r in Resistor:
        if g.joined == "parallel" and r.group == g and absent(g.siemens_disputed):
            g.siemens = sum(r.siemens)

# ── Writing quantities ───────────────────────────────────────────────────────

def unit(name: string, mass: int, length: int, time: int, current: int, temperature: int, factor: real, source: string):
    # One unit, with the source of its value.
    update(Unit(name).mass, mass, source)
    update(Unit(name).length, length, source)
    update(Unit(name).time, time, source)
    update(Unit(name).current, current, source)
    update(Unit(name).temperature, temperature, source)
    update(Unit(name).factor, factor, source)

def si_units():
    # The seven SI base units but the mole and the candela, which no rule
    # here reads, named SI units of mechanics and electricity, and units in
    # common use beside them. The SI values are from the SI Brochure, 9th
    # edition (2019); the others are exact by definition, as NIST Special
    # Publication 811 (2008) gives them. A pound is 0.45359237 kg and a
    # pound-force is a pound under standard gravity, 9.80665 m/s2.
    si = "SI Brochure 2019"
    nist = "NIST SP 811"
    unit("kg", 1, 0, 0, 0, 0, 1.0, si)
    unit("m", 0, 1, 0, 0, 0, 1.0, si)
    unit("s", 0, 0, 1, 0, 0, 1.0, si)
    unit("A", 0, 0, 0, 1, 0, 1.0, si)
    unit("K", 0, 0, 0, 0, 1, 1.0, si)
    unit("one", 0, 0, 0, 0, 0, 1.0, si)
    unit("rad", 0, 0, 0, 0, 0, 1.0, si)
    unit("g", 1, 0, 0, 0, 0, 0.001, si)
    unit("km", 0, 1, 0, 0, 0, 1000.0, si)
    unit("cm", 0, 1, 0, 0, 0, 0.01, si)
    unit("mm", 0, 1, 0, 0, 0, 0.001, si)
    unit("min", 0, 0, 1, 0, 0, 60.0, si)
    unit("h", 0, 0, 1, 0, 0, 3600.0, si)
    unit("Hz", 0, 0, -1, 0, 0, 1.0, si)
    unit("m_per_s", 0, 1, -1, 0, 0, 1.0, si)
    unit("km_per_h", 0, 1, -1, 0, 0, 1000.0 / 3600.0, si)
    unit("m_per_s2", 0, 1, -2, 0, 0, 1.0, si)
    unit("N", 1, 1, -2, 0, 0, 1.0, si)
    unit("N_s", 1, 1, -1, 0, 0, 1.0, si)
    unit("kg_m_per_s", 1, 1, -1, 0, 0, 1.0, si)
    unit("J", 1, 2, -2, 0, 0, 1.0, si)
    unit("W", 1, 2, -3, 0, 0, 1.0, si)
    unit("Pa", 1, -1, -2, 0, 0, 1.0, si)
    unit("C", 0, 0, 1, 1, 0, 1.0, si)
    unit("V", 1, 2, -3, -1, 0, 1.0, si)
    unit("ohm", 1, 2, -3, -2, 0, 1.0, si)
    unit("kWh", 1, 2, -2, 0, 0, 3600000.0, si)
    unit("eV", 1, 2, -2, 0, 0, 1.602176634e-19, si)
    unit("in", 0, 1, 0, 0, 0, 0.0254, nist)
    unit("ft", 0, 1, 0, 0, 0, 0.3048, nist)
    unit("mi", 0, 1, 0, 0, 0, 1609.344, nist)
    unit("mph", 0, 1, -1, 0, 0, 0.44704, nist)
    unit("lb", 1, 0, 0, 0, 0, 0.45359237, nist)
    unit("lbf", 1, 1, -2, 0, 0, 4.4482216152605, nist)
    unit("lbf_s", 1, 1, -1, 0, 0, 4.4482216152605, nist)
    unit("ft_lbf", 1, 2, -2, 0, 0, 1.3558179483314004, nist)
    unit("psi", 1, -1, -2, 0, 0, 6894.757293168362, nist)
    unit("atm", 1, -1, -2, 0, 0, 101325.0, nist)
    unit("cal", 1, 2, -2, 0, 0, 4.184, nist)

def quantity(name: string, value: real, unit: string):
    # A value in a unit.
    update(Quantity(name).value, value, "physics")
    update(Quantity(name).unit, unit, "physics")

def expect(name: string, unit: string):
    # The unit the quantity's reader takes it in.
    update(Quantity(name).expects, unit, "physics")

def factor(equation: string, quantity: string, power: int):
    # One factor of an equation: positive powers on the left, negative on the
    # right.
    f = equation + "-" + quantity
    update(Factor(f).equation, equation, "physics")
    update(Factor(f).quantity, quantity, "physics")
    update(Factor(f).power, power, "physics")

def resistor(name: string, ohms: real, group: string):
    # A resistor and the group it is joined in.
    update(Resistor(name).ohms, ohms, "physics")
    update(Resistor(name).group, group, "physics")

def group(name: string, joined: string):
    # A group, joined "series" or "parallel".
    update(Group(name).joined, joined, "physics")

# ── Reading what the rules concluded ─────────────────────────────────────────
#
# Each reader reads every value in force, so a fact two sources disagree
# about is reported with both values rather than refused; only a reader of a
# disputed sum refuses (above).

def si(name: string) -> real:
    # The quantity's value in SI: the least, where values disagree, and 0.0
    # when there is none.
    vs = list_values(Quantity.si)
    out = 0.0
    found = False
    for v in vs:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == name:
            match v.value:
                case real as x:
                    if not found or x < out:
                        out = x
                        found = True
                case default:
                    pass
    return out

def converted(name: string) -> real:
    # The quantity's value in the unit its reader expects: the least, where
    # values disagree, and 0.0 when there is none.
    vs = list_values(Quantity.converted)
    out = 0.0
    found = False
    for v in vs:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == name:
            match v.value:
                case real as x:
                    if not found or x < out:
                        out = x
                        found = True
                case default:
                    pass
    return out

def resistance(name: string) -> real:
    # A group's resistance in ohms, series or parallel: the least, where
    # values disagree, and 0.0 when there is none. It refuses a disputed group,
    # which has no resistance until the dispute is settled.
    if name in disputed_groups():
        raise ValueError("group " + name + " is disputed: a resistor of it is bound to it twice; forget the claim that is not to stand")
    out = 0.0
    found = False
    for f in ["series", "parallel"]:
        if f == "series":
            vs = list_values(Group.ohms)
        else:
            vs = list_values(Group.parallel)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name:
                match v.value:
                    case real as x:
                        if not found or x < out:
                            out = x
                            found = True
                    case default:
                        pass
    return out

def disputed_equations() -> list<string>:
    # The equations a factor is bound to twice by its share of some dimension,
    # whose sum for that dimension waits until one value is forgotten, sorted.
    out = []
    for d in ["mass", "length", "time", "current", "temperature"]:
        if d == "mass":
            vs = list_values(Equation.mass_disputed)
        elif d == "length":
            vs = list_values(Equation.length_disputed)
        elif d == "time":
            vs = list_values(Equation.time_disputed)
        elif d == "current":
            vs = list_values(Equation.current_disputed)
        else:
            vs = list_values(Equation.temperature_disputed)
        for v in vs:
            parts = v.path.split(".")
            x = parts[len(parts) - 2]
            if x not in out:
                out = out + [x]
    return sorted(out)

def disputed_groups() -> list<string>:
    # The groups a resistor is bound to twice, whose resistance waits until
    # one value is forgotten, sorted.
    out = []
    for f in ["ohms", "siemens"]:
        if f == "ohms":
            vs = list_values(Group.ohms_disputed)
        else:
            vs = list_values(Group.siemens_disputed)
        for v in vs:
            parts = v.path.split(".")
            x = parts[len(parts) - 2]
            if x not in out:
                out = out + [x]
    return sorted(out)

def dimension(name: string) -> string:
    # What is left over of an equation's dimensions, as powers of the base
    # dimensions' symbols, M, L, T, I and Θ: "L1 T-1" for `E = m v`, and
    # "consistent" when nothing is. "unknown unit" when a factor's unit is not
    # recorded, "no factors" when the equation has none, and "disputed" when a
    # dimension's sum is disputed.
    if name in disputed_equations():
        return "disputed"
    sums = {"M": 0, "L": 0, "T": 0, "I": 0, "Θ": 0}
    found = False
    for d in ["M", "L", "T", "I", "Θ"]:
        if d == "M":
            vs = list_values(Equation.mass)
        elif d == "L":
            vs = list_values(Equation.length)
        elif d == "T":
            vs = list_values(Equation.time)
        elif d == "I":
            vs = list_values(Equation.current)
        else:
            vs = list_values(Equation.temperature)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name:
                match v.value:
                    case int as x:
                        if x != 0:
                            sums[d] = x
                        found = True
                    case default:
                        pass
    us = list_values(Factor.unknown)
    es = list_values(Factor.equation)
    for u in us:
        uparts = u.path.split(".")
        for w in es:
            wparts = w.path.split(".")
            if wparts[len(wparts) - 2] == uparts[len(uparts) - 2] and str(w.value) == name:
                return "unknown unit"
    if not found:
        return "no factors"
    left = []
    for d in ["M", "L", "T", "I", "Θ"]:
        if sums[d] != 0:
            left = left + [d + str(sums[d])]
    if len(left) == 0:
        return "consistent"
    return " ".join(left)

def checks() -> list<string>:
    # Everything the rules found: units the table does not hold, equations
    # whose dimensions do not balance, quantities handed over in another unit
    # than their reader's, and disputed equations and groups.
    out = []
    units = {}
    vs = list_values(Quantity.unit)
    for v in vs:
        parts = v.path.split(".")
        units[parts[len(parts) - 2]] = str(v.value)
    wanted = {}
    vs = list_values(Quantity.expects)
    for v in vs:
        parts = v.path.split(".")
        wanted[parts[len(parts) - 2]] = str(v.value)
    vs = list_values(Quantity.unknown)
    for v in vs:
        parts = v.path.split(".")
        q = parts[len(parts) - 2]
        out = out + ["unknown unit: " + q + " in " + units[q]]
    vs = list_values(Quantity.unexpected)
    for v in vs:
        parts = v.path.split(".")
        q = parts[len(parts) - 2]
        out = out + ["unknown unit: " + q + " expected in " + wanted[q]]
    vs = list_values(Equation.unbalanced)
    for v in vs:
        parts = v.path.split(".")
        e = parts[len(parts) - 2]
        left = dimension(e)
        line = "unbalanced: " + e + ", " + left + " left over"
        if left != "disputed" and line not in out:
            out = out + [line]
    for e in disputed_equations():
        out = out + ["disputed: equation " + e]
    for g in disputed_groups():
        out = out + ["disputed: group " + g]
    vs = list_values(Quantity.handoff)
    for v in vs:
        parts = v.path.split(".")
        q = parts[len(parts) - 2]
        line = "handed over: " + q + " in " + units[q] + " where " + wanted[q] + " is expected, " + str(v.value)
        if line not in out:
            out = out + [line]
    return sorted(out)

def unsettled() -> list<string>:
    # Each unit, quantity, factor, group and resistor that holds two different
    # values for a fact you record about it: a unit's five powers and its
    # factor, a quantity's value, unit and expected unit, a factor's equation,
    # quantity and power, a group's joining and a resistor's ohms and group, as
    # "unit x", "quantity x", "factor x", "group x" or "resistor x", sorted.
    out = []
    for f in ["mass", "length", "time", "current", "temperature", "factor", "value", "unit", "expects", "equation", "quantity", "power", "joined", "ohms", "group"]:
        if f == "mass":
            vs = list_values(Unit.mass)
            label = "unit "
        elif f == "length":
            vs = list_values(Unit.length)
            label = "unit "
        elif f == "time":
            vs = list_values(Unit.time)
            label = "unit "
        elif f == "current":
            vs = list_values(Unit.current)
            label = "unit "
        elif f == "temperature":
            vs = list_values(Unit.temperature)
            label = "unit "
        elif f == "factor":
            vs = list_values(Unit.factor)
            label = "unit "
        elif f == "value":
            vs = list_values(Quantity.value)
            label = "quantity "
        elif f == "unit":
            vs = list_values(Quantity.unit)
            label = "quantity "
        elif f == "expects":
            vs = list_values(Quantity.expects)
            label = "quantity "
        elif f == "equation":
            vs = list_values(Factor.equation)
            label = "factor "
        elif f == "quantity":
            vs = list_values(Factor.quantity)
            label = "factor "
        elif f == "power":
            vs = list_values(Factor.power)
            label = "factor "
        elif f == "joined":
            vs = list_values(Group.joined)
            label = "group "
        elif f == "ohms":
            vs = list_values(Resistor.ohms)
            label = "resistor "
        else:
            vs = list_values(Resistor.group)
            label = "resistor "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = label + parts[len(parts) - 2]
            x = str(v.value)
            if s in seen and seen[s] != x and s not in out:
                out = out + [s]
            seen[s] = x
    return sorted(out)

# ── Calculation: the pendulum ────────────────────────────────────────────────

def small_period(length: real, g: real) -> real:
    # The period of a pendulum swinging through a small angle, 2π√(L/g), the
    # same for every amplitude. Where L/g passes the range of a real, or falls
    # below the smallest normal one, while L and g are reals, it is taken as
    # √L/√g.
    if not (length >= 0.0):
        raise ValueError("a length is not below 0")
    if not (g > 0.0):
        raise ValueError("g is above 0")
    q = length / g
    if length == 0.0 or (q >= 2.2250738585072014e-308 and math.isfinite(q)) or not (math.isfinite(length) and math.isfinite(g)):
        return 2.0 * math.pi * math.sqrt(q)
    return 2.0 * math.pi * (math.sqrt(length) / math.sqrt(g))

def agm(a: real, b: real) -> real:
    # The arithmetic–geometric mean of two numbers not below 0: take the mean
    # and the geometric mean of the pair until they agree. Where the sum of
    # the pair passes the range of a real, the mean is taken as the sum of
    # their halves, and where their product passes it or falls below the
    # smallest normal real, the geometric mean as the product of their square
    # roots.
    if not (a >= 0.0 and b >= 0.0):
        raise ValueError("an arithmetic–geometric mean takes two numbers not below 0")
    x = a
    y = b
    for step in range(40):
        nx = (x + y) / 2.0
        if not math.isfinite(nx):
            nx = x / 2.0 + y / 2.0
        p = x * y
        if (p >= 2.2250738585072014e-308 and math.isfinite(p)) or x == 0.0 or y == 0.0 or not (math.isfinite(x) and math.isfinite(y)):
            y = math.sqrt(p)
        else:
            y = math.sqrt(x) * math.sqrt(y)
        x = nx
        if abs(x - y) <= 1e-15 * x:
            break
    return x

def period(theta0: real, length: real, g: real) -> real:
    # The exact period of a pendulum released at rest from `theta0` radians:
    # the small-angle period divided by the arithmetic–geometric mean of 1 and
    # cos(theta0 / 2). An amplitude is from -π to π.
    if not (abs(theta0) <= math.pi):
        raise ValueError("an amplitude is from -π to π")
    return small_period(length, g) / agm(1.0, math.cos(theta0 / 2.0))

def pendulum(theta0: real, length: real, g: real, dt: real) -> real:
    # The period of the same pendulum by integrating its equation of motion,
    # θ'' = −(g/L) sin θ, with the fourth-order Runge–Kutta method at steps of
    # `dt` seconds: four times the time it takes to swing down to the vertical.
    # The amplitude is above 0 and at most π, and the integration runs for up
    # to the exact period; a step too long to reach the vertical in that time
    # is refused.
    if not (theta0 > 0.0 and theta0 <= math.pi):
        raise ValueError("an amplitude is above 0 and at most π")
    if not (length > 0.0):
        raise ValueError("a length is above 0")
    if not (dt > 0.0):
        raise ValueError("a step is above 0")
    w = g / length
    th = theta0
    om = 0.0
    tm = 0.0
    steps = int(period(theta0, length, g) / dt) + 1
    for step in range(steps):
        k1t = om
        k1o = -w * math.sin(th)
        k2t = om + 0.5 * dt * k1o
        k2o = -w * math.sin(th + 0.5 * dt * k1t)
        k3t = om + 0.5 * dt * k2o
        k3o = -w * math.sin(th + 0.5 * dt * k2t)
        k4t = om + dt * k3o
        k4o = -w * math.sin(th + dt * k3t)
        nth = th + dt * (k1t + 2.0 * k2t + 2.0 * k3t + k4t) / 6.0
        nom = om + dt * (k1o + 2.0 * k2o + 2.0 * k3o + k4o) / 6.0
        if nth <= 0.0:
            return 4.0 * (tm + dt * th / (th - nth))
        th = nth
        om = nom
        tm = tm + dt
    raise ValueError("the step is too long for the swing, and the integration does not reach the vertical")

# ── Calculation: an orbit ────────────────────────────────────────────────────

def orbit(e: real, method: string, periods: int, steps: int) -> list<real>:
    # A body on a Kepler orbit of eccentricity `e` round a mass with GM = 1,
    # on a semi-major axis of 1, so a period is 2π, started at its closest
    # point. `method` is "euler", "verlet" or "rk4", taking `steps` steps a
    # period for `periods` periods. Returns the relative error in the energy
    # and in the angular momentum at the end, and the distance from the start,
    # which an exact orbit closes to 0. The eccentricity is from 0 to below 1.
    if not (e >= 0.0 and e < 1.0):
        raise ValueError("an eccentricity is from 0 to below 1")
    x = 1.0 - e
    y = 0.0
    vx = 0.0
    vy = math.sqrt((1.0 + e) / (1.0 - e))
    e0 = 0.5 * (vx * vx + vy * vy) - 1.0 / math.sqrt(x * x + y * y)
    h0 = x * vy - y * vx
    dt = 2.0 * math.pi / steps
    for step in range(periods * steps):
        if method == "euler":
            r3 = (x * x + y * y) * math.sqrt(x * x + y * y)
            ax = -x / r3
            ay = -y / r3
            x = x + dt * vx
            y = y + dt * vy
            vx = vx + dt * ax
            vy = vy + dt * ay
        elif method == "verlet":
            r3 = (x * x + y * y) * math.sqrt(x * x + y * y)
            hx = vx - 0.5 * dt * x / r3
            hy = vy - 0.5 * dt * y / r3
            x = x + dt * hx
            y = y + dt * hy
            r3 = (x * x + y * y) * math.sqrt(x * x + y * y)
            vx = hx - 0.5 * dt * x / r3
            vy = hy - 0.5 * dt * y / r3
        else:
            r3 = (x * x + y * y) * math.sqrt(x * x + y * y)
            k1x = vx
            k1y = vy
            k1u = -x / r3
            k1v = -y / r3
            x2 = x + 0.5 * dt * k1x
            y2 = y + 0.5 * dt * k1y
            r3 = (x2 * x2 + y2 * y2) * math.sqrt(x2 * x2 + y2 * y2)
            k2x = vx + 0.5 * dt * k1u
            k2y = vy + 0.5 * dt * k1v
            k2u = -x2 / r3
            k2v = -y2 / r3
            x3 = x + 0.5 * dt * k2x
            y3 = y + 0.5 * dt * k2y
            r3 = (x3 * x3 + y3 * y3) * math.sqrt(x3 * x3 + y3 * y3)
            k3x = vx + 0.5 * dt * k2u
            k3y = vy + 0.5 * dt * k2v
            k3u = -x3 / r3
            k3v = -y3 / r3
            x4 = x + dt * k3x
            y4 = y + dt * k3y
            r3 = (x4 * x4 + y4 * y4) * math.sqrt(x4 * x4 + y4 * y4)
            k4x = vx + dt * k3u
            k4y = vy + dt * k3v
            k4u = -x4 / r3
            k4v = -y4 / r3
            x = x + dt * (k1x + 2.0 * k2x + 2.0 * k3x + k4x) / 6.0
            y = y + dt * (k1y + 2.0 * k2y + 2.0 * k3y + k4y) / 6.0
            vx = vx + dt * (k1u + 2.0 * k2u + 2.0 * k3u + k4u) / 6.0
            vy = vy + dt * (k1v + 2.0 * k2v + 2.0 * k3v + k4v) / 6.0
    e1 = 0.5 * (vx * vx + vy * vy) - 1.0 / math.sqrt(x * x + y * y)
    h1 = x * vy - y * vx
    dx = x - 1.0 + e
    q = dx * dx + y * y
    gap = math.sqrt(q)
    if q < 2.2250738585072014e-308 and (dx != 0.0 or y != 0.0):
        # The squares fall below the smallest normal real: the distance by
        # hypot, which squares nothing.
        gap = math.hypot(dx, y)
    return [(e1 - e0) / abs(e0), (h1 - h0) / h0, gap]

# ── Calculation: heat ────────────────────────────────────────────────────────

def heat(cells: int, r: real, steps: int) -> real:
    # A rod of unit length and unit diffusivity, held at 0 at both ends and
    # started as a triangle rising to 1 at its middle, stepped by forward
    # differences in time and centred differences in space: each step sets a
    # point to itself plus `r` times its neighbours' sum less twice itself,
    # where `r` is the time step over the square of the spacing. Returns the
    # temperature at the middle after `steps` steps, at time
    # steps · r / cells². `cells` is even.
    u = []
    for j in range(cells + 1):
        x = j / cells
        u = u + [1.0 - abs(2.0 * x - 1.0)]
    nxt = u[:]
    inner = range(1, cells)
    for step in range(steps):
        for j in inner:
            nxt[j] = u[j] + r * (u[j - 1] - 2.0 * u[j] + u[j + 1])
        for j in inner:
            u[j] = nxt[j]
    return u[cells // 2]

def heat_exact(time: real) -> real:
    # The same rod's temperature at its middle at `time`, from Fourier's
    # series for the triangle: the sum over odd n of 8 / (n²π²) e^(−n²π²t).
    total = 0.0
    for j in range(200):
        n = 2 * j + 1
        total = total + 8.0 / (n * n * math.pi * math.pi) * math.exp(-n * n * math.pi * math.pi * time)
    return total

# ── Calculation: linear systems and the bridge ───────────────────────────────

def solve(a: list<list<real>>, b: list<real>) -> list<real>:
    # The solution of the linear system a · x = b, by Gaussian elimination
    # with partial pivoting: at each column the row with the largest entry is
    # swapped up before it is used. An empty list when the system is
    # singular, that is where a pivot is 0. Where x passes the range of a
    # real, or a pivot or a multiplier falls below the smallest normal real
    # and loses the row it scales, while a and b are reals, the system is
    # solved again with each row of a scaled by a power of two that brings its
    # largest entry near 1, and b by another, which changes no x, and x is
    # scaled back.
    n = len(b)
    rows = a
    rhs = b
    shift = 0
    for attempt in range(2):
        m = []
        for j in range(n):
            m = m + [rows[j][:] + [rhs[j]]]
        lost = False
        for c in range(n):
            p = c
            for j in range(c + 1, n):
                if abs(m[j][c]) > abs(m[p][c]):
                    p = j
            if m[p][c] == 0.0:
                return []
            if abs(m[p][c]) < 2.2250738585072014e-308:
                lost = True
            if p != c:
                row = m[c]
                m[c] = m[p]
                m[p] = row
            for j in range(c + 1, n):
                f = m[j][c] / m[c][c]
                if m[j][c] != 0.0 and abs(f) < 2.2250738585072014e-308:
                    lost = True
                row = m[j]
                for col in range(c, n + 1):
                    row[col] = row[col] - f * m[c][col]
                m[j] = row
        x = []
        for j in range(n):
            x = x + [0.0]
        for back in range(n):
            j = n - 1 - back
            s = m[j][n]
            for col in range(j + 1, n):
                s = s - m[j][col] * x[col]
            x[j] = s / m[j][j]
        if attempt == 1:
            return [math.ldexp(v, shift) for v in x]
        if (all(math.isfinite(v) for v in x) and not lost) or not all(math.isfinite(v) for row in a for v in row) or not all(math.isfinite(v) for v in b):
            return x
        rows = []
        rhs = []
        for j in range(n):
            top = 0.0
            for v in a[j]:
                if abs(v) > top:
                    top = abs(v)
            e = 0
            if top > 0.0:
                e = math.frexp(top)[1]
            rows = rows + [[math.ldexp(v, -e) for v in a[j]]]
            rhs = rhs + [math.ldexp(b[j], -e)]
        big = 0.0
        for v in rhs:
            if abs(v) > big:
                big = abs(v)
        shift = n.bit_length() + n + 1
        if big > 1.0:
            shift = shift + math.frexp(big)[1]
        rhs = [math.ldexp(v, -shift) for v in rhs]
    return []

def bridge(volts: real, r1: real, r2: real, r3: real, r4: real, rg: real) -> real:
    # The current in amperes through the meter of a Wheatstone bridge: a
    # source of `volts` across two arms, `r1` above `r2` and `r3` above `r4`,
    # with a meter of `rg` ohms across their middles. The two middle voltages
    # follow from Kirchhoff's current law at each middle node, solved by
    # `solve()`. 0 when the bridge balances, r1 / r2 = r3 / r4. Every
    # resistance is above 0. Where a conductance or a current on the way
    # passes the range of a real while the inputs are reals, each node's
    # equation is taken with its resistances scaled by a power of two that
    # brings the least near 1, which changes neither middle voltage, and a
    # source past 1 volt with its volts scaled by another, which scales both,
    # and the current is scaled back.
    for r in [r1, r2, r3, r4, rg]:
        if not (r > 0.0):
            raise ValueError("a resistance is above 0")
    a = [[1.0 / r1 + 1.0 / r2 + 1.0 / rg, -1.0 / rg], [-1.0 / rg, 1.0 / r3 + 1.0 / r4 + 1.0 / rg]]
    b = [volts / r1, volts / r3]
    v = solve(a, b)
    i = (v[0] - v[1]) / rg
    if math.isfinite(i) or not all(math.isfinite(u) for u in [volts, r1, r2, r3, r4, rg]):
        return i
    s1 = -math.frexp(min(r1, r2, rg))[1]
    s2 = -math.frexp(min(r3, r4, rg))[1]
    q1 = math.ldexp(r1, s1)
    q2 = math.ldexp(r2, s1)
    g1 = math.ldexp(rg, s1)
    q3 = math.ldexp(r3, s2)
    q4 = math.ldexp(r4, s2)
    g2 = math.ldexp(rg, s2)
    sv = 0
    if abs(volts) > 1.0:
        sv = math.frexp(volts)[1]
    u = math.ldexp(volts, -sv)
    a = [[1.0 / q1 + 1.0 / q2 + 1.0 / g1, -1.0 / g1], [-1.0 / g2, 1.0 / q3 + 1.0 / q4 + 1.0 / g2]]
    b = [u / q1, u / q3]
    v = solve(a, b)
    return math.ldexp((v[0] - v[1]) / rg, sv)

# ── Calculation: the Ising model ─────────────────────────────────────────────

def ising(n: int, temperature: real, sweeps: int, seed: int) -> real:
    # A magnet as an n by n square of spins, each up or down, every spin
    # coupled to its four neighbours with the edges wrapped round, at
    # `temperature` in units of the coupling over Boltzmann's constant. Each
    # sweep visits every spin once and flips it by the Metropolis rule: always
    # when that lowers the energy, and otherwise with probability e^(−ΔE/T).
    # Starts with every spin up, and returns the mean of |M| per spin over the
    # second half of the sweeps. The random numbers are Park and Miller's
    # minimal standard generator from `seed`, so a run is the same every time.
    spins = []
    for j in range(n * n):
        spins = spins + [1]
    cells = range(n * n)
    p4 = math.exp(-4.0 / temperature)
    p8 = math.exp(-8.0 / temperature)
    x = seed
    total = 0.0
    counted = 0
    for sweep in range(sweeps):
        for j in cells:
            row = j // n
            col = j % n
            near = spins[((row + 1) % n) * n + col] + spins[((row + n - 1) % n) * n + col] + spins[row * n + (col + 1) % n] + spins[row * n + (col + n - 1) % n]
            de = 2 * spins[j] * near
            flip = de <= 0
            if de == 4:
                x = (16807 * x) % 2147483647
                flip = x / 2147483647.0 < p4
            elif de == 8:
                x = (16807 * x) % 2147483647
                flip = x / 2147483647.0 < p8
            if flip:
                spins[j] = -spins[j]
        if sweep >= sweeps // 2:
            mag = 0
            for j in cells:
                mag = mag + spins[j]
            total = total + abs(mag) / (n * n)
            counted = counted + 1
    return total / counted