# SPDX-License-Identifier: MIT

# biology — read metabolic networks and pedigrees, and calculate inheritance,
# populations and gene regulation.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply two kinds of record as facts.
#
# A metabolic network: compounds (`biology.Compound`), with the ones the growth
# medium supplies marked, and reactions (`biology.Reaction`), each with one or
# two substrates and one or two products. The rules conclude which reactions
# can run and which compounds the network can make from the medium, which is
# its scope. `needs()` asks which one compound, added to the medium, would let
# the network make a compound it cannot, and writes nothing to find out:
#
#     import biology
#     biology.medium("glucose")
#     biology.medium("atp")
#     biology.reaction("hexokinase", "glucose", "atp", "g6p", "adp")
#     biology.reaction("isomerase", "g6p", "", "f6p", "")
#     print(biology.scope())
#     print(biology.needs("x5p", ["xylose", "ribose"]))
#
# A pedigree: people (`biology.Person`), each with a family, a mother and a
# father where known, a sex and whether they are affected, and each family's
# mode of inheritance (`biology.Family`). The rules conclude every person's
# ancestors, a child whose parents are related by descent, the carriers the
# mode requires, and each person whose record the mode cannot explain:
#
#     biology.family("tudor", "X-linked recessive")
#     biology.person("anne", "tudor", "F", False, "", "")
#     biology.person("hal", "tudor", "M", True, "anne", "")
#     print(biology.carriers("tudor"), biology.unexplained("tudor"))
#
# Each conclusion answers `why`, and withdraws itself when a fact it rests on
# is corrected. A compound that only a cycle of its own products would make is
# not reached once the medium stops supplying the cycle. Where two sources
# disagree about a fact, both values are reasoned about, no query refuses, and
# `unsettled()` names what they disagree about. A family counts each member
# the mode cannot explain once, however many rules or values find them so; a
# person whose family two sources disagree on is counted in each family a
# value names.
#
# CALCULATION. Value scripts that return a value and write nothing:
# transcription and translation in the standard genetic code and the
# vertebrate mitochondrial one, a Mendelian cross and the chi-square test of
# its counts, Hardy–Weinberg proportions, logistic growth, predator and prey,
# genetic drift, and a ring of genes that repress one another.

import math

# ── Metabolic networks as facts ──────────────────────────────────────────────

class Compound:
    # A compound, by a name you choose. `medium` is True for one the growth
    # medium supplies. The rules conclude `reached` for each compound the
    # network can make from the medium, the medium's own included.
    medium: bool
    reached: bool

class Reaction:
    # A reaction, by a name you choose: `s1` and `s2` its substrates and `p1`
    # and `p2` its products, each the name of a compound. Leave `s2` or `p2`
    # unwritten for a reaction with one. The rules conclude `runs` when every
    # substrate is reached.
    s1: string
    s2: string
    p1: string
    p2: string
    runs: bool

# ── Pedigrees as facts ───────────────────────────────────────────────────────

class Family:
    # A family, by a name you choose, and `mode`, how the trait its members are
    # recorded for is inherited: "autosomal dominant", "autosomal recessive"
    # or "X-linked recessive". The rules conclude `unexplained`, how many of
    # its members the mode cannot explain, each counted once however many
    # rules or values find it so.
    mode: string
    unexplained: int

class Person:
    # A person, by a name you choose: `family`, `mother` and `father` (each a
    # person's name, where known), `sex` "F" or "M", and `affected`. The rules
    # conclude `inbred` for a child whose parents are related by descent,
    # `carrier` for an unaffected person the mode requires to carry the
    # allele, and `unexplained` for a person whose record the mode cannot
    # explain. Every one of these assumes the trait shows in everyone who has
    # its genotype, and that no new mutation arose in the family: an
    # unexplained person is where one of those, or a record, does not hold.
    family: string
    mother: string
    father: string
    sex: string
    affected: bool
    inbred: bool
    carrier: bool
    unexplained: bool

class Parent:
    # The relation `Parent.child.parent`, True for each parent recorded.
    pass

class Ancestor:
    # The relation `Ancestor.person.ancestor`, True for each parent, each
    # parent's parent, and so on to the first generation recorded.
    pass

# ── Standing rules: a network's scope ────────────────────────────────────────

@rule
def supplied():
    for c in Compound:
        if c.medium == True:
            c.reached = True

@rule
def runs_one():
    for r in Reaction, a in Compound:
        if r.s1 == a and a.reached == True and absent(r.s2):
            r.runs = True

@rule
def runs_two():
    for r in Reaction, a in Compound, b in Compound:
        if r.s1 == a and a.reached == True and r.s2 == b and b.reached == True:
            r.runs = True

@rule
def makes_first():
    for r in Reaction, c in Compound:
        if r.runs == True and r.p1 == c:
            c.reached = True

@rule
def makes_second():
    for r in Reaction, c in Compound:
        if r.runs == True and r.p2 == c:
            c.reached = True

# ── Standing rules: descent ──────────────────────────────────────────────────

@rule
def parent_mother():
    for c in Person, m in Person:
        if c.mother == m:
            Parent.c.m = True

@rule
def parent_father():
    for c in Person, f in Person:
        if c.father == f:
            Parent.c.f = True

@rule
def ancestor_parent():
    for x in Person, y in Person:
        if Parent.x.y == True:
            Ancestor.x.y = True

@rule
def ancestor_further():
    for x in Person, y in Person, z in Person:
        if Ancestor.x.y == True and Parent.y.z == True:
            Ancestor.x.z = True

@rule
def inbred_shared():
    # Parents with an ancestor in common: cousins, siblings, an uncle and a
    # niece.
    for c in Person, m in Person, f in Person, a in Person:
        if c.mother == m and c.father == f and Ancestor.m.a == True and Ancestor.f.a == True:
            c.inbred = True

@rule
def inbred_father_line():
    # A father who is the mother's own ancestor.
    for c in Person, m in Person, f in Person:
        if c.mother == m and c.father == f and Ancestor.m.f == True:
            c.inbred = True

@rule
def inbred_mother_line():
    # A mother who is the father's own ancestor.
    for c in Person, m in Person, f in Person:
        if c.mother == m and c.father == f and Ancestor.f.m == True:
            c.inbred = True

# ── Standing rules: carriers ─────────────────────────────────────────────────

@rule
def recessive_parent():
    # Each unaffected parent of an affected child gave the child one copy.
    for y in Family, c in Person, p in Person:
        if y.mode == "autosomal recessive" and c.family == y and c.affected == True and Parent.c.p == True and p.affected == False:
            p.carrier = True

@rule
def recessive_child():
    # Each unaffected child of an affected parent was given one copy.
    for y in Family, c in Person, p in Person:
        if y.mode == "autosomal recessive" and c.family == y and c.affected == False and Parent.c.p == True and p.affected == True:
            c.carrier = True

@rule
def x_linked_mother():
    # A son's one X is his mother's, so the unaffected mother of an affected
    # son carries the allele.
    for y in Family, c in Person, m in Person:
        if y.mode == "X-linked recessive" and c.family == y and c.sex == "M" and c.affected == True and c.mother == m and m.affected == False:
            m.carrier = True

@rule
def x_linked_daughter():
    # A father gives every daughter his one X.
    for y in Family, c in Person, f in Person:
        if y.mode == "X-linked recessive" and c.family == y and c.sex == "F" and c.affected == False and c.father == f and f.affected == True:
            c.carrier = True

# ── Standing rules: what a mode cannot explain ───────────────────────────────

@rule
def dominant_unexplained():
    # An affected child of two unaffected parents: neither had an allele to
    # give.
    for y in Family, c in Person, m in Person, f in Person:
        if y.mode == "autosomal dominant" and c.family == y and c.affected == True and c.mother == m and m.affected == False and c.father == f and f.affected == False:
            c.unexplained = True

@rule
def recessive_unexplained():
    # An unaffected child of two affected parents: both had only the allele
    # to give.
    for y in Family, c in Person, m in Person, f in Person:
        if y.mode == "autosomal recessive" and c.family == y and c.affected == False and c.mother == m and m.affected == True and c.father == f and f.affected == True:
            c.unexplained = True

@rule
def x_linked_daughter_unexplained():
    # An affected daughter has her father's X, so he must be affected too.
    for y in Family, c in Person, f in Person:
        if y.mode == "X-linked recessive" and c.family == y and c.sex == "F" and c.affected == True and c.father == f and f.affected == False:
            c.unexplained = True

@rule
def x_linked_son_unexplained():
    # An affected mother has the allele on both her X chromosomes, so every
    # son is affected.
    for y in Family, c in Person, m in Person:
        if y.mode == "X-linked recessive" and c.family == y and c.sex == "M" and c.affected == False and c.mother == m and m.affected == True:
            c.unexplained = True

@rule
def family_unexplained():
    for y in Family, c in Person:
        if c.family == y and c.unexplained == True:
            y.unexplained = distinct(c)

# ── Writing facts ────────────────────────────────────────────────────────────

def medium(name: string):
    # A compound the growth medium supplies.
    update(Compound(name).medium, True, "medium")

def reaction(name: string, s1: string, s2: string, p1: string, p2: string):
    # A reaction from `s1` and `s2` to `p1` and `p2`; pass "" for a second
    # substrate or product the reaction does not have.
    update(Reaction(name).s1, s1, "network")
    update(Reaction(name).p1, p1, "network")
    if s2 != "":
        update(Reaction(name).s2, s2, "network")
    if p2 != "":
        update(Reaction(name).p2, p2, "network")

def family(name: string, mode: string):
    # A family and the mode of inheritance its trait is recorded under.
    update(Family(name).mode, mode, "pedigree")

def person(name: string, family: string, sex: string, affected: bool, mother: string, father: string):
    # A person; pass "" for a mother or father not recorded.
    update(Person(name).family, family, "pedigree")
    update(Person(name).sex, sex, "pedigree")
    update(Person(name).affected, affected, "pedigree")
    if mother != "":
        update(Person(name).mother, mother, "pedigree")
    if father != "":
        update(Person(name).father, father, "pedigree")

# ── 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.

def scope() -> list<string>:
    # Every compound the network can make from its medium, sorted.
    out = []
    vs = list_values(Compound.reached)
    for v in vs:
        parts = v.path.split(".")
        c = parts[len(parts) - 2]
        if c not in out:
            out = out + [c]
    return sorted(out)

def blocked() -> list<string>:
    # Every reaction that cannot run on what the medium makes possible,
    # sorted.
    running = {}
    vs = list_values(Reaction.runs)
    for v in vs:
        parts = v.path.split(".")
        running[parts[len(parts) - 2]] = True
    out = []
    vs = list_values(Reaction.s1)
    for v in vs:
        parts = v.path.split(".")
        r = parts[len(parts) - 2]
        if r not in running and r not in out:
            out = out + [r]
    return sorted(out)

def needs(target: string, candidates: list<string>) -> list<string>:
    # Each candidate that, added to the medium alone, would let the network
    # make `target`, in the order given. Each is a hypothesis: nothing is
    # written, and the medium is left as it was.
    out = []
    for c in candidates:
        try:
            v = suppose(Compound.c.medium, True, Compound.target.reached)
            if str(v) == "True":
                out = out + [c]
        except refuse:
            pass
    return out

def people(name: string) -> list<string>:
    # Every person recorded in a family, sorted.
    out = []
    vs = list_values(Person.family)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if str(v.value) == name and p not in out:
            out = out + [p]
    return sorted(out)

def ancestors(name: string) -> list<string>:
    # Every recorded ancestor of a person, sorted.
    everyone = []
    vs = list_values(Person.family)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if p not in everyone:
            everyone = everyone + [p]
    vs = list_values(Person.mother)
    for v in vs:
        p = str(v.value)
        if p not in everyone:
            everyone = everyone + [p]
    vs = list_values(Person.father)
    for v in vs:
        p = str(v.value)
        if p not in everyone:
            everyone = everyone + [p]
    out = []
    for a in everyone:
        r = current(Ancestor.name.a)
        match r:
            case empty:
                pass
            case default:
                out = out + [a]
    return sorted(out)

def inbred() -> list<string>:
    # Every child whose parents are related by descent, sorted.
    out = []
    vs = list_values(Person.inbred)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if p not in out:
            out = out + [p]
    return sorted(out)

def carriers(name: string) -> list<string>:
    # Every member of a family the mode requires to carry the allele, sorted.
    members = people(name)
    out = []
    vs = list_values(Person.carrier)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if p in members and p not in out:
            out = out + [p]
    return sorted(out)

def unexplained(name: string) -> list<string>:
    # Every member of a family whose record its mode cannot explain, sorted.
    members = people(name)
    out = []
    vs = list_values(Person.unexplained)
    for v in vs:
        parts = v.path.split(".")
        p = parts[len(parts) - 2]
        if p in members and p not in out:
            out = out + [p]
    return sorted(out)

def unsettled() -> list<string>:
    # Each compound, reaction, family and person that holds two different
    # values for a fact you write, as "compound x medium", "person x mother"
    # and so on, sorted.
    out = []
    for f in ["medium", "s1", "s2", "p1", "p2", "mode", "family", "mother", "father", "sex", "affected"]:
        if f == "medium":
            vs = list_values(Compound.medium)
            label = "compound "
        elif f == "s1":
            vs = list_values(Reaction.s1)
            label = "reaction "
        elif f == "s2":
            vs = list_values(Reaction.s2)
            label = "reaction "
        elif f == "p1":
            vs = list_values(Reaction.p1)
            label = "reaction "
        elif f == "p2":
            vs = list_values(Reaction.p2)
            label = "reaction "
        elif f == "mode":
            vs = list_values(Family.mode)
            label = "family "
        elif f == "family":
            vs = list_values(Person.family)
            label = "person "
        elif f == "mother":
            vs = list_values(Person.mother)
            label = "person "
        elif f == "father":
            vs = list_values(Person.father)
            label = "person "
        elif f == "sex":
            vs = list_values(Person.sex)
            label = "person "
        else:
            vs = list_values(Person.affected)
            label = "person "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = label + parts[len(parts) - 2] + " " + f
            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 genetic code ────────────────────────────────────────────

def complement(dna: string) -> string:
    # The reverse complement of a DNA strand: the other strand, read in its
    # own direction.
    pair = {"A": "T", "T": "A", "C": "G", "G": "C"}
    out = ""
    n = len(dna)
    for k in range(n):
        out = out + pair[dna[n - 1 - k]]
    return out

def transcribe(dna: string) -> string:
    # The messenger RNA a gene's coding strand is copied into: the same
    # sequence with uracil for thymine.
    return dna.replace("T", "U")

def translate(rna: string, table: int) -> string:
    # The protein an mRNA codes for, in one-letter amino acids, read from the
    # first AUG to the first stop codon, in NCBI translation table 1 (the
    # standard code) or 2 (vertebrate mitochondria). "" when there is no AUG.
    code = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
    if table == 2:
        code = "FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIMMTTTTNNKKSS**VVVVAAAADDEEGGGG"
    base = {"U": 0, "C": 1, "A": 2, "G": 3}
    start = rna.find("AUG")
    if start < 0:
        return ""
    out = ""
    k = start
    for step in range(len(rna)):
        if k + 3 > len(rna):
            break
        aa = code[16 * base[rna[k]] + 4 * base[rna[k + 1]] + base[rna[k + 2]]]
        if aa == "*":
            break
        out = out + aa
        k = k + 3
    return out

def gc(dna: string) -> real:
    # The share of a strand's bases that are G or C.
    if len(dna) == 0:
        return 0.0
    return (dna.count("G") + dna.count("C")) / len(dna)

# ── Calculation: crosses ─────────────────────────────────────────────────────

def gametes(genotype: string) -> list<string>:
    # The gametes a genotype written two letters a gene ("AaBb") makes, one
    # allele of each gene, every combination once.
    out = [""]
    for g in range(len(genotype) // 2):
        next = []
        for part in out:
            next = next + [part + genotype[2 * g], part + genotype[2 * g + 1]]
        out = next
    return out

def zygote(a: string, b: string) -> string:
    # The genotype two gametes make, each gene's dominant allele first.
    out = ""
    for g in range(len(a)):
        x = a[g]
        y = b[g]
        if y < x:
            out = out + y + x
        else:
            out = out + x + y
    return out

def cross(a: string, b: string) -> dict<string, int>:
    # How many of the offspring of every pairing of two parents' gametes have
    # each genotype.
    counts: dict<string, int> = {}
    for x in gametes(a):
        for y in gametes(b):
            z = zygote(x, y)
            if z in counts:
                counts[z] = counts[z] + 1
            else:
                counts[z] = 1
    return counts

def phenotypes(a: string, b: string) -> dict<string, int>:
    # The same offspring by what they show: each gene's dominant letter where
    # the genotype has one, its recessive letter where it has two.
    counts: dict<string, int> = {}
    genotypes = cross(a, b)
    for z in genotypes:
        p = ""
        for g in range(len(z) // 2):
            p = p + z[2 * g]
        if p in counts:
            counts[p] = counts[p] + genotypes[z]
        else:
            counts[p] = genotypes[z]
    return counts

def chi_square(observed: list<real>, ratio: list<real>) -> real:
    # Pearson's statistic for counts against the ratio expected of them: the
    # sum over classes of (observed - expected)^2 / expected. A class the ratio
    # expects none of adds nothing where none were observed in it, and makes
    # the statistic infinite where some were. Where a sum, an expected count
    # or a square on the way passes the range of a real, or falls below the
    # smallest normal one, while the counts and the shares are reals, it is
    # taken again with the counts and the shares each scaled by a power of
    # two that brings the largest near 1, which leaves every expected share
    # as it was and scales the statistic with the counts, and the statistic
    # is scaled back.
    for k in range(len(observed)):
        if not (observed[k] >= 0.0):
            raise ValueError("an observed count is not below 0")
        if not (ratio[k] >= 0.0):
            raise ValueError("a share of the ratio is not below 0")
    for k in range(len(observed)):
        if ratio[k] == 0.0 and observed[k] > 0.0:
            return math.inf
    obs = observed
    shares = ratio
    shift = 0
    out = 0.0
    for attempt in range(2):
        total = 0.0
        parts = 0.0
        for k in range(len(obs)):
            total = total + obs[k]
            parts = parts + shares[k]
        out = 0.0
        lost = False
        for k in range(len(obs)):
            if ratio[k] == 0.0:
                continue
            e = total * shares[k] / parts
            if abs(e) < 2.2250738585072014e-308:
                lost = True
            if attempt == 0 and not lost:
                d = obs[k] - e
                sq = d * d
                if d != 0.0 and abs(sq) < 2.2250738585072014e-308:
                    lost = True
                else:
                    out = out + sq / e
            if attempt == 1:
                d = obs[k] - e
                sq = d * d
                if d != 0.0 and abs(sq) < 2.2250738585072014e-308:
                    out = out + d * (d / e)
                else:
                    out = out + sq / e
        if attempt == 1:
            return math.ldexp(out, shift)
        if (math.isfinite(out) and not lost) or not all(math.isfinite(v) for v in observed) or not all(math.isfinite(v) for v in ratio):
            return out
        top = 0.0
        for v in observed:
            if v > top:
                top = v
        big = 0.0
        for v in ratio:
            if v > big:
                big = v
        shift = 0
        if top > 0.0:
            shift = math.frexp(top)[1] + len(observed).bit_length()
        rshift = math.frexp(big)[1] + len(ratio).bit_length()
        obs = [math.ldexp(v, -shift) for v in observed]
        shares = [math.ldexp(v, -rshift) for v in ratio]
    return out

def hardy_weinberg(aa: int, ab: int, bb: int) -> list<real>:
    # From the counts of the three genotypes of one gene: the frequency of the
    # first allele, and the chi-square of the counts against the proportions
    # random mating gives, p^2 : 2pq : q^2.
    n = aa + ab + bb
    p = (2 * aa + ab) / (2 * n)
    q = 1.0 - p
    return [p, chi_square([aa * 1.0, ab * 1.0, bb * 1.0], [p * p, 2.0 * p * q, q * q])]

# ── Calculation: populations ─────────────────────────────────────────────────

def logistic(r: real, k: real, n0: real, t: real, dt: real) -> real:
    # A population growing at rate r toward a carrying capacity k, dN/dt =
    # rN(1 - N/k), from n0 to time t by fourth-order Runge–Kutta in steps of
    # dt. Where a step passes the range of a real while the inputs are reals,
    # the steps are taken again with n0 and k scaled by a power of two that
    # brings the larger near 1, which scales every step by the same power
    # and changes no rounding, and the population is scaled back.
    if not (k > 0.0):
        raise ValueError("a carrying capacity is above 0")
    if not (n0 >= 0.0):
        raise ValueError("a population is not below 0")
    shift = 0
    n = n0
    for attempt in range(2):
        n = math.ldexp(n0, -shift)
        cap = math.ldexp(k, -shift)
        for step in range(int(t / dt + 0.5)):
            k1 = r * n * (1.0 - n / cap)
            n2 = n + 0.5 * dt * k1
            k2 = r * n2 * (1.0 - n2 / cap)
            n3 = n + 0.5 * dt * k2
            k3 = r * n3 * (1.0 - n3 / cap)
            n4 = n + dt * k3
            k4 = r * n4 * (1.0 - n4 / cap)
            n = n + dt * (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
        if attempt == 1 or math.isfinite(n) or not (math.isfinite(r) and math.isfinite(k) and math.isfinite(n0) and math.isfinite(t) and math.isfinite(dt)):
            return math.ldexp(n, shift)
        shift = math.frexp(max(k, n0))[1] + 2
    return n

def logistic_exact(r: real, k: real, n0: real, t: real) -> real:
    # The same population, from the equation's solution. A population that
    # starts at 0 or at k stays there. Where the factor (k - n0) / n0, the
    # exponential or their product passes the range of a real, or the
    # exponential falls below the smallest normal one, while the inputs are
    # reals, each is carried as a fraction and a power of two: the exponential
    # of an argument past 700 as the square of the exponential of its half,
    # and of one past 4000, where the population is already 0 or k, as of
    # 4000.
    if not (k > 0.0):
        raise ValueError("a carrying capacity is above 0")
    if not (n0 >= 0.0):
        raise ValueError("a population is not below 0")
    if n0 == 0.0:
        return 0.0
    q = (k - n0) / n0
    if q == 0.0:
        return k
    g = math.exp(-r * t)
    if (math.isfinite(q) and math.isfinite(g) and g >= 2.2250738585072014e-308 and math.isfinite(q * g)) or not (math.isfinite(r) and math.isfinite(k) and math.isfinite(n0) and math.isfinite(t)):
        return k / (1.0 + q * g)
    fq = math.frexp(k - n0)
    fn = math.frexp(n0)
    x = max(-4000.0, min(4000.0, -r * t))
    halves = 0
    for i in range(3):
        if abs(x) > 700.0:
            x = x / 2.0
            halves = halves + 1
    fg = math.frexp(math.exp(x))
    gm = fg[0]
    ge = fg[1]
    for i in range(halves):
        f = math.frexp(gm * gm)
        gm = f[0]
        ge = 2 * ge + f[1]
    f = math.frexp(fq[0] / fn[0] * gm)
    pe = fq[1] - fn[1] + ge + f[1]
    if pe > 1000:
        fk = math.frexp(k)
        return math.ldexp(fk[0] / f[0], fk[1] - pe)
    return k / (1.0 + math.ldexp(f[0], pe))

def predator_prey(alpha: real, beta: real, gamma: real, delta: real, x0: real, y0: real, t: real, dt: real) -> list<real>:
    # Prey x and predators y, dx/dt = alpha x - beta x y and dy/dt = delta x y
    # - gamma y, from x0 and y0 to time t by fourth-order Runge–Kutta: the two
    # populations, and the relative change in the quantity the exact motion
    # keeps, delta x - gamma ln x + beta y - alpha ln y. Where a step passes
    # the range of a real while the inputs are reals, the steps are taken
    # again with x0 and y0 scaled by powers of two that bring each near 1,
    # and delta and beta by the same powers the other way, which scales every
    # step and changes no rounding, and the populations are scaled back;
    # populations that pass the range even so are refused.
    v0 = delta * x0 - gamma * math.log(x0) + beta * y0 - alpha * math.log(y0)
    lx0 = math.log(x0)
    ly0 = math.log(y0)
    whole = all(math.isfinite(c) for c in [alpha, beta, gamma, delta, x0, y0, t, dt])
    sx = 0
    sy = 0
    x = x0
    y = y0
    for attempt in range(2):
        x = math.ldexp(x0, -sx)
        y = math.ldexp(y0, -sy)
        b = math.ldexp(beta, sy)
        d = math.ldexp(delta, sx)
        for step in range(int(t / dt + 0.5)):
            ax1 = alpha * x - b * x * y
            ay1 = d * x * y - gamma * y
            x2 = x + 0.5 * dt * ax1
            y2 = y + 0.5 * dt * ay1
            ax2 = alpha * x2 - b * x2 * y2
            ay2 = d * x2 * y2 - gamma * y2
            x3 = x + 0.5 * dt * ax2
            y3 = y + 0.5 * dt * ay2
            ax3 = alpha * x3 - b * x3 * y3
            ay3 = d * x3 * y3 - gamma * y3
            x4 = x + dt * ax3
            y4 = y + dt * ay3
            ax4 = alpha * x4 - b * x4 * y4
            ay4 = d * x4 * y4 - gamma * y4
            x = x + dt * (ax1 + 2.0 * ax2 + 2.0 * ax3 + ax4) / 6.0
            y = y + dt * (ay1 + 2.0 * ay2 + 2.0 * ay3 + ay4) / 6.0
        x = math.ldexp(x, sx)
        y = math.ldexp(y, sy)
        if (math.isfinite(x) and math.isfinite(y)) or not whole:
            break
        if attempt == 1:
            raise ValueError("the populations pass the largest real")
        sx = math.frexp(x0)[1]
        sy = math.frexp(y0)[1]
    v = delta * x - gamma * math.log(x) + beta * y - alpha * math.log(y)
    change = (v - v0) / v0
    if not math.isfinite(change) and whole:
        # The quantity is linear in the four rates, so it is taken again with
        # them scaled by a power of two that keeps each of its four terms an
        # eighth of the largest real or less, which leaves the relative
        # change as it was.
        e = 1021 - math.frexp(max(abs(alpha), abs(beta), abs(gamma), abs(delta)))[1] - math.frexp(max(abs(x0), abs(y0), abs(x), abs(y), 1024.0))[1]
        a = math.ldexp(alpha, e)
        b = math.ldexp(beta, e)
        c = math.ldexp(gamma, e)
        d = math.ldexp(delta, e)
        u0 = d * x0 - c * lx0 + b * y0 - a * ly0
        u = d * x - c * math.log(x) + b * y - a * math.log(y)
        change = (u - u0) / u0
    return [x, y, change]

def drift(n: int, p0: real, generations: int, seed: int) -> list<real>:
    # One population of n diploid individuals under the Wright–Fisher model:
    # each generation draws its 2n gene copies at random from the last one's.
    # The draws come from the Park–Miller generator from `seed`, so a seed
    # gives the same run every time. Returns the allele's frequency at the
    # end, and the generation it was fixed or lost in, -1.0 if neither.
    x = seed
    p = p0
    ended = -1.0
    for g in range(generations):
        k = 0
        for draw in range(2 * n):
            x = 16807 * x % 2147483647
            if x / 2147483647 < p:
                k = k + 1
        p = k / (2 * n)
        if k == 0 or k == 2 * n:
            ended = g + 1.0
            break
    return [p, ended]

def fixation(n: int, p0: real, runs: int, generations: int, seed: int) -> real:
    # The share of `runs` populations, drawn one after another from one
    # Park–Miller stream, in which the allele is fixed within `generations`.
    # Every draw is a value the run builds, so runs, generations and 2n
    # multiply toward a run's limit.
    x = seed
    fixed = 0
    for run in range(runs):
        p = p0
        for g in range(generations):
            k = 0
            for draw in range(2 * n):
                x = 16807 * x % 2147483647
                if x / 2147483647 < p:
                    k = k + 1
            p = k / (2 * n)
            if k == 0 or k == 2 * n:
                break
        if p == 1.0:
            fixed = fixed + 1
    return fixed / runs

# ── Calculation: a ring of repressors ────────────────────────────────────────

def ring(start: string, steps: int) -> list<string>:
    # A ring of genes in which each represses the next, written one character
    # a gene, "1" for on: at each step, all at once, a gene is on exactly when
    # the gene before it in the ring was off. Returns the states from `start`
    # through `steps` steps.
    out = [start]
    state = start
    n = len(start)
    for step in range(steps):
        next = ""
        for g in range(n):
            if state[(g + n - 1) % n] == "1":
                next = next + "0"
            else:
                next = next + "1"
        state = next
        out = out + [state]
    return out

def period(start: string) -> int:
    # How many steps the ring takes to return to a state it has been in, from
    # `start`: 1 for a state that stays as it is.
    states = ring(start, 2 ** len(start))
    seen: dict<string, int> = {}
    for k in range(len(states)):
        s = states[k]
        if s in seen:
            return k - seen[s]
        seen[s] = k
    return 0