# SPDX-License-Identifier: MIT

# mathematics — classify numbers, finite operations and triangles, and
# calculate number theory and numerical methods.
#
# Two halves that share nothing but the subject.
#
# STRUCTURE. You supply three kinds of record as facts.
#
# Whole numbers (`mathematics.Number`), each written with a row for every
# candidate divisor up to half of it (`mathematics.Candidate`). The rules
# conclude which candidates divide it, the sum of its proper divisors, and
# whether it is prime or composite, and perfect, abundant or deficient:
#
#     import mathematics
#     mathematics.numbers(1, 30)
#     print(mathematics.numbers_that_are("prime"))
#     print(mathematics.numbers_that_are("perfect"), mathematics.divisors(28))
#
# A binary operation on a finite set, given as its table (`Operation`,
# `Element`, `Cell`, and a `Triple` for each three elements in order). The
# rules conclude its identity, whether it is closed, commutative and
# associative, which elements have a two-sided inverse, and whether it is a
# group, and an abelian one. Where a property fails, the cells that break it
# are named:
#
#     els = ["0", "1", "2", "3"]
#     add = [[str((i + j) % 4) for j in range(4)] for i in range(4)]
#     sub = [[str((i - j) % 4) for j in range(4)] for i in range(4)]
#     mathematics.operation("z4", els, add)
#     mathematics.operation("minus4", els, sub)
#     print(mathematics.properties("z4"), mathematics.properties("minus4"))
#     print(mathematics.nonassociative_triples("minus4")[0])
#
# Triangles by their three sides in whole units (`Triangle`). The rules
# conclude whether three lengths make a triangle at all, and if they do,
# whether it is acute, right or obtuse, and equilateral, isosceles or
# scalene:
#
#     mathematics.triangle("egypt", 3, 4, 5)
#     mathematics.triangle("none", 1, 2, 5)
#     print(mathematics.kinds("egypt"), mathematics.kinds("none"))
#
# Each conclusion answers `why`, and `why_not` names the premise that failed.
# A corrected fact withdraws what rested on it: a table cell rewritten breaks
# the properties that cell held up. Where two sources disagree about a fact,
# both values are reasoned about, no query refuses, and `unsettled()` names
# what they disagree about.
#
# The divisor sum is the exception, since a sum adds every value it reads: a
# candidate whose `k` a second source writes, even the same `k`, would be added
# twice, and 6 would read as abundant, which no source says. So while a
# divisor is bound to its number twice, as a second `k` binds it, the number
# is `disputed`: no divisor sum, `perfect`, `abundant` or `deficient` is
# concluded for it, `numbers_that_are("disputed")` names it, and `why_not` on
# its divisor sum names the dispute. Forgetting the claim that is not to stand
# settles it. A candidate whose number two sources disagree on is added to each
# number a value names, so each of those divisor sums is one source's.
#
# CALCULATION. Value scripts that return a value and write nothing: Euclid's
# algorithm and its extension, modular inverses and powers, a primality test,
# factorization, Euler's totient, the sieve, binomial coefficients,
# factorials and Fibonacci numbers, Collatz trajectories, continued fractions
# and Pell's equation, roots by Newton's method, integrals from samples by
# the trapezoidal rule and Simpson's, and determinants and linear systems by
# elimination.
#
# An `int` is 64 bits. A result past 9,223,372,036,854,775,807 refuses as
# `arithmetic_overflow` rather than wrapping, so 20! is computed and 21! is
# refused, and each script below says how far it reaches.

import math

# ── Numbers as facts ─────────────────────────────────────────────────────────

class Number:
    # A whole number, by a name you choose, and `n` its value, at least 1.
    # `number()` writes it, named "n" and its value, with its candidates. The
    # rules conclude `divisor_sum`, the sum of its proper divisors (each
    # divisor but itself), `prime` for a number above 1 that nothing between
    # 1 and itself divides and `composite` for one that something does, and
    # `perfect`, `abundant` or `deficient` as its divisor sum equals, exceeds
    # or falls short of it. 1 has no proper divisor, so it has no divisor
    # sum, is neither prime nor composite, and is deficient. The rules also
    # conclude `claims` and `divisor_count`, which `disputed` compares: a
    # divisor bound to the number twice makes the number `disputed`, and then
    # no divisor sum is concluded.
    n: int
    claims: int
    divisor_count: int
    disputed: bool
    divisor_sum: int
    prime: bool
    composite: bool
    perfect: bool
    abundant: bool
    deficient: bool

class Candidate:
    # One candidate divisor `k` of the number named `of`, whose value is `n`,
    # and `remainder`, what is left when `k` divides `n`. `number()` writes
    # one for each `k` from 1 to half of `n`, since no proper divisor is
    # larger. The rules conclude `divides` where the remainder is 0.
    of: string
    n: int
    k: int
    remainder: int = n % k
    divides: bool

@rule
def divides():
    for d in Candidate:
        if d.remainder == 0:
            d.divides = True

@rule
def number_claims():
    for x in Number, d in Candidate:
        if d.of == x and d.divides == True and d.k > -1000000000000000000:
            x.claims = count(d)

@rule
def number_divisors():
    for x in Number, d in Candidate:
        if d.of == x and d.divides == True and d.k > -1000000000000000000:
            x.divisor_count = distinct(d)

@rule
def number_disputed():
    # A sum adds every value it reads, so a divisor bound to the number twice
    # would be added twice, which no source gives.
    for x in Number:
        if x.claims > x.divisor_count:
            x.disputed = True

@rule
def divisor_sum():
    for x in Number, d in Candidate:
        if d.of == x and d.divides == True and absent(x.disputed):
            x.divisor_sum = sum(d.k)

@rule
def composite():
    for x in Number, d in Candidate:
        if d.of == x and d.k > 1 and d.divides == True:
            x.composite = True

@rule
def prime():
    for x in Number:
        if x.n > 1 and none(d in Candidate: d.of == x and d.k > 1 and d.divides == True):
            x.prime = True

@rule
def perfect():
    for x in Number:
        if x.divisor_sum == x.n:
            x.perfect = True

@rule
def abundant():
    for x in Number:
        if x.divisor_sum > x.n:
            x.abundant = True

@rule
def deficient():
    for x in Number:
        if x.divisor_sum < x.n:
            x.deficient = True

@rule
def one_is_deficient():
    for x in Number:
        if x.n == 1:
            x.deficient = True

# ── Finite operations as facts ───────────────────────────────────────────────

class Operation:
    # A binary operation on a finite set, by a name you choose, and `order`,
    # how many elements the set has. `operation()` writes it from its table.
    # The rules conclude `identity`, the element that leaves every element
    # as it is from either side; `closed` when every cell's result is in the
    # set; `commutative` when no two cells a·b and b·a differ; `associative`
    # when no three elements give (x·y)·z and x·(y·z) apart; `group` when it
    # is closed and associative and has an identity that every element has a
    # two-sided inverse for; and `abelian` for a commutative group.
    order: int
    identity: string
    closed: bool
    commutative: bool
    associative: bool
    group: bool
    abelian: bool

class Element:
    # One element `el` of the set the operation named `op` acts on, named
    # `<op>_<el>`. The rules conclude `identity`, `invertible` for an element
    # with a two-sided inverse, and `noninvertible` for one without, where
    # the operation has an identity.
    op: string
    el: string
    identity: bool
    invertible: bool
    noninvertible: bool

class Cell:
    # One cell of the table, a·b = c, for the operation named `op`, named
    # `<op>_<a>_<b>`: `left` and `right` name the elements a and b. The
    # rules conclude `outside` for a result not in the set, `noncommuting`
    # where b·a differs, `to_identity` where the result is the identity, and
    # `twosided` where b·a gives the identity too.
    op: string
    a: string
    b: string
    c: string
    left: string
    right: string
    outside: bool
    noncommuting: bool
    to_identity: bool
    twosided: bool

class Triple:
    # Three elements x, y and z in order, for the operation named `op`, named
    # `<op>_<x>_<y>_<z>`: `xy` names the cell x·y and `yz` the cell y·z. The
    # rules conclude `left`, the result of (x·y)·z, `right`, the result of
    # x·(y·z), and `breaks` where the two differ. A triple whose bracket
    # leaves the set has no result to compare, so it is not judged; that is
    # what `closed` reports.
    op: string
    x: string
    z: string
    xy: string
    yz: string
    left: string
    right: string
    breaks: bool

@rule
def identity():
    for e in Element:
        if none(c in Cell: c.left == e and c.c != c.b) and none(c in Cell: c.right == e and c.c != c.a):
            e.identity = True

@rule
def operation_identity():
    for t in Operation, e in Element:
        if e.op == t and e.identity == True:
            t.identity = e.el

@rule
def outside():
    for c in Cell:
        if none(e in Element: e.op == c.op and e.el == c.c):
            c.outside = True

@rule
def closed():
    for t in Operation:
        if none(c in Cell: c.op == t and c.outside == True):
            t.closed = True

@rule
def noncommuting():
    for p in Cell, q in Cell:
        if q.op == p.op and q.a == p.b and q.b == p.a and q.c != p.c:
            p.noncommuting = True

@rule
def commutative():
    for t in Operation:
        if none(c in Cell: c.op == t and c.noncommuting == True):
            t.commutative = True

@rule
def bracket_left():
    # (x·y)·z: the cell whose left operand is the result of x·y.
    for t in Triple, p in Cell, q in Cell:
        if t.xy == p and q.op == p.op and q.a == p.c and q.b == t.z:
            t.left = q.c

@rule
def bracket_right():
    # x·(y·z): the cell whose right operand is the result of y·z.
    for t in Triple, r in Cell, s in Cell:
        if t.yz == r and s.op == r.op and s.a == t.x and s.b == r.c:
            t.right = s.c

@rule
def breaks():
    for t in Triple:
        if t.left != t.right:
            t.breaks = True

@rule
def associative():
    for o in Operation:
        if none(t in Triple: t.op == o and t.breaks == True):
            o.associative = True

@rule
def to_identity():
    for c in Cell, t in Operation:
        if c.op == t and c.c == t.identity:
            c.to_identity = True

@rule
def twosided():
    for p in Cell, q in Cell:
        if q.op == p.op and q.a == p.b and q.b == p.a and p.to_identity == True and q.to_identity == True:
            p.twosided = True

@rule
def invertible():
    for e in Element, c in Cell:
        if c.left == e and c.twosided == True:
            e.invertible = True

@rule
def noninvertible():
    for e in Element, u in Element:
        if u.op == e.op and u.identity == True and none(c in Cell: c.left == e and c.twosided == True):
            e.noninvertible = True

@rule
def group():
    for t in Operation, u in Element:
        if u.op == t and u.identity == True and t.closed == True and t.associative == True and none(e in Element: e.op == t and e.noninvertible == True):
            t.group = True

@rule
def abelian():
    for t in Operation:
        if t.group == True and t.commutative == True:
            t.abelian = True

# ── Triangles as facts ───────────────────────────────────────────────────────

class Triangle:
    # A triangle, by a name you choose, by its three sides in whole units.
    # `triangle()` writes them in order, so `a <= b <= c`. Computed from
    # them: `perimeter`; `gap`, a² + b² − c², whose sign is the cosine of
    # the angle facing the longest side; `slack`, a + b − c, which is
    # positive exactly when the three lengths close up into a triangle, as
    # `valid` says; and `heron`, sixteen times the square of the area. The
    # rules conclude `acute`, `right` or `obtuse` by `gap`, and
    # `equilateral`, `isosceles` (two sides equal or three) or `scalene`,
    # each only for a valid triangle.
    a: int
    b: int
    c: int
    perimeter: int = a + b + c
    gap: int = a * a + b * b - c * c
    slack: int = a + b - c
    heron: int = (a + b + c) * (b + c - a) * (a + c - b) * (a + b - c)
    valid: bool = slack > 0
    acute: bool
    right: bool
    obtuse: bool
    equilateral: bool
    isosceles: bool
    scalene: bool

@rule
def acute():
    for t in Triangle:
        if t.valid == True and t.gap > 0:
            t.acute = True

@rule
def right():
    for t in Triangle:
        if t.valid == True and t.gap == 0:
            t.right = True

@rule
def obtuse():
    for t in Triangle:
        if t.valid == True and t.gap < 0:
            t.obtuse = True

@rule
def equilateral():
    for t in Triangle:
        if t.valid == True and t.a == t.b and t.b == t.c:
            t.equilateral = True

@rule
def isosceles_base():
    for t in Triangle:
        if t.valid == True and t.a == t.b:
            t.isosceles = True

@rule
def isosceles_top():
    for t in Triangle:
        if t.valid == True and t.b == t.c:
            t.isosceles = True

@rule
def scalene():
    for t in Triangle:
        if t.valid == True and t.a < t.b and t.b < t.c:
            t.scalene = True

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

def number(v: int):
    # The number v, named "n" and its value, with a candidate for each k from
    # 1 to v // 2.
    name = "n" + str(v)
    update(Number(name).n, v, "numbers")
    for k in range(1, v // 2 + 1):
        d = name + "_" + str(k)
        update(Candidate(d).of, name, "numbers")
        update(Candidate(d).n, v, "numbers")
        update(Candidate(d).k, k, "numbers")

def numbers(lo: int, hi: int):
    # Every number from lo to hi. The candidates grow as the square of hi,
    # so this is for numbers in the hundreds, not the millions.
    for v in range(lo, hi + 1):
        number(v)

def operation(name: string, elements: list<string>, table: list<list<string>>):
    # An operation from its table: `table[i][j]` is elements[i] · elements[j].
    # It writes n² cells and n³ triples for n elements, so it is for tables
    # of a few elements.
    update(Operation(name).order, len(elements), "table")
    for x in elements:
        se = name + "_" + x
        update(Element(se).op, name, "table")
        update(Element(se).el, x, "table")
    for i in range(len(elements)):
        for j in range(len(elements)):
            s = name + "_" + elements[i] + "_" + elements[j]
            update(Cell(s).op, name, "table")
            update(Cell(s).a, elements[i], "table")
            update(Cell(s).b, elements[j], "table")
            update(Cell(s).c, table[i][j], "table")
            update(Cell(s).left, name + "_" + elements[i], "table")
            update(Cell(s).right, name + "_" + elements[j], "table")
            for k in range(len(elements)):
                st = s + "_" + elements[k]
                syz = name + "_" + elements[j] + "_" + elements[k]
                update(Triple(st).op, name, "table")
                update(Triple(st).x, elements[i], "table")
                update(Triple(st).z, elements[k], "table")
                update(Triple(st).xy, s, "table")
                update(Triple(st).yz, syz, "table")

def cell(name: string, a: string, b: string, c: string):
    # Rewrite one cell of an operation's table: a · b = c.
    s = name + "_" + a + "_" + b
    update(Cell(s).c, c, "table")

def triangle(name: string, x: int, y: int, z: int):
    # A triangle with sides x, y and z, in any order.
    s = sorted([x, y, z])
    update(Triangle(name).a, s[0], "figures")
    update(Triangle(name).b, s[1], "figures")
    update(Triangle(name).c, s[2], "figures")

# ── 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 numbers_that_are(property: string) -> list<int>:
    # The values of every number that is "prime", "composite", "perfect",
    # "abundant", "deficient" or "disputed", sorted.
    value: dict<string, int> = {}
    vs = list_values(Number.n)
    for v in vs:
        parts = v.path.split(".")
        match v.value:
            case int as k:
                value[parts[len(parts) - 2]] = k
            case default:
                pass
    if property == "prime":
        vs = list_values(Number.prime)
    elif property == "composite":
        vs = list_values(Number.composite)
    elif property == "perfect":
        vs = list_values(Number.perfect)
    elif property == "abundant":
        vs = list_values(Number.abundant)
    elif property == "disputed":
        vs = list_values(Number.disputed)
    else:
        vs = list_values(Number.deficient)
    out = []
    for v in vs:
        parts = v.path.split(".")
        x = parts[len(parts) - 2]
        if x in value and value[x] not in out:
            out = out + [value[x]]
    return sorted(out)

def divisors(v: int) -> list<int>:
    # The proper divisors of the number v, sorted: the candidates the rules
    # found to divide it.
    name = "n" + str(v)
    out = []
    vs = list_values(Candidate.divides)
    for d in vs:
        parts = d.path.split(".")
        c = parts[len(parts) - 2]
        if c.startswith(name + "_"):
            k = int(c[len(name) + 1:])
            if k not in out:
                out = out + [k]
    return sorted(out)

def properties(name: string) -> list<string>:
    # The properties that hold of an operation, from "closed",
    # "commutative", "associative", "group" and "abelian", in that order.
    out = []
    for f in ["closed", "commutative", "associative", "group", "abelian"]:
        if f == "closed":
            vs = list_values(Operation.closed)
        elif f == "commutative":
            vs = list_values(Operation.commutative)
        elif f == "associative":
            vs = list_values(Operation.associative)
        elif f == "group":
            vs = list_values(Operation.group)
        else:
            vs = list_values(Operation.abelian)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name and f not in out:
                out = out + [f]
    return out

def identity_of(name: string) -> string:
    # The identity of an operation, or "" where it has none.
    vs = list_values(Operation.identity)
    for v in vs:
        parts = v.path.split(".")
        if parts[len(parts) - 2] == name:
            return str(v.value)
    return ""

def noncommuting_pairs(name: string) -> list<string>:
    # Every pair "a b" whose cells a·b and b·a differ, sorted.
    out = []
    vs = list_values(Cell.noncommuting)
    for v in vs:
        parts = v.path.split(".")
        c = parts[len(parts) - 2]
        if c.startswith(name + "_"):
            ab = c[len(name) + 1:].split("_")
            s = ab[0] + " " + ab[1]
            if s not in out:
                out = out + [s]
    return sorted(out)

def nonassociative_triples(name: string) -> list<string>:
    # Every three elements "x y z" for which (x·y)·z and x·(y·z) differ,
    # sorted.
    out = []
    vs = list_values(Triple.breaks)
    for v in vs:
        parts = v.path.split(".")
        c = parts[len(parts) - 2]
        if c.startswith(name + "_"):
            xyz = c[len(name) + 1:].split("_")
            s = xyz[0] + " " + xyz[1] + " " + xyz[2]
            if s not in out:
                out = out + [s]
    return sorted(out)

def without_inverse(name: string) -> list<string>:
    # Every element of an operation with an identity that has no two-sided
    # inverse, sorted.
    out = []
    vs = list_values(Element.noninvertible)
    for v in vs:
        parts = v.path.split(".")
        c = parts[len(parts) - 2]
        if c.startswith(name + "_"):
            x = c[len(name) + 1:]
            if x not in out:
                out = out + [x]
    return sorted(out)

def kinds(name: string) -> list<string>:
    # What a triangle is, from "acute", "right", "obtuse", "equilateral",
    # "isosceles" and "scalene", in that order; empty where its sides make
    # no triangle.
    out = []
    for f in ["acute", "right", "obtuse", "equilateral", "isosceles", "scalene"]:
        if f == "acute":
            vs = list_values(Triangle.acute)
        elif f == "right":
            vs = list_values(Triangle.right)
        elif f == "obtuse":
            vs = list_values(Triangle.obtuse)
        elif f == "equilateral":
            vs = list_values(Triangle.equilateral)
        elif f == "isosceles":
            vs = list_values(Triangle.isosceles)
        else:
            vs = list_values(Triangle.scalene)
        for v in vs:
            parts = v.path.split(".")
            if parts[len(parts) - 2] == name and f not in out:
                out = out + [f]
    return out

def area(name: string) -> real:
    # A triangle's area by Heron's formula, the square root of `heron` over
    # 4; 0.0 where its sides make no triangle.
    h = current(Triangle.name.heron)
    match h.value:
        case int as k:
            if k > 0:
                return math.sqrt(k) / 4
        case default:
            pass
    return 0.0

def unsettled() -> list<string>:
    # Each number, candidate, cell and triangle side that holds two different
    # values for a fact you write, as "number x n", "cell x c" and so on,
    # sorted.
    out = []
    for f in ["n", "k", "c", "a", "b", "side"]:
        if f == "n":
            vs = list_values(Number.n)
            label = "number "
        elif f == "k":
            vs = list_values(Candidate.k)
            label = "candidate "
        elif f == "c":
            vs = list_values(Cell.c)
            label = "cell "
        elif f == "a":
            vs = list_values(Triangle.a)
            label = "triangle "
        elif f == "b":
            vs = list_values(Triangle.b)
            label = "triangle "
        else:
            vs = list_values(Triangle.c)
            label = "triangle "
        seen = {}
        for v in vs:
            parts = v.path.split(".")
            s = label + parts[len(parts) - 2] + " " + parts[len(parts) - 1]
            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: divisibility ────────────────────────────────────────────────

def gcd(a: int, b: int) -> int:
    # The greatest common divisor by Euclid's algorithm. No pair of 64-bit
    # numbers needs more than 93 steps.
    x = abs(a)
    y = abs(b)
    for step in range(100):
        if y == 0:
            break
        r = x % y
        x = y
        y = r
    return x

def lcm(a: int, b: int) -> int:
    # The least common multiple, 0 when either is 0.
    if a == 0 or b == 0:
        return 0
    return abs(a) // gcd(a, b) * abs(b)

def extended_gcd(a: int, b: int) -> list<int>:
    # [g, x, y] with a·x + b·y = g, the greatest common divisor, by the
    # extended Euclidean algorithm, for a and b at least 0.
    r0 = a
    r1 = b
    s0 = 1
    s1 = 0
    t0 = 0
    t1 = 1
    for step in range(100):
        if r1 == 0:
            break
        q = r0 // r1
        r = r0 - q * r1
        r0 = r1
        r1 = r
        s = s0 - q * s1
        s0 = s1
        s1 = s
        t = t0 - q * t1
        t0 = t1
        t1 = t
    return [r0, s0, t0]

def modinv(a: int, m: int) -> int:
    # The inverse of a modulo m, from 1 to m − 1, or 0 where a and m share a
    # factor and there is none. m is at least 2.
    e = extended_gcd(a % m, m)
    if e[0] != 1:
        return 0
    return e[1] % m

def modpow(b: int, e: int, m: int) -> int:
    # b to the power e modulo m by repeated squaring, for e at least 0 and m
    # from 1 to 3,037,000,499, so that the product of two remainders fits in
    # an int.
    r = 1 % m
    x = b % m
    k = e
    for step in range(64):
        if k == 0:
            break
        if k % 2 == 1:
            r = r * x % m
        x = x * x % m
        k = k // 2
    return r

def isqrt(n: int) -> int:
    # The largest whole number whose square is at most n, n at least 0, by
    # Newton's method in whole numbers. Its first step halves n rather than
    # adding 1 to it, so n = 2⁶³ − 1 does not refuse.
    if n < 2:
        return n
    x = n
    y = n // 2 + n % 2
    for step in range(100):
        if y >= x:
            break
        x = y
        y = (x + n // x) // 2
    return x

def is_prime(n: int) -> bool:
    # Whether n is prime, by the Miller–Rabin test to the bases 2, 3, 5 and
    # 7, which decide every n below 3,215,031,751. n is at most
    # 3,037,000,499, the bound of `modpow`.
    if n < 2:
        return False
    for p in [2, 3, 5, 7]:
        if n % p == 0:
            return n == p
    d = n - 1
    s = 0
    for step in range(64):
        if d % 2 != 0:
            break
        d = d // 2
        s = s + 1
    for a in [2, 3, 5, 7]:
        x = modpow(a, d, n)
        if x != 1 and x != n - 1:
            witness = True
            for r in range(s - 1):
                x = x * x % n
                if x == n - 1:
                    witness = False
                    break
            if witness:
                return False
    return True

def factorize(n: int) -> list<int>:
    # The prime factors of n, at least 2, smallest first and each as often
    # as it divides n, by trial division. Its steps grow as the square root
    # of n's second largest prime factor, so a prime near 3·10⁹ takes the
    # most a run of the CLI allows.
    out = []
    m = n
    for p in [2, 3]:
        for step in range(64):
            if m % p != 0:
                break
            out = out + [p]
            m = m // p
    d = 5
    for step in range(isqrt(m) // 6 + 2):
        if d * d > m:
            break
        if m % d == 0:
            for e in range(64):
                if m % d != 0:
                    break
                out = out + [d]
                m = m // d
        if m % (d + 2) == 0:
            for e in range(64):
                if m % (d + 2) != 0:
                    break
                out = out + [d + 2]
                m = m // (d + 2)
        d = d + 6
    if m > 1:
        out = out + [m]
    return out

def totient(n: int) -> int:
    # Euler's totient: how many of 1 to n share no factor with n.
    r = n
    last = 0
    for p in factorize(n):
        if p != last:
            r = r // p * (p - 1)
            last = p
    return r

def primes(n: int) -> list<int>:
    # Every prime up to n, by the sieve of Eratosthenes. The sieve holds n
    # values, so n is for the thousands.
    if n < 2:
        return []
    sieve = [True] * (n + 1)
    sieve[0] = False
    sieve[1] = False
    for i in range(2, isqrt(n) + 1):
        if sieve[i]:
            for j in range(i * i, n + 1, i):
                sieve[j] = False
    out = []
    for i in range(n + 1):
        if sieve[i]:
            out = out + [i]
    return out

# ── Calculation: counting ────────────────────────────────────────────────────

def factorial(n: int) -> int:
    # n!, for n from 0 to 20.
    r = 1
    for i in range(2, n + 1):
        r = r * i
    return r

def binomial(n: int, k: int) -> int:
    # n choose k, exactly. Each step divides out a common factor first, so
    # no step is larger than the coefficient it reaches, and only a
    # coefficient past 2⁶³ − 1 refuses.
    if k < 0 or k > n:
        return 0
    j = k
    if n - k < j:
        j = n - k
    r = 1
    for i in range(j):
        g = gcd(r, i + 1)
        r = r // g * ((n - i) // ((i + 1) // g))
    return r

def fibonacci(n: int) -> int:
    # The n-th Fibonacci number, F(0) = 0 and F(1) = 1, for n from 0 to 92.
    # It computes nothing past F(n), so F(92) does not refuse for F(93).
    if n == 0:
        return 0
    a = 0
    b = 1
    for i in range(n - 1):
        t = a + b
        a = b
        b = t
    return b

def collatz(n: int) -> int:
    # How many steps n takes to reach 1, halving an even number and taking
    # 3n + 1 for an odd one; −1 if it has not reached 1 in 1,024 steps.
    # The steps run in blocks of 64, since a run is charged for each value a
    # `range` holds, so a short trajectory costs little.
    x = n
    steps = 0
    for block in range(16):
        for k in range(64):
            if x == 1:
                break
            if x % 2 == 0:
                x = x // 2
            else:
                x = 3 * x + 1
            steps = steps + 1
        if x == 1:
            break
    if x != 1:
        return -1
    return steps

# ── Calculation: continued fractions ─────────────────────────────────────────

def continued_fraction(p: int, q: int) -> list<int>:
    # The terms [a0; a1, a2, …] of the fraction p/q, q at least 1, by
    # Euclid's algorithm.
    out = []
    x = p
    y = q
    for step in range(100):
        if y == 0:
            break
        a = x // y
        out = out + [a]
        r = x - a * y
        x = y
        y = r
    return out

def convergents(terms: list<int>) -> list<list<int>>:
    # The convergents [p, q] of a continued fraction, one for each term.
    out = []
    p0 = 1
    q0 = 0
    p1 = terms[0]
    q1 = 1
    out = out + [[p1, q1]]
    for i in range(1, len(terms)):
        p = terms[i] * p1 + p0
        q = terms[i] * q1 + q0
        p0 = p1
        q0 = q1
        p1 = p
        q1 = q
        out = out + [[p1, q1]]
    return out

def sqrt_period(n: int) -> list<int>:
    # The continued fraction of the square root of n: [a0, a1, …, ar], where
    # a1 to ar repeat forever and ar is 2·a0. [a0] alone for a square.
    a0 = isqrt(n)
    out = [a0]
    if a0 * a0 == n:
        return out
    # The pair (m, d) that carries each step has m from 0 to a0 and d from 1
    # to 2·a0, so no period is longer than 2·a0·(a0 + 1): the steps run in
    # blocks of 64 up to that bound, and the answer is never cut short.
    m = 0
    d = 1
    a = a0
    for block in range(2 * a0 * (a0 + 1) // 64 + 1):
        for k in range(64):
            m = d * a - m
            d = (n - m * m) // d
            a = (a0 + m) // d
            out = out + [a]
            if a == 2 * a0:
                break
        if a == 2 * a0:
            break
    return out

def pell(n: int) -> list<int>:
    # The smallest x and y at least 1 with x² − n·y² = 1, for n not a
    # square: [x, y]. With r the length of the repeating part of the square
    # root's continued fraction, it is the convergent before the period's
    # end when r is even, and the one before the end of the second period
    # when r is odd. It refuses only where x itself passes 2⁶³ − 1.
    period = sqrt_period(n)
    r = len(period) - 1
    terms = [period[k] for k in range(r)]
    if r % 2 == 1:
        terms = terms + [period[r]] + [period[k] for k in range(1, r)]
    c = convergents(terms)
    return c[len(c) - 1]

# ── Calculation: numerical methods ───────────────────────────────────────────

def root(x: real, k: int) -> real:
    # The k-th root of x, x above 0 and k at least 1, by Newton's method
    # from above, stopped when a step no longer brings it down.
    g = x
    if g < 1.0:
        g = 1.0
    moving = True
    for block in range(32):
        for step in range(64):
            p = 1.0
            for i in range(k - 1):
                p = p * g
            h = ((k - 1) * g + x / p) / k
            if h >= g:
                moving = False
                break
            g = h
        if not moving:
            break
    return g

def trapezoid(ys: list<real>, h: real) -> real:
    # The integral of a function sampled at equal steps h, by the
    # trapezoidal rule. Where the sum, or the sum times h, passes the largest
    # real while every sample and h are reals, it is summed again with each
    # sample scaled by a power of two, which changes no digit, and the integral
    # scaled back.
    s = (ys[0] + ys[len(ys) - 1]) / 2
    for i in range(1, len(ys) - 1):
        s = s + ys[i]
    v = s * h
    if math.isfinite(v) or not math.isfinite(h) or not all(math.isfinite(y) for y in ys):
        return v
    k = len(ys).bit_length() + 1
    s = (math.ldexp(ys[0], -k) + math.ldexp(ys[len(ys) - 1], -k)) / 2
    for i in range(1, len(ys) - 1):
        s = s + math.ldexp(ys[i], -k)
    return math.ldexp(s * h, k)

def simpson(ys: list<real>, h: real) -> real:
    # The integral of a function sampled at equal steps h, by Simpson's rule,
    # over an odd number of samples. Where the sum, or the sum times h,
    # passes the largest real while every sample and h are reals, it is taken
    # again as trapezoid does.
    s = ys[0] + ys[len(ys) - 1]
    for i in range(1, len(ys) - 1):
        if i % 2 == 1:
            s = s + 4 * ys[i]
        else:
            s = s + 2 * ys[i]
    v = s * h / 3
    if math.isfinite(v) or not math.isfinite(h) or not all(math.isfinite(y) for y in ys):
        return v
    k = len(ys).bit_length() + 3
    s = math.ldexp(ys[0], -k) + math.ldexp(ys[len(ys) - 1], -k)
    for i in range(1, len(ys) - 1):
        if i % 2 == 1:
            s = s + 4 * math.ldexp(ys[i], -k)
        else:
            s = s + 2 * math.ldexp(ys[i], -k)
    return math.ldexp(s * h / 3, k)

def machin() -> real:
    # π by Machin's formula of 1706, 4·(4·atan(1/5) − atan(1/239)).
    return 4 * (4 * math.atan(1 / 5) - math.atan(1 / 239))

def determinant(m: list<list<real>>) -> real:
    # The determinant of a square matrix by Gaussian elimination with
    # partial pivoting: at each column, the row with the largest entry. The
    # product of the pivots is kept as a mantissa and a power of two, as
    # math.frexp gives them, which changes no digit of a product that stays in
    # range. Where an entry passes the largest real on the way, or a
    # multiplier falls below the smallest normal real and loses the row it
    # scales, while every entry of m is a real, it is taken again with each
    # row scaled by a power of two that brings its largest entry near 1, and
    # the powers added back.
    n = len(m)
    rows = m
    shift = 0
    for attempt in range(2):
        a = [[rows[i][j] for j in range(n)] for i in range(n)]
        dm = 1.0
        de = shift
        lost = False
        for k in range(n):
            p = k
            for i in range(k + 1, n):
                if abs(a[i][k]) > abs(a[p][k]):
                    p = i
            if a[p][k] == 0:
                return 0.0
            if p != k:
                t = a[k]
                a[k] = a[p]
                a[p] = t
                dm = -dm
            fp = math.frexp(a[k][k])
            fd = math.frexp(dm * fp[0])
            dm = fd[0]
            de = de + fp[1] + fd[1]
            for i in range(k + 1, n):
                f = a[i][k] / a[k][k]
                if a[i][k] != 0.0 and abs(f) < 2.2250738585072014e-308:
                    lost = True
                for j in range(k, n):
                    a[i][j] = a[i][j] - f * a[k][j]
        if attempt == 1 or (math.isfinite(dm) and not lost) or not all(math.isfinite(v) for row in m for v in row):
            return math.ldexp(dm, de)
        rows = []
        for i in range(n):
            top = 0.0
            for v in m[i]:
                if abs(v) > top:
                    top = abs(v)
            e = math.frexp(top)[1]
            rows = rows + [[math.ldexp(v, -e) for v in m[i]]]
            shift = shift + e
    return 0.0

def solve(m: list<list<real>>, b: list<real>) -> list<real>:
    # The x with m·x = b, by Gaussian elimination with partial pivoting and
    # back substitution; empty where m is singular. 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 m and b are reals, the system
    # is solved again with each row of m 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(m)
    rows = m
    rhs = b
    c = 0
    for attempt in range(2):
        a = [[rows[i][j] for j in range(n)] + [rhs[i]] for i in range(n)]
        lost = False
        for k in range(n):
            p = k
            for i in range(k + 1, n):
                if abs(a[i][k]) > abs(a[p][k]):
                    p = i
            if a[p][k] == 0:
                return []
            if abs(a[p][k]) < 2.2250738585072014e-308:
                lost = True
            if p != k:
                t = a[k]
                a[k] = a[p]
                a[p] = t
            for i in range(k + 1, n):
                f = a[i][k] / a[k][k]
                if a[i][k] != 0.0 and abs(f) < 2.2250738585072014e-308:
                    lost = True
                for j in range(k, n + 1):
                    a[i][j] = a[i][j] - f * a[k][j]
        x = [0.0] * n
        for step in range(n):
            i = n - 1 - step
            s = a[i][n]
            for j in range(i + 1, n):
                s = s - a[i][j] * x[j]
            x[i] = s / a[i][i]
        if attempt == 1:
            return [math.ldexp(v, c) 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 m for v in row) or not all(math.isfinite(v) for v in b):
            return x
        rows = []
        rhs = []
        for i in range(n):
            top = 0.0
            for v in m[i]:
                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 m[i]]]
            rhs = rhs + [math.ldexp(b[i], -e)]
        big = 0.0
        for v in rhs:
            if abs(v) > big:
                big = abs(v)
        c = n.bit_length() + n + 1
        if big > 1.0:
            c = c + math.frexp(big)[1]
        rhs = [math.ldexp(v, -c) for v in rhs]
    return []

def hilbert(n: int) -> list<list<real>>:
    # The n-by-n Hilbert matrix, 1/(i + j + 1), whose determinant falls
    # toward zero faster than elimination can follow.
    return [[1 / (i + j + 1) for j in range(n)] for i in range(n)]