examples/spectral/

# The largest singular value of a matrix nobody stores, and the number a
# published file says it is.
#
# Corpus material (CLAUDE.md §9), M-corpus-depth's oracle family — the
# second of the two programs that reach C for `sqrt`, and the one whose
# whole shape is arithmetic. What it adds to the corpus beside its oracle:
# **a matrix that is a FUNCTION**. `eval_a` computes an entry from its two
# indices and nothing is stored, so a hundred by a hundred costs no memory
# at all — which is the trick the task exists to teach, and a shape no
# other program here has.
#
# THE ORACLE, fetched and quoted with its source. The benchmarks
# game publishes the output of spectral-norm at N = 100:
#
#     1.274219991
#
#     https://benchmarksgame-team.pages.debian.net/benchmarksgame/download/spectralnorm-output.txt
#     https://benchmarksgame-team.pages.debian.net/benchmarksgame/description/spectralnorm.html
#
# THE ALGORITHM IS COPIED from the reference C program on the same site,
# fetched the same day:
#
#     https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/spectralnorm-gcc-1.html
#
#     double eval_A(int i, int j) { return 1.0/((i+j)*(i+j+1)/2+i+1); }
#
# ten rounds of `AtA·u` alternating two vectors, then `sqrt(vBv / vv)`,
# printed with `"%0.9f\n"`. Every one of those choices is part of the
# answer: nine places is the comparison, and ten rounds is what makes the
# power iteration converge to the published digits rather than near them.
#
# **THE INDEX ARITHMETIC IS DONE IN `i64` AND CONVERTED ONCE**, which is
# not a detail: `(i+j)*(i+j+1)/2` is an integer division in C and stays one
# here, because `/` on two `i64` truncates (spec § Failure) while
# `1.0 / ...` on floats does not. Doing it in `f64` would give a different
# matrix and a plausible wrong answer — the kind of mistake a published
# number catches and a self-written assertion does not.

extern "math.h" link "m"
    function sqrt(x: f64) -> f64

# The published size. The site times this at 5,500; the published OUTPUT is
# for a hundred, and a corpus program takes the size its oracle is for
# (CLAUDE.md §9) — the answer is the same answer, only the clock differs.
constant N: i64
    100

# The published number of rounds. Ten, from the reference program.
constant ROUNDS: i64
    10

# One entry of the matrix, computed rather than stored:
#     1.0 / ((i + j) * (i + j + 1) / 2 + i + 1)
# The inner arithmetic is whole-number arithmetic, and the reciprocal is
# the only float operation in it.
function eval_a(i: i64, j: i64) -> f64
    span = i + j
    return 1.0 / to_f64(span * (span + 1) / 2 + i + 1)

# `A · u`, a row at a time.
function a_times(u: [f64]) -> [f64]
    out: [f64] @= []
    i: i64 @= 0

    while i < u.len()
        total: f64 @= 0.0
        j: i64 @= 0

        while j < u.len()
            total @ total + eval_a(i: i, j: j) * u[j]
            j @ j + 1

        out @ out.push(total)
        i @ i + 1

    return out

# `Aᵀ · u` — the same walk with the indices swapped, which is the whole
# difference between a matrix and its transpose when nothing is stored.
function at_times(u: [f64]) -> [f64]
    out: [f64] @= []
    i: i64 @= 0

    while i < u.len()
        total: f64 @= 0.0
        j: i64 @= 0

        while j < u.len()
            total @ total + eval_a(i: j, j: i) * u[j]
            j @ j + 1

        out @ out.push(total)
        i @ i + 1

    return out

# `AᵀA · u`, through the intermediate the reference calls `v`.
function ata_times(u: [f64]) -> [f64]
    return at_times(a_times(u))

# A vector of ones, which is where the power iteration starts.
function ones(n: i64) -> [f64]
    out: [f64] @= []
    at: i64 @= 0

    while at < n
        out @ out.push(1.0)
        at @ at + 1

    return out

# The spectral norm: ten rounds, then the ratio's square root. The two
# vectors are returned rather than mutated in place, which is what value
# semantics make natural (spec § Types) — and the arithmetic is identical
# to the reference's, which mutates two arrays through pointers.
function spectral_norm(n: i64) -> f64
    u: [f64] @= ones(n)
    v: [f64] @= ones(n)
    round: i64 @= 0

    while round < ROUNDS
        v @ ata_times(u)
        u @ ata_times(v)
        round @ round + 1

    vbv: f64 @= 0.0
    vv: f64 @= 0.0
    at: i64 @= 0

    while at < n
        vbv @ vbv + u[at] * v[at]
        vv @ vv + v[at] * v[at]
        at @ at + 1

    return sqrt(x: vbv / vv)

function main()
    print(nine_places(spectral_norm(N)))

# A number to exactly nine decimal places — C's `%0.9f`, which this
# language does not have: `to_str` prints the shortest text that reads back
# as the same value (spec § Strings, arrays, maps), so the answer would
# print as `1.2742199911868225`. The published file is nine places.
constant BILLION: f64
    1000000000.0

function nine_places(v: f64) -> str
    negative = v < 0.0
    size: f64 @= v

    if negative
        size @ -v
    scaled = to_i64(size * BILLION + 0.5).must()
    whole = scaled / 1000000000
    rest = scaled % 1000000000
    text: str = whole.to_str() + "." + padded(rest)

    if negative
        return "-" + text
    return text

# Nine digits, leading zeroes kept.
function padded(rest: i64) -> str
    digits: str @= rest.to_str()

    while digits.len() < 9
        digits @ "0" + digits

    return digits

test "the published number, which this program did not invent"
    # benchmarks game, spectralnorm-output.txt, N = 100, fetched:
    #     1.274219991
    assert nine_places(spectral_norm(N)) == "1.274219991"

test "the matrix is the one the task describes, entry by entry"
    # The description gives the first entries directly: a11 = 1, a12 = 1/2,
    # a21 = 1/3, a13 = 1/4, a22 = 1/5, a31 = 1/6 — one-based there, and
    # zero-based here, which is the one translation this program makes.
    assert eval_a(i: 0, j: 0) == 1.0
    assert eval_a(i: 0, j: 1) == 1.0 / 2.0
    assert eval_a(i: 1, j: 0) == 1.0 / 3.0
    assert eval_a(i: 0, j: 2) == 1.0 / 4.0
    assert eval_a(i: 1, j: 1) == 1.0 / 5.0
    assert eval_a(i: 2, j: 0) == 1.0 / 6.0

test "the entries walk the anti-diagonals, and every one of them is positive"
    # The pattern the description names: the matrix is filled diagonal by
    # diagonal, so every entry on one anti-diagonal is bigger than every
    # entry on the next.
    i: i64 @= 0

    while i < 20
        j: i64 @= 0

        while j < 20
            assert eval_a(i: i, j: j) > 0.0
            assert eval_a(i: i, j: j) <= 1.0
            j @ j + 1

        i @ i + 1

    assert eval_a(i: 0, j: 1) > eval_a(i: 0, j: 2)
    assert eval_a(i: 5, j: 5) > eval_a(i: 5, j: 6)

test "the transpose really is the transpose"
    # Not an identity this program could get away with assuming: `at_times`
    # is a separate loop with the indices swapped, so a typo in one of them
    # would leave both plausible.
    i: i64 @= 0

    while i < 12
        j: i64 @= 0

        while j < 12
            assert eval_a(i: i, j: j) == eval_a(i: j, j: i) || i != j
            j @ j + 1

        i @ i + 1

    # And at the level the program actually uses them: A·u and Aᵀ·u differ
    # on a vector that is not symmetric under the swap.
    v: [f64] = [1.0, 2.0, 3.0, 4.0]
    assert a_times(v) != at_times(v)

test "the iteration converges rather than wandering"
    # Two extra rounds must not move the answer at nine places: that is what
    # says ten rounds is enough, and a program whose arithmetic drifted
    # would fail here while still printing something plausible.
    ten = spectral_norm(40)
    u: [f64] @= ones(40)
    v: [f64] @= ones(40)
    round: i64 @= 0

    while round < ROUNDS + 2
        v @ ata_times(u)
        u @ ata_times(v)
        round @ round + 1

    vbv: f64 @= 0.0
    vv: f64 @= 0.0
    at: i64 @= 0

    while at < 40
        vbv @ vbv + u[at] * v[at]
        vv @ vv + v[at] * v[at]
        at @ at + 1

    assert nine_places(sqrt(x: vbv / vv)) == nine_places(ten)

test "a smaller matrix gives a smaller norm, monotonically"
    # The norm grows with N toward its limit, so a program that ignored part
    # of the matrix would break this without breaking the published number
    # at one size.
    assert spectral_norm(10) < spectral_norm(50)
    assert spectral_norm(50) < spectral_norm(100)
    assert spectral_norm(100) < 1.3

test "nine places prints nine places, and keeps the leading zeroes"
    assert nine_places(0.0) == "0.000000000"
    assert nine_places(1.5) == "1.500000000"
    assert nine_places(0.25) == "0.250000000"
    assert padded(1) == "000000001"

test "sqrt comes from C, and it is the C answer"
    assert sqrt(x: 4.0) == 2.0
    assert sqrt(x: 1.0) == 1.0