μ

math in μ

Math in μ

Each card shows a mathematical formula on the left and the equivalent μ expression on the right. μ is evaluated strictly right-to-left with no operator precedence; parentheses make order explicit. Verbs apply element-wise over lists, and adverbs / (fold) and \ (scan) reduce or accumulate. Plotting stories and visual narratives live on the plots page.

Finance

Sharpe Ratio

// returns vector r, risk-free rate rf
r: 0.10 0.12 0.08 0.15 0.09 0.11
rf: 0.02
(mean(r) - rf) % stddev r
// → 3.897

The expression computes the ex-post Sharpe ratio S = (r̄ − rᶠ) / σᵣ, where r̄ is the sample mean of realised excess returns and σᵣ their sample standard deviation. Both mean and stddev are single-pass stdlib functions that reduce the full list in one traversal; the optimizer fuses the stddev body into a Welford-style accumulator when it detects the variance pattern, eliminating all intermediate allocations. The result here is a raw period Sharpe — not annualised — so multiply by √T (where T is the number of periods per year) to obtain the conventional annualised figure; the sqrt verb and scalar multiplication extend this trivially.

Bond Price

// F=1000, c=50, r=0.05, n=10
face:1000.0; c:50.0; r:0.05; n:10
coupons: +/ c * exp(-r * 1.0+!n)
coupons + face * exp(-r*n)
// → 990.245

!n generates the integer sequence {0, 1, …, n−1}; adding the scalar 1.0 shifts it to {1, 2, …, n}, giving the coupon payment times t. Multiplying by −r and applying exp element-wise yields the continuous discount factors e^−rt, which are then scaled by the coupon c and reduced by +/ (fold-sum) to give the present value of the coupon stream. Continuous compounding is exact here; for semi-annual compounding with the standard bond convention, replace exp(-r * t) with (1.0 + r%2) pow -2*t. The same structure generalises directly to amortising bonds or floating-rate notes by replacing the scalar c with a coupon vector.

Black-Scholes Call Price

// Normal CDF — A&S rational approx (max err < 1.5e-7)
erf0:{[t] p:1.0%(1.0+0.3275911*t);
  1.0 - exp(-t*t)*(
    (0.254829592*p) + ((-0.284496736)*p*p)
    + (1.421413741*p*p*p)
    + ((-1.453152027)*p*p*p*p)
    + (1.061405429*p*p*p*p*p))}
sq2: sqrt 2.0
sign:{1-2*(0>x)}
ncdf:{[v] 0.5+0.5*(sign v)*(erf0((abs v)%sq2))}

// Black-Scholes call: S K r sig T
bscall:{[S;K;r;sig;T]
  (S*(ncdf (log(S%K)+(r+0.5*sig*sig)*T)%(sig*sqrt T)))
  -(K*exp(-r*T)*(ncdf ((log(S%K)+(r+0.5*sig*sig)*T)%(sig*sqrt T))-sig*sqrt T))}

bscall 100 100 0.05 0.2 1.0
// → 10.450  (exact: 10.4506)

The Abramowitz & Stegun five-term rational approximation to Φ attains a maximum absolute error below 1.5×10⁻⁷, sufficient for option pricing where model uncertainty dominates. The symmetry identity Φ(−x) = 1 − Φ(x) is handled by computing erf0 on |v| and then correcting the sign, keeping the approximation in the positive half-plane where it is most accurate. The values d₁ and d₂ are the log-moneyness adjusted for drift and discounting respectively; d₂ = d₁ − σ√T is the risk-neutral probability that the option expires in-the-money. To sweep C(K) across a strike grid and visualise the call price curve, see the Black-Scholes strike sweep on the plots page.

Statistics

Z-Score Normalisation

zscore: {(x - mean x) % stddev x}
zscore 2 4 4 4 5 5 7 9
// → -1.5 -0.5 -0.5 -0.5 0.0 0.0 1.0 2.0

The lambda standardises the vector x to zero mean and unit variance by computing x̃ᵢ = (xᵢ − x̄) / σₓ element-wise — the implicit argument x receives the entire list, and scalar broadcasting handles the subtraction and division in one pass each. The optimizer recognises the stddev body shape and rewrites it to a single-pass Welford accumulator, so the denominator σₓ is computed without allocating the centred vector (x − x̄) as an intermediate list. The result lies in the standardised space where ‖x̃‖² = n − 1, which is the natural input for PCA, Mahalanobis distance, and ridge regression; apply this transform to each column of a feature matrix before fitting.

Simple Linear Regression

xs: 1.0 2.0 3.0 4.0 5.0
ys: 2.0 4.0 5.0 4.0 5.0
mx: mean xs
my: mean ys
slope: (+/(xs-mx)*(ys-my)) % +/(xs-mx)*(xs-mx)
intcpt: my - slope*mx
slope    // → 0.6
intcpt   // → 2.2

The OLS slope β̂ = Cov(x, y) / Var(x) is computed as the ratio of two inner products: ⟨x − x̄, y − ȳ⟩ / ⟨x − x̄, x − x̄⟩, where +/u*v is the dot product of the centred vectors. The intercept α̂ = ȳ − β̂ x̄ follows from the normal equations' first-order condition that the regression plane passes through the centroid (x̄, ȳ). This two-expression form is exact for simple linear regression; for multiple regression with a design matrix X, use A %. b (matrix solve) on the normal equations X^TX β = X^Ty, or see the predict and residuals functions in the Stat Arb section for ready-made OLS wrappers.

Exponential Moving Average

// weights w_i = alpha*(1-alpha)^i, newest first
alpha: 0.3
prices: 100.0 102.0 101.0 105.0 103.0 107.0
n: #prices
wts: alpha * exp(log(1.0-alpha) * !n)
(+/ wts * prices) % +/wts
// → 103.65  (exponentially weighted mean)

The EMA is equivalently expressed as a weighted mean with geometrically decaying weights wᵢ = α(1−α)ⁱ, newest observation first. Here exp(log(1−α) · i) evaluates to (1−α)ⁱ exactly in IEEE 754 arithmetic and avoids the numerical drift that accumulates when multiplying (1−α) iteratively in floating-point. Normalising by +/wts (the finite geometric sum) rather than assuming an infinite history corrects for edge effects at the start of the series. The result is the EWMA estimate of the current price level; for volatility estimation, replace prices with squared log-returns and set α to correspond to the desired half-life τ via α = 1 − exp(−ln 2 / τ).

Pure Math

Cosine Similarity

a: 3.0 1.0 4.0
b: 1.0 5.0 9.0
(+/ a*b) % (sqrt(+/a*a) * sqrt(+/b*b))
// → 0.834

The expression evaluates cos θ = (a · b) / (‖a‖ · ‖b‖) directly in terms of fold-sums: +/a*b is the inner product, and sqrt(+/a*a) is the Euclidean norm ‖a‖₂. Parentheses around the fold-sum argument to sqrt are required since math keywords bind to a single atom — sqrt(+/a*a) groups the sum first. Cosine similarity lies in [−1, 1] and is invariant to scaling, making it the canonical distance metric for high-dimensional embedding spaces (e.g. word vectors, factor loadings). For pairwise similarity across a matrix of row vectors, the full cosine matrix is D_ij = (A A^T)_ij / (‖aᵢ‖ · ‖aⱼ‖), computable via *. and element-wise normalisation.

Newton's Method — √2

// 8 iterations from seed x₀ = 1.0
8 {[v] 0.5*(v + 2.0%v)}/ 1.0
// → 1.4142135623731

// or converge automatically (stops when stable)
{[v] 0.5*(v + 2.0%v)}/ 1.0
// → 1.4142135623731

The iteration x_n+1 = ½(x_n + 2/x_n) is Newton's method applied to f(x) = x² − 2 = 0, and converges quadratically: the error εₙ₊₁ ≈ εₙ²/2 so the number of correct decimal digits roughly doubles each step. Starting from x₀ = 1, double precision (≈15.9 significant digits) is reached in 5 iterations; 8 are shown here for clarity. n f/ seed applies f exactly n times and returns the final value, equivalent to the fold fⁿ(seed). Use f/ seed (converge) to iterate until the result stabilises automatically. The same template handles arbitrary fixed-point and contraction-mapping iterations: replace the lambda with any f satisfying |f′(x*)| < 1 at the fixed point x*.

Fibonacci via Pair Recurrence

// fibstep advances pair (a b) → (b a+b)
fibstep: {(1_ x) , +/x}
// fib n: apply n steps, take first element
fib: {*(x fibstep/ 0 1)}
fib 10   // → 55
fib 20   // → 6765

fibstep implements one application of the companion matrix [[0,1],[1,1]] acting on the pair (Fₙ, Fₙ₊₁): 1_ x drops the head (tail), +/x sums the pair to give Fₙ₊₂, and , concatenates them into (Fₙ₊₁, Fₙ₊₂). This state-pair encoding avoids both recursion (stack depth O(n)) and the full memoisation table (O(n) space), maintaining only a constant-size state vector. The monadic * takes the first element of the resulting pair. For large n the companion-matrix approach gives O(log n) time via repeated squaring, but μ's n f/ is O(n); the pair encoding is still preferable to naive recursion for clarity and stack safety. n fibstep/ seed uses the iterate-fold adverb: n is the count, fibstep the step function, and 0 1 the initial pair.

Linear Algebra

Matrix Multiplication — Portfolio Variance

// weights w, covariance matrix Σ
c: (cut 2 (0.04 0.01 0.01 0.09))
w: 0.3 0.7
portvar: +/ w * (c *. w)
// → 0.0519 (5.2% volatility)

The portfolio variance σ²ₚ = w^TΣw is the quadratic form evaluated as two successive linear operations: c *. w computes the matrix-vector product Σw ∈ ℝ², then +/ w * (Σw) takes the inner product with w. cut 2 reshapes a flat list into a 2×2 matrix stored row-major. The result σ²ₚ = 0.0519 gives a portfolio volatility σₚ ≈ 22.8% — noticeably less than the weighted average of individual asset vols (√0.04·0.3 + √0.09·0.7 ≈ 27%), reflecting the diversification benefit from the off-diagonal covariance ρσ₁σ₂ = 0.01. To find the minimum-variance portfolio weights, solve the KKT system Σw = λ1 subject to 1^Tw = 1 using %..

Matrix Solve — Linear System

// Solve 2×2 system: 3x + y = 9, x + 2y = 8
A: (cut 2 (3.0 1.0 1.0 2.0))
b: 9.0 8.0
x: A %. b
// → 2.0 3.0

%. solves Ax = b via Gaussian elimination with partial pivoting, which is numerically preferable to forming A⁻¹ explicitly: the condition number of the solve is κ(A), whereas explicit inversion can square the effective condition number in floating-point arithmetic. For the 2×2 system 3x + y = 9, x + 2y = 8, the exact solution is x = 2, y = 3, recoverable analytically as x = (9·2 − 8·1)/(3·2 − 1·1) = 10/5 by Cramer's rule. The same operator handles larger square systems; for overdetermined systems (m > n) use the normal equations X^TX β = X^Ty, forming X^TX via *. before passing to %..

Correlation Matrix Heatmap

// 3-asset correlation matrix
corr: (cut 3 (1.0 0.8 0.3 0.8 1.0 0.5 0.3 0.5 1.0))
heatmap[corr; `Correlations]

heatmap renders the matrix as a colour grid scaled between its minimum and maximum values, mapping each ρᵢⱼ ∈ [−1, 1] to a colour band. The correlation matrix shown — ρ₁₂ = 0.8, ρ₁₃ = 0.3, ρ₂₃ = 0.5 — is positive definite (all eigenvalues positive), a necessary condition for it to be a valid covariance structure. Visual inspection of the heatmap immediately reveals clustering: assets 1 and 2 are highly correlated, while asset 3 is a partial diversifier. In portfolio construction this informs block-diagonal covariance approximations and hierarchical risk parity; in factor modelling it exposes residual correlation that a linear factor structure has not explained.

linalg module

Covariance & Correlation Matrix — covmat / corrmat

// Two return series (each a column of daily log-returns)
r1: 0.01,-0.02,0.03,0.01,-0.01
r2: 0.02,0.01,0.01,-0.02,0.03
cols: [r1; r2]

S: (covmat[cols])    // 2×2 covariance matrix
C: (corrmat[cols])   // 2×2 correlation matrix (diagonal = 1)

// With real CSV data, extract columns first:
// t: csv "prices.csv"
// cols: each[{logrets x}; t@1]
// S: (covmat[cols])

covmat and corrmat accept a list of equal-length return vectors and produce symmetric N×N matrices by applying cov (or cor) to every (i, j) pair via outer. The outer-product map outer[f; xs; ys] evaluates f(xs@i, ys@j) for all i, j, returning a list-of-lists of shape #xs × #ys — the same primitive that drives both matrix functions. To build these matrices from price data, pipe through logrets first: logrets p computes 1 _ log ratios p, dropping the initial reference level and yielding a vector of length #p − 1. With real CSV data, extract each column, apply logrets via each, then pass the resulting list to covmat.

Markowitz Tangency Portfolio — mkw

// Minimal synthetic example (2 assets)
S: cut 2 (0.04 0.01 0.01 0.09)   // covariance matrix
mu: 0.10 0.15                      // expected returns
rf: 0.02                           // risk-free rate

w: (mkw[S; mu; rf])
// → 0.573 0.427  (weights sum to 1)

// Full pipeline from price matrix:
// prices: [(100 103 101 106 104 108); (50 51 53 52 55 54)]
// cols: each[{logrets x}; prices]
// S: (covmat[cols])
// mu: each[{mean x}; cols]
// w: (mkw[S; mu; 0.02])

mkw solves the KKT system for the maximum-Sharpe (tangency) portfolio via two-fund separation: it computes the raw vector z = Σ⁻¹(μ − rf·1) using %. (Gaussian elimination), then normalises by z's sum so weights sum to 1. No quadratic programming is needed — for a long-only portfolio with a risk-free asset, the tangency weights are the unique solution to this linear system. The result w* is the risky portfolio that maximises the Sharpe ratio S = (μₚ − rf) / σₚ over all possible combinations. To run this daily with fresh price targets, replace mu with analyst 1-year target returns and recompute w each day.

Kelly Fraction & Portfolio Sizing — kelly

// Given Markowitz weights w (from mkw above)
S: cut 2 (0.04 0.01 0.01 0.09)
mu: 0.10 0.15
rf: 0.02
w: (mkw[S; mu; rf])

// Portfolio expected return and variance
mup:  (portret[w; mu])          // → 0.121
sig2: (portvar[w; S])           // → 0.034

// Full-Kelly leverage fraction
f: (kelly[mup; sig2; rf])       // → ~2.9×

// Half-Kelly (more conservative in practice)
0.5 * f

The Kelly criterion f* = (μₚ − rf) / σ²ₚ gives the leverage that maximises the long-run expected log-wealth growth rate E[log(1 + f·(rₚ − rf))]. For a continuously rebalanced portfolio, this is exactly the growth-optimal fraction. In practice, full-Kelly is aggressive — f* > 1 means borrowing — so halving (half-Kelly) is common: it preserves 75% of the growth rate while cutting variance by 75%. Here portret[w; mu] evaluates w·μ (the weighted-average expected return) and portvar[w; S] evaluates w^TΣw (the quadratic form), both as O(N) and O(N²) operations respectively. Together with mkw, this pipeline — logrets → covmat → mkw → kelly — is the complete daily portfolio construction loop.

Statistical Arbitrage

Pairs Trading — Cointegration Spread

// x and y price series
xs: 10 11 12 13 14 15
ys: 20 22 24 26 28 30
// Compute spread (residuals from OLS)
spread: (residuals[xs; ys])
// Z-score for entry/exit signals
z: ((spread - mean spread) % stddev spread)

The OLS residuals εₜ = yₜ − (α̂ + β̂ xₜ) represent the deviation of yₜ from its long-run equilibrium relationship with xₜ. If the pair is cointegrated — meaning εₜ is I(0) rather than I(1) — then εₜ is mean-reverting and the spread oscillates around zero. The Z-score z = (εₜ − ε̄) / σ_ε standardises the spread into units of its own volatility; entry signals at |z| > 2 correspond to the spread being more than two standard deviations away from its historical mean, under the assumption that the mean-reversion speed (Ornstein-Uhlenbeck κ) is large enough to close the gap within the holding period. Verify cointegration formally with an Augmented Dickey-Fuller test on the residuals before trading.

OLS Regression — Prediction

// Fit on historical data
xs: 1 2 3 4 5
ys: 2 4 6 8 10
// Predict for new x values
predict[xs; ys; 6 7 8]
// → 12 14 16

predict computes the OLS estimates α̂ and β̂ from the training pair (xs, ys) and evaluates ŷ = α̂ + β̂ x for each new x in the third argument. In this case the data is perfectly linear (y = 2x), so β̂ = 2, α̂ = 0 exactly and out-of-sample predictions are exact within floating-point precision. In practice, prediction error decomposes as Var[ŷ − y] = σ²(1 + 1/n + (x − x̄)²/Σ(xᵢ − x̄)²), growing quadratically as the new x moves away from the training centroid x̄ — a formal statement that extrapolation amplifies uncertainty. For heteroskedastic residuals or time-series data with autocorrelation, replace OLS standard errors with HAC (Newey-West) estimates before constructing prediction intervals.

Volatility

SVI Volatility Surface Fitting

// Fit SVI to market implied volatilities
// Log-moneyness (k) and market variances (w)
ks:(-0.2),(-0.1),0.0,0.1,0.2
ws:0.04,0.035,0.03,0.035,0.045

// Auto-fit with parameter constraints
params:(svi_fit_auto[ks;ws])
// → fitted [a; b; ρ; m; σ]

// Evaluate fitted curve
fitted:(svi_fitted[params;ks])

// Plot market vs fitted (classic approach)
plot[ks;ws;`MarketVol]
plothline fitted

// Or with the layered chart system:
chart[scatter[ks;ws]; line[ks;fitted]; `SVIFit]

The Stochastic Volatility Inspired (SVI) parameterisation represents the total implied variance w(k) = σ²_impl(k) · T as a function of log-moneyness k = log(K/F), where the five parameters {a, b, ρ, m, σ} control the ATM level, wings slope, skew, ATM shift, and curvature respectively. The parametric form is guaranteed to be free of static arbitrage (no butterfly spread violations) when the Gatheral conditions a ≥ 0, b ≥ 0, |ρ| < 1, σ > 0 are satisfied — constraints that svi_fit_auto enforces during optimisation. Fitted parameters can be used to interpolate implied volatilities across strikes for option pricing, or to extract the local volatility surface via the Dupire formula σ²_loc(k, T) = ∂w/∂T / (1 − k/w · ∂w/∂k + ¼(−¼ − 1/w + k²/w²)(∂w/∂k)²).

Data I/O

Load & Inspect a CSV File

// Load CSV → column-oriented table (headers; columns)
t: csv "data.csv"
t@0              // → `ticker `price `volume `change `sector
(t@1)@1          // → 189.5 420.1 195.3 460.2 185.6 110.8 155.4
mean (t@1)@1     // → 245.27  (mean price)
#((t@1)@0)       // → 7  (row count)

csv reads a CSV file and returns a column-oriented table (headers; columns) — a 2-element list where element 0 is a list of symbols (column names) and element 1 is a list of column vectors. Values are auto-detected as integer, float, or string: a column is promoted to float if any field contains a decimal point, and to string if any field is non-numeric. This Q-style columnar layout is natural for vectorised analytics — applying mean, stddev, or any verb to a column requires no reshaping. Index headers with t@0 and column i with (t@1)@i.

Column Arithmetic on Loaded Data

// Load integer CSV
t: csv "test/test_int.csv"
t@0              // → `a `b `c

// Arithmetic on columns
+/(t@1)@0        // → 6  (sum of column a: 1+2+3)
mean (t@1)@1     // → 5  (mean of column b: 4,5,6)

// Element-wise operations across columns
(t@1)@0 + (t@1)@2   // → 8 10 12  (column a + column c)

Because columns are plain μ lists, all standard verbs and adverbs apply directly — +/ folds a column to its sum, mean computes the arithmetic average, and dyadic + adds two columns element-wise. This is the core advantage of columnar storage: no iteration constructs or index manipulation are needed. For more complex analytics, compose stdlib functions: stddev (t@1)@i for volatility, <(t@1)@i for sorted order, or {(x - mean x) % stddev x} (t@1)@i for z-score normalisation of a column.

Write & Round-Trip

// Load, transform, and write back
t: csv "data.csv"
write["output.csv"; t]

// Verify round-trip
t2: csv "output.csv"
t2@0             // headers preserved
(t2@1)@1         // column values preserved

write["path"; table] serialises a table back to CSV format — headers from t@0 become the first row, and each column vector is written as a comma-separated data row. The round-trip property csv(write(csv f)) ≡ csv f holds for numeric and string data (modulo floating-point display precision). write also accepts a plain list of char-list lines for raw text output. This keyword is CLI-only; in the WASM playground, use the JavaScript bridge MuWASM.loadCSV(url) for fetching remote CSV data.

3D Surfaces

Each example below constructs a 32×32 grid, defines a two-argument surface function, applies it across the grid with outer[f;x;y], and renders the result with plot3d. outer[f;xs;ys] evaluates f(xs@i, ys@j) for every pair (i, j), returning a list-of-lists of shape #xs × #ys — the same primitive used by covmat and corrmat.

Water Drop — Circular Wave Interference

x: -10 + 0.625 * !32
y: -10 + 0.625 * !32
ripple: {[a;b] sin(sqrt((a*a)+(b*b)) * 0.3)}
surface: outer[ripple;x;y]
plot3d[x;y;surface;`WaterDrop]

The ripple function computes the radial distance r = √(x²+y²) for each grid point and then evaluates sin(0.3r), producing concentric wave crests that decrease in frequency as r grows — just as a stone dropped in still water creates circular waves of constant wavelength. !32 generates the integer sequence {0,…,31}; scaling by 0.625 and subtracting 10 maps the grid to [−10, +9.375], centring the origin inside the 32-point window. The outer adverb applies the dyadic lambda to every (x, y) pair without an explicit loop, producing the 32×32 z-matrix that plot3d renders with a height-coloured surface and a floor-grid projection.

Gaussian Bell — 2D Probability Kernel

x: -3 + 0.1875 * !32
y: -3 + 0.1875 * !32
gauss: {[a;b] exp(-((a*a)+(b*b))%8)}
z: outer[gauss;x;y]
plot3d[x;y;z;`GaussianBell]

This is the unnormalised isotropic 2D Gaussian kernel K(x,y) = exp(−r²/(2σ²)) with σ² = 4 (bandwidth parameter 8 = 2σ²). The surface is strictly positive, rotationally symmetric around the origin, and decays to near-zero within the grid window [−3, +2.8]. In kernel density estimation, summing a copy of this kernel centred at each data point gives a smooth probability density estimate; the bandwidth controls the trade-off between bias and variance. The peak at (0, 0) has height 1.0; level sets z = c are circles of radius r = √(−8 ln c). Replacing 8 with a smaller value sharpens the peak, while a larger value flattens it — directly controlling the kernel width.

Saddle — Hyperbolic Paraboloid

x: -2 + 0.125 * !32
y: -2 + 0.125 * !32
saddle: {[a;b] a*a-b*b}
z: outer[saddle;x;y]
plot3d[x;y;z;`Saddle]

The surface z = x² − y² is a hyperbolic paraboloid — a saddle point at the origin where the gradient is zero but the Hessian is indefinite (eigenvalues +2 and −2). Cross-sections at constant y are upward-opening parabolas; at constant x they are downward-opening parabolas. This geometry is central to second-order optimisation: gradient descent on a loss surface with a saddle point stalls at the origin because the gradient vanishes, while curvature information (Newton's method, Adam) can escape by following the negative curvature direction. The grid spans [−2, +1.875] in both axes, covering two full wavelengths of the parabolic arms and making the saddle geometry clearly visible in the rendered surface.

Sinc Surface — Diffraction Pattern

x: -8 + 0.5 * !32
y: -8 + 0.5 * !32
sinc: {[a;b] sin(0.01+sqrt((a*a)+(b*b))) % (0.01+sqrt((a*a)+(b*b)))}
z: outer[sinc;x;y]
plot3d[x;y;z;`SincSurface]

The two-dimensional sinc function sinc(r) = sin(r)/r is the radial Fourier transform of a circular top-hat aperture: it describes the Fraunhofer diffraction pattern of a circular lens or pupil, and the point-spread function of a diffraction-limited optical system (the Airy pattern). The central peak at r = 0 has height 1 (using the limiting value lim_{r→0} sin(r)/r = 1); concentric rings of decreasing amplitude surround it, crossing zero at r = kπ. The small offset 0.01 added to r inside both sin and the denominator guards against exact division by zero when a grid point falls on the origin, without visibly distorting the surface. The grid spans [−8, +7.5], covering roughly 2.5 full periods of the oscillation and showing both the central lobe and several surrounding rings.

Internals

Single-Pass Fusion (Optimizer)

// The optimizer detects the variance / stddev body shape
// and rewrites it to a fused C loop automatically.
// No source change needed.
data: ?50000
variance data   // ~0.3 ms  (fused)

// Aliases also fuse — avg is structurally equal to mean:
myvar: {(+/(x-avg x)*(x-avg x))%(#x)}
myvar data      // ~0.3 ms  (fused via inlining)

// Without fusion: ~13 ms (5 passes, 3 x 50k intermediate lists)

The naïve two-pass variance formula Var(x) = (1/n)Σ(xᵢ − x̄)² allocates three intermediate lists of length n (centred values, squared values, and the running sum), touching memory 5 times total. The optimizer rewrites it to the equivalent single-pass identity Var(x) = (Σxᵢ²)/n − ((Σxᵢ)/n)², which requires only two scalar accumulators and zero heap allocation — a ~40× speedup on 50 000 elements under -O3 -ffast-math on Apple M-series. The rewrite operates on the AST: the opt pass first inlines known stdlib definitions into lambda bodies, then pattern-matches the variance/stddev structural shape irrespective of whether mean, avg, or any user-defined alias is used. This means the fusion applies automatically to any expression structurally equivalent to the variance formula, including user-defined rolling statistics or weighted variants, without requiring any source-level change.

Monte Carlo Pricing

See the grammar page for full Monte Carlo examples: a European call option via Box-Muller (mcCall) and an arithmetic-average Asian call via path simulation (mcAsian). Visual plots and convergence stories are on the plots page.

List Comprehensions

Square Numbers

squares: [x*x ~ x <- 1 2 3 4 5 6 7 8 9 10]
squares
// → 1 4 9 16 25 36 49 64 81 100

The expression [x*x ~ x <- 1 2 3] reads as "x squared for each x drawn from the list 1 2 3". It desugars to {[x] x*x}' 1 2 3 — a lambda mapped over the source list with the each adverb. List comprehensions provide declarative syntax for common map operations without explicit lambda syntax.

Filter with Guard — Even Numbers

evens: [x ~ x <- !20; x%2=0]
evens
// → 0 2 4 6 8 10 12 14 16 18

A guard (x%2=0) filters elements — only values where the guard is truthy are included. Multiple guards can be chained with semicolons. The desugared form is a conditional lambda: {[x] x%2=0; x; _}' !20 where _ (underscore) means "drop this element". This compiles to efficient bytecode with no intermediate allocations.

Cartesian Product of Indices

pairs: [[i;j] ~ i <- !3; j <- !3]
pairs
// → [0;0] [0;1] [0;2] [1;0] [1;1] [1;2] [2;0] [2;1] [2;2]

Multiple generators nest left-to-right: for each i in 0..2, iterate j in 0..2. The result is a flat list of vector pairs — useful for iterating over matrix indices, building lookup tables, or enumerating combinatorial spaces. The desugared form nests each' calls: {[i] {[j] [i;j]}' !3}' !3.