#!/usr/bin/env python3
"""
PROVENANCE: SURROGATE

Computation 186 -- The Cauchy slice is a slab with a matched thickness, and the missing
power of mass is that thickness
=========================================================================================
WHAT THIS DOES.  comp 185 reduced the W4 register item to a single residue, the
slice-wise CAR limit, and did not supply it.  This file supplies the construction and the
scaling condition it needs, and measures both, on the i.i.d. surrogate this book already
uses for A2 and A3.

THE OBSTRUCTION, STATED EXACTLY.  sec:spin-statistics-structural asks that the substrate
relations {c_a, c_b^dagger} = delta_ab converge "on a Cauchy slice" to the equal-time
continuum CAR with delta^(3)(x-y), and says outright: "The slice matters and the
construction here does not supply it."  The reason it does not is sharper than a gap in
rigour.  A Cauchy slice Sigma is a THREE-dimensional subset of the four-dimensional M,
hence of measure zero, while rho_n(D_n) is a FINITE point set.  So rho_n^{-1}(Sigma) is
EMPTY with probability one, measured in part 1 and true at every size.  "Pose the
relations on a Cauchy slice" therefore cannot mean "restrict to the properties realised
on Sigma", because there are none.

THE BOOK'S OWN DIAGNOSIS POINTS AT THE FIX.  The same passage records what went wrong in
an earlier version: it "wrote delta^(4)(x-y), summing over sites realised into the
four-dimensional M rather than into a slice", which is "dimensionally wrong by one power
of mass".  THAT MISSING POWER OF MASS IS A THICKNESS.  Take the slab

    S_n(t) = { a in D_n : |tau(rho_n(a)) - t| <= h_n/2 },

tau the sublimation-derived time function.  Its site density PER UNIT 3-VOLUME on Sigma
is nu_n = lambda_n h_n = h_n / a_n^4, carrying [mass]^3 where the 4-density lambda_n
carries [mass]^4.  Smearing with nu_n^{-1/2} then gives

    {psi_n(f), psi_n(g)^dagger} = nu_n^{-1} sum_{a in S_n} <f,g>(rho_n(a))
                               -> integral_Sigma <f,g> d^3x,

which is the distributional delta^(3) the book asks for.  The slab thickness supplies
exactly the one power of mass that was missing, and nothing else has to change.

TWO CONDITIONS, AND EACH KILLS ONE ERROR TERM.  Part 3 separates them and measures both:

    h_n -> 0           kills the BIAS, which is O(h^2): a slab that never shrinks
                       returns the slab average and not the value on the slice
    lambda_n h_n -> oo kills the FLUCTUATION, which is O((lambda_n h_n)^{-1/2}): this is
                       the CLT rate, and it is the Berry-Esseen rate the book says the
                       convergence should match

Both hold on a WIDE window, a_n^4 << h_n << 1, so the construction is not a knife edge.

AND THE TWO TERMS PICK OUT A PREFERRED THICKNESS.  Balancing h^2 against
(lambda h)^{-1/2} gives

    h_n ~ lambda_n^{-1/5} = a_n^{4/5},      total error ~ lambda_n^{-2/5} = a_n^{8/5}.

Measured in part 4: the RMS error over slab thicknesses has its minimum at h = 0.12
against the predicted scaling 0.076 at that size, agreeing to a constant of order one,
and the convergence exponent in lambda comes out -0.36 against the predicted -0.40.
Note that a_n^{4/5} >> a_n: THE OPTIMAL SLAB IS THICKER THAN THE SITE SPACING, so it
contains a diverging number of sites in the time direction too, a_n^{-1/5} of them.

WHAT THIS DOES NOT DO, and it is the same boundary every other matched-scaling result in
this book sits behind.  This is the I.I.D. SURROGATE, the same standing as A2's
verification and A3's corroboration (Computations 70, 83).  For the
emergence-determined rho_n the equidistribution that part 3 leans on is exactly A3, which
P3's maximum-entropy sublimation supplies rather than this file.  So the item is NOT
closed here unilaterally: it is reduced to the typicality condition the Mosco chain
already carries, and it needs the sublimation-derived time function tau that P3 supplies.
Whether the register may then call it closed is a question about A3's reach, not about
this construction.
"""

import sys

import numpy as np


def rule():
    print("=" * 100)


# --- the surrogate: an equidistributed cloud in a unit 4-box ---------------------
# f and g carry EVEN time dependence about the slice, so the slab bias does not
# cancel by symmetry and can be measured rather than hidden.
def fg(t, x):
    d = t - 0.5
    f = 1.0 + 4.0 * d * d + np.sin(2 * np.pi * x[:, 0])
    g = 2.0 + 1.0 * d * d + np.cos(2 * np.pi * x[:, 1])
    return f * g


TARGET = 2.0   # on the slice d = 0: int (1 + sin)(2 + cos) d^3x over the unit cube


def slab_estimate(N, h, seed, norm="slab", chunk=200_000):
    """nu^{-1} sum over the slab, with nu = lambda h (correct) or lambda (the old
    delta^4 reading).

    DRAWN IN CHUNKS, and that is not a detail.  This book publishes its computations
    to be re-run, and the website runs them in the browser under Pyodide, whose WASM
    heap is far smaller than a workstation's.  Allocating the whole cloud at once
    needs 4N floats: 128 MB at the largest size used here, which is the kind of single
    allocation that fails in a browser while passing silently on a laptop.  A chunked
    draw holds 6 MB instead.  The numbers are UNCHANGED: Generator.random fills in C
    order from one bit stream, so chunk-wise draws concatenate to exactly the array a
    single draw would have produced, which part 0 checks rather than assumes."""
    r = np.random.default_rng(seed)
    nu = N * h if norm == "slab" else N
    total, hits, done = 0.0, 0, 0
    while done < N:
        k = min(chunk, N - done)
        done += k
        p = r.random((k, 4))
        t, x = p[:, 0], p[:, 1:]
        m = np.abs(t - 0.5) <= h / 2
        if m.any():
            hits += int(m.sum())
            total += float(np.sum(fg(t[m], x[m])))
    if hits == 0:
        return np.nan
    return total / nu


def part0():
    rule()
    print("  0.  THE CHUNKED DRAW IS NUMERICALLY INERT")
    rule()
    print("      The cloud is drawn in chunks so the peak allocation stays small enough")
    print("      for the browser runtime this book publishes to.  That is a change to the")
    print("      code, so it is checked rather than asserted, and in two steps because")
    print("      the first attempt claimed more than is true.")
    print()
    print("      (a) THE POINTS ARE THE SAME, BIT FOR BIT.  Generator.random fills in C")
    print("      order from one bit stream, so chunk-wise draws concatenate to exactly")
    print("      the array a single draw gives.  Compared directly:")
    print()
    print(f"      {'N':>9}{'chunks':>9}{'identical cloud?':>19}")
    ok = True
    for N, ch in ((200_000, 50_000), (200_000, 7_919)):
        a = np.random.default_rng(0).random((N, 4))
        r = np.random.default_rng(0)
        parts, done = [], 0
        while done < N:
            k = min(ch, N - done)
            done += k
            parts.append(r.random((k, 4)))
        same = bool(np.array_equal(a, np.concatenate(parts)))
        ok &= same
        print(f"      {N:>9}{ch:>9}{str(same):>19}")
    print()
    print("      (b) THE SUM AGREES TO ONE ULP, NOT BIT FOR BIT, and the difference is")
    print("      SUMMATION ORDER rather than sampling.  Floating-point addition is not")
    print("      associative, and accumulating per chunk adds the same numbers in a")
    print("      different order; where the chunk divides N the partial sums line up and")
    print("      the result is identical, where it does not they differ in the last bit.")
    print()
    print(f"      {'N':>9}{'chunk':>9}{'chunked':>15}{'single draw':>15}"
          f"{'relative diff':>15}")
    for N, ch in ((200_000, 50_000), (200_000, 7_919), (400_000, 200_000)):
        h = N ** -0.2
        got = slab_estimate(N, h, 0, chunk=ch)
        ref = slab_estimate(N, h, 0, chunk=N)
        rel = abs(got - ref) / abs(ref)
        ok &= rel < 1e-14
        print(f"      {N:>9}{ch:>9}{got:>15.10f}{ref:>15.10f}{rel:>15.1e}")
    print()
    print("      One part in 1e16 on a quantity whose own statistical error is one part")
    print("      in 1e2, so nothing in this file moves.  Peak memory falls from about")
    print("      210 MB to about 50 MB, which is the point of the change.")
    print()
    print(f"      -> same points, same answer to one ULP, a quarter the memory : {ok}")
    print()
    return ok


def part1():
    rule()
    print("  1.  WHY THE CONSTRUCTION CANNOT SUPPLY THE SLICE AS IT STANDS")
    rule()
    print("      sec:spin-statistics-structural asks the substrate relations to converge")
    print("      'on a Cauchy slice' to the equal-time CAR with delta^(3), and says")
    print("      outright that the construction does not supply the slice.  The reason is")
    print("      sharper than a gap in rigour and is worth stating as a count.")
    print()
    print("      A Cauchy slice Sigma is THREE-dimensional inside the four-dimensional M,")
    print("      hence of measure zero.  rho_n(D_n) is a FINITE point set.  So")
    print("      rho_n^{-1}(Sigma) is EMPTY with probability one: there are no properties")
    print("      realised ON the slice to pose anything on.")
    print()
    print(f"      {'|D_n|':>10}{'hits on the exact slice t = 1/2':>34}")
    ok = True
    for N in (10 ** 4, 10 ** 5, 10 ** 6):
        r = np.random.default_rng(N)
        t = r.random(N)
        hits = int(np.sum(t == 0.5))
        ok &= (hits == 0)
        print(f"      {N:>10}{hits:>34}")
    print()
    print("      So 'pose the relations slice-wise' cannot mean 'restrict to the sites")
    print("      realised on Sigma'.  It has to mean something else, and the book's own")
    print("      diagnosis of the earlier error says what.")
    print()
    print(f"      -> the slice carries no sites, so a slab is forced     : {ok}")
    print()
    return ok


def part2():
    rule()
    print("  2.  THE MISSING POWER OF MASS IS A THICKNESS")
    rule()
    print("      The book records the earlier error precisely: it 'wrote delta^(4)(x-y),")
    print("      summing over sites realised into the four-dimensional M rather than into")
    print("      a slice', which is 'dimensionally wrong by one power of mass'.")
    print()
    print("      Take the slab S_n(t) = { a : |tau(rho_n(a)) - t| <= h_n/2 }.  Counting:")
    print()
    print("          4-density of sites in M          lambda_n = 1/a_n^4   [mass]^4")
    print("          sites in the slab per 3-volume   nu_n = lambda_n h_n  [mass]^3")
    print()
    print("      and [mass]^3 is exactly what {psi, psi^dagger} carries in four")
    print("      dimensions.  Smearing with nu_n^{-1/2} per mode,")
    print()
    print("          {psi_n(f), psi_n(g)^+} = nu_n^{-1} sum_{a in S_n} <f,g>(rho_n(a))")
    print("                                -> integral_Sigma <f,g> d^3x,")
    print()
    print("      which is the distributional delta^(3).  THE SLAB THICKNESS SUPPLIES THE")
    print("      MISSING POWER AND NOTHING ELSE HAS TO CHANGE.  Measured, against the")
    print("      exact slice integral, with the old reading beside it:")
    print()
    N, h = 10 ** 6, (10 ** 6) ** -0.25
    slab = float(np.mean([slab_estimate(N, h, s) for s in range(8)]))
    four = float(np.mean([slab_estimate(N, h, s, "four") for s in range(8)]))
    print(f"      {'normalisation':>34}{'estimate':>12}{'target':>10}")
    print(f"      {'nu = lambda h   (slab, [mass]^3)':>34}{slab:>12.5f}{TARGET:>10.3f}")
    print(f"      {'nu = lambda     (the old delta^4)':>34}{four:>12.5f}{TARGET:>10.3f}")
    print()
    print(f"      The old reading is low by a factor {slab / four:.1f}, and h = {h:.4f},")
    print(f"      so 1/h = {1/h:.1f}: the discrepancy IS the thickness, made visible.")
    ok = abs(slab - TARGET) < 0.05 and abs(four / slab - h) < 0.05 * h + 1e-3
    print()
    print(f"      -> the slab normalisation lands on the slice integral  : {ok}")
    print()
    return ok


def part3():
    rule()
    print("  3.  TWO CONDITIONS, AND EACH KILLS ONE ERROR TERM")
    rule()
    print("      The error splits, and the split is the whole content of the scaling.")
    print()
    print("          h_n -> 0            kills the BIAS, O(h^2): a slab that never")
    print("                              shrinks returns the slab average, not the slice")
    print("          lambda_n h_n -> oo  kills the FLUCTUATION, O((lambda_n h_n)^{-1/2}),")
    print("                              which is the CLT rate")
    print()
    print("      Measured separately at |D_n| = 4e5 over 60 clouds, against the predicted")
    print("      coefficients rather than against a fit:")
    print()
    N, S = 400000, 60
    print(f"      {'h':>9}{'|S_n|':>9}{'bias':>11}{'predicted':>11}{'ratio':>8}"
          f"{'s.e.m.':>10}{'resolved?':>11}{'std':>10}{'predicted':>11}")
    bias_ok = fluc_ok = True
    for h in (0.40, 0.20, 0.10, 0.05, 0.025):
        v = np.array([slab_estimate(N, h, s) for s in range(S)])
        bias, std = float(np.mean(v) - TARGET), float(np.std(v))
        sem = std / np.sqrt(S)
        pb, pf = h * h / 3 * 2, (N * h) ** -0.5 * 2.3
        resolved = bias > 2 * sem
        # the bias claim is only testable where the bias clears the noise on the mean
        if resolved:
            bias_ok &= (0.5 * pb < bias < 2.0 * pb)
        fluc_ok &= (0.5 * pf < std < 2.0 * pf)
        print(f"      {h:>9.3f}{int(N*h):>9}{bias:>11.5f}{pb:>11.5f}{bias/pb:>8.2f}"
              f"{sem:>10.5f}{('yes' if resolved else 'no'):>11}{std:>10.5f}{pf:>11.5f}")
    print()
    print("      A RESOLUTION FLOOR, stated rather than tuned around.  The bias is the")
    print("      shift of a MEAN over 60 clouds, so it is only measurable while it")
    print("      exceeds the standard error of that mean.  It does for the first three")
    print("      rows, where it tracks the predicted h^2 coefficient to a constant of")
    print("      about 1.2; below that the rows carry no information about the bias and")
    print("      are marked so.  Thinning the slab shrinks the bias and inflates the very")
    print("      noise that measures it, which is the trade-off this part is about.")
    print()
    print("      The bias falls as h^2 and tracks its predicted coefficient; the spread")
    print("      grows as the slab thins, exactly as 1/sqrt of its population.  So the")
    print("      two conditions are not one condition stated twice: thinning the slab")
    print("      helps one term and hurts the other.")
    print()
    print("      THE ADMISSIBLE WINDOW IS WIDE.  Both hold whenever a_n^4 << h_n << 1,")
    print("      and with lambda_n = a_n^{-4} that is satisfied by h_n = a_n itself, or")
    print("      by any power a_n^p with 0 < p < 4.  This is not a knife edge.")
    print()
    ok = bias_ok and fluc_ok
    print(f"      -> bias is O(h^2) and fluctuation is O(nu^-1/2)        : {ok}")
    print()
    return ok


def part4():
    rule()
    print("  4.  THE TWO TERMS PICK OUT A PREFERRED THICKNESS")
    rule()
    print("      Balancing h^2 against (lambda h)^{-1/2}:")
    print()
    print("          h_n ~ lambda_n^{-1/5} = a_n^{4/5},   error ~ lambda_n^{-2/5} = a_n^{8/5}")
    print()
    print("      Note a_n^{4/5} >> a_n: THE OPTIMAL SLAB IS THICKER THAN THE SITE")
    print("      SPACING, holding a_n^{-1/5} sites across its thickness, so it is a slab")
    print("      in the time direction too and not a thickened slice of one site.")
    print()
    N, S = 400000, 60
    print(f"      RMS error over slab thickness at |D_n| = {N}:")
    print()
    print(f"      {'h':>9}{'RMS':>11}")
    best = None
    for h in (0.30, 0.20, 0.12, 0.08, 0.05, 0.03, 0.02, 0.012, 0.008):
        v = np.array([slab_estimate(N, h, s) for s in range(S)])
        rms = float(np.sqrt(np.mean((v - TARGET) ** 2)))
        if best is None or rms < best[1]:
            best = (h, rms)
        print(f"      {h:>9.3f}{rms:>11.5f}")
    pred = N ** -0.2
    print()
    print(f"      minimum at h = {best[0]:.3f}, predicted scaling lambda^(-1/5) = "
          f"{pred:.4f}")
    print("      The prediction is a SCALING and not a constant, so agreeing to a factor")
    print("      of order one is what it asserts; the constant is not derived here.")
    print()
    print("      Convergence with h = lambda^(-1/5):")
    print()
    print(f"      {'|D_n|':>10}{'h':>10}{'|S_n|':>10}{'RMS error':>13}")
    Ns, es = (10 ** 4, 10 ** 5, 10 ** 6, 4 * 10 ** 6), []
    for n in Ns:
        h = n ** -0.2
        v = np.array([slab_estimate(n, h, s) for s in range(30)])
        e = float(np.sqrt(np.mean((v - TARGET) ** 2)))
        es.append(e)
        print(f"      {n:>10}{h:>10.4f}{int(n*h):>10}{e:>13.5f}")
    slope = float(np.polyfit(np.log(Ns), np.log(es), 1)[0])
    ok = (0.5 * pred < best[0] < 3.0 * pred) and (-0.55 < slope < -0.28)
    print()
    print(f"      fitted exponent in lambda: {slope:+.2f}, predicted -0.40.  The shortfall")
    print("      is the same order-one constant the minimum shows and is reported rather")
    print("      than tuned away.")
    print()
    print(f"      -> the preferred thickness and rate are as predicted    : {ok}")
    print()
    return ok


def part5():
    rule()
    print("  5.  TAMPER TESTS: DOES THE MEASUREMENT DETECT A WRONG SLICING?")
    rule()
    print("      Three ways to get it wrong, each changing one thing:")
    print()
    N = 10 ** 6
    h = N ** -0.25
    good = float(np.mean([slab_estimate(N, h, s) for s in range(8)]))
    four = float(np.mean([slab_estimate(N, h, s, "four") for s in range(8)]))
    fixed = float(np.mean([slab_estimate(N, 0.40, s) for s in range(8)]))
    thin = np.array([slab_estimate(N, 4.0 / N, s) for s in range(40)])
    print(f"      {'what is changed':>40}{'estimate':>12}{'verdict':>22}")
    print(f"      {'nothing (h = a_n, correct)':>40}{good:>12.5f}{'lands on 2.000':>22}")
    print(f"      {'nu = lambda, the old delta^4 reading':>40}{four:>12.5f}"
          f"{'low by 1/h':>22}")
    print(f"      {'h fixed at 0.40, slab never shrinks':>40}{fixed:>12.5f}"
          f"{'slab average, biased':>22}")
    print(f"      {'h = 4/lambda, so nu = O(1)':>40}{np.nanmean(thin):>12.5f}"
          f"{'spread ' + format(float(np.nanstd(thin)), '.2f'):>22}")
    print()
    print("      The third is the informative one.  A slab that does not shrink still")
    print("      CONVERGES, and converges to a perfectly well-defined number; it is just")
    print("      the wrong number, the average over the slab rather than the value on the")
    print("      slice.  So the test distinguishes a slice from a slab, which is the")
    print("      distinction the whole item is about, rather than merely detecting noise.")
    print()
    ok = (abs(good - TARGET) < 0.05 and abs(four - TARGET) > 1.0
          and abs(fixed - TARGET) > 0.05 and float(np.nanstd(thin)) > 0.5)
    print(f"      -> each wrong slicing is detected, and differently     : {ok}")
    print()
    return ok


def part6():
    rule()
    print("  6.  WHAT THIS SUPPLIES, AND THE BOUNDARY IT DOES NOT CROSS")
    rule()
    print("      SUPPLIED.  A construction for the slice-wise CAR limit, the scaling")
    print("      condition it needs (a_n^4 << h_n << 1), the preferred thickness inside")
    print("      that window (h_n ~ a_n^{4/5}), and the rate (lambda_n^{-2/5}), with the")
    print("      fluctuation term being the CLT rate the book says the convergence should")
    print("      match.  The condition is of the same KIND as A1: a matched scaling")
    print("      between two vanishing quantities, stated and then used.")
    print()
    print("      NOT SUPPLIED, and this is the same boundary every matched-scaling result")
    print("      in this book sits behind.  Everything above is the I.I.D. SURROGATE, the")
    print("      standing of A2's verification and of A3's corroboration (comps 70, 83).")
    print("      Part 3's limit leans on the slab population being equidistributed over")
    print("      Sigma, and for the emergence-determined rho_n that is EXACTLY A3, which")
    print("      P3's maximum-entropy sublimation supplies rather than this file.  The")
    print("      construction also needs the sublimation-derived time function tau to")
    print("      exist, which P3 supplies and which is what makes the slicing non-trivial")
    print("      for a pre-causal substrate in the first place.")
    print()
    print("      SO THE ITEM IS NOT CLOSED HERE UNILATERALLY.  What changes is its")
    print("      character: it was 'the construction does not supply the slice', an")
    print("      absence; it is now 'the slab construction supplies it on the surrogate,")
    print("      and the emergence-determined case rests on A3', which is the condition")
    print("      the Mosco chain already carries.  Whether the register may call it closed")
    print("      is a question about A3's reach and not about this construction, and it is")
    print("      left to the register rather than answered here.")
    print()
    print("      -> supplied on the surrogate, reduced to A3 beyond it   : True")
    print()
    return True


def main():
    rule()
    print("  COMPUTATION 186 -- THE CAUCHY SLICE IS A SLAB WITH A MATCHED THICKNESS")
    rule()
    print()
    r = [part0(), part1(), part2(), part3(), part4(), part5(), part6()]
    rule()
    print("  SUMMARY")
    rule()
    print(f"  0. same points bit for bit, same answer to one ULP, 1/4 memory : {r[0]}")
    print(f"  1. the slice carries no sites, so a slab is forced            : {r[1]}")
    print(f"  2. the missing power of mass is the slab thickness            : {r[2]}")
    print(f"  3. bias is O(h^2), fluctuation is O(nu^-1/2), window is wide  : {r[3]}")
    print(f"  4. preferred thickness a_n^(4/5), error a_n^(8/5), measured   : {r[4]}")
    print(f"  5. tamper: a non-shrinking slab converges to the WRONG number : {r[5]}")
    print(f"  6. supplied on the surrogate; beyond it this is A3            : {r[6]}")
    print()
    print("""  RESULT.  comp 185 reduced the W4 register item to one residue, the slice-wise
  CAR limit.  This supplies the construction it was missing, and the construction is
  forced by a count the book had already made.

  A Cauchy slice is three-dimensional inside a four-dimensional M and the realised
  cloud is finite, so no property is realised ON the slice: the preimage is empty at
  every size.  'Pose the relations slice-wise' therefore cannot mean restricting to
  sites on Sigma.  The book's own diagnosis of the earlier delta^(4) error says what
  it does mean: that error was 'dimensionally wrong by one power of mass', and THE
  MISSING POWER IS A THICKNESS.  A slab of 4-thickness h_n has site density
  nu_n = lambda_n h_n per unit 3-volume, which carries [mass]^3, and smearing with
  nu_n^{-1/2} lands on the slice integral, which is the distributional delta^(3).

  The error then splits in two and each condition kills one term.  h_n -> 0 kills a
  bias measured here as O(h^2) against its predicted coefficient; lambda_n h_n -> oo
  kills a fluctuation measured as O(nu^{-1/2}), which is the CLT rate and therefore
  the Berry-Esseen rate the book says the convergence should match.  Thinning the
  slab helps one and hurts the other, so these are two conditions and not one stated
  twice, and both hold on the wide window a_n^4 << h_n << 1.

  Balancing them picks out h_n ~ a_n^{4/5} with total error a_n^{8/5}: the optimal
  slab is THICKER than the site spacing and holds a_n^{-1/5} sites across its
  thickness.  The measured RMS minimum sits within a factor of order one of the
  predicted scaling and the convergence exponent at -0.36 against -0.40, both
  reported rather than tuned.  The tamper test that matters is the slab that never
  shrinks: it converges perfectly well, to the slab average instead of the slice
  value, so the measurement distinguishes a slice from a slab rather than detecting
  noise.

  This is the i.i.d. surrogate, the standing of A2's verification and A3's
  corroboration.  For the emergence-determined realisation the equidistribution this
  leans on IS A3, supplied by P3 rather than here, and the construction needs P3's
  time function besides.  So the item is not closed here: it has gone from 'the
  construction does not supply the slice' to 'the slab construction supplies it on
  the surrogate, and beyond the surrogate this is A3', which is the condition the
  Mosco chain already carries.""")
    print()
    return all(r)


if __name__ == "__main__":
    # Exit only on FAILURE.  The website runs these in the browser under Pyodide,
    # where SystemExit is fatal even at status 0: it tears the runtime down, and the
    # page prints ERROR after output that was in fact correct.  Falling off the end
    # hands the build gates the same status 0, and a failure still exits 1.
    if not main():
        sys.exit(1)
