# SPDX-License-Identifier: MIT

# chemistry — read molecular structure and reaction enthalpy, and calculate
# reactions.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply molecules as facts: their atoms (`chemistry.Atom`),
# each naming its element and molecule, and the bonds between them
# (`chemistry.Bond`), each with its order. Hydrogens are not written: every
# atom carries the hydrogens its valence leaves over, as a chemist's skeletal
# drawing does. `elements()` writes the table of elements the rules read, and
# `bond()` writes a bond:
#
#     import chemistry
#     chemistry.elements()
#     chemistry.atom("c1", "C", "acid")
#     chemistry.atom("c2", "C", "acid")
#     chemistry.atom("o1", "O", "acid")
#     chemistry.atom("o2", "O", "acid")
#     chemistry.bond("c1", "c2", 1)
#     chemistry.bond("c2", "o1", 2)
#     chemistry.bond("c2", "o2", 1)
#     print(chemistry.formula("acid"), chemistry.weight("acid"))
#     print(chemistry.groups("acid"))
#
# The standing rules then conclude every atom's valence, the hydrogens it
# carries, an atom bonded past its valence, each molecule's weight and its
# hydrogens, and the functional groups on its atoms. Over reactions written as
# species with coefficients, they conclude each reaction's enthalpy by Hess's
# law from the species' enthalpies of formation. 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 atom's bonds, a molecule's weight and hydrogens and a reaction's enthalpy
# are sums, and a sum adds every value it reads: a bond whose order a second
# source writes, even the same order, would be counted twice, and its atom
# would carry a hydrogen fewer than any source gives it. 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. An atom whose bonds are
# disputed has no hydrogens or mass, so its molecule is disputed too rather
# than weighed without it. `weight()`, `formula()` and `enthalpy()` refuse
# what is disputed, `faults()` names each disputed atom, molecule and
# reaction, 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. An
# atom whose molecule, a bond whose atom, or a term whose reaction 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
# formula's molar mass, the balancing of an equation in exact integers, the pH
# of an acid and a titration curve, and the kinetics of consecutive and of
# reversible reactions.

import math

# ── Molecules as facts ───────────────────────────────────────────────────────

class Element:
    # An element, by its symbol: its atomic number `z`, its usual `valence` in
    # an organic molecule, the number of bonds it forms there, 0 where it has
    # none here, and its standard atomic `weight`.
    z: int
    valence: int
    weight: real

class Atom:
    # An atom, by a name you choose: its element's symbol and the molecule it
    # belongs to. The rules conclude its `valence` and `weight` from the
    # element, `bonds` as the sum of its bonds' orders, and `hweight` as
    # hydrogen's weight; `hydrogens` is what its valence leaves over and `mass`
    # its weight with those hydrogens. An element with no valence here, such as
    # a metal, carries no hydrogens. The rules also conclude how many halves of
    # its bonds carry an order (`bond_halves`) and how many claims of order
    # they hold (`bond_claims`); a half bound to the atom twice makes
    # `bonds_disputed`, and then no `bonds` is concluded, nor its hydrogens or
    # mass, and its molecule is disputed.
    element: string
    molecule: string
    valence: int
    weight: real
    bond_halves: int
    bond_claims: int
    bonds_disputed: bool
    bonds: int
    hweight: real
    hydrogens: int = valence - bonds
    mass: real = weight + hydrogens * hweight
    # What the rules conclude about the atom: bonded past its valence, a
    # carbon double-bonded to an oxygen, and the functional group it carries.
    overbonded: bool
    carbonyl: bool
    group: string

class Bond:
    # Half a bond, from atom `a` to atom `b`, of `order` 1, 2 or 3. `bond()`
    # writes both halves, so each atom reaches the other. The rules conclude
    # the element at the far end, its hydrogens, and whether it is a carbonyl
    # carbon, which the functional groups read.
    a: string
    b: string
    order: int
    belement: string
    bhydrogens: int
    bcarbonyl: bool

class Molecule:
    # A molecule, by the name its atoms give. The rules conclude its weight,
    # the sum of its atoms' masses, its hydrogens and its heavy atoms, every
    # atom that is not a hydrogen, each counted once. For each sum they also
    # conclude how many atoms hold a value of it (`weight_atoms`,
    # `hydrogen_atoms`) and how many claims they hold (`weight_claims`,
    # `hydrogen_claims`); an atom bound to the molecule twice by a sum, or an
    # atom whose bonds are disputed, disputes that sum (`weight_disputed`,
    # `hydrogens_disputed`), and it is not concluded.
    weight_atoms: int
    weight_claims: int
    weight_disputed: bool
    hydrogen_atoms: int
    hydrogen_claims: int
    hydrogens_disputed: bool
    weight: real
    hydrogens: int
    heavy: int

# ── Reactions as facts ───────────────────────────────────────────────────────

class Species:
    # A substance in a state, by a name you choose, such as "water_l" for
    # liquid water: its standard enthalpy of formation, in kJ/mol at 298.15 K.
    formation: real

class Term:
    # One species of a reaction and its coefficient: positive for a product,
    # negative for a reactant. The rules conclude its `formation` from the
    # species, and `contribution` is the coefficient times it.
    reaction: string
    species: string
    coefficient: int
    formation: real
    contribution: real = coefficient * formation
    missing: bool

class Reaction:
    # A reaction, by a name you choose. The rules sum its terms'
    # contributions into `enthalpy`, in kJ, and conclude `incomplete` when a
    # term's species has no enthalpy of formation recorded. They also conclude
    # `terms` and `claims`, which `disputed` compares: a term bound to the
    # reaction twice, as two values of its coefficient or its species'
    # formation bind it, makes the reaction `disputed`, and then no enthalpy
    # is concluded.
    terms: int
    claims: int
    disputed: bool
    enthalpy: real
    incomplete: bool

# ── Standing rules: each atom ────────────────────────────────────────────────

@rule
def valence_of():
    for a in Atom, e in Element:
        if a.element == e and e.valence > 0:
            a.valence = e.valence

@rule
def metal_valence():
    # An element with no valence here, a metal, carries no hydrogens: its
    # atom's valence is the bonds it has, so a salt written with its bond
    # weighs what its formula does.
    for a in Atom, e in Element:
        if a.element == e and e.valence == 0 and a.bonds >= 0:
            a.valence = a.bonds

@rule
def weight_of():
    for a in Atom, e in Element:
        if a.element == e:
            a.weight = e.weight

@rule
def hydrogen_weight():
    # Hydrogen is the element of atomic number 1.
    for a in Atom, h in Element:
        if a.element != "" and h.z == 1:
            a.hweight = h.weight

@rule
def bond_halves():
    for a in Atom, x in Bond:
        if x.a == a and x.order > -1000000000000000000:
            a.bond_halves = distinct(x)

@rule
def bond_claims():
    for a in Atom, x in Bond:
        if x.a == a and x.order > -1000000000000000000:
            a.bond_claims = count(x)

@rule
def bonds_disputed():
    # A sum adds every value it reads, so a bond bound to the atom twice would
    # be counted twice, which no source gives.
    for a in Atom:
        if a.bond_claims > a.bond_halves:
            a.bonds_disputed = True

@rule
def bonds_of():
    for a in Atom, x in Bond:
        if x.a == a and absent(a.bonds_disputed):
            a.bonds = sum(x.order)

@rule
def unbonded():
    # An atom with no bond carries its whole valence in hydrogens, as the
    # oxygen of water and the carbon of methane do.
    for a in Atom:
        if a.element != "" and none(x in Bond: x.a == a):
            a.bonds = 0

@rule
def overbonded():
    # More bonds than the element's valence: a structure written wrong.
    for a in Atom:
        if a.hydrogens < 0:
            a.overbonded = True

# ── Standing rules: across a bond ────────────────────────────────────────────

@rule
def far_element():
    for x in Bond, a in Atom:
        if x.b == a:
            x.belement = a.element

@rule
def far_hydrogens():
    for x in Bond, a in Atom:
        if x.b == a:
            x.bhydrogens = a.hydrogens

@rule
def far_carbonyl():
    for x in Bond, a in Atom:
        if x.b == a and a.carbonyl == True:
            x.bcarbonyl = True

# ── Standing rules: each molecule ────────────────────────────────────────────

@rule
def molecule_weight_atoms():
    for m in Molecule, a in Atom:
        if a.molecule == m and a.mass > -1.0e300:
            m.weight_atoms = distinct(a)

@rule
def molecule_weight_claims():
    for m in Molecule, a in Atom:
        if a.molecule == m and a.mass > -1.0e300:
            m.weight_claims = count(a)

@rule
def molecule_weight_disputed():
    # A sum adds every value it reads, so an atom bound to the molecule twice
    # would be counted twice, which no source gives.
    for m in Molecule:
        if m.weight_claims > m.weight_atoms:
            m.weight_disputed = True

@rule
def molecule_hydrogen_atoms():
    for m in Molecule, a in Atom:
        if a.molecule == m and a.hydrogens > -1000000000000000000:
            m.hydrogen_atoms = distinct(a)

@rule
def molecule_hydrogen_claims():
    for m in Molecule, a in Atom:
        if a.molecule == m and a.hydrogens > -1000000000000000000:
            m.hydrogen_claims = count(a)

@rule
def molecule_hydrogens_disputed():
    for m in Molecule:
        if m.hydrogen_claims > m.hydrogen_atoms:
            m.hydrogens_disputed = True

@rule
def molecule_atom_weight_disputed():
    # An atom whose bonds are disputed has no mass, and would be left out of
    # its molecule's weight, which no source gives.
    for m in Molecule, a in Atom:
        if a.molecule == m and a.bonds_disputed == True:
            m.weight_disputed = True

@rule
def molecule_atom_hydrogens_disputed():
    for m in Molecule, a in Atom:
        if a.molecule == m and a.bonds_disputed == True:
            m.hydrogens_disputed = True

@rule
def molecule_weight():
    for m in Molecule, a in Atom:
        if a.molecule == m and absent(m.weight_disputed):
            m.weight = sum(a.mass)

@rule
def molecule_hydrogens():
    for m in Molecule, a in Atom:
        if a.molecule == m and absent(m.hydrogens_disputed):
            m.hydrogens = sum(a.hydrogens)

@rule
def molecule_heavy():
    # Each atom once, however many claims place it in the molecule.
    for m in Molecule, a in Atom:
        if a.molecule == m:
            m.heavy = distinct(a)

# ── Standing rules: functional groups ────────────────────────────────────────
#
# Each is concluded on the atom that carries it: the carbonyl carbon of an
# acid, an ester, an amide, an aldehyde and a ketone, the oxygen of an alcohol
# and the nitrogen of an amine.

@rule
def carbonyl():
    # A carbon double-bonded to an oxygen.
    for c in Atom, x in Bond:
        if c.element == "C" and x.a == c and x.order == 2 and x.belement == "O":
            c.carbonyl = True

@rule
def carboxylic_acid():
    # A carbonyl carbon with an OH: acetic acid's carbon.
    for c in Atom, y in Bond:
        if c.carbonyl == True and y.a == c and y.order == 1 and y.belement == "O" and y.bhydrogens == 1:
            c.group = "carboxylic acid"

@rule
def ester():
    # A carbonyl carbon single-bonded to an oxygen that carries no hydrogen
    # and is bonded on to another carbon: methyl acetate's carbon.
    for c in Atom, o in Atom, y in Bond, z in Bond:
        if c.carbonyl == True and y.a == c and y.order == 1 and y.b == o and o.element == "O" and o.hydrogens == 0 and z.a == o and z.belement == "C" and z.b != y.a:
            c.group = "ester"

@rule
def amide():
    # A carbonyl carbon single-bonded to a nitrogen: acetamide's carbon.
    for c in Atom, y in Bond:
        if c.carbonyl == True and y.a == c and y.order == 1 and y.belement == "N":
            c.group = "amide"

@rule
def aldehyde():
    # A carbonyl carbon that carries a hydrogen, with no single bond to an
    # oxygen or a nitrogen: acetaldehyde's carbon.
    for c in Atom:
        if c.carbonyl == True and c.hydrogens >= 1 and none(y in Bond: y.a == c and y.order == 1 and y.belement == "O") and none(y in Bond: y.a == c and y.order == 1 and y.belement == "N"):
            c.group = "aldehyde"

@rule
def ketone():
    # A carbonyl carbon between two carbons: acetone's middle carbon.
    for c in Atom, y in Bond, z in Bond:
        if c.carbonyl == True and y.a == c and z.a == c and y.order == 1 and z.order == 1 and y.belement == "C" and z.belement == "C" and y.b != z.b:
            c.group = "ketone"

@rule
def alcohol():
    # An oxygen with one hydrogen, bonded to a carbon that is not a carbonyl
    # carbon: ethanol's oxygen. On a carbonyl carbon the same OH is an acid's.
    for o in Atom, y in Bond:
        if o.element == "O" and o.hydrogens == 1 and y.a == o and y.belement == "C" and absent(y.bcarbonyl):
            o.group = "alcohol"

@rule
def amine():
    # A nitrogen bonded to carbon by single bonds only, none of them to a
    # carbonyl carbon: methylamine's nitrogen. Bonded to a carbonyl carbon, it
    # is an amide's; with a double or triple bond, an imine's or a nitrile's,
    # which are not named here.
    for n in Atom, y in Bond:
        if n.element == "N" and y.a == n and y.belement == "C" and none(z in Bond: z.a == n and z.bcarbonyl == True) and none(w in Bond: w.a == n and w.order > 1):
            n.group = "amine"

# ── Standing rules: Hess's law ───────────────────────────────────────────────

@rule
def term_formation():
    for t in Term, s in Species:
        if t.species == s:
            t.formation = s.formation

@rule
def term_missing():
    # A term whose species has no enthalpy of formation recorded.
    for t in Term:
        if t.species != "" and absent(t.formation):
            t.missing = True

@rule
def reaction_terms():
    for r in Reaction, t in Term:
        if t.reaction == r and t.contribution > -1.0e300:
            r.terms = distinct(t)

@rule
def reaction_claims():
    for r in Reaction, t in Term:
        if t.reaction == r and t.contribution > -1.0e300:
            r.claims = count(t)

@rule
def reaction_disputed():
    # A sum adds every value it reads, so a term bound to the reaction twice
    # would be counted twice, which no source gives.
    for r in Reaction:
        if r.claims > r.terms:
            r.disputed = True

@rule
def reaction_enthalpy():
    # The enthalpy of a reaction is the sum, over its species, of each one's
    # coefficient times its enthalpy of formation, whatever path the reaction
    # takes: Hess's law.
    for r in Reaction, t in Term:
        if t.reaction == r and absent(r.disputed):
            r.enthalpy = sum(t.contribution)

@rule
def reaction_incomplete():
    for r in Reaction, t in Term:
        if t.reaction == r and t.missing == True:
            r.incomplete = True

# ── Writing molecules and reactions ──────────────────────────────────────────

def element(symbol: string, z: int, valence: int, weight: real):
    update(Element(symbol).z, z, "CIAAW")
    update(Element(symbol).valence, valence, "CIAAW")
    update(Element(symbol).weight, weight, "CIAAW")

def elements():
    # The elements the rules and `molar_mass()` read: the abridged standard
    # atomic weights of IUPAC's Commission on Isotopic Abundances and Atomic
    # Weights, and the usual valence in an organic molecule. Write an `Element`
    # yourself for any other.
    element("H", 1, 1, 1.008)
    element("C", 6, 4, 12.011)
    element("N", 7, 3, 14.007)
    element("O", 8, 2, 15.999)
    element("Na", 11, 0, 22.990)
    element("P", 15, 3, 30.974)
    element("S", 16, 2, 32.06)
    element("Cl", 17, 1, 35.45)
    element("K", 19, 0, 39.098)
    element("Ca", 20, 0, 40.078)
    element("Cr", 24, 0, 51.996)
    element("Mn", 25, 0, 54.938)
    element("Fe", 26, 0, 55.845)

def atom(name: string, symbol: string, molecule: string):
    # An atom of an element in a molecule.
    update(Atom(name).element, symbol, "chemistry")
    update(Atom(name).molecule, molecule, "chemistry")

def bond(a: string, b: string, order: int):
    # A bond between atoms `a` and `b`, as its two halves.
    ab = a + "-" + b
    ba = b + "-" + a
    update(Bond(ab).a, a, "chemistry")
    update(Bond(ab).b, b, "chemistry")
    update(Bond(ab).order, order, "chemistry")
    update(Bond(ba).a, b, "chemistry")
    update(Bond(ba).b, a, "chemistry")
    update(Bond(ba).order, order, "chemistry")

def species(name: string, formation: real):
    # A species and its standard enthalpy of formation, in kJ/mol.
    update(Species(name).formation, formation, "chemistry")

def term(reaction: string, species: string, coefficient: int):
    # A species in a reaction: a positive coefficient for a product, negative
    # for a reactant.
    t = reaction + "-" + species
    update(Term(t).reaction, reaction, "chemistry")
    update(Term(t).species, species, "chemistry")
    update(Term(t).coefficient, coefficient, "chemistry")

# ── 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 weight(name: string) -> real:
    # A molecule's weight, in g/mol: the least, where values disagree, and 0.0
    # when it has none. It refuses a molecule whose weight is disputed, which
    # has none until the dispute is settled.
    ws = list_values(Molecule.weight_disputed)
    for w in ws:
        parts = w.path.split(".")
        if parts[len(parts) - 2] == name:
            raise ValueError("molecule " + name + " is disputed: an atom of it, or a bond of one, is bound to it twice; forget the claim that is not to stand")
    vs = list_values(Molecule.weight)
    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 enthalpy(reaction: string) -> real:
    # A reaction's enthalpy, in kJ: the least, where values disagree, and 0.0
    # when it has none. It refuses a disputed reaction, which has none until
    # the dispute is settled.
    if reaction in disputed_reactions():
        raise ValueError("reaction " + reaction + " is disputed: a term of it is bound to it twice; forget the claim that is not to stand")
    vs = list_values(Reaction.enthalpy)
    out = 0.0
    found = False
    for v in vs:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == reaction:
            match v.value:
                case real as x:
                    if not found or x < out:
                        out = x
                        found = True
                case default:
                    pass
    return out

def formula(name: string) -> string:
    # A molecule's formula in the Hill order: carbon first, then hydrogen, then
    # the other elements alphabetically, or every element alphabetically when
    # there is no carbon. The hydrogens are the ones the rules concluded. Each
    # atom is counted once, however many sources give its element. It refuses
    # a molecule whose hydrogens are disputed, or one of whose atoms holds two
    # elements.
    hd = list_values(Molecule.hydrogens_disputed)
    for w in hd:
        parts = w.path.split(".")
        if parts[len(parts) - 2] == name:
            raise ValueError("molecule " + name + " is disputed: an atom of it, or a bond of one, is bound to it twice; forget the claim that is not to stand")
    names = {}
    vs = list_values(Atom.molecule)
    for v in vs:
        parts = v.path.split(".")
        if str(v.value) == name:
            names[parts[len(parts) - 2]] = True
    elems = {}
    vs = list_values(Atom.element)
    for v in vs:
        parts = v.path.split(".")
        at = parts[len(parts) - 2]
        if at in names:
            s = str(v.value)
            if at in elems and elems[at] != s:
                raise ValueError("molecule " + name + " is disputed: atom " + at + " holds two elements; forget the claim that is not to stand")
            elems[at] = s
    counts = {}
    for at in elems:
        s = elems[at]
        if s in counts:
            counts[s] = counts[s] + 1
        else:
            counts[s] = 1
    hs = 0
    vs = list_values(Molecule.hydrogens)
    for v in vs:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == name:
            match v.value:
                case int as x:
                    hs = x
                case default:
                    pass
    if hs > 0:
        if "H" in counts:
            counts["H"] = counts["H"] + hs
        else:
            counts["H"] = hs
    order = []
    if "C" in counts:
        order = ["C"]
        if "H" in counts:
            order = order + ["H"]
    rest = []
    for s in counts:
        if s not in order:
            rest = rest + [s]
    out = ""
    for s in order + sorted(rest):
        out = out + s
        if counts[s] != 1:
            out = out + str(counts[s])
    return out

def groups(name: string) -> list<string>:
    # The functional groups on a molecule's atoms, each named once, sorted.
    names = {}
    vs = list_values(Atom.molecule)
    for v in vs:
        parts = v.path.split(".")
        if str(v.value) == name:
            names[parts[len(parts) - 2]] = True
    out = []
    vs = list_values(Atom.group)
    for v in vs:
        parts = v.path.split(".")
        g = str(v.value)
        if parts[len(parts) - 2] in names and g not in out:
            out = out + [g]
    return sorted(out)

def disputed_molecules() -> list<string>:
    # The molecules whose weight or hydrogens are disputed, sorted.
    out = []
    for f in ["weight", "hydrogens"]:
        if f == "weight":
            vs = list_values(Molecule.weight_disputed)
        else:
            vs = list_values(Molecule.hydrogens_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_reactions() -> list<string>:
    # The reactions a term is bound to twice, sorted.
    out = []
    vs = list_values(Reaction.disputed)
    for v in vs:
        parts = v.path.split(".")
        out = out + [parts[len(parts) - 2]]
    return sorted(out)

def faults() -> list<string>:
    # Every atom bonded past its valence, every reaction with a species whose
    # enthalpy of formation is not recorded, and every disputed atom, molecule
    # and reaction, sorted.
    out = []
    vs = list_values(Atom.bonds_disputed)
    for v in vs:
        parts = v.path.split(".")
        out = out + ["disputed: atom " + parts[len(parts) - 2]]
    for m in disputed_molecules():
        out = out + ["disputed: molecule " + m]
    for r in disputed_reactions():
        out = out + ["disputed: reaction " + r]
    vs = list_values(Atom.overbonded)
    for v in vs:
        parts = v.path.split(".")
        out = out + ["overbonded: " + parts[len(parts) - 2]]
    vs = list_values(Reaction.incomplete)
    for v in vs:
        parts = v.path.split(".")
        line = "incomplete: " + parts[len(parts) - 2]
        if line not in out:
            out = out + [line]
    return sorted(out)

def unsettled() -> list<string>:
    # Each element, atom, bond, species and term that holds two different
    # values for a fact you record about it: an element's atomic number,
    # valence and weight, an atom's element and molecule, a bond's two atoms
    # and order, a species' formation and a term's reaction, species and
    # coefficient, as "element x", "atom x" and so on, sorted.
    out = []
    for f in ["z", "valence", "weight", "element", "molecule", "a", "b", "order", "formation", "reaction", "species", "coefficient"]:
        if f == "z":
            vs = list_values(Element.z)
            label = "element "
        elif f == "valence":
            vs = list_values(Element.valence)
            label = "element "
        elif f == "weight":
            vs = list_values(Element.weight)
            label = "element "
        elif f == "element":
            vs = list_values(Atom.element)
            label = "atom "
        elif f == "molecule":
            vs = list_values(Atom.molecule)
            label = "atom "
        elif f == "a":
            vs = list_values(Bond.a)
            label = "bond "
        elif f == "b":
            vs = list_values(Bond.b)
            label = "bond "
        elif f == "order":
            vs = list_values(Bond.order)
            label = "bond "
        elif f == "formation":
            vs = list_values(Species.formation)
            label = "species "
        elif f == "reaction":
            vs = list_values(Term.reaction)
            label = "term "
        elif f == "species":
            vs = list_values(Term.species)
            label = "term "
        else:
            vs = list_values(Term.coefficient)
            label = "term "
        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: formulas ────────────────────────────────────────────────────

def parse(formula: string) -> dict<string, int>:
    # How many atoms of each element a formula holds: "Ca(OH)2" is one
    # calcium, two oxygens and two hydrogens. Brackets nest. A charge written
    # after a caret, as in "MnO4^-" or "Fe^3+", is left out here and read by
    # `charge()`.
    text = formula.split("^")[0]
    n = len(text)
    kinds = []
    nums = []
    j = 0
    for step in range(n + 1):
        if j >= n:
            break
        ch = text[j]
        j = j + 1
        if ch == "(":
            kinds = kinds + ["("]
            nums = nums + [1]
            continue
        s = ch
        if ch != ")":
            for k in range(n):
                if j < n and text[j] in "abcdefghijklmnopqrstuvwxyz":
                    s = s + text[j]
                    j = j + 1
                else:
                    break
        num = 0
        for k in range(n):
            if j < n and text[j] in "0123456789":
                num = num * 10 + int(text[j])
                j = j + 1
            else:
                break
        if num == 0:
            num = 1
        kinds = kinds + [s]
        nums = nums + [num]
    # Each element's multiplier is its own count times the count after every
    # bracket that closes around it, found from the right.
    mult = []
    for k in range(len(kinds)):
        mult = mult + [1]
    outer = [1]
    for back in range(len(kinds)):
        k = len(kinds) - 1 - back
        if kinds[k] == ")":
            outer = outer + [outer[len(outer) - 1] * nums[k]]
        elif kinds[k] == "(":
            outer = outer[:len(outer) - 1]
        else:
            mult[k] = nums[k] * outer[len(outer) - 1]
    counts: dict<string, int> = {}
    for k in range(len(kinds)):
        s = kinds[k]
        if s != "(" and s != ")":
            if s in counts:
                counts[s] = counts[s] + mult[k]
            else:
                counts[s] = mult[k]
    return counts

def charge(formula: string) -> int:
    # The charge written after a caret: "^2-" is -2, "^+" is +1, and no caret
    # is 0.
    parts = formula.split("^")
    if len(parts) < 2:
        return 0
    c = parts[1]
    digits = ""
    sign = 1
    for ch in c:
        if ch in "0123456789":
            digits = digits + ch
        elif ch == "-":
            sign = -1
    if digits == "":
        return sign
    return sign * int(digits)

def molar_mass(formula: string) -> real:
    # A formula's molar mass in g/mol, from the elements' weights in the
    # store. 0.0 when an element in it has no weight recorded.
    counts = parse(formula)
    total = 0.0
    for s in counts:
        w = current(Element(s).weight)
        match w:
            case empty:
                return 0.0
            case default:
                pass
        match w.value:
            case real as x:
                total = total + x * counts[s]
            case default:
                return 0.0
    return total

# ── Calculation: balancing an equation ───────────────────────────────────────

def gcd(a: int, b: int) -> int:
    x = abs(a)
    y = abs(b)
    for step in range(200):
        if y == 0:
            break
        r = x % y
        x = y
        y = r
    return x

def balance(reactants: list<string>, products: list<string>) -> list<int>:
    # The least whole-number coefficients that balance an equation, reactants
    # first: every element, and the charge, the same on both sides. The
    # element matrix, with reactants counted positive and products negative, is
    # reduced by fraction-free Gauss–Jordan elimination, which keeps every
    # entry a whole number. An empty list when no single set of coefficients
    # balances it.
    formulas = reactants + products
    cols = len(formulas)
    symbols = []
    counts = []
    for f in formulas:
        c = parse(f)
        counts = counts + [c]
        for s in c:
            if s not in symbols:
                symbols = symbols + [s]
    rows = []
    for s in symbols:
        row = []
        for j in range(cols):
            v = 0
            if s in counts[j]:
                v = counts[j][s]
            if j >= len(reactants):
                v = -v
            row = row + [v]
        rows = rows + [row]
    crow = []
    nonzero = False
    for j in range(cols):
        q = charge(formulas[j])
        if j >= len(reactants):
            q = -q
        if q != 0:
            nonzero = True
        crow = crow + [q]
    if nonzero:
        rows = rows + [crow]
    m = len(rows)
    prev = 1
    r = 0
    pivots = []
    for c in range(cols):
        if r >= m:
            break
        p = -1
        for i in range(r, m):
            if p < 0 and rows[i][c] != 0:
                p = i
        if p < 0:
            continue
        if p != r:
            t = rows[r]
            rows[r] = rows[p]
            rows[p] = t
        piv = rows[r][c]
        for i in range(m):
            if i != r:
                fi = rows[i][c]
                row = rows[i]
                for j in range(cols):
                    row[j] = (piv * row[j] - fi * rows[r][j]) // prev
                rows[i] = row
        prev = piv
        pivots = pivots + [c]
        r = r + 1
    free = []
    for c in range(cols):
        if c not in pivots:
            free = free + [c]
    if len(free) != 1:
        return []
    fc = free[0]
    x = []
    for c in range(cols):
        x = x + [0]
    x[fc] = prev
    for i in range(len(pivots)):
        x[pivots[i]] = -rows[i][fc]
    g = 0
    for v in x:
        g = gcd(g, v)
    if g == 0:
        return []
    neg = x[fc] < 0
    out = []
    for v in x:
        w = v // g
        if neg:
            w = -w
        if w <= 0:
            return []
        out = out + [w]
    return out

# ── Calculation: acids ───────────────────────────────────────────────────────

def charge_ph(acid: real, base: real, pkas: list<real>) -> real:
    # The pH at which the charges in a solution balance: a weak acid of total
    # concentration `acid`, in mol/L, with the acid constants `pkas`, one for
    # each proton it gives up, beside `base` mol/L of a strong base's cation.
    # Water's ion product is 1.0e-14, its value at 25 °C. Found by bisection
    # on the pH between -2 and 16. A concentration below 0 refuses. Where an
    # acid constant is so large that a term of the balance passes the largest
    # real, the share of protons given up is taken from the terms' common
    # logarithms less the largest of them, scaled by a power of two so that no
    # sum of them passes the range either.
    if not (acid >= 0.0 and base >= 0.0):
        raise ValueError("a concentration is not below 0")
    kw = 1.0e-14
    kas = []
    for pk in pkas:
        kas = kas + [math.exp(-pk * math.log(10.0))]
    lo = -2.0
    hi = 16.0
    for step in range(100):
        mid = (lo + hi) / 2.0
        h = math.exp(-mid * math.log(10.0))
        terms = [1.0]
        prod = 1.0
        for j in range(len(kas)):
            prod = prod * kas[j] / h
            terms = terms + [prod]
        denom = 0.0
        given = 0.0
        for j in range(len(terms)):
            denom = denom + terms[j]
            given = given + j * terms[j]
        f = h + base - kw / h - acid * given / denom
        if not math.isfinite(f):
            sc = math.ldexp(1.0, -(len(pkas).bit_length() + 1))
            lt = [0.0]
            acc = 0.0
            for pk in pkas:
                acc = acc + (mid * sc - pk * sc)
                lt = lt + [acc]
            top = max(lt)
            denom = 0.0
            given = 0.0
            for j in range(len(lt)):
                w = math.exp((lt[j] - top) / sc * math.log(10.0))
                denom = denom + w
                given = given + j * w
            f = h + base - kw / h - acid * (given / denom)
        if f > 0.0:
            lo = mid
        else:
            hi = mid
    return (lo + hi) / 2.0

def ph(acid: real, pkas: list<real>) -> real:
    # The pH of a weak acid of concentration `acid` mol/L in water.
    return charge_ph(acid, 0.0, pkas)

def titration(acid: real, va: real, base: real, vb: real, pkas: list<real>) -> real:
    # The pH after `vb` mL of a strong base of concentration `base` mol/L is
    # added to `va` mL of a weak acid of concentration `acid`. A volume below
    # 0 refuses. Where the total volume or a concentration after dilution
    # passes the largest real, it is the concentration times the share of the
    # volume, of halves.
    if not (va >= 0.0 and vb >= 0.0):
        raise ValueError("a volume is not below 0")
    total = va + vb
    ca = acid * va / total
    cb = base * vb / total
    if not (math.isfinite(total) and math.isfinite(ca) and math.isfinite(cb)):
        half = 0.5 * va + 0.5 * vb
        ca = acid * (0.5 * va / half)
        cb = base * (0.5 * vb / half)
    return charge_ph(ca, cb, pkas)

# ── Calculation: kinetics ────────────────────────────────────────────────────

def consecutive(k1: real, k2: real, t: real, dt: real) -> list<real>:
    # A → B → C, first order with rate constants `k1` and `k2`, from pure A,
    # integrated by the fourth-order Runge–Kutta method at steps of `dt` to
    # time `t`. Returns the fractions of A, B and C. A rate constant or a time
    # below 0, or a step not above 0, refuses, and so does a step so long for
    # the rates that the integration diverges past the range of a real.
    if not (k1 >= 0.0 and k2 >= 0.0):
        raise ValueError("a rate constant is not below 0")
    if not (t >= 0.0 and dt > 0.0):
        raise ValueError("a time is not below 0 and a step is above 0")
    a = 1.0
    b = 0.0
    steps = int(t / dt + 0.5)
    for step in range(steps):
        ka1 = -k1 * a
        kb1 = k1 * a - k2 * b
        a2 = a + 0.5 * dt * ka1
        b2 = b + 0.5 * dt * kb1
        ka2 = -k1 * a2
        kb2 = k1 * a2 - k2 * b2
        a3 = a + 0.5 * dt * ka2
        b3 = b + 0.5 * dt * kb2
        ka3 = -k1 * a3
        kb3 = k1 * a3 - k2 * b3
        a4 = a + dt * ka3
        b4 = b + dt * kb3
        ka4 = -k1 * a4
        kb4 = k1 * a4 - k2 * b4
        a = a + dt * (ka1 + 2.0 * ka2 + 2.0 * ka3 + ka4) / 6.0
        b = b + dt * (kb1 + 2.0 * kb2 + 2.0 * kb3 + kb4) / 6.0
        if not (math.isfinite(a) and math.isfinite(b)):
            raise ValueError("the step is too long for the rates, and the integration diverges")
    return [a, b, 1.0 - a - b]

def consecutive_exact(k1: real, k2: real, t: real) -> list<real>:
    # The same fractions in closed form, for `k1` different from `k2`. A rate
    # constant or a time below 0 refuses.
    if not (k1 >= 0.0 and k2 >= 0.0):
        raise ValueError("a rate constant is not below 0")
    if not (t >= 0.0):
        raise ValueError("a time is not below 0")
    a = math.exp(-k1 * t)
    b = k1 / (k2 - k1) * (math.exp(-k1 * t) - math.exp(-k2 * t))
    return [a, b, 1.0 - a - b]

def peak_time(k1: real, k2: real) -> real:
    # When B is at its most: ln(k1 / k2) / (k1 - k2), and 1 / k1 when the two
    # are equal. Where k1 / k2 is not a normal real, its logarithm is
    # ln k1 - ln k2. A rate constant below 0 refuses.
    if k1 < 0.0 or k2 < 0.0:
        raise ValueError("a rate constant is not below 0")
    if k1 == k2:
        return 1.0 / k1
    q = k1 / k2
    if not (math.isfinite(q) and q >= 2.2250738585072014e-308) and k1 > 0.0 and k2 > 0.0:
        return (math.log(k1) - math.log(k2)) / (k1 - k2)
    return math.log(q) / (k1 - k2)

def equilibrium(kf: real, kr: real, a0: real, b0: real, t: real, dt: real) -> real:
    # A + B ⇌ C, forward rate constant `kf` and reverse `kr`, from `a0` and
    # `b0` mol/L of A and B and no C, integrated by the fourth-order
    # Runge–Kutta method to time `t`. Returns the concentration of C.
    c = 0.0
    steps = int(t / dt + 0.5)
    for step in range(steps):
        k1 = kf * (a0 - c) * (b0 - c) - kr * c
        c2 = c + 0.5 * dt * k1
        k2 = kf * (a0 - c2) * (b0 - c2) - kr * c2
        c3 = c + 0.5 * dt * k2
        k3 = kf * (a0 - c3) * (b0 - c3) - kr * c3
        c4 = c + dt * k3
        k4 = kf * (a0 - c4) * (b0 - c4) - kr * c4
        c = c + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
    return c

def equilibrium_exact(k: real, a0: real, b0: real) -> real:
    # The concentration of C at equilibrium, where C / (A B) equals the
    # equilibrium constant `k` = kf / kr: the smaller root of
    # k x² − (k (a0 + b0) + 1) x + k a0 b0 = 0. It is taken as
    # 2 k a0 b0 / (p + r), p = k a0 + k b0 + 1 and r the root of
    # (k a0 − k b0)² + 2 (k a0 + k b0) + 1, which is the same root and
    # subtracts nothing, where the textbook (p − r) / 2k lost every digit of
    # a small one; and where k times the larger concentration is above 1, the
    # same divided through by k and by 8, so no step passes the range of a
    # real. A constant or a concentration below 0 refuses; a constant of
    # 0 gives 0.
    if not (k >= 0.0):
        raise ValueError("an equilibrium constant is not below 0")
    if not (a0 >= 0.0 and b0 >= 0.0):
        raise ValueError("a concentration is not below 0")
    big = max(a0, b0)
    small = min(a0, b0)
    if k * big <= 1.0:
        ka = k * a0
        kb = k * b0
        p = ka + kb + 1.0
        r = math.hypot(ka - kb, math.sqrt(2.0 * (ka + kb) + 1.0))
        return 2.0 * (k * big) * small / (p + r)
    h = 0.5 * a0 + 0.5 * b0
    ik = 1.0 / k
    r = math.hypot(0.125 * a0 - 0.125 * b0, 0.25 * (math.sqrt(h) / math.sqrt(k)), 0.125 * ik)
    return small * (0.25 * big / (0.25 * h + 0.125 * ik + r))