# Viv prompt pack

> Everything a model needs to write working Viv, in one file. Generated from https://viv-lang.org/reference, the builtin table and the documented examples; every example's output below is what the interpreter printed.

Viv is a small language for numerical programs. A number is a 64-bit
float (`Num`) unless it is declared otherwise: `n: Int` is a 64-bit integer,
`Float32` an IEEE binary32 number, and `2.0 + 3.0i` a complex one. Arrays are
rank 1 to 4, and programs end in `test` blocks that check their claims. One
program runs in three engines: `viv run f.viv` interprets it (the site runs the
same interpreter compiled to wasm), `viv f.viv -o f.c` compiles it to C99, and
`viv f.viv -o f.wasm` compiles it to wasm. All three are meant to print the
same bytes. Figures are records on stdout, so a plot is part of the output and
needs no library.

How to use this pack: write Viv, not Python. The grammar is short and differs
from Python in the ways the "Python habits" table lists; the examples show
every feature with the exact output it prints.

## 1. Grammar and types

### Lexical structure

- **Numbers**: A number is a 64-bit float. `10` lexes as `10.0`, scientific notation works (`1e-10`), and a 64-bit integer is an `Int`, declared with an annotation (`n: Int`). `%` and `//` follow Python's rules for floats (see operators (#operators)); `mod(a, b)` is C's fmod, whose result takes the sign of `a`.

- **No implicit multiplication**: A name cannot follow a number directly: `2x` is an error that says to write `2 * x`. The `i` of an imaginary literal (`2.5i`) and the `e` of an exponent (`1e-10`) are part of the number.

- **Names**: A name starts with a letter or `_` and continues with letters, digits and `_`. Letters are Unicode letters, so `Δx` and `σ` are names. Names that start with `viv__` are used by the parser for the hidden variables of `a[i] += x` and the longer `for` forms.

- **The dot rule**: A `.` is part of a number only with digits on both sides. So `1.5` is a float, `x.y` is field access, and `1.` is neither.

- **Statements end at the newline**: No semicolons. Inside parentheses or brackets the lexer suppresses newlines, so long expressions wrap freely.

- **Comments**: `#` to end of line, dropped by the lexer. The annotation markers — `#|` what, `#||` how, `#|||` why, `#~` tags, `#:` variable glossary, `#[`/`#]` regions — are comments the compiler ignores and `viv doc` extracts. Above a `def` they form its structured doc.

- **Strings**: `"text"`, with the escapes `\n` `\t` `\r` `\0` `\\` `\"` `\'` and `\u{1F600}` (any other backslash is an error), and f-strings, `f"x = {x:.3f}"`. A string is a value: see Strings (#strings).

- **Docstrings**: A bare string as the first statement of a `def` documents it and compiles to nothing.

### Statements

| Form | Notes |
|---|---|
| `x = expr` | Assignment. At top level it creates a global; inside a `def`, a local hoisted to function scope. |
| `x += e` `-=` `*=` `/=` | Compound assignment; desugars to the plain form, `x = x + e`. |
| `a[i, j] += e` `-=` `*=` `/=` | Compound assignment on an array element: `a[i, j] = a[i, j] + e`, with each index evaluated once. An index that is not a number or a variable is computed into a hidden variable first, so `c[floor(10 * rand())] += 1` draws one random number, not two. |
| `a[i, j] = expr` | Indexed store, up to rank 4. The base may be any array-valued expression, including a record field. |
| `a[lo:hi] = v` `A[:, j] = v` | A write into a slice (#slices): `v` is a number, which every selected element takes, or an array of exactly the slice's shape. The value is copied before the write, so `s[1:5] = s[0:4]` shifts. `a[lo:hi] += v` and the other compound forms read the slice, combine and write it back, with each bound evaluated once. |
| `return expr` | A `def` that ends without a `return` gives nothing: it is called as a statement, and the compiler refuses a use of its value (`show f(x)`, `f(x) + 1.0`). |
| `return a, b` | Several results. Every `return` of the def gives the same number of values, each with one type (a number, a boolean, a complex number, an array, a record); `return g(x)` passes on the values of a def `g` that returns as many. The caller unpacks them, `a, b = f(x)`, on a line of its own, or calls `f(x)` alone and drops them. `f(x) + 1`, `show f(x)` and `f` as a function value are errors that name `f`. The def returns a record the compiler declares for it (Records (#types)), so every engine runs it. |
| `a, b = e` | Unpacking. `e` is a call of a def that returns as many values, a list of values (`a, b = b, a` swaps: every value on the right is computed first, then the targets are assigned left to right), or an array with exactly one element per name, such as `B, C, H, W = shape(X)`; an array of another length stops the run, located at the line. A target is a name, an element `v[i]` or a slice `v[lo:hi]`, and its indices are computed after the right side, as in Python: `v[i], v[j] = v[j], v[i]` swaps two elements. `loss, g = value_and_grad(...)` is this form too. |
| `if c { } else { }` | With `else if` chains. |
| `while c { }` | With `break` and `continue`. |
| `for i to n { }` | `i` runs 0 … n−1. |
| `for i from a to b { }` | `i` runs a, a+1, … while it is below b. As in `for i to n`, `to` names the end and the end is left out. |
| `for i from a to b by s { }` | `i` runs a, a+s, a+2s, … and stops before reaching b: below b for a positive s, above b for a negative one (`for i from 10 to 0 by -3` is 10, 7, 4, 1). The loop makes ceil((b − a) / s) passes, none when that is 0 or less, and pass k has `i = a + k·s`. a, b and s are evaluated once, in that order. A step of 0 is an error: at compile time when it is a number, at run time otherwise. |
| `for i to n reversed { }` | The same values in the opposite order: `for t to T reversed` is T−1, …, 0 and `for i from 1 to n reversed` is n−1, …, 1. `reversed` goes last, after any `from` and `by`. These longer forms are rewritten by the parser into `for k to count` over a hidden counter, with `i` assigned at the top of each pass, so every engine runs them. |
| `parallel for i to n { }` | Sequential in the interpreter; the C backend emits an OpenMP pragma. Only this form: no `from`, `by` or `reversed`. Each pass may write only the elements of a shared array indexed by `i` (`a[i]`, `m[i, j]`, `m[:, i] = col`, `setslab(W, i, mine)`, or a def given `i` as the index it writes at), and may read a shared array it writes only at those elements; a write two passes could share, a write of a global through a def, a `break` or a `return` is refused before the run, in every engine. Arrays a pass makes are its own, and are freed when the pass ends. |
| `show expr` | Prints the value. A variable, or `fmt` of one, is labeled with its name (`x: 2`); any other expression prints bare, as `print` would. |
| `show a, b, c` | Several labeled values on one line, two spaces apart: `a: 1.5 b: 0.3 v: [1, 2]`. A variable, or `fmt` of one, is labeled with its name, so `show fmt(x, 17)` prints `x: 0.30000000000000004`. Arrays print on one line here. It is `print("a: ", a, " b: ", b, ..., sep="")`; see Output (#output). |
| `print(...)` | Text and values on one line, a space between them; see Output (#output). |
| `for x in e { }` | Each element of `e` in order: a string's characters, a list's strings, a vector's numbers. `for i, x in enumerate(e)` gives the index too, and `enumerate(e, 1)` (or `start=1`) counts from 1. `for a, b in zip(u, v)` walks several sequences together and stops at the end of the shortest. `for row in rows(M)` gives the rows of a matrix, `M[k, :]` (an array's rank is known only when the program runs, so `for x in M` reads numbers and stops on a matrix). `e` is evaluated once; the parser writes the loop as `for k to len(e)` with `x = e[k]` at the top of each pass. |
| `stop("text")` | Ends the run with `error: text` and the line it is on, exit status 1, in every engine: the error a failed check should give. The argument is one string, a literal or a value (`stop(f"no atom {name}")`). Output printed before it stays printed (the wasm engine reports the error alone, as for every run-time error). |
| `narrate "text"` | A progress line. A non-string argument is commentary and compiles to nothing. |
| `test "name" { }` | Runs the body and reports under the name: one `PASS name` or `FAIL name` line, however many asserts are inside. A failure makes the process exit nonzero. |
| `def f(a, b, c = 1.0) { }` | A parameter with a default comes after those without and is passed by name only: `f(1, 2, c=0.5)`; a default may read the parameters before it (`scale = 1.0 / sqrt(rows)`) and is computed at each call that leaves it out. The front end fills every call out, so it costs nothing at run time; a call through a function value passes its arguments in order. Optional annotations: `def f(a: Num, v: Arr) -> Num`. Recursion works to a depth of 2500 nested calls, the same in every engine; a call past that stops the run with `recursion deeper than 2500 nested calls`, located at the def. A def and a top-level variable cannot share a name, even when the def comes from a `#! use` module: the compiled program has one namespace for both, so every engine refuses it and names the two lines. |
| `step name "text" { }` | A named phase of a def (or of the top level), with optional text. A step is a block, not a scope: its variables are the def's. The engines record entering and leaving it, and the explorer, the story and the intent and debug lenses address it as `def/name`, so its name is unique in the def. Not inside a `parallel for`. |
| `step name: statement` | A step of one statement, on that statement's line: `step attn_residual: X1 = X + Y` is `step attn_residual { X1 = X + Y }` to every engine and tool, hashes included. The statement starts after the `:` on the same line and ends where a statement ends, so `step bump: for i to n {` runs to the loop's closing brace, and an `else` on the next line continues an `if`. `;` is a line break, so `step a: x = 1; y = 2` holds only `x = 1`. With text: `step name "text": statement`. |
| `### … ###` | A fenced block of text: everything between two lines that are exactly `###` is skipped, so prose needs no `#` in front of it. A post (https://viv-lang.org/posts) keeps its title and its notes in one. |
| `rung "title" { }` | One step of a ladder (https://viv-lang.org/ladders): defs, tests and statements that build on the rungs before it. A def in a later rung replaces an earlier one of the same name. `viv run` runs the top rung; `viv ladder` runs every rung and checks the ladder's rules. |
| `record P { x: Num, v: Arr }` | Record type (a user-defined struct), in every engine. Build one with every field named, `p = P { x: 1, v: zeros(3) }`; read a field with `p.x`; assign one with `p.x = 2`; write into an array field with `p.v[i] = e`. A record is a value: `q = p` gives `q` its own, and `q.x = 3` leaves `p` as it was (an array inside is shared until it is written into, as arrays are). A write is into a variable of the def it is in: a def does not change a record it was given, or a global. |
| `record M { blocks: [B] }` | An array of records: `[b1, b2]`, `append(bs, b)`, `bs[i]`, `len(bs)`, `bs[i] = b`, and a place through it, `m.blocks[l].mlp.up = W` or `m.blocks[l].mlp.up[:] = W`. A value, as each record in it is. |
| `a + b`, `zeros_like(r)`, `for_leaves(f, r1, r2)` | Over a record of arrays (every field an array, a record of arrays, or an array of them): `+ - * /` with another of its type or a number, `-r`, `sqrt(r)` and the other elementwise functions, act on each array; `zeros_like(r)` is one of the same shape; `sum(r)` adds its arrays' sums; `for_leaves(f, r1, ..., x)` calls `f` on each array of the records at once, with `x` after them, and `f` may write them in place. `loss, g = value_and_grad(f, model, ...)` differentiates a record of arrays and gives `g` of its type. A number field is refused, with its name: keep sizes in a record of their own. |
| `fn(x) -> expr` | A lambda, and a value: assign it, pass it, call it through the parameter. |

### Output with `print` and `fmt`

`print` writes its arguments on one line with a space between them, as Python's does, then a newline. `sep=", "` puts other text between them and `sep=""` none; `sep` is a string literal or a variable. There is no `end=`: every `print` ends its line. An argument is a string (written as is), `fmt(x, ...)` (a value in a format), or any other value (as `show` writes it: six significant digits, the way C's `%g` does). Every value is computed before anything is written. The interpreter, C and wasm print the same bytes, checked over 10^6 doubles (`tools/print-fmt-sweep/`).

| Form | Writes | Python |
|---|---|---|
| `print("loss", loss)` | `loss 0.693147` | `print(f"loss {loss:g}")` |
| `print("loss", fmt(loss))` | `loss 0.6931471805599453`: the shortest digits that read back as the same double | `print("loss", loss)` (`repr`) |
| `print(a, b, c, sep=", ")` | `1, 2, 3` | the same |
| `print("[", x, "]", sep="")` | `[0.5]` | the same |
| `fmt(x, 17)` | 17 significant digits, C's `%.17g`; any count 1 to 17 | `{x:.17g}` |
| `fmt(x, ".6g")`, `".3e"`, `".2f"` | C's `%.6g`, `%.3e`, `%.2f` | the same spec |
| `fmt(x, "10.4f")` | at least 10 characters, padded on the left | the same spec |
| `fmt(n, "d")`, `"5d"` | a whole number, no decimal point (`%.0f`) | `{n}`, `{n:5d}` for an int |
| `print(v)`, `fmt(v, ".2f")` | an array on one line, `[1, 2.5]`, `[[1, 2], [3, 4]]`, each element in the format | — |
| `print(flag, z)` | `true 1+2i`: booleans and complex numbers, bare only | — |
| `print(x)` | one number: the one-argument `print`, unchanged | — |
| `print()` | an empty line | `print()` |

A format is `[width][.precision]` then `e`, `f`, `g` or `d`, or empty for `repr`; precision and width are at most 40 and 100. The rest of Python's spec (fill, alignment, `+`, `,`, `%`, a precision with no type) is refused with its location and what to write instead. The digit count and the format are written in the source, not computed. `nan` prints as `nan` whatever its sign, and `fmt` keeps the sign of a negative zero (`-0.0`, `-0`) where a bare value writes `0`. A program's own `def fmt` or `def print` replaces the builtin, and no separator is added to its calls.

Until 2026-09-30 `print` put nothing between its arguments, and programs wrote `print(a, " ", b)`. `tools/migrate_print_sep.py` rewrote every program in the repository to the new form, checked call by call to reach the engines unchanged, so each prints the bytes it printed before.

### Operator precedence

| Level | Operators | Notes |
|---|---|---|
| 1 | `x ^ y` | power, the same as `pow(x, y)` (so an array base works elementwise). Right-associative, and tighter than a minus on its left: `-x ^ 2` is `-(x ^ 2)`, `2 ^ 3 ^ 2` is 512, `2 ^ -1` is 0.5 |
| 2 | `-x` `not x` | unary |
| 3 | `*` `/` `@` `%` `//` | `*` on two arrays is elementwise; `@` is always the matrix product (NumPy's `matmul`: a stack of matrices multiplies matrix by matrix, under Axes (#arrays)). `a % b` is the floored remainder and `a // b` the floor of the quotient, Python's rules for floats: the remainder takes the sign of `b`, so `-7 % 3` is 2 and `(0 - 1) // 2` is -1, and `a == (a // b) * b + a % b`. They take numbers; a zero divisor stops the program with `division by zero in a % b` |
| 4 | `+` `-` | elementwise on arrays; all four arithmetic ops broadcast between an array and a scalar, either side, and between two (not complex) arrays by NumPy's rule (below (#arrays)) |
| 5 | `<` `<=` `>` `>=` `==` `!=` `in` `not in` | `x in a`: a number in a real array of any rank (some element `==` x, so a NaN is in nothing, as in NumPy), a key in a map, a substring in a string, a string in a list of strings. On numbers, the comparisons give a Bool. With a real array on either side they compare elementwise and give an array of 1.0 where the comparison holds and 0.0 where it does not, so `dY * (Y > 0.0)` is ReLU's adjoint on any rank. `all`, `any` and `where` read such a mask; `if a == b` on arrays is refused, with a note that says to write `all(a == b)` |
| 6 | `and` `or` | short-circuit in every engine: the right side runs only when the left does not decide, so `j > 0 and a[j - 1] > v` never reads `a[-1]` |
| 7 | `a if c else b` | a when c is true, else b; only the one chosen is computed. Loosest of all, as in Python: `a if c else b if d else e` is `a if c else (b if d else e)` |

Measured: `2.0 + 3.0 * 4.0` is 14, `-7 % 3` is 2 and `7 // -2` is -4. Parentheses override everything.

### Comprehensions, `if` expressions and keys

A comprehension builds an array, or a list of strings, from loops, as Python's does. Its `for` takes any head a `for` statement takes, and its `if` keeps the elements it holds for:

```viv
xs = [3, 1, 4, 1, 5, 9]
evens = [x for x in xs if x % 2 == 0]           # [4]
squares = [i * i for i to 5]                    # [0, 1, 4, 9, 16]
table = [[i * j for j from 1 to 4] for i from 1 to 4]
words = [upper(w) for w in split("a bb ccc")]   # ['A', 'BB', 'CCC']
print(" ".join(f"{x:.2f}" for x in xs[:3]))     # 3.00 1.00 4.00
print(all(x > 0 for x in xs), max(words, key=len))
sign = "neg" if xs[-1] < 0 else "pos"
print(sorted(words, key=len, reverse=true))
```

| Form | Gives |
|---|---|
| `[e for x in xs]`, `[e for i to n if c]` | a vector of numbers, a matrix when every `e` is a vector of one length (a row each), or a list of strings. More `for`s nest, the first outermost; `xs` may be an array, a string, a list of strings, `enumerate(...)`, `zip(...)` or `rows(M)`. An empty result is `[]` |
| `{k: v for x in xs}` | a map (#maps), Python's dict comprehension: `stoi = {c: i for i, c in enumerate(vocab)}`. A key seen again keeps its first place and takes the later value. The same heads and filters as a list's. A brace in an expression is only this: a map literal is `map()` and stores, and a set (`{x, y}`) is refused with a pointer to `unique` |
| `f(e for x in xs)` | a generator as the only argument of a call: `f([e for x in xs])`, so `sum`, `max` and `" ".join` take one |
| `all(c for x in xs)`, `any(...)` | a Bool, and the loop stops at the first element that decides it, as Python's does |
| `a if c else b` | a when c holds, else b; only the one chosen is computed |
| `max(xs, key=f)`, `min(xs, key=f)` | the first element whose key is largest (smallest), as Python's; a NaN key is never chosen after the first element; an empty `xs` stops with `max() arg is an empty sequence` |
| `sorted(xs, key=f, reverse=r)`, `argsort(xs, key=f, reverse=r)` | the elements (or their positions) in the order of their keys, numbers or strings: stable, so equal keys keep their order, with `reverse=true` too, as Python's; a NaN key goes after every number |

A key is a def's name, a builtin's (`len`, `abs`), a variable holding a function, or a lambda `fn(w) -> ...`, which may read the program's variables. A comprehension's variables are its own: `x` outside it is neither read nor changed, as in Python 3. The parser writes each of these forms as the loop it stands for, computed just before the statement that uses it, so every engine runs code it already runs; what a statement computes before it (a call of a def, a `print`) still runs first, and in a `while` condition it is computed again each time. For that reason one cannot be the body of a `fn(...) -> ...` (except a key's): write a def.

### Axes, broadcasting, and agreement with NumPy

The array builtins take NumPy's names and give NumPy's results: reshapes and transposes, `eye`, `diag`, `arange`, reductions along an axis, gathers and masks, `cumsum`, `argsort`, `einsum`, and solvers that return what they compute. Their tables are under Builtins (#builtins); the rules they share are here.

```viv
X = reshape(arange(12), 3, 4)
Xc = X - mean(X, 0)                 # the column means, broadcast down the rows
C = transpose(Xc) @ Xc / 3.0
vals, vecs = eig_sym(C)
w = lin_solve(C + eye(4), ones(4))
top = argsort(vals)
picked = take(X, [2.0, 0.0], 0)     # X[[2, 0], :]
big = extract(X > 5.0, X)           # X[X > 5]
```

- **Broadcasting**: `+ - * /`, the comparisons and two-array `max` and `min` between two arrays of different shapes (not complex ones) line up their last axes; each pair of lengths is equal, or one of them is 1 and repeats. So a row vector against a matrix repeats down its rows, `reshape(v, n, 1)` is a column that repeats across, and shapes that do not fit stop with the words they always did: `array shape mismatch: [2, 3] vs [3, 2]`. Complex arrays do not broadcast. `grad` sums a broadcast operand's gradient back to its shape.

- **Matrix products**: `A @ B` is NumPy's `matmul`. Two matrices give their product; a vector on the left is a row and on the right a column, and the axis it adds is dropped. An array of rank 3 or 4 is a stack of matrices in its last two axes, and the axes before them broadcast, so `X @ W` with `X` of shape (4, 8, 16) and `W` (16, 8) multiplies each of the four matrices by `W`. Two vectors are refused (`dot(u, v)` is their product). Each entry adds its terms in order of the inner index from 0.0, skipping a term whose left factor is zero, with a separate multiply and add (no fused multiply-add); the loop is blocked for the cache without changing that order, so every engine gives the same bits (`lang/viv/src/interp/matmul.rs`). NumPy calls BLAS, which sums in its own order and agrees to rounding.

- **Axes**: An axis is 0 to rank − 1, or −rank to −1 counted from the end, as in NumPy and as an index is. A reduction along an axis (`axis=`) drops it; with `keepdims=true` it keeps it with length 1, as NumPy's does, which is what lets `A - mean(A, axis=1, keepdims=true)` subtract each row's mean. A builtin's options are always given by name, with NumPy's names (`axis=`, `keepdims=`, `k=` for `triu`/`tril`/`diag`, `offset=` for `trace`, `endpoint=`, `axes=`); written by position they are refused. A reduction along the only axis of a vector would give a number, and is refused with the whole-array form to write instead (`mean(v)`).

- **Orders of summation**: `sum`, `mean`, `var` and `std`, along an axis or of a whole array, and `trace` add in NumPy's order: from 0.0, a contiguous run in blocks of 8192 entries, each block by NumPy's pairwise sum; along an axis that is not the last, the slabs in order. They agree with NumPy to the bit (`tools/py2viv/probe_numpy_order.py` states the order), so `mean(a)` is `sum(a) / numel(a)` to the bit. The order binds the exact engines; fast mode may add in another. `prod`, `cumsum` and `cumprod` go left to right, as NumPy does. `dot` and `norm` go left to right; NumPy calls BLAS, which agrees to rounding. `einsum` loops over the output entries in row-major order and over the summed indices in their order of first appearance, adding each product to 0.0; NumPy's einsum has its own loops, so the two agree to rounding.

- **NaN and signed zeros**: `amax`, `amin`, `argmax` and `argmin` are NumPy's: a NaN wins, and `amax` puts 0.0 above −0.0. (`max(a)` of one array is older and skips NaN, as C's `fmax` does.) `argsort` puts NaN last and keeps −0.0 and 0.0 in their order. `clip` keeps a NaN entry and treats a NaN bound as no bound.

- **Masks and indices**: A mask is an array of 1.0 and 0.0, as a comparison gives, so an array cannot say whether it holds truth values or indices. NumPy's `a[mask]` is `extract(mask, a)`, `a[mask] = v` is `putmask(a, mask, v)`, `a[idx]` along an axis is `take(a, idx, axis=k)`, and `np.flatnonzero` is `flatnonzero`. The indices `take` and `put` read are not negative: `a[-1]` counts from the end, `take(a, [-1.0], 0)` stops.

- **Solvers that return**: `x = lin_solve(A, b)` (a matrix `B` solves for every column), `inv(A)`, `L = cholesky(A)`, `vals, vecs = eig_sym(A)`, `w, V = eig_herm(A)` and `U, s, Vt = svd(A)` run the same code as the forms that fill arrays you pass, which still work, and give the same bits. The compiler writes an unpacking of `eig_sym`, `eig_herm` or `svd` as that call with arrays it allocates. `svd` is one-sided Jacobi (`lib/svd.viv`'s algorithm); its vectors can have other signs than LAPACK's.

`tools/arrays/numpy_cases.py` writes 4,812 cases whose results NumPy computes in an order Viv follows, and 89 that agree with NumPy within a tolerance (`einsum`, the solvers, `dot`); `cargo test --test arrays` runs them in the interpreter, C and wasm, which print NumPy's results byte for byte.

### Types

- **Num**: f64, and the default: a literal is a Num. The other scalars are declared: `Int` (a 64-bit integer, `n: Int`), `Float32` (IEEE binary32) and `Cplx` (below). A double holds every integer up to 2^53 exactly, so counters, indices and `n == 1.0` are exact.

- **Bool**: `true` or `false`: what a comparison yields and what `if` and `while` take. It can be stored, passed and returned (`def f(x: Num) -> Bool`), combined with `and`, `or`, `not`, and compared with `==` and `!=`. It is **not a number**: `true + 1.0` and `if count {` are refused, identically, by the interpreter and the compiler. `show` prints it as `true` or `false`. Arrays hold numbers, so a mask is still an array of 1.0 and 0.0: a comparison with an array on either side makes one, and `all(m)` and `any(m)` turn one back into a Bool. See booleans (https://viv-lang.org/docs/booleans).

- **Str**: An immutable UTF-8 string, indexed by code point (Strings (#strings)). A list of them is `[Str]`.

- **Arr**: Rank 1 to 4, row-major, f64 elements. Arrays pass into functions and return from them, and scalars broadcast against them from either side, and two arrays (not complex ones) broadcast as NumPy's do (Axes, broadcasting, NumPy (#arrays)). **Assignment and parameter passing alias:** `b = a` shares storage, and a function that writes into an array parameter writes into the caller's array. A real copy is `copy(a)`, and a slice (#slices) is a copy too — see scope & mutation (https://viv-lang.org/docs/scope). Indices count from 0, and an index outside its dimension stops the run in the interpreter and in C with the same message, `index 2 out of bounds for dimension 0 (size 2)`, and the line it happened on. A negative index counts from the end, as in Python: `a[-1]` is the last element, and `a[-3]` of two elements stops the run. A fractional index is cut toward zero, so `a[1.7]` is `a[1]`.

- **Map**: Numbers under number keys, in the order the keys were first stored, made by `map()`. A handle: assignment, parameters and record fields share one map. A parameter may be annotated `m: Map`. See Maps (#maps).

- **Cplx**: A first-class complex scalar: the `i` suffix on a numeric literal (`2.0 + 3.0i`), full arithmetic with real promotion, and complex overloads of `abs/exp/log/sqrt/sin/cos/tan`. Compiled to C99 `double complex`. Every engine computes these the same way, to the bit, with Viv's own formulas rather than C's operators or `<complex.h>`: A real operand is lifted to `(x, 0.0)` first, so `2.0 * z` is the complex product with its zero terms. `+` and `-` work part by part. `a * b` is `(ar br − ai bi) + i (ar bi + ai br)`, each product rounded on its own (no fused multiply-add). `a / b` is Smith's algorithm (1962): divide through by the larger part of `b`. If `|br| ≥ |bi|`, then with `r = bi/br` and `d = br + bi r` the quotient is `((ar + ai r)/d, (ai − ar r)/d)`; otherwise `r = br/bi`, `d = bi + br r`, and the quotient is `((ar r + ai)/d, (ai r − ar)/d)`. It never squares `b`, so it overflows only where the quotient does, not once `|b|` passes 1e154 as the textbook `a conj(b) / |b|²` does. A real divisor divides each part once, so `(1 + 2i) / 3` is exactly `1/3 + (2/3)i`. Dividing by zero gives `nan + nan i`. `abs(z)` is `hypot(re, im)`. `exp` is `e^re (cos im + i sin im)`. `log` is `ln|z| + i atan2(im, re)`. `sqrt` is `√|z|` at half the argument. `sin` and `cos` use the `cosh`/`sinh` formulas, and `tan` is `sin z / cos z` by the division above. All of them use Viv's libm. A complex array divided by a scalar divides each element the same way.

- **CArr**: An array of complex numbers, rank 1 to 4: `czeros(n)`, or a literal with a complex element in it (`[1.0 + 2.0i, 3.0]` — the reals promote). Indexing gives a `Cplx`; storing a real promotes it. It is **packed**: interleaved (re, im) doubles, which is exactly how C99 lays out `double complex`, so both engines walk the same memory and `fft` agrees between them to the byte. A complex and a real array are different types and do not mix silently — `sin(a)` on a complex array is an error. Annotate a parameter as `CArr` when you want to say so; otherwise it is inferred from the call site.

- **Records**: User-defined structs with typed fields, copied by value: the C backend compiles them to C structs, the wasm backend to a block of memory, and the interpreter holds the fields. Array fields refer to the original array, so bundling a mesh into one value copies no array data. A def that takes or returns one annotates it (`def f(p: P) -> P`). Records have no operators and no `show`: compare and print their fields. The several results of `return a, b` are a record the compiler declares for the def, with one field per value and the field types inferred from the returns; the program never names it, because the caller unpacks it.

- **Function values**: A `def` name or a lambda is a value. The compiled backends implement higher-order calls by monomorphization — `run_it(sq, x)` specializes to `run_it__sq` — with `viv_fn` function pointers for the cases where a function is held in a variable. Arrays of functions work in the interpreter and are a compile error in C.

### What a test prints

A `test` block is the unit that passes or fails, so it prints exactly one line under the name you gave it — not one line per assertion:

```viv
test "simpson integrates exp correctly" {
    assert_near(simpson(fn(x) -> exp(x), 0.0, 1.0, 60), exp(1.0) - 1.0, 1.0e-9)
}

$ viv run integrate.viv
PASS simpson integrates exp correctly
```

Several asserts in one block still give one line. A failure names the test and carries the first failing assertion's detail, and a block with no assertions says so instead of printing PASS:

```
PASS two asserts, one verdict
FAIL a failing one -- assert_near: |1 - 2| = 1 > 1e-09
PASS nothing asserted (no assertions)
assert_near PASS
1 assertion FAILED
```

That last `assert_near PASS` is an assertion at top level, outside any `test` — there is no name to report under, so it prints its own line as before. The trailing count is the summary, and it is what makes the process exit nonzero, so a build script can fail on the exit code. It goes to stderr, after everything the program printed, in the interpreter and in a compiled program alike, so stdout holds only the program's own output.

## 2. Python habits that do not work

Each message below is what the interpreter printed for the input on the left; each fix on the right runs.

| Habit | Python-style input | What Viv says | Write instead |
|---|---|---|---|
| power | `x = 2.0 ** 3.0` | error: expected an expression | `x = 2.0 ^ 3.0` |
| range loop | `for i in range(3) { ⏎ print(i) ⏎ }` | error: expected `to` | `for i to 3 { ⏎ print(i) ⏎ }` |
| colon blocks | `def f(x): ⏎ return x * x` | error: expected `{` | `def f(x) { ⏎ return x * x ⏎ }` |
| elif | `x = 1.0 ⏎ if x > 2.0 { ⏎ print(1.0) ⏎ } elif x > 0.0 { ⏎ print(2.0) ⏎ }` | error: a brace in an expression is a map comprehension, `{k: v for x in xs}` | `x = 1.0 ⏎ if x > 2.0 { ⏎ print(1.0) ⏎ } else if x > 0.0 { ⏎ print(2.0) ⏎ }` |
| imports | `import numpy as np` | error: line 1: `numpy`: no module by that name (looked for ./numpy.viv, lib/numpy.viv, numpy.viv) | `# no numpy to import: every builtin is always in scope; a module of lib/ is imported as Python does it, from nn import rmsnorm` |
| np. prefix | `a = np.zeros(3)` | error: undefined variable `np` | `a = zeros(3)` |
| chained indexing | `A = zeros(2, 2) ⏎ show A[1][1]` | error: array has rank 2 but was indexed with 1 index/indices | `A = zeros(2, 2) ⏎ show A[1, 1]` |
| .T and methods | `A = zeros(2, 2) ⏎ B = A.T` | error: cannot access field `T` of a array | `A = zeros(2, 2) ⏎ B = transpose(A)` |
| list.append | `a = [] ⏎ a.append(1.0)` | error: `append` gives back the longer array, so it is written `a = append(a, x)`, with the same name on both sides | `a = [] ⏎ a = append(a, 1.0)` |
| dict | `d = {} ⏎ d["a"] = 1.0` | error: a brace in an expression is a map comprehension, `{k: v for x in xs}` | `d = map() ⏎ d["a"] = 1.0` |
| a def name to simpson | `def sq(x) { ⏎ return x * x ⏎ } ⏎ show simpson(sq, 0.0, 1.0, 10)` | error: `simpson`'s first argument must be a `fn(x) -> ...` lambda, not the def `sq` (wrap it: `fn(x) -> sq(x)`) | `def sq(x) { ⏎ return x * x ⏎ } ⏎ show simpson(fn(x) -> sq(x), 0.0, 1.0, 10)` |
| lambda | `f = lambda x: x * x` | error: expected `=` after the type | `f = fn(x) -> x * x` |
| True / False / None | `x = True` | error: undefined variable `True` | `x = true` |

How to write it so it runs:

- `show name` prints `name: value`; `show` of any other expression prints the bare value, and so does `print`. Numbers print with 6 significant digits.
- Unary minus works anywhere: `-1.0`, `exp(-t)`, `[-1.0, 2.0]`. (Some examples write `0.0 - x`; both forms are the same.)
- The solvers return their results: `x = lin_solve(A, b)`, `L = cholesky(A)`, `vals, vecs = eig_sym(A)`, `U, s, Vt = svd(A)`. (The older forms that fill arrays you allocate, `eig_sym(A, vals, vecs)`, still work.)
- NumPy's array functions have their NumPy names and results: `reshape(a, 2, 3)`, `sum(A, 0)`, `mean(A, 1, true)` (keepdims), `argmax`, `eye`, `outer`, `einsum("ij,jk->ik", A, B)`, and two arrays broadcast. `A.T` is `transpose(A)`, `a[mask]` is `extract(mask, a)`, `a[idx]` is `take(a, idx, 0)`, `np.max(a, axis=0)` is `amax(a, 0)` (`max(a, b)` is elementwise).
- Set a test's tolerance from an error estimate, or run the program once and read the printed error first. In a trial of this pack (five programs written from it) the one failure was a tolerance guessed too tight (Simpson's rule with 8 panels on exp over [0, 1] is off by 2.3e-6, and the test asked for 1e-6).
- A program should pass in both engines: `viv run f.viv` and `viv f.viv -o f.c && cc -O2 f.c -lm && ./a.out` print the same lines. The interpreter accepts a few things the compiler refuses; section 3 lists them.

## 3. The interpreter runs these; the C compiler refuses them

**deriv of a function parameter.** deriv, simpson and trapezoid take a lambda written at the call. Write the function inside the deriv call, or pass the derivative in as a second function parameter. (deriv cannot see inside a call to a def either.)

```viv
def newton(f, x0) {
    x = x0
    for k to 20 {
        x = x - f(x) / deriv(f, x)
    }
    return x
}
show newton(fn(x) -> x * x - 2.0, 1.0)
```

`viv f.viv -o f.c` says: `` `deriv`'s first argument must be a `fn(x) -> ...` ``

Instead:

```viv
def newton(x0) {
    x = x0
    for k to 20 {
        x = x - (x * x - 2.0) / deriv(fn(t) -> t * t - 2.0, x)
    }
    return x
}
show newton(1.0)
```

## 4. Builtins

All 267 builtins run in both the interpreter and the C backend. Figures need `figure(...)` before the first series.

### Real math

Every one-argument function here also maps elementwise over an array, so `y = sin(x)` needs no loop.

- `exp(x)`: e^x
- `log(x)`: natural log
- `ln(x)`: alias for log
- `sqrt(x)`: square root
- `pow(x, y)`: x to the y
- `abs(x)`: magnitude
- `fabs(x)`: alias for abs
- `sin(x)`
- `cos(x)`
- `tan(x)`
- `asin(x)`
- `acos(x)`
- `atan(x)`
- `atan2(y, x)`: quadrant-correct angle
- `sinh(x)`
- `cosh(x)`
- `tanh(x)`
- `erf(x)`: error function
- `gamma(x)`: Γ, C’s tgamma
- `floor(x)`
- `ceil(x)`
- `round(x)`
- `min(a, b) · min(a, b, c, ...)`: the smallest; three or more arguments are min(min(a, b), c)
- `max(a, b) · max(a, b, c, ...)`: the largest, as min
- `hypot(x, y)`: sqrt(x^2 + y^2) with no overflow on the way (musl's, as abs of a complex number)
- `copysign(x, y)`: |x| with the sign of y; -0.0 is negative
- `fma(a, b, c)`: a*b + c rounded once (IEEE 754 fusedMultiplyAdd), the same bits in every engine and every build
- `fsum(xs)`: the correctly rounded sum, as Python's math.fsum (Shewchuk's partials)
- `mod(a, b)`: C’s fmod: the result takes the sign of a. The operator `a % b` is the floored remainder, with the sign of b

### Ints

Int is a 64-bit integer. Overflow, `//` and `%` are checked and stop the run with the same words in every engine; these builtins are the unchecked forms an algorithm on 64-bit patterns needs (a random generator, a hash), with the pattern read as unsigned where the name says so.

- `int_zeros(n)`: an [Int] array of n zeros
- `wrap_add(a, b) · wrap_sub(a, b) · wrap_mul(a, b)`: the result mod 2^64, as a pattern: no overflow check
- `umulhi(a, b)`: the high 64 bits of the 128-bit unsigned product of the two patterns
- `ult(a, b)`: a < b with both read as unsigned
- `ustr(a)`: the pattern as an unsigned decimal string

### Binary data

A file's bytes are an [Int], one element per byte, 0 to 255. The readers take that array and an offset and decode the bytes there, little-endian, as Python's struct with `<` does; an offset past the end, or an element that is not a byte, stops the run. In the browser `read_bytes` reads the files the page attaches to the run.

- `read_bytes(path) · read_bytes(path, offset, count)`: a file's bytes as an [Int]; with offset and count, only those (count left out: to the end)
- `le_u8(b, i) · le_i8(b, i) · le_u16(b, i) · le_i16(b, i) · le_u32(b, i) · le_i32(b, i)`: the Int that the 1, 2 or 4 bytes at offset i encode, unsigned or two's complement
- `le_u64(b, i) · le_i64(b, i)`: the Int that the 8 bytes at i encode; le_u64 stops when the value is above 2^63 - 1
- `le_f16(b, i) · le_bf16(b, i) · le_f32(b, i) · le_f64(b, i)`: the number an IEEE half, a bfloat16, a single or a double at i holds, exactly (subnormals, infinities and NaN included)
- `fround(x)`: x rounded to the nearest binary32 (ties to even), as a number: JavaScript's Math.fround, NumPy's float64(float32(x)); a format that computes in single precision rounds each step with it

### Arrays

- `zeros(d1 … d4)`: allocate, rank 1 to 4
- `len(a)`: the leading dimension
- `numel(a)`: the number of elements, every axis multiplied
- `shape(a)`: the lengths of the axes, as a vector
- `linspace(lo, hi, n)`: n evenly spaced points
- `fill(a, v)`: set every element in place
- `sort(a)`: in place, ascending, stable; -0.0 before 0.0, NaN last; no value (copy first for a sorted copy)
- `transpose(A)`
- `sum(a)`: the sum of every element of a real array, in NumPy's order in every exact engine: the entries in row-major order, from 0.0, in blocks of 8192, each block by NumPy's pairwise sum (below 8 entries left to right; to 128, eight running sums combined as a tree; above, two halves, the first a multiple of 8 long). Fast mode is not bound by this order
- `histogram(a, lo, hi, nbins)`: counts of a's entries in nbins equal bins over [lo, hi], the last bin closed on the right, entries outside not counted. Draw it with bars()
- `roll(a, k) · roll(A, k, axis)`: circular shift: element i moves to i + k, wrapping. A periodic stencil in one line: roll(u, 1) - 2.0 * u + roll(u, -1)
- `tridiag_solve(sub, diag, super, rhs)`: the Thomas algorithm: solves a tridiagonal system in O(n), for the diagonally dominant systems that implicit PDE schemes produce. sub[0] and super[n-1] are unused
- `max(a, s) · max(a, b)`: elementwise, against a scalar (either side) or another array; `max(rho, 1.0e-30)` floors a density
- `min(a, s) · min(a, b)`: elementwise, as max
- `max(a) · min(a)`: with ONE array: its largest / smallest entry. `max(abs(u - exact))` is the max-norm error of a solve
- `max(xs, key=f) · min(xs, key=f)`: the first element whose key is largest (smallest), as Python's; xs an array, a list of strings, a string or rows(M); f a def, a builtin or fn(x) -> ...
- `sorted(xs, key=f, reverse=r)`: a new array or list in the order of the keys (numbers or strings), stable, as Python's; a NaN key last
- `argsort(xs, key=f, reverse=r)`: the positions of xs in the order sorted(xs, key=f) puts them
- `all(m)`: true when every entry of a real array is nonzero, as in the mask `a == b` gives; true for an empty array
- `any(m)`: true when some entry is nonzero; false for an empty array
- `where(m, a, b)`: a where m is nonzero, b where it is 0. a and b are each an array or a number, and the three broadcast to one shape, as NumPy's `where`: `where(x > 0.0, x, 0.0)`
- `pow(a, s)`: every entry of a real array to the same scalar power: `pow(rho, 1.0 / 3.0)`
- `slab(A, k)`: a COPY of the k-th slab along the leading axis: row k of a matrix, plane k of a rank-3 array. Real or complex
- `setslab(A, k, S)`: overwrite that slab with S, in place
- `take(A, idx)`: gather: A viewed flat (row-major), idx a real array of flat indices; the result has idx's shape
- `put(A, idx, V)`: scatter, in place: the inverse of take
- `a[lo:hi:s] · A[i, :] · A[:, j]`: a COPY of a block: lo to hi - 1 in steps of s, any part left out; a plain index drops its axis. Real or complex. See Slices above
- `a[lo:hi] = v · A[:, j] = v`: write a number, or an array of the slice's shape, into the block, in place
- `copy(a)`: a new array with a's shape and elements (assignment aliases; this does not)
- `concat(a, b, ...)`: vectors joined end to end
- `vcat(A, B, ...)`: rows stacked; a vector is one row
- `hcat(A, B, ...)`: matrices side by side; a vector is one column
- `[]`: an empty vector, the start of an array grown by append
- `a = append(a, x)`: a with x on the end: a number onto a vector, or a vector as a new last row of a matrix ([] takes a first row). Written rebinding a; amortized O(1). See Growable arrays above
- `a = extend(a, b)`: vector a followed by vector b, grown like an append

### Shapes and constructors

NumPy's names and results (interp/arrays.rs; each case checked against NumPy in tools/arrays/numpy_cases.py). Every result is a new array. An axis may be negative, as in NumPy: -1 is the last.

- `reshape(a, d1 … d4)`: the same entries in row-major order under a new shape; one dimension may be -1 and is worked out
- `ravel(a) · flatten(a)`: every entry as a vector, row-major
- `transpose(a) · transpose(a, axes=axes)`: the axes reversed (a vector is copied); with axes, axis d of the result is axis axes[d] of a
- `swapaxes(a, i, j)`: two axes exchanged
- `eye(n) · eye(n, m) · eye(n, m, k)`: ones on diagonal k (k > 0 above the main one)
- `diag(v) · diag(v, k=k) · diag(A) · diag(A, k=k)`: a vector onto diagonal k of a square matrix, or diagonal k of a matrix as a vector
- `ones(d1 … d4) · zeros_like(a) · ones_like(a)`: filled arrays
- `arange(stop) · arange(start, stop) · arange(start, stop, step)`: NumPy's arange: ceil((stop - start) / step) entries, entry i start + i * ((start + step) - start)
- `linspace(lo, hi, n, endpoint=endpoint)`: NumPy's linspace: i * step + lo, with the end point included or left out; n may be 0 or 1
- `rand(d1 … d4) · randn() · randn(d1 … d4)`: draws from rand()'s generator in row-major order; randn is Box-Muller from two draws, lib/random's formula

### Reductions, products and broadcasting

Arithmetic and comparisons between two arrays (not complex ones) broadcast as NumPy's do: the last axes line up, and each pair of lengths is equal or one of them is 1, so `A - mean(A, axis=0)` subtracts the column means and `A * transpose([w])` scales rows. `sum`, `mean`, `var`, `std` along an axis and `trace` add in NumPy's order (pairwise, in blocks of 8192), so they agree with NumPy to the bit; so does `sum(a)` of a whole array. `keepdims=true` keeps the reduced axis with length 1 (NumPy's keepdims).

- `sum(a, ordered=true)`: the sum of every element left to right from 0.0 (row-major), the order of a Fortran or C loop, in every engine; for a port that must match a reference's last digits. `ordered=false` is `sum(a)`. `dot` and `prod` always go left to right and take `ordered=true` to say so
- `sum(a, axis=axis) · sum(a, axis=axis, keepdims=keepdims)`: along one axis, in NumPy's order
- `mean(a) · mean(a, axis=axis) · mean(a, axis=axis, keepdims=keepdims)`: NumPy's pairwise sum divided by the count
- `var(a) · std(a) · var(a, axis=axis) · std(a, axis=axis)`: population variance (divides by n) and its square root, NumPy's order
- `prod(a) · prod(a, axis=axis)`: left to right
- `amax(a) · amin(a) · amax(a, axis=axis) · amin(a, axis=axis)`: NumPy's max and min: a NaN wins, -0.0 is below 0.0 (max(a) of one array skips NaN, as C's fmax)
- `argmax(a) · argmin(a) · argmax(a, axis=axis) · argmin(a, axis=axis)`: the index of the first largest (smallest) entry, flat or along an axis; the first NaN if there is one
- `any(a, axis=axis) · all(a, axis=axis)`: masks of 1.0 and 0.0 along an axis
- `count_nonzero(a)`: how many entries are not 0
- `outer(u, v)`: u[i] * v[j], each argument flattened
- `dot(u, v) · dot(u, v, ordered=true)`: two vectors of one length, summed left to right (A @ B is the matrix product); `ordered=true` states that order and changes nothing
- `norm(a)`: sqrt(dot(a, a)) over every entry (the Frobenius norm of a matrix)
- `trace(A) · trace(A, offset=k)`: the sum of diagonal k, NumPy's order
- `triu(A) · triu(A, k=k) · tril(A) · tril(A, k=k)`: the entries on and above (below) diagonal k; the rest 0
- `einsum("ij,jk->ik", A, B)`: a contraction of one to three arrays written in NumPy's subscripts (a string literal; an implicit output is the letters written once); a number when the output has no index. Summed left to right over the indices in their order of first appearance, so it agrees with NumPy to rounding

### Gathers, masks, sorting and scans

A mask is an array of 1.0 and 0.0, as a comparison gives, so an array cannot say whether it is a mask or a list of indices: selecting by a mask is `extract`, and by indices `take`.

- `take(A, idx, axis=axis)`: NumPy's take along an axis: the result has A's axes before it, idx's shape, then A's axes after it. take(A, idx) with two arguments is flat
- `extract(mask, a)`: the entries of a where mask is nonzero, as a vector (NumPy's a[mask])
- `putmask(a, mask, v)`: a[mask] = v in place: v a number, or an array of a's shape whose entries go where mask is nonzero
- `flatnonzero(a)`: the flat indices of the nonzero entries
- `argsort(a)`: indices that sort each row (the last axis), stable, NaN last, -0.0 equal to 0.0: NumPy's argsort(kind="stable")
- `cumsum(a) · cumsum(a, axis=axis) · cumprod(a) · cumprod(a, axis=axis)`: running sums and products, left to right; without an axis over every entry, flat
- `diff(a) · diff(a, n=n) · diff(a, n=n, axis=axis)`: a[i + 1] - a[i] along the last axis (or axis), n times
- `clip(x, lo, hi)`: a number or an array held to [lo, hi]; a NaN entry stays NaN and a NaN bound is no bound, as NumPy's clip

### Strings

Python's semantics, by code point; `a.f(b)` is `f(a, b)`. See Strings above.

- `len(s)`: code points (a list: its strings)
- `str(x)`: a number as Python's repr, a boolean as true/false
- `repr(x)`: a string in quotes, as Python writes it
- `float(s)`: CPython's float(str), correctly rounded
- `int(s[, base])`: CPython's int(str, base), to the nearest double
- `chr(n)`: the character with code point n
- `ord(c)`: the code point of one character
- `upper(s)`: full Unicode mapping (ß is SS)
- `lower(s)`: full mapping, final sigma
- `strip(s[, chars])`: whitespace (or chars) off both ends; lstrip, rstrip one end
- `lstrip(s[, chars])`
- `rstrip(s[, chars])`
- `split(s[, sep[, maxsplit]])`: a list of strings; no sep: runs of whitespace
- `rsplit(s[, sep[, maxsplit]])`: from the right
- `splitlines(s)`: Python's line boundaries
- `join(sep, words)`: Python's sep.join(words)
- `find(s, sub[, start[, end]])`: first position, or -1; rfind the last
- `rfind(s, sub[, start[, end]])`
- `index(s, sub)`: find, but a miss stops; also a list's index
- `rindex(s, sub)`
- `count(s, sub)`: non-overlapping occurrences
- `replace(s, old, new[, count])`
- `startswith(s, p)`
- `endswith(s, p)`
- `contains(s, sub)`: the `in` operator: sub in s
- `ljust(s, width[, fill])`: rjust and center likewise
- `rjust(s, width[, fill])`
- `center(s, width[, fill])`
- `zfill(s, width)`: zeros after a sign
- `isdigit(s)`: also isspace isalpha isalnum isdecimal isnumeric isupper islower isascii isprintable
- `isalpha(s) · isalnum(s) · isdecimal(s) · isnumeric(s)`: Python's character classes, by CPython 3.11's Unicode tables: true when s is not empty and every character is in the class
- `isspace(s) · isupper(s) · islower(s) · isascii(s) · isprintable(s)`: Python's, character for character (isascii and isprintable are true for an empty string)
- `unique(s) · unique(words) · unique(xs)`: each character of a string, or each string of a list, once, in code point order: Python's sorted(set(s)); of numbers, the distinct values ascending, as np.unique
- `sorted(words) · sorted(xs)`: a new list, by code point; of numbers, a new array, ascending and stable, NaN last
- `reversed(x)`: x[::-1]: an array, a string or a list of strings, back to front
- `pop(xs) · pop(xs, i)`: a statement: `y = pop(xs, i)` takes element i (the last when left out) out of the list xs, as Python's
- `read_text(path) · read_lines(path) · write_text(path, s)`: a UTF-8 file; read_lines(path) its lines; write_text(path, s)

### Databases (SQLite)

SQLite (third_party/sqlite, public domain) is built into viv and compiled beside a C program that uses these (`viv export --with-sqlite`); the browser does not have them yet. A database is a number, its handle. Parameters are `?` in the SQL, bound from a list of strings (as text) or an array of numbers (NaN as NULL), never pasted into it. `db_exec` runs inside a savepoint: SQL that fails changes nothing. A query is one statement giving one column. NULL reads as NaN in `db_col` and is an error in `db_col_int` and `db_col_str` (select COALESCE(x, k)). SQLite's errors are VivError; a misuse (a bad handle, a wrong count of parameters, a NULL in an Int column) is ValueError.

- `db_open(path) · db_close(db)`: open or create the database file; its handle. db_close closes it
- `db_exec(db, sql) · db_exec(db, sql, params)`: run SQL (several statements, or one with `?` parameters) inside a savepoint; the rows it changed
- `db_col(db, sql, params) · db_col_int(db, sql, params) · db_col_str(db, sql, params)`: the query's one column as an array of numbers (NULL is NaN), an [Int], or a list of strings; params may be left out
- `db_stage(db, name, column) · db_insert(db, table)`: set a column ([Num], [Int] or [Str]) aside, then write every staged column as rows in one transaction (NaN as NULL); the rows written

### Maps

A map holds numbers or strings under number keys or string keys (one kind of each per map), in the order the keys were first stored. It is one handle: `b = m`, a def's parameter and a record field all change the same map. See Maps above.

- `map()`: a new empty map
- `m[k] = v · m[k]`: store v under k (a key already there keeps its place), or read it; a key that is not there stops the run
- `m["key"] = v`: a string key: printed in its repr, `{'the': 2}`; keys(m) is then a list of strings
- `has(m, k)`: true when k is a key of m
- `get(m, k, d)`: m[k] when k is there, else d
- `remove(m, k) · remove(m, k, d)`: take k out and give its value; without d a key that is not there stops the run
- `keys(m) · values(m)`: a new vector of the keys, or the values, in insertion order: loop over a map with `ks = keys(m)` then `for i to len(ks)`
- `copy(m)`: a new map with the same keys, values and order
- `clear(m)`: take every key out
- `len(m)`: the number of keys

### Group-by on Int keys

Written in Viv (lang/viv/std/groupby.viv) and added to a program that calls them, so every engine runs the same code. The keys stay Ints (no float conversion), hashed into an open-addressing table whose slot holds the key beside its group id. A group-by is factorize, then a dense loop: `s[ids[i]] += v[i]`.

- `uniq, ids = factorize(keys)`: the distinct Ints of an [Int] column in first-appearance order, and each row's position in them (`uniq[ids[i]] == keys[i]`), as pandas.factorize(sort=False); any Int is a key
- `ka, kb, ids = factorize2(a, b)`: factorize on two key columns at once: the distinct pairs in first-appearance order, as two columns, and each row's pair id; three columns: factorize2 its ids with the third
- `match_rows(build, probe)`: a hash join's probe: for each Int of probe, the first row of build with that key, or -1 (pandas.Index.get_indexer)
- `lo, hi = key_range(keys)`: the smallest and largest Int of a column, in one pass (0, 0 for no keys)
- `uniq, s, c = group_sum(keys, v)`: per-group sum and row count of an [Int] column in one pass, the groups in first-appearance order as factorize's; keys spanning few values (found inside the call) are slotted by key - lo, others hashed with the sums beside the table; group_sum_num sums a number column
- `ka, kb, s, c = group_sum2(a, b, v)`: group_sum keyed by two [Int] columns: one pass over pair slots when both span few values, else factorize2 and a loop

### Linear algebra

In debug mode each solve checks its own result (a residual, a pivot check or a trace invariant) and stops with an error instead of returning a wrong answer.

- `eig_sym(A, vals, vecs)`: symmetric spectrum by cyclic Jacobi; values ascending, vectors in columns
- `spectrum_sym(A, vals)`: the same eigenvalues eig_sym computes; in debug mode it runs every available method and asserts they agree
- `lin_solve(A, b, x)`: LU with partial pivoting; checks the residual in debug mode
- `lin_solve_multi(A, B, X)`: factor once, many right-hand sides
- `cholesky(A, L)`: factors a symmetric positive-definite matrix; checks every pivot
- `x = lin_solve(A, b) · X = lin_solve(A, B) · inv(A)`: the same LU as lin_solve and lin_solve_multi, into a new array: the same bits
- `L = cholesky(A)`: the same factorization, into a new array
- `vals, vecs = eig_sym(A) · w, V = eig_herm(A)`: the eigensolvers, unpacked into new arrays (written as the call that fills them)
- `U, s, Vt = svd(A) · svd(A, U, s, Vt)`: A = U diag(s) Vt with k = min(m, n) singular values, descending: one-sided Jacobi, lib/svd.viv's algorithm. The signs of the vectors can differ from LAPACK's
- `eig_general(A, re, im)`: balancing + Hessenberg + Francis double-shift QR; complex pairs split across re/im
- `eig_near(A, sigma, k, re, im)`: k eigenvalues near a real shift by shift-invert Arnoldi — one LU instead of a full QR

### Higher-order and contracts

A contract is named for what it computes, not for the method. In debug mode it runs every method that applies and asserts they agree; `viv run --release` runs the auto-selected one. The methods are ordinary Viv code, in lang/viv/std/contracts.viv.

- `simpson(f, a, b, n)`: composite Simpson
- `trapezoid(f, a, b, n)`
- `deriv(f, x0)`: the exact derivative, computed by forward-mode autodiff when the program is compiled
- `apply(f, args…)`: call a function held in a variable or an array
- `release_mode()`: true under --release, else false, so a program can choose what to run in each mode
- `clock()`: wall time in seconds from a monotonic clock; only the difference of two readings means anything, and only a program with a `#! timing` line may call it (every engine refuses it elsewhere; three-engines lists a timing program and does not compare it)

### Derivatives and records

grad and value_and_grad are reverse mode over a whole def (lang/viv/src/grad.rs): the compiler writes the gradient def before any engine runs, so the three engines run the same generated code.

- `grad(f, x, args…)`: the gradient of the number def f returns with respect to its first argument x (a number, an array, or a record of arrays), with the shape of x; the other arguments are constants
- `l, g = value_and_grad(f, x, args…)`: the value and the gradient in one forward and one backward pass, unpacked on a line of its own
- `for_leaves(f, r1, r2, …, x)`: call def f on each array of one or more records of the same type at once (the records' arrays first, then x), in field order; f may write the arrays in place: an optimizer step over a model

### Complex scalars

The `i` suffix works only on a number, so a loop variable named `i` is still a variable. For arrays of them, see the next group.

- `complex(re, im)`
- `re(z)`: real part
- `im(z)`: imaginary part
- `conj(z)`

### Complex arrays and the Fourier transform

A complex array is packed — interleaved (re, im) doubles, which is exactly how C99 lays out `double complex`, so both engines walk the same memory and agree to the byte. `fft` is radix-2 at a power of two, mixed radix (with dedicated radix-3 and radix-5 butterflies) when the primes of the length are all ≤ 7, and Bluestein everywhere else, so no length is a special case and 2·3·5-smooth ones are fast. The twiddle factors and factorization for each length are computed once and reused. Whole-array arithmetic works on complex arrays: `a + b`, `a * b`, a real array or a real or complex scalar on either side of `*`, `/` by a scalar, and `@`. Convention: `fft` is unscaled with exp(−2πikn/N); `ifft` carries the 1/N — the same as numpy and FFTW.

- `czeros(d1 … d4)`: a complex array of zeros, rank 1 to 4
- `fft(a)`: forward transform, in place, any length
- `ifft(a)`: inverse transform, in place; ifft(fft(x)) returns x to rounding error
- `re(a)`: real parts, as a real array of the same shape
- `im(a)`: imaginary parts, as a real array
- `abs(a)`: magnitudes — what you plot as a spectrum
- `conj(a)`: the conjugate, as another complex array
- `complex(re, im)`: a complex array from two real ones of the same shape
- `fft3(a)`: forward transform of a rank-3 complex array along all three axes, in place; all-zero lines are skipped
- `ifft3(a)`: its inverse, carrying 1/N
- `eig_herm(A, w, V)`: the Hermitian eigenproblem: ascending eigenvalues into w, eigenvectors into the COLUMNS of V. Householder to real tridiagonal, then implicit QL
- `transpose(A)`: plain transpose, NO conjugation: the adjoint is conj(transpose(A)), spelled out where it is meant
- `slab(A, k) · setslab · take · put`: slab, setslab, take and put work on complex arrays too; source and target must both be real or both complex

### Randomness

Park–Miller MINSTD in pure f64, so the stream is bit-identical between the interpreter and compiled C — a seeded Monte Carlo gives the same answer either way. `#! use lib/random` adds randn, randn_matrix, randint, shuffle, randperm, choice, arange and pi, written in Viv over rand(), each documented as the formula it computes so Python can replay it (tools/random_replay.py); `#! use lib/stats` adds mean and std.

- `srand(seed)`: seed the stream
- `rand()`: uniform on [0, 1)

### Figures and geometry

Each of these prints a structured record on stdout, and the page draws the chart from it. Compiled C prints the same record, so a figure is part of the program's output, not a call into a plotting library.

- `plot(x, y, "label")`: a line series; consecutive calls sharing x overlay
- `scatter(x, y, "label")`: points
- `bars(x, y, "label")`: a bar chart: a histogram's counts at its bin centers. A plot() on the same x afterwards draws a curve over the bars
- `plot2d(Z, "label")`: a rank-2 array as a viridis heatmap with a colorbar
- `plot3d(x, y, z, "label")`: a curve in 3-space, rotatable
- `quiver(x, y, u, v, "label")`: a vector field
- `scene("title")`: open a 3D scene; wire/verts output belongs to it
- `wire(x, y, z, "group")`: a 3D polyline, optionally colormapped by a 4th array
- `verts(x, y, z, "group")`: point markers — joints, atoms, nodes
- `figure(title, xlabel, ylabel[, xlo, xhi, ylo, yhi[, style]])`: REQUIRED before the first plot. The program sets the axis labels and limits; the page adds none of its own. On a plot2d the window is the picture's physical extent and gives it real axes; style is "contour" (overlay level lines) or "pixels" (no smoothing: the cells are categories or counts)
- `xtick("text", x)`: inside a figure: a tick labeled with the text at x (a high-symmetry point, a category). A chart with any xtick shows only those x ticks.
- `xscale("log") · yscale("log")`: a logarithmic x or y axis for the open figure: plot the values, not their logs. Ticks at powers of ten; a value of zero or less is a gap, and the chart says how many it left out. "linear" is the default
- `aspect("equal")`: one unit of x as long on the screen as one unit of y, for an orbit or a shape in a plane: the plot takes the height that makes a circle round, and widens a range only if that height would be too tall or too flat. "auto" is the default
- `fill_between("label", x, lo, hi)`: a band shaded from lo to hi over x: a confidence band, a gap, a tolerance. Grouped like a plot; takes its label's color, so a band labeled like its curve is drawn in the curve's color
- `rangebars("label", x, lo, hi)`: a capped bar from lo to hi at each x: error bars, spreads, ranges. It takes the ends, not a center and an error, so mean ± err is rangebars("m", x, m - err, m + err) and an asymmetric interval is no harder; the center is a plot or scatter of its own
- `grid("title", rows, cols) · share("x")`: the next rows * cols figures become the panels of one picture, row by row, each drawn as it would be alone. share("x") gives each column one x range and labels it under the lowest panel; share("y") does the same along each row; share("xy") both
- `yright("series label", "axis label")`: the series with that label read against a second, linear y axis on the right, with its own ticks and label: two quantities in different units over one x. The legend marks the entry "(right axis)"
- `style("label", "dash s2 ring thick")`: how the open figure draws the series with that label: a line (solid, dash, dot, noline), a palette slot (s1 to s6), a marker (circle, ring, square, cross, nomarker) and a width (thin, thick), at most one of each. Any other word is a compile error. With no style, series take the palette in order, solid
- `ytick("text", y)`: inside a figure: a y tick labeled with the text; a chart with any shows only those y ticks
- `origin("lower") · colormap("diverging")`: heatmap settings for the open figure: row 0 at the bottom instead of the top; a blue-white-red map centered on zero for a signed field
- `clim(lo, hi) · levels(values)`: heatmap color limits (beyond them, the end color), and contour lines at exactly these values
- `box("id", "style", x0, y0, x1, y1)`: in a figure with the style "diagram": a box by two corners. Style words, at most one of each: a color (s1 to s6, ink, muted), a fill (plain, tint, clear), a line (solid, dash, dot, noline), a width (thin, thick), corners (round, square)
- `arrow("id", "style", xs, ys)`: in a diagram: a line through the points with a head at the last. Words: a color, a line (solid, dash, dot), a width, head or nohead, straight or curved
- `note("id", "text", "style", x, y)`: in a diagram: text at a point. Words: a color, left, center or right, top, middle or bottom (which part of the text is at the point), small, normal or large (11, 13, 15 units), bold, mono
- `checkpoint()`: save the figures so far as one frame; a sequence of frames plays as an animation

### Output and testing

A failed assert prints what failed (for assert_near, the two values, their difference and the tolerance), keeps running, and makes the process exit nonzero at the end, so a build can stop on it while you still see every failure in one run.

- `print(x)`: the bare value
- `print("text", x, ..., sep=" ")`: text and values on one line, a space (or sep) between them, as Python's print
- `fmt(x)`: a print argument: the shortest digits that read back as x, as Python's repr
- `fmt(x, 17)`: a print argument: 1 to 17 significant digits, as %.17g
- `fmt(x, ".3f")`: a print argument: a Python format, [[fill]align][sign][0][width][,][.precision] then e, f, g, d, % (or E, F, G)
- `narrate "text"`: prints a progress line while the program runs, separate from data output
- `text(codes)`: writes an array of Unicode code points as one line of text (the counterpart of chars)
- `assert_near(a, b, tol)`
- `assert_true(cond)`
- `assert_gt(a, b)`
- `assert_lt(a, b)`
- `assert_allclose(a, b, tol)`: elementwise, over arrays

Left out here: the bits and animation builtins (in https://viv-lang.org/llms-full.txt).

## 5. Examples, by topic

### Numbers & arithmetic

Numbers are 64-bit floats, Int is a 64-bit integer, and Bool holds truth values.

**Integer literals are floats.** 10 is read as 10.0, so division never truncates.

```viv
n = 10
show n / 4.0
```

Prints:
```text
2.5
```

**print vs show.** print gives the bare value. show labels a variable with its name, and prints any other expression bare, as print would.

```viv
x = 2.0
print(x)
show x
show x + 1.0
```

Prints:
```text
2
x: 2
3
```

**Scientific notation.**

```viv
tiny = 2.5e-3
show tiny * 4.0e3
```

Prints:
```text
10
```

**Compound assignment.** x += 1 means x = x + 1, and the same holds for -=, *= and /=.

```viv
x = 9.0
x += 1.0
x -= 2.0
x *= 3.0
x /= 4.0
show x
```

Prints:
```text
x: 6
```

**The math builtins.** 29 functions. pow, atan2, hypot, copysign and mod take two arguments, fma takes three, min and max take two or more, and fsum gives the correctly rounded sum of an array. mod is C's fmod.

```viv
show pow(2.0, 10.0)
show atan2(1.0, 1.0)
show mod(7.5, 2.0)
show round(2.5)
```

Prints:
```text
1024
0.785398
1.5
3
```

**erf and gamma.** The engines agree to double precision; gamma(5) is 4! and gamma(1/2) is √π.

```viv
show erf(1.0)
show gamma(5.0)
show gamma(0.5) * gamma(0.5)
```

Prints:
```text
0.842701
24
3.14159
```

**Precedence.** Multiplication binds tighter than addition; parentheses override.

```viv
show 2.0 + 3.0 * 4.0
show (2.0 + 3.0) * 4.0
```

Prints:
```text
14
20
```

### Booleans

true and false are their own type: what a comparison yields, what if and while take, and a value you can store, pass and return.

**A comparison is a Bool.** It prints as true or false, in every engine. It is not the number 1: there is no arithmetic on it.

```viv
show 2.0 > 1.0
show 2.0 == 3.0
```

Prints:
```text
true
false
```

**Store one in a variable.** A flag is a Bool: a comparison gives one, and `and`, `or` and `not` combine them.

```viv
found = false
n = 7.0
small = n < 10.0
show small
both = small and not found
show both
```

Prints:
```text
small: true
both: true
```

**while not done.** not flips a boolean, so while not done keeps looping until done becomes true.

```viv
done = false
count = 0.0
while not done {
    count = count + 1.0
    done = count >= 4.0
}
show count
```

Prints:
```text
count: 4
```

### Control flow

if/else, while, the counted for loop, and the loop-control statements.

**if / else if / else.**

```viv
x = 2.0
if x < 1.0 {
    show 1.0
} else if x < 3.0 {
    show 2.0
} else {
    show 3.0
}
```

Prints:
```text
2
```

**The counted loop.** for i to n runs i from 0 to n−1.

```viv
total = 0.0
for i to 5 {
    total += i
}
show total
```

Prints:
```text
total: 10
```

**while and break.**

```viv
p = 1.0
while 1.0 < 2.0 {
    p = p * 2.0
    if p > 100.0 { break }
}
show p
```

Prints:
```text
p: 128
```

**continue.** Sum only the loop indices above 2.

```viv
s = 0.0
for i to 6 {
    if i < 2.5 { continue }
    s += i
}
show s
```

Prints:
```text
s: 12
```

**Conditions compose with and / or / not.**

```viv
x = 5.0
if x > 0.0 and not (x > 10.0) {
    show x
}
```

Prints:
```text
x: 5
```

**Nested loops.** break leaves the innermost loop only.

```viv
hits = 0.0
for i to 3 {
    for j to 5 {
        if j > 1.5 { break }
        hits += 1.0
    }
}
show hits
```

Prints:
```text
hits: 6
```

**Computed and fractional bounds.** The bound is any expression; i counts 0, 1, 2, … while i < bound, so a bound of 2.5 runs three times.

```viv
n = 2.0 * 3.0
c1 = 0.0
for i to n { c1 += 1.0 }
show c1
c2 = 0.0
for k to 2.5 { c2 += 1.0 }
show c2
```

Prints:
```text
c1: 6
c2: 3
```

### Functions

def, return, docstrings, recursion, and optional type annotations.

**def and return.**

```viv
def hyp(a, b) {
    return sqrt(a * a + b * b)
}
show hyp(3.0, 4.0)
```

Prints:
```text
5
```

**Docstrings.** A bare string as the first statement documents the function and compiles to nothing.

```viv
def double(x) {
    "twice the input"
    return 2.0 * x
}
show double(21.0)
```

Prints:
```text
42
```

**Recursion.** A double holds every integer up to 2^53 exactly, so a counter can stop the recursion with equality; a comparison such as `n < 1.5` also stops it when n is not a whole number.

```viv
def fact(n) {
    if n < 1.5 {
        return 1.0
    }
    return n * fact(n - 1.0)
}
show fact(6.0)
```

Prints:
```text
720
```

**Type annotations.** Optional; the compiler checks them before the program runs.

```viv
def scale(v: Arr, s: Num) -> Arr {
    return v * s
}
show scale([1.0, 2.0], 10.0)
```

Prints:
```text
[10, 20]
```

**A def with no return gives nothing.** A def that ends without a `return` is called as a statement. Its value is nothing: the compiler refuses `show f(x)` or `f(x) + 1.0` for it.

```viv
def report(x) {
    show x + 1.0
}
report(5.0)
```

Prints:
```text
6
```

**Early return.**

```viv
def clamp01(x) {
    if x < 0.0 { return 0.0 }
    if x > 1.0 { return 1.0 }
    return x
}
show clamp01(0.0 - 2.0)
show clamp01(0.4)
show clamp01(9.0)
```

Prints:
```text
0
0.4
1
```

**Call before definition.** Top-level defs are collected before anything runs, so source order does not matter.

```viv
show later(4.0)
def later(x) {
    return x * 10.0
}
```

Prints:
```text
40
```

### Arrays

Literals, zeros, indexing, elementwise arithmetic, and passing arrays through functions.

**Elementwise math on an array.** Every 1-argument math builtin maps elementwise over an array. No loop.

```viv
x = linspace(0.0, 3.14159265, 5)
y = sin(x)
show y
```

Prints:
```text
y: [0, 0.707107, 1, 0.707107, 3.58979e-09]
```

**Broadcasting an array and a scalar.** +, -, *, / broadcast between arrays and scalars, in either order.

```viv
x = linspace(0.0, 2.0, 3)
show 2.0 * x + 1.0
show 1.0 - x
show x / 2.0
```

Prints:
```text
[1, 3, 5]
[1, 0, -1]
[0, 0.5, 1]
```

**Resetting an array in place with fill.** fill(a, v) sets every element. Rebinding (a = zeros(n)) makes a new array, and a caller that passed a in never sees the change; fill changes the array in place.

```viv
a = zeros(3)
a[0] = 7.0
fill(a, 0.0)
show a
fill(a, 1.5)
show a
```

Prints:
```text
a: [0, 0, 0]
a: [1.5, 1.5, 1.5]
```

**Literals and indexing.** Comma-separated subscripts, zero-based.

```viv
v = [10.0, 20.0, 30.0]
M = [[1.0, 2.0], [3.0, 4.0]]
show v[1]
show M[1, 0]
```

Prints:
```text
20
3
```

**zeros and indexed stores.** zeros allocates rank 1 through 4.

```viv
Z = zeros(2, 2)
Z[0, 1] = 7.0
show Z
```

Prints:
```text
Z: [[0, 7],
    [0, 0]]
```

**Elementwise arithmetic.** * between arrays is elementwise (the Hadamard product); array times scalar scales.

```viv
a = [1.0, 2.0]
b = [10.0, 20.0]
show a + b
show a * b
show a * 5.0
```

Prints:
```text
[11, 22]
[10, 40]
[5, 10]
```

**len and linspace.** len reads the leading dimension.

```viv
g = linspace(0.0, 1.0, 5)
show g
show len(g)
```

Prints:
```text
g: [0, 0.25, 0.5, 0.75, 1]
5
```

**Functions take and return arrays.** head_tail takes an array and returns a new one.

```viv
def head_tail(v) {
    w = zeros(2)
    w[0] = v[0]
    w[1] = v[len(v) - 1.0]
    return w
}
show head_tail([3.0, 9.0, 27.0])
```

Prints:
```text
[3, 27]
```

**Scalars broadcast.** A scalar combines elementwise with an array from either side, for every arithmetic operator.

```viv
v = [1.0, 2.0, 4.0]
show v + 10.0
show 1.0 + v
show 2.0 / v
show v / 2.0
```

Prints:
```text
[11, 12, 14]
[2, 3, 5]
[2, 1, 0.5]
[0.5, 1, 2]
```

**Matrices combine elementwise too.**

```viv
A = [[1.0, 2.0], [3.0, 4.0]]
B = [[10.0, 20.0], [30.0, 40.0]]
show A + B
show A * B
```

Prints:
```text
[[11, 22],
 [33, 44]]
[[10, 40],
 [90, 160]]
```

**Rank 3 and 4 exist.** Allocate with zeros; index with one subscript per axis. show prints the flat data for rank 3+.

```viv
z = zeros(2, 2, 2)
z[1, 0, 1] = 5.0
show z[1, 0, 1]
```

Prints:
```text
5
```

**linspace endpoints are inclusive.** Both endpoints are in the array, as in NumPy. With n = 1 it is the start alone, and with n = 0 it is empty.

```viv
show linspace(0.0, 10.0, 5)
show linspace(3.0, 9.0, 1)
```

Prints:
```text
[0, 2.5, 5, 7.5, 10]
[3]
```

**dot is the loop you would write.** dot(a, b) is a builtin. The def below is the loop it runs, left to right, and gives the same sum.

```viv
def mydot(a, b) {
    s = 0.0
    for i to len(a) {
        s += a[i] * b[i]
    }
    return s
}
u = [1.0, 2.0, 3.0]
v = [4.0, 5.0, 6.0]
show mydot(u, v)
show dot(u, v)
```

Prints:
```text
32
32
```

**Negative indices count from the end.** As in Python: a[-1] is the last element, a[:-1] all but it, A[-1, :] the last row.

```viv
a = [10.0, 11.0, 12.0]
show a[-1]
show a[:-1]
A = [[1.0, 2.0], [3.0, 4.0]]
show A[-1, :]
```

Prints:
```text
12
[10, 11]
[3, 4]
```

**A comprehension builds an array.** Its for takes any loop head, and its if keeps the elements it holds for.

```viv
xs = [3.0, 1.0, 4.0, 1.0, 5.0]
show [x * x for x in xs if x > 1.0]
show [i * 2.0 for i to 4]
```

Prints:
```text
[9, 16, 25]
[0, 2, 4, 6]
```

**in, and a key for max and sorted.** x in a is NumPy's: some element equals x. A key is a def, a builtin or fn(x) -> ....

```viv
a = [3.0, -7.0, 4.0]
show 4.0 in a
show max(a, key=abs)
show sorted(a, key=fn(v) -> v * v)
```

Prints:
```text
true
-7
[3, 4, -7]
```

**norm is the square root of dot.** norm(v) is a builtin too, the square root of dot(v, v).

```viv
v = [3.0, 4.0]
show norm(v)
show sqrt(dot(v, v))
```

Prints:
```text
5
5
```

### Complex numbers

A complex number type: the `i` suffix on a number makes it imaginary, arithmetic mixes complex and real values, and every engine (the interpreter, C and wasm) computes the complex functions with the same formulas, so they print the same digits.

**Imaginary literals and arithmetic.** `2.5i` is an imaginary literal (a bare identifier `i` is still a variable — only digits take the suffix). +, -, *, / work between complex and real values; reals promote.

```viv
z = 2.0 + 3.0i
w = complex(1.0, -1.0)
show z
show z * w
show z / w
show z + 1.0
```

Prints:
```text
z: 2+3i
5+1i
-0.5+2.5i
3+3i
```

**re, im, conj, abs.** abs of a complex number is its modulus; z times its conjugate is exactly |z| squared.

```viv
z = 3.0 + 4.0i
show re(z)
show im(z)
show conj(z)
show abs(z)
show re(z * conj(z))
```

Prints:
```text
3
4
3-4i
5
25
```

**Euler's identity, computed.** exp, log, sqrt, sin, cos, tan all work on complex arguments. e^{i pi} comes out as -1 with an imaginary part of 1.2e-16, which is sin of the double nearest pi. sqrt takes the half-angle formulas, so the square root of -1 is exactly i.

```viv
pi = 3.141592653589793
show exp(pi * 1.0i)
show abs(exp(1.0i))
show sqrt(-1.0 + 0.0i)
```

Prints:
```text
-1+1.22465e-16i
1
0+1i
```

**A rotation is a complex multiply.** Multiplying by e^{i theta} rotates the plane by theta. A function that takes and returns complex values needs no annotations; the compiler infers the types.

```viv
def rotate(z, theta) {
    return z * exp(complex(0.0, theta))
}
q = rotate(1.0 + 0.0i, 1.5707963267948966)
show im(q)
test "quarter turn sends 1 to i" {
    assert_near(im(q), 1.0, 1.0e-12)
    assert_near(re(q), 0.0, 1.0e-12)
}
```

Prints:
```text
1
PASS quarter turn sends 1 to i
```

### Complex arrays & the FFT

Packed complex arrays, and a Fourier transform that works at every length and gives identical output in every engine.

**A complex array.** czeros allocates; indexing gives a complex number, and storing a real promotes it.

```viv
a = czeros(3)
a[0] = 1.0 + 2.0i
a[1] = 5.0
show a
show len(a)
```

Prints:
```text
a: [1+2i, 5+0i, 0+0i]
3
```

**A literal with a complex in it.** One complex element makes the whole literal a packed complex array; the reals promote.

```viv
z = [1.0 + 2.0i, 3.0, 0.0 - 1.5i]
show z
```

Prints:
```text
z: [1+2i, 3+0i, 0-1.5i]
```

**Real, imaginary, magnitude.** re, im and abs over a complex array give a REAL array of the same shape. abs is what you plot as a spectrum.

```viv
z = [3.0 + 4.0i, 0.0 - 2.0i]
show re(z)
show im(z)
show abs(z)
```

Prints:
```text
[3, 0]
[4, -2]
[5, 2]
```

**A delta transforms to all ones.** A unit impulse at index 0 transforms to all ones, which also confirms that fft does not scale its output.

```viv
d = czeros(8)
d[0] = 1.0
fft(d)
show d
```

Prints:
```text
d: [1+0i, 1+0i, 1+0i, 1+0i, 1+0i, 1+0i, 1+0i, 1+0i]
```

**A pure tone shows up in two bins.** A cosine at frequency 3 on 16 points puts N/2 = 8 into bins 3 and 13. The test checks that each of the other 14 bins is below 1e-12.

```viv
PI = 3.141592653589793
n = 16
t = czeros(n)
for j to n {
    t[j] = cos(2.0 * PI * 3.0 * j / n)
}
fft(t)
m = abs(t)
show m[3]
show m[13]
test "every other bin is rounding noise" {
    for k to n {
        if k != 3.0 and k != 13.0 {
            assert_lt(m[k], 1.0e-12)
        }
    }
}
```

Prints:
```text
8
8
PASS every other bin is rounding noise
```

**ifft inverts fft.** ifft carries the 1/N, so the round trip returns what you started with — the same convention as numpy and FFTW.

```viv
x = [1.0 + 0.5i, 2.0, 0.0 - 1.0i, 4.0]
y = [1.0 + 0.5i, 2.0, 0.0 - 1.0i, 4.0]
fft(y)
ifft(y)
test "the round trip is the identity" {
    for j to 4 {
        assert_lt(abs(y[j] - x[j]), 1.0e-14)
    }
}
```

Prints:
```text
PASS the round trip is the identity
```

**Any length, not just powers of two.** Seven is prime. Lengths that are not a power of two go through Bluestein's algorithm, so they cost N log N too and nothing about the call changes.

```viv
PI = 3.141592653589793
n = 7
x = czeros(n)
for j to n {
    x[j] = cos(2.0 * PI * 2.0 * j / n)
}
fft(x)
m = abs(x)
show m[2]
test "a prime-length tone still lands in its bin" {
    assert_near(m[2], 3.5, 1.0e-12)
    assert_lt(m[1], 1.0e-12)
}
```

Prints:
```text
3.5
PASS a prime-length tone still lands in its bin
```

**A real array is not a complex one.** fft refuses a real array and says what to do instead. The two are different types and never mix silently.

```viv
x = zeros(8)
fft(x)
```

Stops with an error containing: `` `fft` needs a complex array ``

### Linear algebra

The matrix product, transpose, and the symmetric eigensolver.

**The matrix product.** @ always means matmul, never elementwise.

```viv
A = [[1.0, 2.0], [3.0, 4.0]]
B = [[0.0, 1.0], [1.0, 0.0]]
show A @ B
```

Prints:
```text
[[2, 1],
 [4, 3]]
```

**transpose.**

```viv
A = [[1.0, 2.0], [3.0, 4.0]]
show transpose(A)
```

Prints:
```text
[[1, 3],
 [2, 4]]
```

**The symmetric eigensolver.** eig_sym fills its arguments: values ascending, eigenvectors in columns. Trace = sum of eigenvalues.

```viv
A = [[2.0, 1.0], [1.0, 3.0]]
vals = zeros(2)
vecs = zeros(2, 2)
eig_sym(A, vals, vecs)
show vals
show vals[0] + vals[1]
```

Prints:
```text
vals: [1.38197, 3.61803]
5
```

**Matrix-vector products via matmul.** A 2×2 times a length-2 vector.

```viv
A = [[0.0, 1.0], [1.0, 0.0]]
v = [3.0, 7.0]
show A @ v
```

Prints:
```text
[7, 3]
```

**Chained products.**

```viv
R = [[0.0, 0.0 - 1.0], [1.0, 0.0]]
show R @ R @ R @ R
```

Prints:
```text
[[1, 0],
 [0, 1]]
```

**A 3x3 eigenproblem.** The eigenvalues come out in ascending order and sum to the trace, 9.

```viv
A = [[2.0, 1.0, 0.0], [1.0, 3.0, 1.0], [0.0, 1.0, 4.0]]
vals = zeros(3)
vecs = zeros(3, 3)
eig_sym(A, vals, vecs)
show vals
show vals[0] + vals[1] + vals[2]
```

Prints:
```text
vals: [1.26795, 3, 4.73205]
9
```

### Lambdas & higher-order builtins

fn values, quadrature, and exact differentiation.

**fn is a value.**

```viv
square = fn(x) -> x * x
show square(7.0)
```

Prints:
```text
49
```

**Quadrature.** Integral of x² over [0,1]; exact answer 1/3.

```viv
show simpson(fn(x) -> x * x, 0.0, 1.0, 40)
```

Prints:
```text
0.333333
```

**Exact derivatives.** deriv uses forward-mode autodiff, so the slope is exact to rounding: d/dx[x sin x] at 2 against the product rule.

```viv
show deriv(fn(x) -> x * sin(x), 2.0)
show sin(2.0) + 2.0 * cos(2.0)
```

Prints:
```text
0.0770038
0.0770038
```

**deriv knows log.** d/dx log(x) at 3 is 1/3.

```viv
show deriv(fn(x) -> log(x), 3.0)
```

Prints:
```text
0.333333
```

**Two-argument lambdas.**

```viv
add = fn(a, b) -> a + b
show add(3.0, 4.0)
```

Prints:
```text
7
```

**simpson accepts an odd panel count.** Simpson's rule needs an even number of panels, so simpson rounds an odd n up to the next even number: 7 panels give exactly the 8-panel answer. The integrand is exp, which Simpson's rule does not integrate exactly, so the two would differ otherwise.

```viv
odd = simpson(fn(x) -> exp(x), 0.0, 1.0, 7)
even = simpson(fn(x) -> exp(x), 0.0, 1.0, 8)
show odd == even
show trapezoid(fn(x) -> x * x, 0.0, 1.0, 1000)
```

Prints:
```text
true
0.333333
```

**deriv composes through the chain rule.** d/dx exp(sin x) at 1 = exp(sin 1) · cos 1.

```viv
show deriv(fn(x) -> exp(sin(x)), 1.0)
show exp(sin(1.0)) * cos(1.0)
```

Prints:
```text
1.25338
1.25338
```

**The functions deriv can differentiate.** Besides + - * and /, deriv knows exp, sqrt, pow, sin, cos, log and tan. The last line uses the quotient rule.

```viv
show deriv(fn(x) -> sqrt(x), 4.0)
show deriv(fn(x) -> pow(x, 3.0), 2.0)
show deriv(fn(x) -> tan(x), 1.0)
show deriv(fn(x) -> sin(x) / x, 2.0)
```

Prints:
```text
0.25
12
3.42552
-0.435398
```

**deriv wants a constant exponent.** pow(x, c) differentiates only for literal c.

```viv
show deriv(fn(x) -> pow(x, x), 2.0)
```

Stops with an error containing: `constant exponent`

**Quadrature error shrinks with panels.** For the integral of sin on [0, pi], exactly 2, Simpson's error falls from 0.0046 with 4 panels to 6.5e-8 with 64.

```viv
pi_v = 3.141592653589793
show abs(simpson(fn(x) -> sin(x), 0.0, pi_v, 4) - 2.0)
show abs(simpson(fn(x) -> sin(x), 0.0, pi_v, 64) - 2.0)
```

Prints:
```text
0.00455975
6.453e-08
```

### Scope & mutation

Where names live, what functions can see, and which values are shared.

**Functions read globals.** A top-level scalar is visible inside every def.

```viv
g = 5.0
def f() {
    return g + 1.0
}
show f()
```

Prints:
```text
6
```

**Assigning inside a def never touches a global scalar.** Assignment in a function body creates a local variable; the global g keeps its value, 5.

```viv
g = 5.0
def f() {
    g = 99.0
    return g
}
show f()
show g
```

Prints:
```text
99
g: 5
```

**A def can change a global array.** A def that stores into an element of a global array changes the global array; no return is needed.

```viv
tape = zeros(3)
def store(v) {
    tape[0] = v
    return 0.0
}
x = store(42.0)
show tape
```

Prints:
```text
tape: [42, 0, 0]
```

**Assignment aliases arrays.** b = a makes b refer to the same array; it does not copy the data. Changing b changes a.

```viv
a = [1.0, 2.0]
b = a
b[0] = 99.0
show a
```

Prints:
```text
a: [99, 2]
```

**Array parameters alias too.** A function mutating its array parameter mutates the caller's array.

```viv
def poke(v) {
    v[0] = 77.0
    return 0.0
}
a = [1.0, 2.0]
x = poke(a)
show a
```

Prints:
```text
a: [77, 2]
```

**Closures see the current global.** Capture is by reference for globals: rebinding y changes what f computes.

```viv
y = 2.0
f = fn(x) -> x * y
y = 10.0
show f(3.0)
```

Prints:
```text
30
```

### Prose in programs

Comments, docstrings, narration, and annotation markers.

**Comments are ignored by the compiler.**

```viv
# this line compiles to nothing
show 1.0  # and so does this tail
```

Prints:
```text
1
```

### Plotting

plot prints a structured record, and the page draws the chart from it.

**plot emits a structured record.** figure declares the chart (title, x label, y label) and plot adds a series. Both print records on stdout; the page draws the chart from them and adds nothing of its own.

```viv
x = linspace(0.0, 3.0, 7)
y = zeros(7)
for i to 7 {
    y[i] = x[i] * x[i]
}
figure("A parabola", "x", "x squared")
plot(x, y, "x squared")
```

Prints:
```text
figure: A parabola ; x ; x squared
plot: x squared ; [0, 0.5, 1, 1.5, 2, 2.5, 3] ; [0, 0.25, 1, 2.25, 4, 6.25, 9]
```

**A matrix as a heatmap with plot2d.** plot2d(Z) or plot2d(Z, "label") renders a rank-2 array as a color map (viridis) with a colorbar. The format of the printed record is fixed, and this example's test fails if it changes.

```viv
Z = zeros(2, 3)
for i to 2 {
    for j to 3 {
        Z[i, j] = i * 3.0 + j
    }
}
figure("A 2 x 3 field", "column", "row")
plot2d(Z, "a tiny field")
```

Prints:
```text
figure: A 2 x 3 field ; column ; row
plot2d: a tiny field ; 2 3 ; [0, 1, 2, 3, 4, 5]
```

**A curve in three dimensions with plot3d.** plot3d(x, y, z, "label") draws a rotatable 3-D curve — drag the rendered figure to spin it. Three equal-length rank-1 arrays.

```viv
n = 3
xs = linspace(0.0, 1.0, n)
ys = zeros(n)
zs = zeros(n)
for i to n {
    ys[i] = xs[i] * xs[i]
    zs[i] = 1.0 - xs[i]
}
figure("An arc in space", "x", "y")
plot3d(xs, ys, zs, "a short arc")
```

Prints:
```text
figure: An arc in space ; x ; y
plot3d: a short arc ; [0, 0.5, 1] ; [0, 0.25, 1] ; [1, 0.5, 0]
```

**A histogram.** histogram counts into equal bins; bars draws the counts at the bin centers. A plot on the same x afterwards draws a curve over the bars, which is how a measured distribution is shown against a theoretical one.

```viv
v = [0.1, 0.4, 0.45, 0.8, 0.9, 0.95, 1.7, 2.2]
h = histogram(v, 0.0, 2.0, 4)
c = linspace(0.25, 1.75, 4)
figure("Counts per bin", "value", "count")
bars(c, h, "counts")
plot(c, 4.0 - 2.0 * c, "a line over them")
```

Prints:
```text
figure: Counts per bin ; value ; count
bars: counts ; [0.25, 0.75, 1.25, 1.75] ; [3, 3, 0, 1]
plot: a line over them ; [0.25, 0.75, 1.25, 1.75] ; [3.5, 2.5, 1.5, 0.5]
```

**Same x merges into one chart.** Consecutive plots sharing an x-array become labeled series on one figure.

```viv
n = 41
x = linspace(0.0, 6.28, n)
s1 = zeros(n)
s2 = zeros(n)
for i to n {
    s1[i] = sin(x[i])
    s2[i] = cos(x[i])
}
figure("Sine and cosine", "x", "value")
plot(x, s1, "sin")
plot(x, s2, "cos")
```

Prints:
```text
figure: Sine and cosine ; x ; value
plot: sin ; [0, 0.157, 0.314, 0.471, 0.628, 0.785, 0.942, 1.099, 1.256, 1.413, 1.57, 1.727, 1.884, 2.041, 2.198, 2.355, 2.512, 2.669, 2.826, 2.983, 3.14, 3.297, …
plot: cos ; [0, 0.157, 0.314, 0.471, 0.628, 0.785, 0.942, 1.099, 1.256, 1.413, 1.57, 1.727, 1.884, 2.041, 2.198, 2.355, 2.512, 2.669, 2.826, 2.983, 3.14, 3.297, …
```

**Different x means different charts.** Each chart gets its own figure.

```viv
a_x = linspace(0.0, 1.0, 5)
b_x = linspace(0.0, 2.0, 9)
a_y = zeros(5)
b_y = zeros(9)
for i to 5 { a_y[i] = a_x[i] }
for i to 9 { b_y[i] = b_x[i] * b_x[i] }
figure("A line", "x", "y")
plot(a_x, a_y, "line")
figure("A parabola", "x", "y")
plot(b_x, b_y, "parabola")
```

Prints:
```text
figure: A line ; x ; y
plot: line ; [0, 0.25, 0.5, 0.75, 1] ; [0, 0.25, 0.5, 0.75, 1]
figure: A parabola ; x ; y
plot: parabola ; [0, 0.25, 0.5, 0.75, 1, 1.25, 1.5, 1.75, 2] ; [0, 0.0625, 0.25, 0.5625, 1, 1.5625, 2.25, 3.0625, 4]
```

**figure is required.** A chart's title and axis labels are part of the program. Nothing outside it can add them, so a plot with no figure declared does not compile — in any engine.

```viv
x = linspace(0.0, 1.0, 5)
plot(x, x, "a line")
```

Stops with an error containing: `` needs a `figure(title, xlabel, ylabel)` before it ``

### Testing

test blocks and the assert family: a program checks its own results and prints which tests passed and which failed.

**A passing test.** assert_near(a, b, tol) inside a named test block. What prints is the name — one line per test, however many assertions are inside.

```viv
test "pythagoras" {
    assert_near(sqrt(3.0 * 3.0 + 4.0 * 4.0), 5.0, 1.0e-12)
}
```

Prints:
```text
PASS pythagoras
```

**A failing test says what failed.** The FAIL line names the test and shows the failing assertion's values, so you see which check failed and by how much.

```viv
test "wrong on purpose" {
    assert_near(2.0 + 2.0, 5.0, 1.0e-9)
}
```

Prints:
```text
FAIL wrong on purpose -- assert_near: |4 - 5| = 1 > 1e-09
```

**Test a function against a property.** Check a function against a mathematical fact it must satisfy (here, the partial sums of 1/k! converge to e), not only against fixed values.

```viv
def esum(n) {
    s = 0.0
    for i to n {
        s += 1.0 / gamma(i + 1.0)
    }
    return s
}
test "partial sums of 1/k! approach e" {
    assert_near(esum(20.0), exp(1.0), 1.0e-12)
}
```

Prints:
```text
PASS partial sums of 1/k! approach e
```

### Errors, exactly as the compiler says them

The error messages are part of the language. Each one below is what the compiler prints, and the site's tests fail if a message changes.

**A string is not changed in place.** Strings are immutable, as in Python; build a new one instead.

```viv
s = "hi"
s[0] = "H"
```

Stops with an error containing: `'str' object does not support item assignment`

**deriv cannot differentiate a call to a def.** Write the expression inline in the fn instead of calling f.

```viv
def f(x) {
    return x * x
}
show deriv(fn(a) -> f(a), 1.0)
```

Stops with an error containing: `cannot differentiate through`

**A def does not change a record it is given.** A record is a value: change a copy and return it. An array inside it can still be written in place.

```viv
record P { x: Num }
def f(p: P) {
    p.x = 2.0
    return p.x
}
show f(P { x: 1.0 })
```

Stops with an error containing: `` `p` is a parameter, and a def does not change the values it is given ``

**Unknown function.**

```viv
show frobnicate(1.0)
```

Stops with an error containing: `unknown function`

**Arity is checked.**

```viv
def f(a, b) {
    return a + b
}
show f(1.0)
```

Stops with an error containing: `expects 2 argument`

**Out of bounds says where.**

```viv
a = [1.0, 2.0]
show a[5]
```

Stops with an error containing: `index 5 out of bounds for dimension 0 (size 2)`

**An index before the start stops.** A negative index counts from the end, as in Python: -1 is the last element. Two elements have no element -3.

```viv
a = [1.0, 2.0]
show a[-1]
show a[0.0 - 3.0]
```

Stops with an error containing: `index -3 out of bounds for dimension 0 (size 2)`

**Shape mismatch.**

```viv
show [1.0, 2.0] + [1.0, 2.0, 3.0]
```

Stops with an error containing: `array shape mismatch: [2] vs [3]`

**break outside a loop.**

```viv
break
show 1.0
```

Stops with an error containing: `` `break`/`continue` outside a loop ``

**Undefined variable.**

```viv
show nonexistent + 1.0
```

Stops with an error containing: `undefined variable`

**plot wants arrays.**

```viv
figure("Nothing to see", "x", "y")
plot(1.0, 2.0)
```

Stops with an error containing: `` `plot` takes two arrays ``

**Keywords cannot name functions.** record, test, show and the other keywords are reserved.

```viv
def record(v) {
    return v
}
show record(1.0)
```

Stops with an error containing: `expected function name`

**A number is not a condition.** if and while take a Bool. Write the comparison you mean.

```viv
count = 3.0
if count {
    show count
}
```

Stops with an error containing: `condition must be a boolean`

**No arithmetic on a Bool.** true is not 1. Every engine refuses this.

```viv
x = true + 1.0
```

Stops with an error containing: `boolean`

**An array is not a condition.** A comparison of two arrays is a mask, not a Bool; the note says to write all(a == b).

```viv
a = [1.0, 2.0]
if a == a {
    show 1.0
}
```

Stops with an error containing: `condition must be a boolean, and this is an array`

**% by zero stops the program.** `/` by zero gives inf, as IEEE floats do, but `%` and `//` by zero have no answer and stop.

```viv
show 7.0 % 0.0
```

Stops with an error containing: `` division by zero in `a % b` ``

### Gotchas

Float facts and Viv-specific surprises, each shown by a program that runs.

**Float equality is a trap.** 0.1 + 0.2 is not 0.3 in binary floating point; compare with a tolerance instead. (Integers are another matter: a double holds them exactly up to 2^53, so `n == 1.0` on a counter is an exact test.)

```viv
show (0.1 + 0.2) == 0.3
show abs((0.1 + 0.2) - 0.3) < 1.0e-9
```

Prints:
```text
false
true
```

**The printer rounds to 6 significant digits.** show uses %g formatting; the stored value is full double precision even when the print looks tidy.

```viv
show 0.1 + 0.2
show (0.1 + 0.2) - 0.3
```

Prints:
```text
0.3
5.55112e-17
```

**eig_sym does not check symmetry.** Feed it a non-symmetric matrix and it computes something — for [[1,5],[0,2]] the true eigenvalues are 1 and 2, and this is not them. Symmetrize first.

```viv
A = [[1.0, 5.0], [0.0, 2.0]]
v = zeros(2)
w = zeros(2, 2)
eig_sym(A, v, w)
show v
```

Prints:
```text
v: [-3.52494, 6.52494]
```

**Copying an array means allocating one.** Because assignment aliases, copying the data takes copy(a), which allocates a new array with the same entries.

```viv
a = [1.0, 2.0]
b = copy(a)
b[0] = 99.0
show a
show b
```

Prints:
```text
a: [1, 2]
b: [99, 2]
```

**Adding 1 to 1e16 loses the 1.** In exact arithmetic both lines would print 1; here both print 0. A double carries about 16 significant digits, so the 1 is rounded away when it meets 1e16, whichever way the sum is grouped.

```viv
big = 1.0e16
show (big + 1.0) - big
show big + (1.0 - big)
```

Prints:
```text
0
0
```

**Catastrophic cancellation.** At x = 1e12, sqrt(x+1) - sqrt(x) subtracts two nearly equal numbers and keeps only 5 correct digits; the rewrite 1 / (sqrt(x+1) + sqrt(x)) has no subtraction and gets them all.

```viv
x = 1.0e12
naive = sqrt(x + 1.0) - sqrt(x)
rewritten = 1.0 / (sqrt(x + 1.0) + sqrt(x))
show naive
show rewritten
```

Prints:
```text
naive: 5.00004e-07
rewritten: 5e-07
```

### Stencils and PDE steps

A finite-difference stencil is arithmetic on shifted copies of an array, and an implicit time step is a tridiagonal solve. These are the pieces the PDE solvers are built from.

**A periodic second difference.** roll(u, k) moves element i to i + k, wrapping around, so roll(u, 1)[i] is u[i - 1]. The second difference of i squared is exactly 2; the two ends show the wrap.

```viv
u = [0.0, 1.0, 4.0, 9.0, 16.0]
d2 = roll(u, 1.0) - 2.0 * u + roll(u, -1.0)
show d2
```

Prints:
```text
d2: [17, 2, 2, 2, -23]
```

**roll along one axis of a grid.** A third argument picks the axis: 0 shifts rows, 1 shifts along each row. A five-point Laplacian is four of these.

```viv
A = [[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]]
B = roll(A, 1.0, 1.0)
C = roll(A, 1.0, 0.0)
show B[0, 0]
show C[0, 0]
```

Prints:
```text
3
4
```

**A tridiagonal solve.** tridiag_solve(sub, diag, super, rhs) is the Thomas algorithm, O(n). Here the matrix is the 1-D Laplacian (-1, 2, -1), and the right-hand side is chosen so the answer is all ones. sub[0] and super[n-1] are never read.

```viv
n = 4
sub = zeros(n) - 1.0
diag = zeros(n) + 2.0
sup = zeros(n) - 1.0
x = tridiag_solve(sub, diag, sup, [1.0, 0.0, 0.0, 1.0])
show x
```

Prints:
```text
x: [1, 1, 1, 1]
```

**One implicit heat step.** Backward Euler for u_t = u_xx with dt / h^2 = 1: solve (I - r D2) u_new = u_old. The spike spreads symmetrically, and no time step is too large for it.

```viv
n = 5
r = 1.0
u = [0.0, 0.0, 1.0, 0.0, 0.0]
off = zeros(n) - r
u = tridiag_solve(off, zeros(n) + (1.0 + 2.0 * r), off, u)
show u
```

Prints:
```text
u: [0.0555556, 0.166667, 0.444444, 0.166667, 0.0555556]
```

**The max-norm error.** max and min of ONE array reduce it to its largest or smallest entry. With two arguments they stay elementwise.

```viv
x = linspace(0.0, 1.0, 5)
err = max(abs(x * x - x))
show err
show min(x * x - x)
```

Prints:
```text
err: 0.25
-0.25
```

**An array to a power.** pow(a, s) raises every entry to the same scalar power.

```viv
c = pow([1.0, 8.0, 27.0], 1.0 / 3.0)
show c[1]
show sum(pow([1.0, 2.0, 3.0], 2.0))
```

Prints:
```text
2
14
```

**Negating an array.** Unary minus works on arrays, real or complex.

```viv
a = [1.0, -2.0, 0.5]
b = -a
show b
```

Prints:
```text
b: [-1, 2, -0.5]
```

### Small algorithms

Classic numerical algorithms, each a few lines long.

**Newton's square root.** Repeat x = (x + a/x) / 2; each step roughly doubles the number of correct digits.

```viv
def newton_sqrt(a) {
    x = a
    for i to 20 {
        x = 0.5 * (x + a / x)
    }
    return x
}
show newton_sqrt(2.0)
show sqrt(2.0)
```

Prints:
```text
1.41421
1.41421
```

**Bisection.** Root of cos(x) on [1, 2]. Each step halves the interval, so 40 halvings leave it 2^-40 ≈ 9e-13 wide and give pi/2 to about twelve digits; show prints the first six.

```viv
def bisect_cos(lo, hi) {
    for i to 40 {
        mid = 0.5 * (lo + hi)
        if cos(lo) * cos(mid) < 0.0 {
            hi = mid
        } else {
            lo = mid
        }
    }
    return 0.5 * (lo + hi)
}
show bisect_cos(1.0, 2.0)
```

Prints:
```text
1.5708
```

**gcd via mod.** Euclid's algorithm. The floats hold whole numbers exactly, so mod is exact.

```viv
def gcd(a, b) {
    if b < 0.5 {
        return a
    }
    return gcd(b, mod(a, b))
}
show gcd(1071.0, 462.0)
```

Prints:
```text
21
```

**Horner evaluation.** A polynomial as an array of coefficients, highest power first.

```viv
def horner(coeffs, x) {
    acc = 0.0
    for i to len(coeffs) {
        acc = acc * x + coeffs[i]
    }
    return acc
}
show horner([2.0, 0.0 - 3.0, 5.0], 2.0)
```

Prints:
```text
7
```

**The power method.** The dominant eigenvalue of a symmetric 2x2, found by repeated matrix-vector products; compare eig_sym.

```viv
def normalize(v) {
    n = sqrt(v[0] * v[0] + v[1] * v[1])
    return v / n
}
def rayleigh(A, v) {
    w = A @ v
    return v[0] * w[0] + v[1] * w[1]
}
A = [[2.0, 1.0], [1.0, 3.0]]
v = [1.0, 0.0]
for i to 30 {
    v = normalize(A @ v)
}
show rayleigh(A, v)
```

Prints:
```text
3.61803
```

**Compensated summation.** Kahan's method carries each addition's rounding error into the next one. Adding 0.1 a million times this way gives exactly 100000.

```viv
def kahan_sum(n, val) {
    s = 0.0
    comp = 0.0
    for i to n {
        y = val - comp
        t = s + y
        comp = (t - s) - y
        s = t
    }
    return s
}
show kahan_sum(1000000.0, 0.1) - 100000.0
```

Prints:
```text
0
```

**Fibonacci two ways.** Naive recursion, whose number of calls grows exponentially with n, against a loop that takes n steps.

```viv
def fib_rec(n) {
    if n < 1.5 {
        return n
    }
    return fib_rec(n - 1.0) + fib_rec(n - 2.0)
}
def fib_loop(n) {
    a = 0.0
    b = 1.0
    for i to n {
        t = a + b
        a = b
        b = t
    }
    return a
}
show fib_rec(15.0)
show fib_loop(15.0)
```

Prints:
```text
610
610
```

## 6. Writing a post

A post is one `.viv` file whose first lines are a fenced header: a line that
is exactly `###`, then `key: value` fields, a blank line, the prose notes, and
`###` again. The lexer skips the block. The site build runs every post and
refuses it unless:

- the header has `post:` (title), `date:` (YYYY-MM-DD) and `tags:` (comma separated); `summary:` is at most 600 characters;
- at least one `test "..." { ... }` block runs, and none fails (a `kind: problem` post is the exception: it must have a failing test);
- a demo (the default kind) or an `kind: article` draws a figure: `figure(...)` then `plot(...)` or another series, or `stage(...)` for an animation;
- it runs in the browser in 10 s or less, and the figure records are at most 160 KB;
- a program written by a model says so with `agent: <model name>`. Such a post starts in the AI lane, and only other people's reviews move it up.

Other kinds: `kind: ladder` (built from `rung "title" { ... }` blocks), `kind: problem`
(stubs and tests that fail), `kind: solution` with `solves: <problem slug>` (the
problem's tests, word for word, all passing). `from: <slug>` with `edge:` (fork,
fix, extend, rewrite, adversarial-review, validate) derives a post from another,
and the build checks the claim by running both.

```viv
###
post: The trapezoid rule's error falls as h squared
date: 2026-09-24
tags: quadrature, calculus
summary: Halving the step divides the trapezoid rule's error on sin by four.
agent: claude-opus-5.5

The integral of sin from 0 to pi is exactly 2. The trapezoid rule's error is
proportional to h^2, so doubling the number of panels should cut it by four.
###

def trap(a, b, n) {
    h = (b - a) / n
    s = 0.5 * (sin(a) + sin(b))
    for k to n - 1 {
        s = s + sin(a + (k + 1) * h)
    }
    return s * h
}

ns = [4.0, 8.0, 16.0, 32.0, 64.0]
errs = zeros(5)
for k to 5 {
    errs[k] = abs(trap(0.0, 3.141592653589793, ns[k]) - 2.0)
}
figure("trapezoid error on sin", "panels n", "|error|")
plot(ns, errs, "error")
ratio = errs[3] / errs[4]
show ratio

test "halving h divides the error by four" {
    assert_near(ratio, 4.0, 1.0e-2)
}
```