diff --git a/.gitignore b/.gitignore
index 998edbd..fc474d4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,5 @@
+.DS_Store
+
# Generated by Cargo
# will have compiled files and executables
debug/
diff --git a/DEMO-EXAMPLES.md b/DEMO-EXAMPLES.md
new file mode 100644
index 0000000..4a45f9b
--- /dev/null
+++ b/DEMO-EXAMPLES.md
@@ -0,0 +1,309 @@
+# Goth Demo Examples - Planned Simulation & Numerical Methods
+
+This document outlines a comprehensive list of example/demo code that would showcase Goth's capabilities for scientific computing, numerical simulation, and computational mathematics.
+
+## Category 1: Classic Numerical Methods
+
+### 1.1 Newton-Raphson Method
+- **Purpose**: Root-finding algorithm showcase
+- **Key Features**:
+ - First-order derivative computation
+ - Iterative convergence criteria
+ - Float precision handling
+ - Error tolerance checking
+- **Example Use Cases**:
+ - Finding cube roots
+ - Solving nonlinear equations
+ - Demonstrating uncertainty propagation in numerical methods
+
+### 1.2 Runge-Kutta Integration (RK4)
+- **Purpose**: Ordinary Differential Equation (ODE) solver
+- **Key Features**:
+ - Four-stage explicit integration
+ - Step size adaptive behavior
+ - State vector management
+ - Temporal evolution
+- **Example Use Cases**:
+ - Projectile motion with drag
+ - Pendulum dynamics
+ - Population growth models
+ - Spring-mass systems
+
+### 1.3 Monte Carlo Methods
+- **Purpose**: Probabilistic simulation and numerical integration
+- **Key Features**:
+ - Random sampling
+ - Statistical aggregation
+ - Uncertainty quantification
+ - Convergence monitoring
+- **Example Use Cases**:
+ - Integration by random sampling
+ - Pi estimation
+ - Portfolio risk analysis
+ - Particle transport simulation
+
+### 1.4 Quadrature Methods (Numerical Integration)
+- **Purpose**: Computing definite integrals numerically
+- **Key Features**:
+ - Simpson's rule
+ - Gaussian quadrature
+ - Adaptive refinement
+- **Example Use Cases**:
+ - Arc length calculation
+ - Work integrals
+ - Probability distributions
+
+## Category 2: Partial Differential Equations (PDEs)
+
+### 2.1 Finite Difference Methods
+- **Purpose**: Discretizing and solving PDEs
+- **Key Features**:
+ - Grid-based computation
+ - Boundary condition handling
+ - Stencil operations
+ - Iterative solvers
+
+#### 2.1.1 Laplace Equation (Electrostatics/Heat Diffusion)
+- **Purpose**: Classic boundary value problem
+- **Key Features**:
+ - 2D grid iteration
+ - Relaxation methods (Jacobi, Gauss-Seidel)
+ - Convergence criteria
+- **Physics**: Electric potential, steady-state heat distribution
+
+#### 2.1.2 Poisson Equation (Inhomogeneous)
+- **Purpose**: Extension of Laplace with source terms
+- **Key Features**:
+ - Source term evaluation
+ - Iterative refinement
+ - Boundary conditions
+- **Physics**: Gravitational potential, steady-state heat with sources
+
+#### 2.1.3 Advection Equation
+- **Purpose**: Transport phenomena
+- **Key Features**:
+ - First-order hyperbolic PDE
+ - Characteristic-based methods
+ - Flux limiters for stability
+- **Physics**: Scalar transport, pollutant diffusion, flow advection
+
+#### 2.1.4 Wave Equation
+- **Purpose**: Hyperbolic wave propagation
+- **Key Features**:
+ - Second-order temporal integration
+ - Stability conditions (CFL)
+ - Boundary reflections
+- **Physics**: Acoustic waves, electromagnetic waves, vibrating strings
+
+#### 2.1.5 Burgers' Equation
+- **Purpose**: Nonlinear convection + diffusion
+- **Key Features**:
+ - Shock formation
+ - Cole-Hopf transformation
+ - Artificial diffusion
+- **Physics**: Shock waves, traffic flow, turbulence modeling
+
+### 2.2 Finite Element Method (FEA)
+- **Purpose**: Variational approach to PDEs
+- **Key Features**:
+ - Basis function construction
+ - Assembly of stiffness matrix
+ - Sparse linear system solving
+ - Post-processing (stress/strain)
+- **Example Problems**:
+ - 1D bar under load (tension/compression)
+ - 2D plate bending
+ - Thermal analysis
+ - Modal analysis (eigenvalues)
+
+### 2.3 Computational Fluid Dynamics (CFD)
+- **Purpose**: Simulating fluid flow
+- **Key Features**:
+ - Navier-Stokes equation discretization
+ - Pressure-velocity coupling (SIMPLE algorithm)
+ - Turbulence modeling (k-ε models)
+ - Incompressible flow assumptions
+
+#### 2.3.1 Incompressible Flow Examples
+- Lid-driven cavity flow
+- Channel flow with obstacles
+- Flow around a cylinder
+- Thermal convection
+
+#### 2.3.2 Shallow Water Equations
+- **Purpose**: Simplified 2D fluid equations
+- **Features**:
+ - Surface height tracking
+ - Velocity fields
+ - Reduced dimensionality
+- **Applications**: Dam break simulation, tsunami modeling
+
+## Category 3: Linear Algebra & Eigenvalue Problems
+
+### 3.1 Eigenvalue/Eigenvector Computation
+- **Purpose**: Finding characteristic modes and values
+- **Key Features**:
+ - Power iteration method
+ - Inverse iteration
+ - QR algorithm
+ - Rayleigh quotient
+- **Example Use Cases**:
+ - Structural natural frequencies
+ - Stability analysis
+ - Principal component analysis
+ - Graph spectral methods
+
+### 3.2 Linear System Solvers
+- **Purpose**: Solving Ax = b
+- **Key Features**:
+ - Gaussian elimination
+ - LU decomposition
+ - Iterative methods (CG, GMRES)
+ - Sparse matrix handling
+- **Example Use Cases**:
+ - FEA system assembly
+ - CFD pressure correction
+ - Circuit analysis
+
+## Category 4: Optimization
+
+### 4.1 Gradient Descent Variants
+- **Purpose**: Unconstrained optimization
+- **Key Features**:
+ - Steepest descent
+ - Conjugate gradient
+ - BFGS quasi-Newton method
+ - Learning rate scheduling
+- **Example Use Cases**:
+ - Function minimization
+ - Machine learning parameter fitting
+ - Inverse problems
+
+### 4.2 Constrained Optimization
+- **Purpose**: Optimization with constraints
+- **Key Features**:
+ - Penalty methods
+ - Lagrange multipliers
+ - Sequential quadratic programming
+- **Example Use Cases**:
+ - Portfolio optimization
+ - Shape optimization
+ - Control problems
+
+## Category 5: Interpolation & Approximation
+
+### 5.1 Polynomial Interpolation
+- **Purpose**: Function reconstruction from points
+- **Key Features**:
+ - Lagrange interpolation
+ - Newton divided differences
+ - Spline interpolation (cubic, B-splines)
+- **Example Use Cases**:
+ - Data smoothing
+ - Function approximation
+ - Mesh refinement
+
+### 5.2 Least Squares Fitting
+- **Purpose**: Finding best-fit functions
+- **Key Features**:
+ - Linear least squares
+ - Nonlinear fitting
+ - Regularization (Ridge, LASSO)
+- **Example Use Cases**:
+ - Curve fitting
+ - Inverse problems
+ - Data analysis
+
+## Category 6: Statistical Methods
+
+### 6.1 Bayesian Inference
+- **Purpose**: Probabilistic inference
+- **Key Features**:
+ - Prior/likelihood/posterior
+ - MCMC sampling
+ - Variational inference
+- **Example Use Cases**:
+ - Parameter estimation
+ - Model selection
+ - Uncertainty quantification
+
+### 6.2 Uncertainty Propagation
+- **Purpose**: Tracking error through computations
+- **Key Features**:
+ - Dual number arithmetic
+ - Adjoint methods
+ - Sensitivity analysis
+ - Forward/reverse mode AD
+- **Example Use Cases**:
+ - Measurement error analysis
+ - Model parameter sensitivity
+ - Risk assessment
+
+## Category 7: Specialized Simulation
+
+### 7.1 N-Body Simulations
+- **Purpose**: Particle dynamics
+- **Key Features**:
+ - Pairwise force computation
+ - Time stepping
+ - Spatial acceleration structures
+ - Energy conservation monitoring
+- **Example Use Cases**:
+ - Gravitational N-body
+ - Molecular dynamics
+ - Particle systems
+
+### 7.2 Finite Difference Time Domain (FDTD)
+- **Purpose**: Electromagnetic wave simulation
+- **Key Features**:
+ - Yee grid discretization
+ - Staggered grids
+ - PML boundary conditions
+- **Example Use Cases**:
+ - Antenna design
+ - Photonic devices
+ - Waveguide analysis
+
+### 7.3 Lattice Boltzmann Method (LBM)
+- **Purpose**: Alternative to direct PDE solving
+- **Key Features**:
+ - Distribution functions
+ - Collision operators
+ - Streaming step
+- **Example Use Cases**:
+ - Complex flow geometry
+ - Multiphase flows
+ - Reactive transport
+
+## Suggested Priority Examples (Tier 1)
+
+These would be the best starting points to showcase Goth's capabilities:
+
+1. **Newton-Raphson for Square Root** - Simple, demonstrates iterative methods and error handling
+2. **Simple ODE with Runge-Kutta** - Projectile motion or pendulum
+3. **2D Laplace Solver** - Heat diffusion on a 2D grid
+4. **Monte Carlo Pi Estimation** - Shows randomness and aggregation
+5. **Simple Linear System Solver** - 3x3 or 5x5 system with Gaussian elimination
+6. **1D Wave Equation** - Shows spatiotemporal discretization
+7. **Eigenvalue Power Iteration** - Simple symmetric matrix case
+8. **Gradient Descent Optimization** - Minimize a simple 2D function
+
+## Key Technical Requirements for Examples
+
+- Clear variable naming showing Goth's type system
+- Uncertainty/dual number demonstrations (AD support)
+- Performance characteristics (iteration counts, convergence rates)
+- Comparison with analytical solutions where available
+- Comments explaining the numerical method
+- Output formatting for visualization-ready data
+- Memory-efficient handling of arrays/vectors
+- Floating-point precision considerations
+
+## Platform Considerations
+
+- Performance benchmarks (C vs interpreted)
+- Memory usage patterns
+- Compilation time
+- Ease of visualization (CSV/JSON output)
+- Integration with Python/plotting tools (matplotlib, Paraview)
+
diff --git a/DEMO-OUTPUT-FORMATS.md b/DEMO-OUTPUT-FORMATS.md
new file mode 100644
index 0000000..e7f1220
--- /dev/null
+++ b/DEMO-OUTPUT-FORMATS.md
@@ -0,0 +1,821 @@
+# Goth Demo Examples - CSV and SVG Output Formats
+
+This document specifies how Goth demo examples should output results in CSV and SVG formats for visualization and analysis.
+
+## Part 1: CSV Output Formats
+
+### 1.1 General CSV Structure
+
+All CSV outputs should follow these standards:
+- UTF-8 encoding
+- Standard comma delimiter
+- Headers in first row
+- Quoted fields containing special characters
+- ISO 8601 timestamps where applicable
+- Consistent precision (e.g., 6 decimal places for floats)
+
+### 1.2 CSV Format by Example Type
+
+#### 1.2.1 Iterative Methods (Newton-Raphson)
+```
+iteration,x_current,x_previous,f_x,df_x,error,converged
+0,2.0,NaN,7.0,12.0,NaN,false
+1,1.41667,2.0,0.00694,4.83333,0.58333,false
+2,1.41421,1.41667,0.0000061,4.828,0.00246,true
+```
+
+**Fields:**
+- `iteration`: Iteration counter (int)
+- `x_current`: Current estimate (float)
+- `x_previous`: Previous estimate (float)
+- `f_x`: Function value at current x
+- `df_x`: Derivative value at current x
+- `error`: Absolute error |x_current - x_previous|
+- `converged`: Boolean convergence flag
+
+#### 1.2.2 ODE Integration (Runge-Kutta)
+```
+time,x,y,vx,vy,ke,pe,total_energy,step_size
+0.0,1.0,0.0,0.0,1.0,0.5,0.5,1.0,0.01
+0.01,1.00995,0.00995,0.0995,0.99975,0.499875,0.500125,1.0,0.01
+0.02,1.01980,0.01980,0.198,0.99901,0.499505,0.500495,1.0,0.01
+```
+
+**Fields:**
+- `time`: Current time value
+- `x`, `y`: Position coordinates
+- `vx`, `vy`: Velocity components
+- `ke`: Kinetic energy
+- `pe`: Potential energy
+- `total_energy`: Total energy (for conservation checking)
+- `step_size`: Adaptive step size used
+
+#### 1.2.3 PDE Solutions (Laplace/Poisson 2D Grid)
+```
+x,y,phi,source,iteration
+0.0,0.0,1.0,0.0,500
+0.1,0.0,0.95432,0.0,500
+0.2,0.0,0.87654,0.0,500
+0.0,0.1,0.96789,0.1,500
+```
+
+**Fields:**
+- `x`, `y`: Grid coordinates
+- `phi`: Solution value (potential, temperature, etc.)
+- `source`: Source term value (for Poisson)
+- `iteration`: Iteration count at which output occurred
+
+#### 1.2.4 Monte Carlo Simulation (Pi Estimation)
+```
+sample_count,pi_estimate,error,std_dev,ci_lower,ci_upper
+1000,3.144,0.00242,0.0512,3.044,3.244
+10000,3.14159,0.0000107,0.0161,3.109,3.173
+100000,3.141592,0.0000011,0.00512,3.130,3.153
+1000000,3.1415926,0.00000011,0.00162,3.138,3.144
+```
+
+**Fields:**
+- `sample_count`: Number of samples
+- `pi_estimate`: Current estimate of pi
+- `error`: Error from true value (3.14159265...)
+- `std_dev`: Standard deviation
+- `ci_lower`, `ci_upper`: 95% confidence interval bounds
+
+#### 1.2.5 Linear System Solution (Gaussian Elimination)
+```
+iteration,residual_norm,solution_norm,condition_number
+0,15.234,NaN,NaN
+1,2.1456,3.4521,12.45
+2,0.00234,3.4589,12.45
+3,0.00001,3.4590,12.45
+```
+
+**Fields:**
+- `iteration`: Iteration number
+- `residual_norm`: ||Ax - b|| (measure of solution quality)
+- `solution_norm`: ||x||
+- `condition_number`: Condition number of matrix
+
+#### 1.2.6 Wave Equation (1D or 2D)
+```
+time,x,u,u_t,u_x,u_xx
+0.0,0.0,1.0,0.0,0.0,0.0
+0.0,0.1,0.99499,0.0,−0.05,−0.1
+0.01,0.0,1.0,0.0,0.0,0.0
+0.01,0.1,0.98994,−0.05,−0.05,−0.1
+```
+
+**Fields:**
+- `time`: Current time
+- `x`: Position (or x,y for 2D)
+- `u`: Solution value
+- `u_t`: Temporal derivative
+- `u_x`: Spatial derivative (1st)
+- `u_xx`: Spatial derivative (2nd)
+
+#### 1.2.7 Eigenvalue Computation
+```
+iteration,lambda_1,lambda_2,lambda_3,residual_1,residual_2,residual_3,converged
+0,1.5,NaN,NaN,2.34,NaN,NaN,false
+1,4.875,NaN,NaN,0.234,NaN,NaN,false
+2,5.0,NaN,NaN,0.00123,NaN,NaN,true
+```
+
+**Fields:**
+- `iteration`: Iteration number
+- `lambda_i`: i-th eigenvalue estimate
+- `residual_i`: Residual for i-th eigenvalue
+- `converged`: Overall convergence flag
+
+#### 1.2.8 Optimization (Gradient Descent)
+```
+iteration,x1,x2,f_value,gradient_norm,step_size,line_search_evals
+0,0.0,0.0,4.0,4.472,NaN,NaN
+1,−0.1,−0.1,3.96,4.412,0.1,3
+2,−0.198,−0.198,3.84,4.234,0.1,2
+3,−0.291,−0.291,3.65,3.98,0.1,2
+```
+
+**Fields:**
+- `iteration`: Iteration number
+- `x1`, `x2`, ...: Decision variables
+- `f_value`: Objective function value
+- `gradient_norm`: ||∇f||
+- `step_size`: Step size used
+- `line_search_evals`: Number of function evaluations in line search
+
+#### 1.2.9 Interpolation/Fitting Results
+```
+x_data,y_data,y_fit,residual,local_error
+0.0,0.0,0.001,−0.001,0.0001
+0.1,0.0995,0.0998,−0.0003,0.00001
+0.2,0.1987,0.1985,0.0002,0.00002
+```
+
+**Fields:**
+- `x_data`: Original x values
+- `y_data`: Original y values
+- `y_fit`: Fitted/interpolated values
+- `residual`: Difference (y_data - y_fit)
+- `local_error`: Local truncation error estimate
+
+### 1.3 CSV Export Function Template (Pseudo-Goth)
+
+```goth
+fn export_csv_header(fields: Array[String]) -> String {
+ join(fields, ",")
+}
+
+fn export_csv_row(values: Array[Float | Int | Bool | String]) -> String {
+ mapped = map(values, fn(v: Float | Int | Bool | String) -> String {
+ match v {
+ String(s) => quote_if_needed(s)
+ Float(f) => format_float(f, 6)
+ Int(i) => to_string(i)
+ Bool(b) => b ? "true" : "false"
+ }
+ })
+ join(mapped, ",")
+}
+
+fn write_csv_file(filename: String, headers: Array[String],
+ rows: Array[Array[Float | Int | Bool | String]]) -> IO[Unit] {
+ handle = open_file(filename, "w")
+ write_line(handle, export_csv_header(headers))
+ foreach(rows, fn(row) {
+ write_line(handle, export_csv_row(row))
+ })
+ close_file(handle)
+}
+```
+
+---
+
+## Part 2: SVG Output Formats
+
+### 2.1 General SVG Structure
+
+All SVG outputs should:
+- Use standard SVG 1.1 specification
+- Include viewBox for proper scaling
+- Include descriptive titles and comments
+- Use consistent styling
+- Include legends where applicable
+- Provide 600x400 minimum resolution (scalable via viewBox)
+- Use readable fonts (Arial, sans-serif)
+
+### 2.2 SVG Output Formats by Example Type
+
+#### 2.2.1 Convergence Plot (Iterative Methods)
+
+**File: `newton_raphson_convergence.svg`**
+
+```xml
+
+
+```
+
+#### 2.2.2 Phase Space Plot (ODE Solutions)
+
+**File: `trajectory_phase_space.svg`**
+
+```xml
+
+
+```
+
+#### 2.2.3 Heat Map (2D PDE Solution)
+
+**File: `laplace_solution_heatmap.svg`**
+
+```xml
+
+
+```
+
+#### 2.2.4 Multi-Series Line Plot (Multiple Solutions)
+
+**File: `multi_comparison.svg`**
+
+```xml
+
+
+```
+
+#### 2.2.5 Convergence Rate Plot (Log Scale)
+
+**File: `convergence_log_scale.svg`**
+
+```xml
+
+
+```
+
+#### 2.2.6 Error Distribution Plot (Histogram)
+
+**File: `error_distribution.svg`**
+
+```xml
+
+
+```
+
+### 2.3 SVG Export Function Template (Pseudo-Goth)
+
+```goth
+fn svg_header(width: Int, height: Int, title: String) -> String {
+ format_string(
+ "\n",
+ "\n"
+}
+
+fn generate_convergence_svg(iterations: Array[Int],
+ errors: Array[Float],
+ filename: String) -> IO[Unit] {
+ handle = create_file(filename)
+ write(handle, svg_header(800, 600, "Convergence Plot"))
+
+ // Scale data to fit SVG coordinates
+ max_iter = max(iterations)
+ max_error = max(errors)
+ scale_x = 700.0 / float(max_iter)
+ scale_y = 450.0 / log10(max_error)
+
+ points = map(zip(iterations, errors), fn(p) -> (Float, Float) {
+ (50.0 + float(p.0) * scale_x, 500.0 - log10(p.1) * scale_y)
+ })
+
+ write(handle, svg_polyline(points, "#2196F3", 2.0, "none"))
+ write(handle, svg_footer())
+ close_file(handle)
+}
+```
+
+---
+
+## Part 3: Integration Strategy
+
+### 3.1 Automatic Output Generation
+
+Each example should include a function that:
+1. Generates intermediate results
+2. Collects data into arrays
+3. Exports to CSV automatically after computation
+4. Generates SVG visualization from CSV data
+
+### 3.2 Example Output Structure
+
+For a single example, expect these files:
+
+```
+example_newton_raphson/
+├── newton_raphson.goth # Main Goth source
+├── output/
+│ ├── convergence_data.csv # Raw numerical data
+│ ├── convergence_plot.svg # Visualization
+│ └── analysis.txt # Convergence metrics
+```
+
+### 3.3 Standardized Metadata in Outputs
+
+Both CSV headers and SVG metadata should include:
+- Algorithm/method name
+- Date and time generated
+- Version of Goth compiler
+- Parameters used
+- Domain/problem description
+- Convergence criteria
+- Execution time
+
diff --git a/LANGUAGE.md b/LANGUAGE.md
index 9580b18..8148afa 100644
--- a/LANGUAGE.md
+++ b/LANGUAGE.md
@@ -55,10 +55,16 @@ Goth uses de Bruijn indices instead of named variables. This eliminates shadowin
λ→ λ→ ₀ + ₁ # add: ₀ is second (inner) arg, ₁ is first (outer) arg
```
-In function declarations, arguments bind in order (`₀` = first arg):
+In multi-arg function declarations, `₀` = **last** argument (most recently bound):
```goth
╭─ sub : I64 → I64 → I64
-╰─ ₀ - ₁ # ₀ = first arg, ₁ = second arg (sub 10 3 = 7)
+╰─ ₁ - ₀ # ₁ = first arg, ₀ = second arg (sub 10 3 = 7)
+```
+
+For single-arg functions, `₀` is the sole argument:
+```goth
+╭─ square : I64 → I64
+╰─ ₀ × ₀ # ₀ = the argument
```
Let bindings shift indices:
@@ -202,6 +208,7 @@ let v : [5]F64 ← [1.0, 2.0, 3.0] in v -- Error: shape mismatch
| `/` | `/` | Division |
| `%` | `%` | Modulo |
| `^` | `^` | Power |
+| `±` | `+-` | Uncertainty |
### Comparison
@@ -234,12 +241,14 @@ let v : [5]F64 ← [1.0, 2.0, 3.0] in v -- Error: shape mismatch
### Indexing
```goth
-arr[0] -- First element
+arr[0] -- First element (bracket must be adjacent, no space)
arr[i] -- Element at index i
tuple.0 -- First tuple field
record.x -- Named field x
```
+**Note:** The `[` must be directly adjacent to the expression (no space) to be parsed as indexing. With a space, `f [1,2]` is function application (passing the array `[1,2]` to `f`).
+
---
## Enums and Pattern Matching
@@ -293,6 +302,8 @@ enum Option α where Some α | None
## Effects
+> **Note:** Effect annotations are parsed into the AST but not yet enforced. Functions with I/O side effects work with or without annotations. See `docs/EFFECT-SYSTEM-ROADMAP.md` for the enforcement plan.
+
Effects are declared with `◇` annotations:
| Effect | Description |
@@ -307,7 +318,7 @@ Effects are declared with `◇` annotations:
╰─ print "Hello!"
```
-Pure functions (no effects) need no annotation.
+Pure functions (no effects) need no annotation. Currently, effect annotations serve as documentation — the runtime does not reject effectful operations in unannotated functions.
---
@@ -316,16 +327,13 @@ Pure functions (no effects) need no annotation.
### Imports
```goth
-use stdlib.prelude
-use stdlib.option
-use mylib.utils
+use "prelude.goth"
+use "../stdlib/option.goth"
```
### Module Files
-Each `.goth` file is a module. The module path corresponds to the file path:
-- `stdlib/prelude.goth` → `stdlib.prelude`
-- `stdlib/option.goth` → `stdlib.option`
+Each `.goth` file is a module. The `use` declaration takes a string path (relative to the importing file) and inlines all declarations from that file into the current namespace.
---
@@ -344,7 +352,7 @@ Each `.goth` file is a module. The module path corresponds to the file path:
|------|-----------|-------------|
| `Σ`, `sum` | `[n]α → α` | Sum elements |
| `Π`, `prod` | `[n]α → α` | Product elements |
-| `length` | `[n]α → I64` | Array length |
+| `len` | `[n]α → I64` | Array length |
### Transformations
@@ -355,7 +363,8 @@ Each `.goth` file is a module. The module path corresponds to the file path:
| `reverse` | `[n]α → [n]α` | Reverse order |
| `take` | `I64 → [n]α → [m]α` | Take first k elements |
| `drop` | `I64 → [n]α → [m]α` | Drop first k elements |
-| `⧺`, `++` | `[n]α → [m]α → [p]α` | Concatenate arrays |
+| `⧺` | `String → String → String` | Concatenate strings |
+| `⊕` | `[n]α → [m]α → [p]α` | Concatenate arrays |
### Linear Algebra
@@ -377,6 +386,25 @@ Each `.goth` file is a module. The module path corresponds to the file path:
| `floor`, `ceil`, `round` | `F64 → F64` |
| `abs` | Numeric → Numeric |
+### Uncertainty Propagation
+
+Goth supports first-class uncertain values using the `±` operator. When uncertain values flow through arithmetic and math functions, uncertainty propagates automatically.
+
+**Creating uncertain values:**
+```goth
+10.5 ± 0.3 -- 10.5 with uncertainty 0.3
+```
+
+**Automatic propagation through operations:**
+```goth
+(10.0 ± 0.3) + (20.0 ± 0.4) -- 30 ± 0.5 (additive: δ = √(δa² + δb²))
+(5.0 ± 0.1) × (3.0 ± 0.2) -- 15 ± 1.04 (relative: quadrature sum)
+√(9.0 ± 0.3) -- 3 ± 0.05 (derivative: δ/(2√x))
+sin (1.0 ± 0.1) -- 0.841 ± 0.054 (|cos x| × δx)
+```
+
+Supported functions: `+`, `-`, `×`, `/`, `^`, `√`, `exp`, `ln`, `log10`, `log2`, `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `sinh`, `cosh`, `tanh`, `abs`, `floor`, `ceil`, `round`, `Γ`.
+
### Type Conversions
| Name | Signature | Description |
@@ -391,8 +419,13 @@ Each `.goth` file is a module. The module path corresponds to the file path:
| Name | Signature | Description |
|------|-----------|-------------|
-| `print` | `α → ()` | Print to stdout |
+| `print` | `α → ()` | Print to stdout (with newline) |
| `readLine` | `() → String` | Read line from stdin |
+| `readFile` | `String → String` | Read file contents |
+| `writeFile` | `String → String → ()` | Write content to file path |
+| `▷` | `String → String → ()` | Write operator: `"content" ▷ "path"` |
+| `toString` | `α → String` | Convert value to string |
+| `strConcat`, `⧺` | `String → String → String` | Concatenate strings |
---
diff --git a/LLM_SPEC.md b/LLM_SPEC.md
index e5bac25..d9eaa29 100644
--- a/LLM_SPEC.md
+++ b/LLM_SPEC.md
@@ -16,14 +16,14 @@ If: if cond then a else b
## De Bruijn Indices (Critical)
Variables are numbered by binding depth, not named:
-- `₀` (or `_0`) = first parameter / most recent binding
-- `₁` (or `_1`) = second parameter / one binding out
-- `₂` (or `_2`) = third parameter / two bindings out
+- `₀` (or `_0`) = most recent binding / last parameter
+- `₁` (or `_1`) = one binding out / second-to-last parameter
+- `₂` (or `_2`) = two bindings out / third-to-last parameter
```goth
-# Function args: ₀ = first arg, ₁ = second arg
+# Function args: ₀ = last arg, ₁ = second-to-last arg
╭─ sub : I64 → I64 → I64
-╰─ ₀ - ₁ # sub 10 3 = 7
+╰─ ₁ - ₀ # sub 10 3 = 7 (₁=10, ₀=3)
# Let shifts indices
╭─ example : I64 → I64
@@ -50,7 +50,7 @@ I64 → I64 Function
## Common Operators
```
-Arithmetic: + - × / % ^ (% also written as mod)
+Arithmetic: + - × / % ^ ± (% also written as mod, ± also +-)
Comparison: = ≠ < > ≤ ≥
Equality: = (value) ≡/== (structural) ≣/=== (referential, reserved)
Logical: ∧ ∨ ¬ (or && || !)
@@ -160,6 +160,19 @@ let v : [5]F64 ← [1.0, 2.0, 3.0] in v # Error: shape mismatch
1 + 2 × x + 3 × x × x
```
+## Uncertainty Propagation
+
+Goth has first-class uncertain values. Create with `±`, and uncertainty propagates automatically through arithmetic and math functions.
+
+```goth
+# Create uncertain value
+╭─ main : F64 → F64 → (F64 ± F64)
+╰─ let a ← ₁ ± ₀ in # a = value ± uncertainty
+ √₀ + (2.0 ± 0.1) # propagates through √ and +
+```
+
+Propagation rules: additive (√(δa²+δb²)) for `+`/`-`, relative error for `×`/`/`, derivative-based for math functions (`√`, `sin`, `cos`, `exp`, `ln`, etc.).
+
## Common Mistakes
1. **Wrong index after let**: Each `let` shifts indices by 1
@@ -178,7 +191,13 @@ let v : [5]F64 ← [1.0, 2.0, 3.0] in v # Error: shape mismatch
# Argument would be ₂ inside the lambda
```
-3. **Type mismatch I vs I64**: Use `I64` in signatures
+3. **Indexing vs application**: `arr[i]` (no space) is indexing; `f [1,2]` (with space) is function application
+ ```goth
+ let arr = [10, 20, 30] in arr[1] # Indexing → 20
+ let f = (λ→ Σ ₀) in f [10, 20, 30] # Application → 60
+ ```
+
+4. **Type mismatch I vs I64**: Use `I64` in signatures
```goth
# WRONG
╭─ main : () → I
@@ -219,9 +238,18 @@ let v : [5]F64 ← [1.0, 2.0, 3.0] in v # Error: shape mismatch
| `⊤` | `true` |
| `⊥` | `false` |
| `ι` | `iota` |
-| `⧺` | `++` |
+| `⧺` | `strConcat` |
+| `⊕` | `concat` |
| `⍉` | `transpose` |
| `·` | `dot` |
+| `±` | `+-` |
+
+## Enforcement Notes
+
+- **Contracts** (`⊢` preconditions, `⊨` postconditions): enforced at runtime
+- **Effect annotations** (`◇io`, `◇mut`, `◇rand`): parsed but **not enforced** — `print` works without `◇io`. Use effect annotations as documentation only; do not rely on them for correctness.
+- **Type annotations**: parsed but type checking is partial
+- **Refinement types** (`{x : F64 | x > 0}`): parsed but predicates not solved
## Running Code
diff --git a/benchmark/tests/uncertainty.json b/benchmark/tests/uncertainty.json
index fa9c234..fcc1cab 100644
--- a/benchmark/tests/uncertainty.json
+++ b/benchmark/tests/uncertainty.json
@@ -12,6 +12,56 @@
{"input": [0.0, 0.1], "expected": "0±0.1"},
{"input": [3.14159, 0.00001], "expected": "3.14159±0.00001"}
]
+ },
+ {
+ "name": "add_uncertain",
+ "file": "examples/uncertainty/add_uncertain.goth",
+ "description": "Add two uncertain values with additive propagation: δ = √(δa² + δb²)",
+ "cases": [
+ {"input": [10.0, 0.3, 20.0, 0.4], "expected": "30±0.5"},
+ {"input": [100.0, 1.0, 200.0, 2.0], "expected": "300±2.23606797749979"},
+ {"input": [5.0, 0.5, -3.0, 0.5], "expected": "2±0.7071067811865476"}
+ ]
+ },
+ {
+ "name": "mul_uncertain",
+ "file": "examples/uncertainty/mul_uncertain.goth",
+ "description": "Multiply two uncertain values with relative error propagation",
+ "cases": [
+ {"input": [5.0, 0.1, 3.0, 0.2], "expected": "15±1.044030650891055"},
+ {"input": [10.0, 0.5, 4.0, 0.2], "expected": "40±2.8284271247461907"},
+ {"input": [2.0, 0.1, 8.0, 0.4], "expected": "16±1.1313708498984762"}
+ ]
+ },
+ {
+ "name": "sqrt_uncertain",
+ "file": "examples/uncertainty/sqrt_uncertain.goth",
+ "description": "Square root with derivative-based propagation: δ = δx / (2√x)",
+ "cases": [
+ {"input": [9.0, 0.3], "expected": "3±0.049999999999999996"},
+ {"input": [25.0, 1.0], "expected": "5±0.1"},
+ {"input": [4.0, 0.2], "expected": "2±0.05"}
+ ]
+ },
+ {
+ "name": "sin_uncertain",
+ "file": "examples/uncertainty/sin_uncertain.goth",
+ "description": "Sine with trig propagation: δ = |cos(x)| × δx",
+ "cases": [
+ {"input": [1.0, 0.1], "expected": "0.8414709848078965±0.05403023058681398"},
+ {"input": [0.0, 0.05], "expected": "0±0.05"},
+ {"input": [1.5708, 0.01], "expected": "0.9999999999932537±0.00000003673205103346574"}
+ ]
+ },
+ {
+ "name": "chained_uncertain",
+ "file": "examples/uncertainty/chained_uncertain.goth",
+ "description": "Multi-step propagation: sin(√(a ± δa) + (b ± δb))",
+ "cases": [
+ {"input": [4.0, 0.2, 1.0, 0.1], "expected": "0.1411200080598672±0.11068452598066628"},
+ {"input": [9.0, 0.5, 0.0, 0.1], "expected": "0.1411200080598672±0.1288681429287275"},
+ {"input": [1.0, 0.1, 0.5, 0.05], "expected": "0.9974949866040544±0.0050018754981393096"}
+ ]
}
]
}
diff --git a/crates/goth-eval/src/eval.rs b/crates/goth-eval/src/eval.rs
index 33191cc..be78551 100644
--- a/crates/goth-eval/src/eval.rs
+++ b/crates/goth-eval/src/eval.rs
@@ -132,7 +132,7 @@ impl Evaluator {
Expr::ArrayFill { shape, value } => { let shape_vals: Vec = shape.iter().map(|e| { let v = self.eval_with_env(e, env)?; v.as_int().map(|n| n as usize).ok_or_else(|| EvalError::type_error("Int", &v)) }).collect::>()?; let fill_val = self.eval_with_env(value, env)?; let size: usize = shape_vals.iter().product(); let data = vec![fill_val; size]; Ok(Value::Tensor(Tensor::from_values(shape_vals, data))) }
Expr::Variant { constructor, payload } => { let payload_val = match payload { Some(p) => Some(self.eval_with_env(p, env)?), None => None }; Ok(Value::variant(constructor.to_string(), payload_val)) }
Expr::Field(base, access) => { let val = self.eval_with_env(base, env)?; self.access_field(val, access) }
- Expr::Index(base, indices) => { let arr = self.eval_with_env(base, env)?; let idx_vals: Vec = indices.iter().map(|e| { let v = self.eval_with_env(e, env)?; v.as_int().map(|n| n as usize).ok_or_else(|| EvalError::type_error("Int", &v)) }).collect::>()?; self.index_value(arr, &idx_vals) }
+ Expr::Index(base, indices) => { let arr = self.eval_with_env(base, env)?; let idx_vals: Vec = indices.iter().map(|e| { let v = self.eval_with_env(e, env)?; v.as_index().ok_or_else(|| EvalError::type_error("Int", &v)) }).collect::>()?; self.index_value(arr, &idx_vals) }
Expr::Slice { array, start, end } => { let arr = self.eval_with_env(array, env)?; let start_idx = match start { Some(e) => { let v = self.eval_with_env(e, env)?; v.as_int().map(|n| n as usize).unwrap_or(0) } None => 0 }; let end_idx = match end { Some(e) => { let v = self.eval_with_env(e, env)?; v.as_int().map(|n| n as usize) } None => None }; self.slice_value(arr, start_idx, end_idx) }
Expr::Annot(inner, _ty) => self.eval_with_env(inner, env),
Expr::Cast { expr, target: _, kind } => { let val = self.eval_with_env(expr, env)?; match kind { CastKind::Static => Ok(val), CastKind::Try => Ok(Value::variant("Some", Some(val))), CastKind::Force => Ok(val) } }
@@ -389,16 +389,20 @@ impl Evaluator {
fn bind_pattern(&self, pattern: &Pattern, val: Value, env: &mut Env) -> EvalResult<()> { if self.match_pattern(pattern, &val, env)? { Ok(()) } else { Err(EvalError::MatchFailed) } }
fn values_to_tensor(&self, values: Vec) -> Value {
+ self.values_to_tensor_shaped(vec![values.len()], values)
+ }
+
+ fn values_to_tensor_shaped(&self, shape: Vec, values: Vec) -> Value {
if values.is_empty() { return Value::Tensor(Tensor::from_ints(vec![])); }
let all_int = values.iter().all(|v| matches!(v, Value::Int(_)));
let all_float = values.iter().all(|v| matches!(v, Value::Float(_) | Value::Int(_)));
let all_bool = values.iter().all(|v| matches!(v, Value::Bool(_)));
let all_char = values.iter().all(|v| matches!(v, Value::Char(_)));
- if all_int { Value::Tensor(Tensor::from_ints(values.iter().map(|v| v.as_int().unwrap()).collect())) }
- else if all_float { Value::Tensor(Tensor::from_floats(values.iter().map(|v| v.coerce_float().unwrap()).collect())) }
- else if all_bool { Value::Tensor(Tensor::from_bools(values.iter().map(|v| v.as_bool().unwrap()).collect())) }
- else if all_char { let chars: Vec = values.iter().map(|v| v.as_char().unwrap()).collect(); let len = chars.len(); Value::Tensor(Tensor { shape: vec![len], data: crate::value::TensorData::Char(chars) }) }
- else { Value::Tensor(Tensor::from_values(vec![values.len()], values)) }
+ if all_int { Value::Tensor(Tensor { shape, data: crate::value::TensorData::Int(values.iter().map(|v| v.as_int().unwrap()).collect()) }) }
+ else if all_float { Value::Tensor(Tensor { shape, data: crate::value::TensorData::Float(values.iter().map(|v| ordered_float::OrderedFloat(v.coerce_float().unwrap())).collect()) }) }
+ else if all_bool { Value::Tensor(Tensor { shape, data: crate::value::TensorData::Bool(values.iter().map(|v| v.as_bool().unwrap()).collect()) }) }
+ else if all_char { Value::Tensor(Tensor { shape, data: crate::value::TensorData::Char(values.iter().map(|v| v.as_char().unwrap()).collect()) }) }
+ else { Value::Tensor(Tensor::from_values(shape, values)) }
}
fn access_field(&self, val: Value, access: &FieldAccess) -> EvalResult {
@@ -418,14 +422,14 @@ impl Evaluator {
fn slice_value(&self, val: Value, start: usize, end: Option) -> EvalResult {
match val {
- Value::Tensor(t) => { if t.rank() != 1 { return Err(EvalError::not_implemented("slicing rank > 1")); } let end = end.unwrap_or(t.len()); if start > end || end > t.len() { return Err(EvalError::IndexOutOfBounds { index: end, size: t.len() }); } let data: Vec = (start..end).map(|i| t.get_flat(i).unwrap()).collect(); Ok(Value::Tensor(Tensor::from_values(vec![data.len()], data))) }
+ Value::Tensor(t) => { if t.rank() != 1 { return Err(EvalError::not_implemented("slicing rank > 1")); } let end = end.unwrap_or(t.len()); if start > end || end > t.len() { return Err(EvalError::IndexOutOfBounds { index: end, size: t.len() }); } let data: Vec = (start..end).map(|i| t.get_flat(i).unwrap()).collect(); Ok(self.values_to_tensor_shaped(vec![data.len()], data)) }
_ => Err(EvalError::type_error("Tensor", &val)),
}
}
fn eval_map(&mut self, arr: Value, func: Value) -> EvalResult {
match arr {
- Value::Tensor(t) => { let results: Vec = t.iter().map(|elem| self.apply(func.clone(), elem)).collect::>()?; Ok(Value::Tensor(Tensor::from_values(t.shape.clone(), results))) }
+ Value::Tensor(t) => { let shape = t.shape.clone(); let results: Vec = t.iter().map(|elem| self.apply(func.clone(), elem)).collect::>()?; Ok(self.values_to_tensor_shaped(shape, results)) }
Value::Tuple(vs) => { let results: Vec = vs.into_iter().map(|elem| self.apply(func.clone(), elem)).collect::>()?; Ok(Value::Tuple(results)) }
_ => Err(EvalError::type_error("Tensor or Tuple", &arr)),
}
@@ -433,7 +437,7 @@ impl Evaluator {
fn eval_filter(&mut self, arr: Value, pred: Value) -> EvalResult {
match arr {
- Value::Tensor(t) => { let results: Vec = t.iter().filter_map(|elem| { let keep = self.apply(pred.clone(), elem.clone()).ok()?; match keep { Value::Bool(true) => Some(elem), _ => None } }).collect(); Ok(Value::Tensor(Tensor::from_values(vec![results.len()], results))) }
+ Value::Tensor(t) => { let results: Vec = t.iter().filter_map(|elem| { let keep = self.apply(pred.clone(), elem.clone()).ok()?; match keep { Value::Bool(true) => Some(elem), _ => None } }).collect(); Ok(self.values_to_tensor_shaped(vec![results.len()], results)) }
Value::Tuple(vs) => { let results: Vec = vs.into_iter().filter_map(|elem| { let keep = self.apply(pred.clone(), elem.clone()).ok()?; match keep { Value::Bool(true) => Some(elem), _ => None } }).collect(); Ok(Value::Tuple(results)) }
_ => Err(EvalError::type_error("Tensor or Tuple", &arr)),
}
@@ -441,7 +445,7 @@ impl Evaluator {
fn eval_bind(&mut self, arr: Value, func: Value) -> EvalResult {
match arr {
- Value::Tensor(t) => { let mut results = Vec::new(); for elem in t.iter() { let mapped = self.apply(func.clone(), elem)?; match mapped { Value::Tensor(inner) => results.extend(inner.iter()), Value::Tuple(inner) => results.extend(inner), other => results.push(other) } } Ok(Value::Tensor(Tensor::from_values(vec![results.len()], results))) }
+ Value::Tensor(t) => { let mut results = Vec::new(); for elem in t.iter() { let mapped = self.apply(func.clone(), elem)?; match mapped { Value::Tensor(inner) => results.extend(inner.iter()), Value::Tuple(inner) => results.extend(inner), other => results.push(other) } } Ok(self.values_to_tensor_shaped(vec![results.len()], results)) }
_ => Err(EvalError::type_error("Tensor", &arr)),
}
}
diff --git a/crates/goth-eval/src/value.rs b/crates/goth-eval/src/value.rs
index e772b2a..fd457a3 100644
--- a/crates/goth-eval/src/value.rs
+++ b/crates/goth-eval/src/value.rs
@@ -127,6 +127,8 @@ impl Value {
pub fn as_tensor(&self) -> Option<&Tensor> { match self { Value::Tensor(t) => Some(t), _ => None } }
pub fn as_tuple(&self) -> Option<&[Value]> { match self { Value::Tuple(vs) => Some(vs), Value::Unit => Some(&[]), _ => None } }
pub fn coerce_float(&self) -> Option { match self { Value::Float(f) => Some(f.0), Value::Int(n) => Some(*n as f64), _ => None } }
+ /// Coerce to a usize index: accepts Int directly, or Float if it's a whole number.
+ pub fn as_index(&self) -> Option { match self { Value::Int(n) => Some(*n as usize), Value::Float(f) => { let v = f.0; if v.fract() == 0.0 && v >= 0.0 { Some(v as usize) } else { None } }, _ => None } }
pub fn type_name(&self) -> &'static str {
match self {
diff --git a/crates/goth-parse/src/parser.rs b/crates/goth-parse/src/parser.rs
index 48f2a2c..a4a6383 100644
--- a/crates/goth-parse/src/parser.rs
+++ b/crates/goth-parse/src/parser.rs
@@ -31,11 +31,13 @@ pub type ParseResult = Result;
/// Parser
pub struct Parser<'a> {
lexer: Lexer<'a>,
+ /// End byte offset of the last consumed token (for adjacency checks)
+ last_token_end: usize,
}
impl<'a> Parser<'a> {
pub fn new(source: &'a str) -> Self {
- Parser { lexer: Lexer::new(source) }
+ Parser { lexer: Lexer::new(source), last_token_end: 0 }
}
// ============ Utilities ============
@@ -44,8 +46,17 @@ impl<'a> Parser<'a> {
self.lexer.peek()
}
+ /// Peek at the start byte offset of the next token (for adjacency checks)
+ fn peek_start(&mut self) -> Option {
+ self.lexer.peek_spanned().map(|s| s.loc.start)
+ }
+
fn next(&mut self) -> Option {
- self.lexer.next()
+ let spanned = self.lexer.next_spanned();
+ if let Some(ref s) = spanned {
+ self.last_token_end = s.loc.end;
+ }
+ spanned.map(|s| s.value)
}
fn expect(&mut self, expected: Token) -> ParseResult<()> {
@@ -379,6 +390,12 @@ impl<'a> Parser<'a> {
if !can_index {
break; // Treat [..] as a separate array, not indexing
}
+ // Whitespace-sensitive: `expr[i]` (adjacent) = indexing,
+ // `expr [1,2]` (space) = function application with array arg.
+ let adjacent = self.peek_start().map_or(false, |s| s == self.last_token_end);
+ if !adjacent {
+ break; // Treat [..] as a separate array expression
+ }
self.next();
let mut indices = vec![self.parse_expr()?];
while self.eat(&Token::Comma) {
diff --git a/docs/BENCHMARK-TESTS.md b/docs/BENCHMARK-TESTS.md
index 24b4b80..5b3f959 100644
--- a/docs/BENCHMARK-TESTS.md
+++ b/docs/BENCHMARK-TESTS.md
@@ -87,6 +87,19 @@ Classic computer science algorithms.
| `isqrt` | ⌊√n⌋ | `50 → 7` | integer sqrt |
| `modpow` | b^e mod m | `2 10 1000 → 24` | fast exponentiation |
+### 6. Uncertainty (6 tests) - `examples/uncertainty/`
+
+Uncertain value creation and error propagation through operations.
+
+| Test | Description | Input → Output | Contract |
+|------|-------------|----------------|----------|
+| `measure` | Create uncertain value | `10.5 0.3 → 10.5±0.3` | value ± uncertainty |
+| `add_uncertain` | Additive propagation | `10 0.3 20 0.4 → 30±0.5` | δ = √(δa² + δb²) |
+| `mul_uncertain` | Multiplicative propagation | `5 0.1 3 0.2 → 15±1.04` | relative error quadrature |
+| `sqrt_uncertain` | Sqrt propagation | `9 0.3 → 3±0.05` | δ = δx / (2√x) |
+| `sin_uncertain` | Trig propagation | `1.0 0.1 → 0.841±0.054` | δ = \|cos x\| × δx |
+| `chained_uncertain` | Multi-step propagation | `4 0.2 1 0.1 → 0.141±0.111` | sin(√(a±δa)+(b±δb)) |
+
## Running Tests
```bash
@@ -142,7 +155,8 @@ For multi-argument functions, indices count from most-recent binding:
| Higher-Order | 10 | ✓ |
| Numeric | 8 | ✓ |
| Algorithms | 6 | ✓ |
-| **Total** | **48** | |
+| Uncertainty | 6 | ✓ |
+| **Total** | **54** | |
## Future Additions
diff --git a/docs/EFFECT-SYSTEM-ROADMAP.md b/docs/EFFECT-SYSTEM-ROADMAP.md
new file mode 100644
index 0000000..e225534
--- /dev/null
+++ b/docs/EFFECT-SYSTEM-ROADMAP.md
@@ -0,0 +1,99 @@
+# Effect System Roadmap
+
+**Status:** Aspirational — parsed but not enforced (as of v0.1.0)
+
+## Current State
+
+Goth's effect system has complete AST infrastructure but no enforcement. Effect annotations are parsed into the syntax tree and serialized, but ignored by the type checker, evaluator, and codegen backends.
+
+### What Works
+
+- **Lexer/Parser**: The `◇` (Diamond) token and effect names (`io`, `mut`, `rand`, `div`) are recognized in function box middle lines (`│ ◇io`)
+- **AST**: `Effect` enum (Pure, Io, Mut, Rand, Div, Exn, Ffi, Custom) with `Effects` set algebra (union, subset, containment) in `goth-ast/src/effect.rs`
+- **Type system**: `Type::Effectful(Box, Effects)` variant exists in `goth-ast/src/types.rs`
+- **Function declarations**: `FnDecl.effects: Effects` field stores parsed annotations in `goth-ast/src/decl.rs`
+- **Serialization**: Effects serialize correctly via serde
+
+### What Doesn't Work
+
+- **Type checker**: Primitives typed without effects; no effect inference or checking (`goth-check/src/infer.rs`, `builtins.rs`)
+- **Evaluator**: I/O primitives dispatch unconditionally; `Closure` doesn't store effects; `EffectNotAllowed` error defined but never raised (`goth-eval/src/prim.rs`, `value.rs`)
+- **Parser limitations**: `◇` only works in function box middle lines, not in type signatures
+- **Codegen**: Both MLIR and LLVM backends explicitly erase effects during lowering
+
+## Implementation Plan
+
+### Phase 1: Type Checker (core enforcement)
+
+**1.1 Annotate primitive types with effects** (~50 lines)
+- File: `goth-check/src/builtins.rs`
+- Wrap I/O functions (`print`, `readLine`, `readFile`, `writeFile`, `write`, `flush`, `readKey`, `rawModeEnter`, `rawModeExit`, `sleep`) with `Type::Effectful(..., Effects::single(Effect::Io))`
+
+**1.2 Add effect context to type checker** (~100 lines)
+- File: `goth-check/src/context.rs` (or equivalent)
+- Add `allowed_effects: Effects` to `Context` struct
+- Methods: `with_effect_context()`, `check_effect_allowed()`
+
+**1.3 Effect inference** (~800 lines)
+- File: `goth-check/src/infer.rs`
+- Modify `infer()` to return `(Type, Effects)` pairs
+- Accumulate effects through: function application, let bindings, match arms, if/else branches, binary operators (Write/Read)
+- Pure functions: reject effectful sub-computations
+- Effectful functions: propagate effects to return type
+
+### Phase 2: Parser Enhancement (~200 lines)
+
+- Support effect annotations in type signatures: `(F64 → F64) ◇io`
+- Handle effect composition in type syntax: `String ◇io ◇mut`
+
+### Phase 3: Runtime Support (~130 lines)
+
+**3.1 Store effects on closures** (~30 lines)
+- File: `goth-eval/src/value.rs`
+- Add `effects: Effects` field to `Closure`
+
+**3.2 Runtime effect checking** (~100 lines)
+- File: `goth-eval/src/eval.rs`
+- Check allowed effects before dispatching I/O primitives
+- Raise `EffectNotAllowed` when violated
+
+### Phase 4: Codegen (optional)
+
+Effects are a static property and can remain erased at codegen. Optionally:
+- Emit MLIR function attributes for effects
+- Preserve effects through MIR lowering for tooling/analysis
+
+### Phase 5: Testing (~300 lines)
+
+- Pure function calling `print` → type error
+- Effectful function calling `print` → allowed
+- Effect inference through let/match/if chains
+- Effect union across branches
+- Polymorphic effects (advanced): `∀ε. (α →ε β) → α →ε β`
+
+## Key Files
+
+| File | Role | Work Needed |
+|------|------|-------------|
+| `goth-ast/src/effect.rs` | Effect data model | Complete — no changes needed |
+| `goth-ast/src/types.rs` | `Type::Effectful` | Complete — no changes needed |
+| `goth-ast/src/decl.rs` | `FnDecl.effects` | Complete — no changes needed |
+| `goth-check/src/builtins.rs` | Primitive types | Add `Type::Effectful` wrappers |
+| `goth-check/src/infer.rs` | Type inference | Add effect inference throughout |
+| `goth-check/src/check.rs` | Type checking | Add effect validation |
+| `goth-eval/src/value.rs` | Closure definition | Add `effects` field |
+| `goth-eval/src/eval.rs` | Evaluation | Add runtime effect checks |
+| `goth-eval/src/prim.rs` | Primitive dispatch | Guard I/O behind effect checks |
+| `goth-parse/src/parser.rs` | Parsing | Extend effect syntax to types |
+
+## Design Decisions (pending)
+
+1. **Strictness**: Should missing `◇io` on a function that calls `print` be a hard error or a warning?
+2. **Inference vs annotation**: Should effects be inferred automatically, or must they be declared?
+3. **Polymorphic effects**: Support `∀ε` effect variables? (Significant complexity increase)
+4. **Effect handlers**: Algebraic effect handlers (like Koka/Eff) or just tracking/checking?
+5. **Backwards compatibility**: How to handle existing `.goth` files that use `print` without `◇io`?
+
+## Estimated Scope
+
+~1,200 lines of new/modified code for core enforcement (Phases 1 + 3 + 5). Parser enhancement (Phase 2) adds ~200 lines. Codegen (Phase 4) is optional.
diff --git a/docs/GOTH-LANGUAGE-REFERENCE-v0.1.md b/docs/GOTH-LANGUAGE-REFERENCE-v0.1.md
index 121c387..ec55a43 100644
--- a/docs/GOTH-LANGUAGE-REFERENCE-v0.1.md
+++ b/docs/GOTH-LANGUAGE-REFERENCE-v0.1.md
@@ -496,6 +496,8 @@ let arr = [10, 20, 30, 40] in arr[2]
# Result: 30 (0-indexed)
```
+**Note:** The `[` must be directly adjacent (no space) to be parsed as indexing. With a space, `f [1,2]` is function application, passing the array `[1,2]` as an argument to `f`.
+
**Multi-dimensional indexing:**
```goth
let matrix = [[1, 2, 3], [4, 5, 6]] in
@@ -856,6 +858,36 @@ F ± F # Float with uncertainty
I ± I # Int with uncertainty
```
+**Creating uncertain values at runtime:**
+```goth
+10.5 ± 0.3 # Value 10.5 with uncertainty 0.3
+```
+
+**Automatic uncertainty propagation:**
+
+When uncertain values flow through arithmetic operators and math functions, uncertainty propagates automatically using standard error propagation rules:
+
+| Operation | Propagation Rule |
+|-----------|-----------------|
+| `(a±δa) + (b±δb)` | δ = √(δa² + δb²) |
+| `(a±δa) - (b±δb)` | δ = √(δa² + δb²) |
+| `(a±δa) × (b±δb)` | δ = \|a×b\| × √((δa/a)² + (δb/b)²) |
+| `(a±δa) / (b±δb)` | δ = \|a/b\| × √((δa/a)² + (δb/b)²) |
+| `√(x±δx)` | δ = δx / (2√x) |
+| `sin(x±δx)` | δ = \|cos(x)\| × δx |
+| `cos(x±δx)` | δ = \|sin(x)\| × δx |
+| `exp(x±δx)` | δ = exp(x) × δx |
+| `ln(x±δx)` | δ = δx / \|x\| |
+
+**Supported functions:** `+`, `-`, `×`, `/`, `^`, `√`, `exp`, `ln`, `log10`, `log2`, `sin`, `cos`, `tan`, `asin`, `acos`, `atan`, `sinh`, `cosh`, `tanh`, `abs`, `floor`, `ceil`, `round`, `Γ`.
+
+**Example — chained propagation:**
+```goth
+╭─ main : F64 → F64 → F64 → F64 → (F64 ± F64)
+╰─ sin (√(₃ ± ₂) + (₁ ± ₀))
+# With inputs 4.0 0.2 1.0 0.1 → 0.1411±0.1107
+```
+
### Refinement Types
**Constrained types:**
@@ -872,6 +904,8 @@ I ± I # Int with uncertainty
### Effect Types
+> **Aspirational:** Effect annotations are parsed and stored in the AST but not enforced by the type checker or evaluator. They currently serve as documentation. See `docs/EFFECT-SYSTEM-ROADMAP.md`.
+
```goth
□ # Pure (no effects)
◇io # I/O effects
diff --git a/docs/tests.md b/docs/tests.md
index 76fc2f4..a79a8df 100644
--- a/docs/tests.md
+++ b/docs/tests.md
@@ -53,6 +53,13 @@ else 0
# Uncertain types (parse test with :ast)
:ast let temp : F64 ± F64 = ⟨μ: 20.0, σ: 0.5⟩
+# Uncertainty propagation (runtime tests)
+10.5 ± 0.3 # → 10.5±0.3
+(10.0 ± 0.3) + (20.0 ± 0.4) # → 30±0.5
+(5.0 ± 0.1) × (3.0 ± 0.2) # → 15±1.044...
+√(9.0 ± 0.3) # → 3±0.05
+sin (1.0 ± 0.1) # → 0.841...±0.054...
+
# Refinement types (parse test)
:ast let positive : {x : F64 | x > 0} = 5.0
diff --git a/examples/README.md b/examples/README.md
new file mode 100644
index 0000000..6effdc0
--- /dev/null
+++ b/examples/README.md
@@ -0,0 +1,260 @@
+# Goth Examples
+
+A collection of examples demonstrating the Goth programming language, organized by concept. All examples run under the eval interpreter.
+
+## Running Examples
+
+From the `crates/` directory:
+
+```sh
+# Run an example with an integer argument
+cargo run --quiet --package goth-cli -- ../examples/basic/square.goth 7
+# Output: 49
+
+# Run with no argument (for () → T signatures)
+cargo run --quiet --package goth-cli -- ../examples/numeric/product_range.goth
+
+# Evaluate an inline expression
+cargo run --quiet --package goth-cli -- -e 'Σ [1, 2, 3, 4, 5]'
+# Output: 15
+```
+
+Most examples define `main : I64 → ...` and expect a single integer argument. Multi-argument functions return a partial application when given fewer arguments than needed — pass additional arguments as space-separated values.
+
+```sh
+# Two-argument function
+cargo run --quiet --package goth-cli -- ../examples/recursion/gcd.goth 12 8
+# Output: 4
+```
+
+### Running All Examples
+
+```sh
+cd crates
+for f in ../examples/basic/*.goth; do
+ echo "$(basename $f): $(cargo run --quiet --package goth-cli -- "$f" 5 2>&1)"
+done
+```
+
+---
+
+## basic/
+
+Elementary single-function programs. Each takes one integer and returns a result.
+
+| File | Description | Example |
+|------|-------------|---------|
+| `identity.goth` | Returns its argument unchanged | `5 → 5` |
+| `add_one.goth` | Adds one | `5 → 6` |
+| `double.goth` | Doubles the input | `5 → 10` |
+| `square.goth` | Squares the input | `5 → 25` |
+| `abs.goth` | Absolute value | `-3 → 3` |
+| `sign.goth` | Sign function (−1, 0, or 1) | `5 → 1` |
+| `is_even.goth` | Evenness check | `4 → ⊤` |
+| `is_positive.goth` | Positivity check | `5 → ⊤` |
+| `max_two.goth` | Maximum of two integers (curried) | `3 7 → 7` |
+| `min_two.goth` | Minimum of two integers (curried) | `3 7 → 3` |
+
+**Demonstrates:** Function declarations, De Bruijn indices, conditionals, boolean returns.
+
+---
+
+## recursion/
+
+Classic recursive algorithms. Demonstrates self-referencing function calls, base cases, and termination.
+
+| File | Description | Example |
+|------|-------------|---------|
+| `factorial.goth` | n! | `5 → 120` |
+| `fibonacci.goth` | Fibonacci numbers | `10 → 55` |
+| `sum_to_n.goth` | Triangular numbers: 1 + 2 + ... + n | `5 → 15` |
+| `power.goth` | base^exp (two args) | `2 10 → 1024` |
+| `gcd.goth` | Euclidean GCD (two args) | `12 8 → 4` |
+| `lcm.goth` | Least common multiple (two args) | `4 6 → 12` |
+| `ackermann.goth` | Ackermann function (two args) — non-primitive recursive | `2 3 → 9` |
+| `sudan.goth` | Sudan function (three args) — non-primitive recursive | `1 1 1 → 4` |
+| `collatz_len.goth` | Length of Collatz (3n+1) sequence | `7 → 16` |
+| `digit_sum.goth` | Sum of decimal digits | `1234 → 10` |
+| `reverse_num.goth` | Reverse digits of a number | `1234 → 4321` |
+| `hyperop.goth` | Hyperoperation / Knuth's up-arrow (three args) | `2 3 4 → 65536` |
+| `tak.goth` | Takeuchi function (three args) — classic call-overhead benchmark | `7 4 2 → 4` |
+| `mccarthy91.goth` | McCarthy 91 function — nested recursion | `85 → 91` |
+
+**Demonstrates:** Recursive self-calls, multi-argument functions, De Bruijn index ordering in curried functions.
+
+---
+
+## higher-order/
+
+Functions that take or return other functions. Demonstrates map, filter, fold, and function composition.
+
+| File | Description | Example |
+|------|-------------|---------|
+| `map_double.goth` | Map: double each element of a range | `5 → [2 4 6 8 10]` |
+| `filter_positive.goth` | Filter: keep only positive numbers | `5 → [1 2 3 4 5]` |
+| `compose.goth` | Function composition: square(double(x)) | `5 → 100` |
+| `apply_twice.goth` | Apply a function twice: f(f(x)) | `5 → 20` |
+| `all_positive.goth` | Check if all elements satisfy a predicate | `5 → ⊤` |
+| `any_negative.goth` | Check if any element satisfies a predicate | `5 → ⊤` |
+| `count_if.goth` | Count elements matching a predicate | `5 → 2` |
+| `fold_sum.goth` | Fold/reduce to compute sum | `5 → 15` |
+| `fold_product.goth` | Fold/reduce to compute product (factorial) | `5 → 120` |
+| `pipeline.goth` | Chained operations: range → filter → map → sum | `6 → 56` |
+| `concat.goth` | Array concatenation with `⊕` | `5 → [1 2 3 4 5 1 2 3 4 5]` |
+| `zip_sum.goth` | Zip two arrays into pairs | `5 → [⟨1,0⟩ ⟨2,1⟩ ...]` |
+
+**Demonstrates:** `↦` (map), `▸` (filter), `Σ`/`Π` (reduce), lambda expressions, closures, `⊕` (concat).
+
+---
+
+## numeric/
+
+Numerical computing: series, special functions, and iterative methods.
+
+| File | Description | Example |
+|------|-------------|---------|
+| `sum_squares.goth` | Sum of squares: 1² + 2² + ... + n² | `5 → 55` |
+| `product_range.goth` | Product of 1..n via `Π` (factorial) | `5 → 120` |
+| `harmonic.goth` | Harmonic number H(n) = 1 + 1/2 + ... + 1/n | `5 → 2.28...` |
+| `pi_leibniz.goth` | π via Leibniz series (slow convergence) | `20 → 2.976...` |
+| `exp_taylor.goth` | e^x via Taylor series | `5.0 → 148.41...` |
+| `sqrt_newton.goth` | √x via Newton-Raphson iteration | `5.0 → 2.236...` |
+| `gamma_fact.goth` | Factorial via the Gamma function: Γ(n+1) = n! | `5.0 → 120.0` |
+| `gamma_half.goth` | Gamma at half-integers: Γ(n+½) | `5.0 → 24.0` |
+
+**Demonstrates:** `Σ`, `Π`, `ι` (iota), `↦` (map), floating-point arithmetic, built-in math functions (`Γ`, `√`, `exp`, `ln`).
+
+---
+
+## algorithms/
+
+Classical algorithms implemented in a functional style.
+
+| File | Description | Example |
+|------|-------------|---------|
+| `isPrime.goth` | Primality test by trial division | `7 → ⊤` |
+| `count_primes.goth` | Count primes ≤ n (prime counting function π(n)) | `10 → 4` |
+| `nth_prime.goth` | Find the nth prime number | `5 → 11` |
+| `isqrt.goth` | Integer square root ⌊√n⌋ | `10 → 3` |
+| `binary_search.goth` | Binary search in a range (three args) | `7 0 10 → 7` |
+| `modpow.goth` | Modular exponentiation base^exp mod m (three args) | `2 10 1000 → 24` |
+
+**Demonstrates:** Multi-argument curried functions, conditional recursion, helper functions.
+
+---
+
+## contracts/
+
+Functions with runtime-checked preconditions (`⊢`) and postconditions (`⊨`).
+
+| File | Description | Contract |
+|------|-------------|----------|
+| `abs_post.goth` | Absolute value | Postcondition: result ≥ 0 |
+| `div_safe.goth` | Safe division (two args) | Precondition: divisor ≠ 0 |
+| `factorial_contract.goth` | Factorial | Pre: n ≥ 0, Post: result ≥ 1 |
+| `log_safe.goth` | Natural logarithm | Precondition: x > 0 |
+| `sqrt_safe.goth` | Square root | Precondition: x ≥ 0 |
+
+**Demonstrates:** `⊢` (preconditions), `⊨` (postconditions), `│` (contract lines in function boxes).
+
+---
+
+## tco/
+
+Paired naive vs. tail-call-optimized versions of the same algorithm. Each pair shows the transformation from non-tail-recursive to accumulator-based tail recursion.
+
+| Pair | Naive | TCO |
+|------|-------|-----|
+| Factorial | `factorial_naive.goth` | `factorial_tco.goth` |
+| Fibonacci | `fibonacci_naive.goth` | `fibonacci_tco.goth` |
+| Sum 1..n | `sum_naive.goth` | `sum_tco.goth` |
+| Collatz length | `collatz_naive.goth` | `collatz_tco.goth` |
+| List length | `length_naive.goth` | `length_tco.goth` |
+
+Both versions produce identical results. The TCO versions use accumulator parameters and place the recursive call in tail position.
+
+**Demonstrates:** Tail recursion, accumulator pattern, helper functions with extra parameters.
+
+---
+
+## io/
+
+File and stream I/O using the `▷` (write) operator.
+
+| File | Description |
+|------|-------------|
+| `write_stdout.goth` | Write `"hello"` to stdout (no newline) |
+| `write_stderr.goth` | Write `"error"` to stderr |
+| `write_file.goth` | Write `"hello world"` to `/tmp/goth_write_test.txt` |
+
+**Demonstrates:** `▷ stdout`, `▷ stderr`, `▷ "filepath"` — the write operator dispatches on the right-hand side.
+
+---
+
+## uncertainty/
+
+First-class uncertain values with automatic error propagation.
+
+| File | Description | Example |
+|------|-------------|---------|
+| `measure.goth` | Create an uncertain value | `10.0 0.5 → 10 ± 0.5` |
+| `add_uncertain.goth` | Addition with error propagation | `10 0.3 20 0.4 → 30 ± 0.5` |
+| `mul_uncertain.goth` | Multiplication with relative error | `5 0.1 3 0.2 → 15 ± 1.04...` |
+| `sqrt_uncertain.goth` | Square root with derivative propagation | `9 0.3 → 3 ± 0.05` |
+| `sin_uncertain.goth` | Sine with derivative propagation | `1.0 0.1 → 0.841... ± 0.054...` |
+| `chained_uncertain.goth` | Multi-step propagation: sin(√a + b) | `4 0.1 1 0.05 → ...` |
+
+**Demonstrates:** The `±` operator, automatic propagation through `+`, `-`, `×`, `/`, `√`, `sin`, `cos`, `exp`, `ln`, etc.
+
+---
+
+## simulation/
+
+Numerical simulations that produce SVG visualizations and CSV data files. These examples write output files to the working directory (typically `crates/`).
+
+| File | Description | Output Files |
+|------|-------------|-------------|
+| `heat1d.goth` | 1D heat diffusion (20 pts, 50 steps, r=0.4) | `heat1d.svg` |
+| `heat2d.goth` | 2D heat diffusion (15×15 grid, 30 steps, r=0.2) | `heat2d.svg` |
+| `heat1d_tunable.goth` | 1D heat with imported config | `heat1d.svg` |
+| `heat2d_tunable.goth` | 2D heat with imported config | `heat2d.svg` |
+| `newton_raphson.goth` | Newton-Raphson cube root finder | `newton.svg`, `newton.csv` |
+| `wave1d.goth` | 1D wave equation (leapfrog scheme) | `wave1d.svg`, `wave1d.csv` |
+| `laplace2d.goth` | 2D Laplace equation (Jacobi iteration) | `laplace2d.svg`, `laplace2d.csv` |
+| `power_iter.goth` | Eigenvalue power iteration | `power_iter.svg`, `power_iter.csv` |
+
+Config files (`heat1d_config.goth`, `heat2d_config.goth`) are imported via `use "file.goth"` and contain bare `let` bindings for tunable parameters.
+
+```sh
+# Run a simulation (argument is unused, pass 0)
+cargo run --quiet --package goth-cli -- ../examples/simulation/heat2d.goth 0
+
+# Open the SVG output
+open heat2d.svg
+```
+
+**Demonstrates:** Recursive time-stepping, flat 2D array indexing, string building for SVG/CSV, `writeFile`/`▷` for file output, `use` imports for configuration, `toString`, `strConcat`/`⧺`.
+
+---
+
+## Language Features Covered
+
+| Feature | Examples |
+|---------|---------|
+| De Bruijn indices (`₀`, `₁`, `₂`) | All |
+| Function declarations (`╭─`/`╰─`) | All |
+| Lambda expressions (`λ→`) | higher-order, numeric, simulation |
+| Let bindings | numeric, simulation, algorithms |
+| Recursion | recursion, tco, algorithms, simulation |
+| Map (`↦`) | higher-order, numeric, simulation |
+| Filter (`▸`) | higher-order |
+| Sum/Product (`Σ`/`Π`) | numeric, higher-order, simulation |
+| Iota (`ι`) | numeric, simulation |
+| Array concat (`⊕`) | higher-order |
+| String concat (`⧺`) | simulation |
+| Conditionals (`if/then/else`) | basic, recursion, algorithms |
+| Contracts (`⊢`/`⊨`) | contracts |
+| Uncertainty (`±`) | uncertainty |
+| File I/O (`▷`, `writeFile`) | io, simulation |
+| Module imports (`use`) | simulation (tunable variants) |
+| Tail-call optimization | tco |
diff --git a/examples/higher-order/all_positive.goth b/examples/higher-order/all_positive.goth
index 7e028a0..ccc89d0 100644
--- a/examples/higher-order/all_positive.goth
+++ b/examples/higher-order/all_positive.goth
@@ -8,4 +8,5 @@
╰─ ₀ > 0
╭─ main : I64 → Bool
-╰─ all isPositive (range 1 (₀ + 1))
+╰─ let arr = range 1 (₀ + 1) in
+ Σ (arr ↦ (λ→ toInt (isPositive ₀))) == len arr
diff --git a/examples/higher-order/any_negative.goth b/examples/higher-order/any_negative.goth
index fb3c67f..be8245d 100644
--- a/examples/higher-order/any_negative.goth
+++ b/examples/higher-order/any_negative.goth
@@ -8,4 +8,4 @@
╰─ ₀ < 0
╭─ main : I64 → Bool
-╰─ any isNegative (range (0 - ₀) (₀ + 1))
+╰─ Σ (range (0 - ₀) (₀ + 1) ↦ (λ→ toInt (isNegative ₀))) > 0
diff --git a/examples/higher-order/map_double.goth b/examples/higher-order/map_double.goth
index d8716ee..9d71821 100644
--- a/examples/higher-order/map_double.goth
+++ b/examples/higher-order/map_double.goth
@@ -10,4 +10,4 @@
# Return doubled list for range [1,n]
╭─ main : I64 → T
-╰─ (range 1 (₀ + 1)) ▷ double
+╰─ (range 1 (₀ + 1)) ↦ double
diff --git a/examples/simulation/heat1d.goth b/examples/simulation/heat1d.goth
new file mode 100644
index 0000000..f20dbf4
--- /dev/null
+++ b/examples/simulation/heat1d.goth
@@ -0,0 +1,68 @@
+# 1D Heat Diffusion with SVG Visualization
+# Explicit finite difference: u[i]' = u[i] + r*(u[i+1] - 2*u[i] + u[i-1])
+# Grid: 20 points, 50 timesteps, r = 0.4
+# Output: heat1d.svg — row of colored rectangles (blue→red)
+
+# Clamp value to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Temperature → red channel (0–255)
+╭─ colorR : F64 → I64
+╰─ toInt (clamp ₀ × 255.0)
+
+# Temperature → blue channel (0–255)
+╭─ colorB : F64 → I64
+╰─ toInt ((1.0 - clamp ₀) × 255.0)
+
+# Build "rgb(R,0,B)" color string from temperature
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat ",0," (strConcat (toString (colorB ₀)) ")")))
+
+# Initial condition: hot center (cells 8–11), cold elsewhere
+╭─ initCell : I64 → F64
+╰─ if ₀ > 7 then
+ if ₀ < 12 then 1.0 else 0.0
+ else 0.0
+
+# Compute new value at index i using neighbors
+# ₁ = grid, ₀ = index i; r = 0.4
+╭─ stencil : [n]F64 → I64 → F64
+╰─ if ₀ = 0 then 0.0
+ else if ₀ = 19 then 0.0
+ else ₁[₀] + 0.4 × (₁[₀ - 1] + ₁[₀ + 1] - 2.0 × ₁[₀])
+
+# One timestep: map stencil over all grid indices
+╭─ stepGrid : [n]F64 → [n]F64
+╰─ let grid = ₀ in
+ ι 20 ↦ (λ→ stencil grid ₀)
+
+# Recursive time-stepping: evolve grid for t steps
+# ₁ = grid, ₀ = remaining steps
+╭─ evolve : [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₁
+ else evolve (stepGrid ₁) (₀ - 1)
+
+# SVG rectangle for one cell
+# ₁ = cell index, ₀ = temperature value
+╭─ cellRect : I64 → F64 → String
+╰─ strConcat "\n")))
+
+# Recursively build SVG rect strings for cells [i..20)
+# ₂ = grid, ₁ = current index i, ₀ = accumulator string
+╭─ buildRects : [n]F64 → I64 → String → String
+╰─ if ₁ = 20 then ₀
+ else buildRects ₂ (₁ + 1) (strConcat ₀ (cellRect ₁ ₂[₁]))
+
+# Build complete SVG from final grid
+╭─ buildSvg : [n]F64 → String
+╰─ let rects = buildRects ₀ 0 "" in
+ strConcat "\n")
+
+# Main: run simulation, write SVG file
+╭─ main : I64 → ()
+╰─ let grid0 = ι 20 ↦ (λ→ initCell ₀) in
+ let final = evolve grid0 50 in
+ buildSvg final ▷ "heat1d.svg"
diff --git a/examples/simulation/heat1d_config.goth b/examples/simulation/heat1d_config.goth
new file mode 100644
index 0000000..ad16f54
--- /dev/null
+++ b/examples/simulation/heat1d_config.goth
@@ -0,0 +1,15 @@
+# 1D Heat Diffusion Configuration
+# Edit these values to tune the simulation
+
+# Grid points along the bar
+let nx = 20
+
+# Number of timesteps
+let nt = 50
+
+# Stability ratio r = α·Δt/Δx² (must be < 0.5 for stability)
+let r = 0.4
+
+# Initial hot region: cells from icLo to icHi (exclusive) are 1.0
+let icLo = 8
+let icHi = 12
diff --git a/examples/simulation/heat1d_tunable.goth b/examples/simulation/heat1d_tunable.goth
new file mode 100644
index 0000000..310e0df
--- /dev/null
+++ b/examples/simulation/heat1d_tunable.goth
@@ -0,0 +1,70 @@
+# 1D Heat Diffusion with Tunable Parameters
+# Reads configuration from heat1d_config.goth via use import
+# Edit that file to change grid size, timesteps, diffusion rate, and IC
+# Output: heat1d.svg — row of colored rectangles (blue→red)
+
+use "heat1d_config.goth"
+
+# Clamp value to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Temperature → red channel (0–255)
+╭─ colorR : F64 → I64
+╰─ toInt (clamp ₀ × 255.0)
+
+# Temperature → blue channel (0–255)
+╭─ colorB : F64 → I64
+╰─ toInt ((1.0 - clamp ₀) × 255.0)
+
+# Build "rgb(R,0,B)" color string from temperature
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat ",0," (strConcat (toString (colorB ₀)) ")")))
+
+# Initial condition: hot region [icLo, icHi), cold elsewhere
+╭─ initCell : I64 → F64
+╰─ if ₀ < icLo then 0.0
+ else if ₀ < icHi then 1.0
+ else 0.0
+
+# Compute new value at index i using neighbors
+# ₁ = grid, ₀ = index i
+╭─ stencil : [n]F64 → I64 → F64
+╰─ if ₀ = 0 then 0.0
+ else if ₀ = (nx - 1) then 0.0
+ else ₁[₀] + r × (₁[₀ - 1] + ₁[₀ + 1] - 2.0 × ₁[₀])
+
+# One timestep: map stencil over all grid indices
+╭─ stepGrid : [n]F64 → [n]F64
+╰─ let grid = ₀ in
+ ι nx ↦ (λ→ stencil grid ₀)
+
+# Recursive time-stepping: evolve grid for t steps
+# ₁ = grid, ₀ = remaining steps
+╭─ evolve : [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₁
+ else evolve (stepGrid ₁) (₀ - 1)
+
+# SVG rectangle for one cell
+# ₁ = cell index, ₀ = temperature value
+╭─ cellRect : I64 → F64 → String
+╰─ strConcat "\n")))
+
+# Recursively build SVG rect strings for cells [i..nx)
+# ₂ = grid, ₁ = current index i, ₀ = accumulator string
+╭─ buildRects : [n]F64 → I64 → String → String
+╰─ if ₁ = nx then ₀
+ else buildRects ₂ (₁ + 1) (strConcat ₀ (cellRect ₁ ₂[₁]))
+
+# Build complete SVG from final grid
+╭─ buildSvg : [n]F64 → String
+╰─ let rects = buildRects ₀ 0 "" in
+ strConcat "\n")))
+
+# Main: run simulation, write SVG file
+╭─ main : I64 → ()
+╰─ let grid0 = ι nx ↦ (λ→ initCell ₀) in
+ let final = evolve grid0 nt in
+ buildSvg final ▷ "heat1d.svg"
diff --git a/examples/simulation/heat2d.goth b/examples/simulation/heat2d.goth
new file mode 100644
index 0000000..b219fe5
--- /dev/null
+++ b/examples/simulation/heat2d.goth
@@ -0,0 +1,95 @@
+# 2D Heat Diffusion with SVG Visualization
+# 5-point stencil: u' = u + r*(u_N + u_S + u_E + u_W - 4*u)
+# Grid: 15×15 = 225 cells (flat), 30 timesteps, r = 0.2
+# Output: heat2d.svg — 15×15 heatmap grid
+
+# Clamp value to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Temperature → red channel (0–255)
+╭─ colorR : F64 → I64
+╰─ toInt (clamp ₀ × 255.0)
+
+# Temperature → blue channel (0–255)
+╭─ colorB : F64 → I64
+╰─ toInt ((1.0 - clamp ₀) × 255.0)
+
+# Build "rgb(R,0,B)" color string from temperature
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat ",0," (strConcat (toString (colorB ₀)) ")")))
+
+# Initial condition: hot 5×5 center block (rows 5–9, cols 5–9)
+# ₀ = flat index
+╭─ initCell : I64 → F64
+╰─ let row = ₀ / 15 in
+ let col = ₁ % 15 in
+ if row > 4 then
+ if row < 10 then
+ if col > 4 then
+ if col < 10 then 1.0 else 0.0
+ else 0.0
+ else 0.0
+ else 0.0
+
+# Boundary-aware grid access: returns 0.0 for out-of-bounds
+# ₂ = grid, ₁ = row, ₀ = col
+╭─ getVal : [n]F64 → I64 → I64 → F64
+╰─ if ₁ < 0 then 0.0
+ else if ₁ > 14 then 0.0
+ else if ₀ < 0 then 0.0
+ else if ₀ > 14 then 0.0
+ else ₂[₁ × 15 + ₀]
+
+# 5-point stencil at flat index k
+# ₁ = grid, ₀ = flat index k
+╭─ stencil2d : [n]F64 → I64 → F64
+╰─ let grid = ₁ in
+ let k = ₁ in
+ let row = k / 15 in
+ let col = k % 15 in
+ let center = grid[k] in
+ let north = getVal grid (row - 1) col in
+ let south = getVal grid (row + 1) col in
+ let west = getVal grid row (col - 1) in
+ let east = getVal grid row (col + 1) in
+ center + 0.2 × (north + south + east + west - 4.0 × center)
+
+# One timestep: map stencil over all grid cells
+╭─ stepGrid : [n]F64 → [n]F64
+╰─ let grid = ₀ in
+ ι 225 ↦ (λ→ stencil2d grid ₀)
+
+# Recursive time-stepping: evolve grid for t steps
+# ₁ = grid, ₀ = remaining steps
+╭─ evolve : [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₁
+ else evolve (stepGrid ₁) (₀ - 1)
+
+# SVG rectangle for one cell at flat index k
+# ₁ = flat index, ₀ = temperature value
+╭─ cellRect : I64 → F64 → String
+╰─ let k = ₁ in
+ let temp = ₁ in
+ let row = k / 15 in
+ let col = k % 15 in
+ strConcat "\n")))))
+
+# Recursively build SVG rect strings for cells [i..225)
+# ₂ = grid, ₁ = current index, ₀ = accumulator string
+╭─ buildRects : [n]F64 → I64 → String → String
+╰─ if ₁ = 225 then ₀
+ else buildRects ₂ (₁ + 1) (strConcat ₀ (cellRect ₁ ₂[₁]))
+
+# Build complete SVG from final grid
+╭─ buildSvg : [n]F64 → String
+╰─ let rects = buildRects ₀ 0 "" in
+ strConcat "\n")
+
+# Main: run simulation, write SVG file
+╭─ main : I64 → ()
+╰─ let grid0 = ι 225 ↦ (λ→ initCell ₀) in
+ let final = evolve grid0 30 in
+ buildSvg final ▷ "heat2d.svg"
diff --git a/examples/simulation/heat2d_config.goth b/examples/simulation/heat2d_config.goth
new file mode 100644
index 0000000..3749273
--- /dev/null
+++ b/examples/simulation/heat2d_config.goth
@@ -0,0 +1,18 @@
+# 2D Heat Diffusion Configuration
+# Edit these values to tune the simulation
+
+# Grid side length (total cells = n × n)
+let n = 15
+
+# Number of timesteps
+let nt = 30
+
+# Stability ratio r = α·Δt/Δx² (must be < 0.25 for 2D stability)
+let r = 0.2
+
+# Initial hot block: rows [icLo, icHi) and cols [icLo, icHi)
+let icLo = 5
+let icHi = 10
+
+# SVG cell size in pixels
+let cellSz = 30
diff --git a/examples/simulation/heat2d_tunable.goth b/examples/simulation/heat2d_tunable.goth
new file mode 100644
index 0000000..6b8cdb8
--- /dev/null
+++ b/examples/simulation/heat2d_tunable.goth
@@ -0,0 +1,97 @@
+# 2D Heat Diffusion with Tunable Parameters
+# Reads configuration from heat2d_config.goth via use import
+# Edit that file to change grid size, timesteps, diffusion rate, and IC
+# Output: heat2d.svg — n×n heatmap grid
+
+use "heat2d_config.goth"
+
+# Clamp value to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Temperature → red channel (0–255)
+╭─ colorR : F64 → I64
+╰─ toInt (clamp ₀ × 255.0)
+
+# Temperature → blue channel (0–255)
+╭─ colorB : F64 → I64
+╰─ toInt ((1.0 - clamp ₀) × 255.0)
+
+# Build "rgb(R,0,B)" color string from temperature
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat ",0," (strConcat (toString (colorB ₀)) ")")))
+
+# Initial condition: hot block [icLo, icHi) in both row and col
+# ₀ = flat index
+╭─ initCell : I64 → F64
+╰─ let row = ₀ / n in
+ let col = ₁ % n in
+ if row < icLo then 0.0
+ else if row < icHi then
+ if col < icLo then 0.0
+ else if col < icHi then 1.0
+ else 0.0
+ else 0.0
+
+# Boundary-aware grid access: returns 0.0 for out-of-bounds
+# ₂ = grid, ₁ = row, ₀ = col
+╭─ getVal : [n]F64 → I64 → I64 → F64
+╰─ if ₁ < 0 then 0.0
+ else if ₁ > (n - 1) then 0.0
+ else if ₀ < 0 then 0.0
+ else if ₀ > (n - 1) then 0.0
+ else ₂[₁ × n + ₀]
+
+# 5-point stencil at flat index k
+# ₁ = grid, ₀ = flat index k
+╭─ stencil2d : [n]F64 → I64 → F64
+╰─ let grid = ₁ in
+ let k = ₁ in
+ let row = k / n in
+ let col = k % n in
+ let center = grid[k] in
+ let north = getVal grid (row - 1) col in
+ let south = getVal grid (row + 1) col in
+ let west = getVal grid row (col - 1) in
+ let east = getVal grid row (col + 1) in
+ center + r × (north + south + east + west - 4.0 × center)
+
+# One timestep: map stencil over all grid cells
+╭─ stepGrid : [n]F64 → [n]F64
+╰─ let grid = ₀ in
+ ι (n × n) ↦ (λ→ stencil2d grid ₀)
+
+# Recursive time-stepping: evolve grid for t steps
+# ₁ = grid, ₀ = remaining steps
+╭─ evolve : [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₁
+ else evolve (stepGrid ₁) (₀ - 1)
+
+# SVG rectangle for one cell at flat index k
+# ₁ = flat index, ₀ = temperature value
+╭─ cellRect : I64 → F64 → String
+╰─ let k = ₁ in
+ let temp = ₁ in
+ let row = k / n in
+ let col = k % n in
+ strConcat "\n")))))))))
+
+# Recursively build SVG rect strings for cells [i..(n × n))
+# ₂ = grid, ₁ = current index, ₀ = accumulator string
+╭─ buildRects : [n]F64 → I64 → String → String
+╰─ if ₁ = (n × n) then ₀
+ else buildRects ₂ (₁ + 1) (strConcat ₀ (cellRect ₁ ₂[₁]))
+
+# Build complete SVG from final grid
+╭─ buildSvg : [n]F64 → String
+╰─ let rects = buildRects ₀ 0 "" in
+ let dim = n × cellSz in
+ strConcat "\n")))))
+
+# Main: run simulation, write SVG file
+╭─ main : I64 → ()
+╰─ let grid0 = ι (n × n) ↦ (λ→ initCell ₀) in
+ let final = evolve grid0 nt in
+ buildSvg final ▷ "heat2d.svg"
diff --git a/examples/simulation/laplace2d.goth b/examples/simulation/laplace2d.goth
new file mode 100644
index 0000000..9d92431
--- /dev/null
+++ b/examples/simulation/laplace2d.goth
@@ -0,0 +1,102 @@
+# 2D Laplace Equation Solver (Jacobi Iteration)
+# Solves ∇²u = 0 on a square domain with Dirichlet BCs
+# Jacobi: u_new[i,j] = (u[i-1,j] + u[i+1,j] + u[i,j-1] + u[i,j+1]) / 4
+# Grid: 15×15, 200 iterations
+# BCs: top edge = 1.0, other edges = 0.0
+# Output: laplace2d.svg — heatmap of steady-state solution
+
+# Grid side length
+let n = 15
+
+# Jacobi iterations
+let nIter = 200
+
+# Clamp to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Temperature → red channel (0–255)
+╭─ colorR : F64 → I64
+╰─ toInt (clamp ₀ × 255.0)
+
+# Temperature → blue channel (0–255)
+╭─ colorB : F64 → I64
+╰─ toInt ((1.0 - clamp ₀) × 255.0)
+
+# Build "rgb(R,0,B)" color string
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat ",0," (strConcat (toString (colorB ₀)) ")")))
+
+# Initial/boundary condition: top edge (row 0) = 1.0, else 0.0
+# ₀ = flat index
+╭─ initCell : I64 → F64
+╰─ if ₀ / n = 0 then 1.0 else 0.0
+
+# Is this cell on a boundary?
+# ₀ = flat index
+╭─ isBoundary : I64 → I64
+╰─ let row = ₀ / n in
+ let col = ₁ % n in
+ if row = 0 then 1
+ else if row = (n - 1) then 1
+ else if col = 0 then 1
+ else if col = (n - 1) then 1
+ else 0
+
+# Boundary-aware grid access
+# ₂ = grid, ₁ = row, ₀ = col
+╭─ getVal : [n]F64 → I64 → I64 → F64
+╰─ if ₁ < 0 then 0.0
+ else if ₁ > (n - 1) then 0.0
+ else if ₀ < 0 then 0.0
+ else if ₀ > (n - 1) then 0.0
+ else ₂[₁ × n + ₀]
+
+# Jacobi stencil at flat index k
+# ₁ = grid, ₀ = flat index k
+╭─ jacobiStep : [n]F64 → I64 → F64
+╰─ if isBoundary ₀ = 1 then ₁[₀]
+ else let row = ₀ / n in
+ let col = ₁ % n in
+ (getVal ₃ (row - 1) col + getVal ₃ (row + 1) col + getVal ₃ row (col - 1) + getVal ₃ row (col + 1)) / 4.0
+
+# One full Jacobi sweep
+╭─ sweep : [n]F64 → [n]F64
+╰─ let grid = ₀ in
+ ι (n × n) ↦ (λ→ jacobiStep grid ₀)
+
+# Recursive iteration
+# ₁ = grid, ₀ = remaining iterations
+╭─ iterate : [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₁
+ else iterate (sweep ₁) (₀ - 1)
+
+# SVG rectangle for one cell
+# ₁ = flat index, ₀ = temperature value
+╭─ cellRect : I64 → F64 → String
+╰─ strConcat "\n")))))
+
+# Recursively build SVG rects
+# ₂ = grid, ₁ = index, ₀ = accumulator
+╭─ buildRects : [n]F64 → I64 → String → String
+╰─ if ₁ = (n × n) then ₀
+ else buildRects ₂ (₁ + 1) (strConcat ₀ (cellRect ₁ ₂[₁]))
+
+# Build full SVG
+╭─ buildSvg : [n]F64 → String
+╰─ let rects = buildRects ₀ 0 "" in
+ strConcat "\n")
+
+# Build CSV: row,col,value
+╭─ buildCsv : [n]F64 → I64 → String → String
+╰─ if ₁ = (n × n) then ₀
+ else buildCsv ₂ (₁ + 1) (strConcat ₀ (strConcat (toString (₁ / n)) (strConcat "," (strConcat (toString (₁ % n)) (strConcat "," (strConcat (toString ₂[₁]) "\n"))))))
+
+# Main
+╭─ main : I64 → ()
+╰─ let grid0 = ι (n × n) ↦ (λ→ initCell ₀) in
+ let final = iterate grid0 nIter in
+ let _ = buildSvg final ▷ "laplace2d.svg" in
+ strConcat "row,col,value\n" (buildCsv final 0 "") ▷ "laplace2d.csv"
diff --git a/examples/simulation/newton_raphson.goth b/examples/simulation/newton_raphson.goth
new file mode 100644
index 0000000..ef533da
--- /dev/null
+++ b/examples/simulation/newton_raphson.goth
@@ -0,0 +1,101 @@
+# Newton-Raphson Root Finder
+# Finds cube root of a target value: solve x³ - target = 0
+# f(x) = x³ - target, f'(x) = 3x²
+# Iteration: x' = x - f(x)/f'(x) = x - (x³ - target)/(3x²)
+# Outputs: newton.svg — bar chart of error magnitude per iteration
+# Also prints CSV to stdout: iteration, estimate, error
+
+# Target value whose cube root we seek
+let target = 27.0
+
+# Number of iterations to run
+let nIter = 10
+
+# Initial guess
+let x0 = 1.0
+
+# True answer for error computation
+let answer = 3.0
+
+# f(x) = x³ - target
+╭─ f : F64 → F64
+╰─ ₀ × ₀ × ₀ - target
+
+# f'(x) = 3x²
+╭─ fprime : F64 → F64
+╰─ 3.0 × ₀ × ₀
+
+# One Newton step: x - f(x)/f'(x)
+╭─ step : F64 → F64
+╰─ ₀ - f ₀ / fprime ₀
+
+# Clamp value to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Absolute value
+╭─ absF : F64 → F64
+╰─ if ₀ < 0.0 then 0.0 - ₀ else ₀
+
+# Collect iteration history: returns array of x values at each step
+# ₂ = current x, ₁ = iteration count, ₀ = max iterations
+# Returns array of [x0, x1, ..., x_n]
+╭─ iterate : F64 → I64 → I64 → [n]F64
+╰─ if ₁ = ₀ then [₂]
+ else [₂] ⊕ iterate (step ₂) (₁ + 1) ₀
+
+# Error at each step (absolute)
+╭─ errorAt : F64 → F64
+╰─ absF (₀ - answer)
+
+# Map error for bar height: scale log10(error) into pixel range
+# SVG chart is 400px tall; we map error [1e-15, 10] → [0, 400]
+# log10(error) ranges from about -15 to 1; map to pixels
+╭─ barHeight : F64 → F64
+╰─ let err = errorAt ₀ in
+ if err < 1.0e-15 then 5.0
+ else let logErr = ln err / ln 10.0 in
+ let frac = (logErr + 16.0) / 17.0 in
+ clamp frac × 380.0 + 10.0
+
+# Temperature-style color for error magnitude (high error = red, low = blue)
+╭─ colorR : F64 → I64
+╰─ toInt (clamp (barHeight ₀ / 390.0) × 255.0)
+
+╭─ colorB : F64 → I64
+╰─ toInt ((1.0 - clamp (barHeight ₀ / 390.0)) × 255.0)
+
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat ",0," (strConcat (toString (colorB ₀)) ")")))
+
+# SVG bar for one iteration
+# ₁ = iteration index, ₀ = x value at that iteration
+╭─ bar : I64 → F64 → String
+╰─ let h = barHeight ₀ in
+ strConcat "\n")))))))
+
+# Recursively build bar SVG strings
+# ₂ = xs array, ₁ = current index, ₀ = accumulator
+╭─ buildBars : [n]F64 → I64 → String → String
+╰─ if ₁ = nIter then ₀
+ else buildBars ₂ (₁ + 1) (strConcat ₀ (bar ₁ ₂[₁]))
+
+# Build CSV output: "iteration,estimate,error\n..."
+# ₂ = xs array, ₁ = current index, ₀ = accumulator
+╭─ buildCsv : [n]F64 → I64 → String → String
+╰─ if ₁ = nIter then ₀
+ else buildCsv ₂ (₁ + 1) (strConcat ₀ (strConcat (toString ₁) (strConcat "," (strConcat (toString ₂[₁]) (strConcat "," (strConcat (toString (errorAt ₂[₁])) "\n"))))))
+
+# Build full SVG document
+╭─ buildSvg : [n]F64 → String
+╰─ let bars = buildBars ₀ 0 "" in
+ strConcat "\n"))
+
+# Main: run Newton-Raphson, write SVG and CSV
+╭─ main : I64 → ()
+╰─ let xs = iterate x0 0 nIter in
+ let _ = buildSvg xs ▷ "newton.svg" in
+ let csv = strConcat "iteration,estimate,error\n" (buildCsv xs 0 "") in
+ csv ▷ "newton.csv"
diff --git a/examples/simulation/power_iter.goth b/examples/simulation/power_iter.goth
new file mode 100644
index 0000000..c1e5bbf
--- /dev/null
+++ b/examples/simulation/power_iter.goth
@@ -0,0 +1,113 @@
+# Eigenvalue Power Iteration
+# Finds dominant eigenvalue and eigenvector of a symmetric matrix
+# Algorithm: v' = A*v / ||A*v||, λ = v'·(A*v')
+# Matrix: 4×4 symmetric (flat storage, row-major)
+# Output: power_iter.svg — bar chart of eigenvector components
+# power_iter.csv — iteration, eigenvalue estimate
+
+# Matrix dimension
+let dim = 4
+
+# Number of iterations
+let nIter = 20
+
+# 4×4 symmetric matrix (flat, row-major):
+# [ 4 1 2 0 ]
+# [ 1 3 1 1 ]
+# [ 2 1 5 2 ]
+# [ 0 1 2 3 ]
+# Dominant eigenvalue ≈ 7.57
+╭─ matA : I64 → F64
+╰─ [4.0, 1.0, 2.0, 0.0, 1.0, 3.0, 1.0, 1.0, 2.0, 1.0, 5.0, 2.0, 0.0, 1.0, 2.0, 3.0][₀]
+
+# Matrix-vector multiply: y = A*x
+# ₁ = A (as function I64→F64), but we use matA directly
+# Computes row i of A·x by dotting row i with x
+# ₁ = x vector, ₀ = row index i
+╭─ matvecRow : [n]F64 → I64 → F64
+╰─ Σ (ι dim ↦ (λ→ matA (₁ × dim + ₀) × ₂[₀]))
+
+# Full matrix-vector product: returns new vector
+╭─ matvec : [n]F64 → [n]F64
+╰─ let x = ₀ in
+ ι dim ↦ (λ→ matvecRow x ₀)
+
+# Vector L2 norm
+╭─ vecNorm : [n]F64 → F64
+╰─ √(Σ (₀ ↦ (λ→ ₀ × ₀)))
+
+# Scale vector by scalar: s * v
+# ₁ = scalar, ₀ = vector
+╭─ scaleVec : F64 → [n]F64 → [n]F64
+╰─ ₀ ↦ (λ→ ₀ × ₂)
+
+# Dot product of two vectors
+╭─ dotVec : [n]F64 → [n]F64 → F64
+╰─ Σ (ι dim ↦ (λ→ ₂[₀] × ₁[₀]))
+
+# One power iteration step: v -> A*v / ||A*v||
+╭─ powerStep : [n]F64 → [n]F64
+╰─ let av = matvec ₀ in
+ let nrm = vecNorm av in
+ scaleVec (1.0 / nrm) av
+
+# Compute eigenvalue estimate: λ = v · (A*v)
+╭─ eigenvalue : [n]F64 → F64
+╰─ dotVec ₀ (matvec ₀)
+
+# Collect eigenvalue history over iterations
+# ₂ = current v, ₁ = iteration, ₀ = max iterations
+╭─ eigenHist : [n]F64 → I64 → I64 → [n]F64
+╰─ if ₁ = ₀ then [eigenvalue ₂]
+ else [eigenvalue ₂] ⊕ eigenHist (powerStep ₂) (₁ + 1) ₀
+
+# Run power iteration to get final eigenvector
+# ₁ = v, ₀ = remaining iterations
+╭─ iterateVec : [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₁
+ else iterateVec (powerStep ₁) (₀ - 1)
+
+# Absolute value
+╭─ absF : F64 → F64
+╰─ if ₀ < 0.0 then 0.0 - ₀ else ₀
+
+# Clamp to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Color for eigenvector component (positive=green, negative=orange)
+╭─ toColorStr : F64 → String
+╰─ if ₀ > 0.0 then strConcat "rgb(0," (strConcat (toString (toInt (clamp ₀ × 255.0))) ",0)")
+ else strConcat "rgb(255," (strConcat (toString (toInt (clamp (1.0 - absF ₀) × 180.0))) ",0)")
+
+# SVG bar for eigenvector component
+# ₁ = index, ₀ = component value
+╭─ vecBar : I64 → F64 → String
+╰─ strConcat " 0.0 then 200 - toInt (absF ₀ × 180.0) else 200)) (strConcat "\" width=\"80\" height=\"" (strConcat (toString (toInt (absF ₀ × 180.0))) (strConcat "\" fill=\"" (strConcat (toColorStr ₀) "\"/>\n")))))))
+
+# Build eigenvector bar chart
+# ₂ = vec, ₁ = index, ₀ = accumulator
+╭─ buildVecBars : [n]F64 → I64 → String → String
+╰─ if ₁ = dim then ₀
+ else buildVecBars ₂ (₁ + 1) (strConcat ₀ (vecBar ₁ ₂[₁]))
+
+# Build eigenvalue convergence CSV
+# ₂ = eigenvalue history, ₁ = index, ₀ = accumulator
+╭─ buildCsv : [n]F64 → I64 → String → String
+╰─ if ₁ = nIter then ₀
+ else buildCsv ₂ (₁ + 1) (strConcat ₀ (strConcat (toString ₁) (strConcat "," (strConcat (toString ₂[₁]) "\n"))))
+
+# Build full SVG
+╭─ buildSvg : [n]F64 → String
+╰─ let bars = buildVecBars ₀ 0 "" in
+ strConcat "\n")))
+
+# Main: run power iteration, output SVG and CSV
+╭─ main : I64 → ()
+╰─ let v0 = ι dim ↦ (λ→ 1.0) in
+ let finalVec = iterateVec v0 nIter in
+ let hist = eigenHist v0 0 nIter in
+ let _ = buildSvg finalVec ▷ "power_iter.svg" in
+ strConcat "iteration,eigenvalue\n" (buildCsv hist 0 "") ▷ "power_iter.csv"
diff --git a/examples/simulation/wave1d.goth b/examples/simulation/wave1d.goth
new file mode 100644
index 0000000..cf2161e
--- /dev/null
+++ b/examples/simulation/wave1d.goth
@@ -0,0 +1,101 @@
+# 1D Wave Equation Simulation
+# Leapfrog scheme: u(t+1,i) = 2*u(t,i) - u(t-1,i) + C²*(u(t,i+1) - 2*u(t,i) + u(t,i-1))
+# where C = c*dt/dx is the Courant number
+# Grid: 40 points, 100 timesteps, C = 0.8 (stable when C ≤ 1)
+# Initial: Gaussian pulse, zero velocity
+# Output: wave1d.svg — wave profile at final time
+
+# Grid points
+let nx = 40
+
+# Timesteps
+let nt = 15
+
+# Courant number C = c*dt/dx (must be ≤ 1 for stability)
+let courant = 0.8
+
+# C² precomputed
+let c2 = 0.64
+
+# Clamp to [0, 1]
+╭─ clamp : F64 → F64
+╰─ if ₀ < 0.0 then 0.0
+ else if ₀ > 1.0 then 1.0
+ else ₀
+
+# Absolute value
+╭─ absF : F64 → F64
+╰─ if ₀ < 0.0 then 0.0 - ₀ else ₀
+
+# Initial condition: Gaussian pulse centered at x = 0.5
+# Grid spans [0, 1], dx = 1/(nx-1)
+# ₀ = cell index
+╭─ initCell : I64 → F64
+╰─ let x = toFloat ₀ / toFloat (nx - 1) in
+ let d = x - 0.5 in
+ exp (0.0 - d × d / 0.005)
+
+# Wave stencil: compute u_new[i] from u_cur and u_prev
+# ₂ = u_cur, ₁ = u_prev, ₀ = index i
+╭─ waveStencil : [n]F64 → [n]F64 → I64 → F64
+╰─ if ₀ = 0 then 0.0
+ else if ₀ = (nx - 1) then 0.0
+ else 2.0 × ₂[₀] - ₁[₀] + c2 × (₂[₀ - 1] + ₂[₀ + 1] - 2.0 × ₂[₀])
+
+# One timestep: produce u_new from u_cur and u_prev
+# ₁ = u_cur, ₀ = u_prev
+╭─ stepWave : [n]F64 → [n]F64 → [n]F64
+╰─ let cur = ₁ in
+ let prev = ₁ in
+ ι nx ↦ (λ→ waveStencil cur prev ₀)
+
+# Recursive time-stepping
+# ₂ = u_cur, ₁ = u_prev, ₀ = remaining steps
+╭─ evolve : [n]F64 → [n]F64 → I64 → [n]F64
+╰─ if ₀ = 0 then ₂
+ else let unew = stepWave ₂ ₁ in
+ evolve unew ₃ (₁ - 1)
+
+# Color mapping: blue for negative, white for zero, red for positive
+╭─ colorR : F64 → I64
+╰─ if ₀ > 0.0 then toInt (clamp ₀ × 255.0)
+ else 0
+
+╭─ colorG : F64 → I64
+╰─ 0
+
+╭─ colorB : F64 → I64
+╰─ if ₀ < 0.0 then toInt (clamp (0.0 - ₀) × 255.0)
+ else 0
+
+╭─ toColorStr : F64 → String
+╰─ strConcat "rgb(" (strConcat (toString (colorR ₀)) (strConcat "," (strConcat (toString (colorG ₀)) (strConcat "," (strConcat (toString (colorB ₀)) ")")))))
+
+# SVG bar for one cell — vertical bar chart of wave amplitude
+# Bar extends up (positive) or down (negative) from centerline at y=200
+# ₁ = index, ₀ = amplitude value
+╭─ cellBar : I64 → F64 → String
+╰─ strConcat " 0.0 then 200 - toInt (absF ₀ × 180.0) else 200)) (strConcat "\" width=\"18\" height=\"" (strConcat (toString (toInt (absF ₀ × 180.0))) (strConcat "\" fill=\"" (strConcat (toColorStr ₀) "\"/>\n")))))))
+
+# Recursively build SVG bars
+# ₂ = grid, ₁ = current index, ₀ = accumulator
+╭─ buildBars : [n]F64 → I64 → String → String
+╰─ if ₁ = nx then ₀
+ else buildBars ₂ (₁ + 1) (strConcat ₀ (cellBar ₁ ₂[₁]))
+
+# Build CSV: index,x,amplitude
+╭─ buildCsv : [n]F64 → I64 → String → String
+╰─ if ₁ = nx then ₀
+ else buildCsv ₂ (₁ + 1) (strConcat ₀ (strConcat (toString ₁) (strConcat "," (strConcat (toString (toFloat ₁ / toFloat (nx - 1))) (strConcat "," (strConcat (toString ₂[₁]) "\n"))))))
+
+# Build full SVG
+╭─ buildSvg : [n]F64 → String
+╰─ let bars = buildBars ₀ 0 "" in
+ strConcat "\n")))
+
+# Main: run simulation, write outputs
+╭─ main : I64 → ()
+╰─ let u0 = ι nx ↦ (λ→ initCell ₀) in
+ let final = evolve u0 u0 nt in
+ let _ = buildSvg final ▷ "wave1d.svg" in
+ strConcat "index,x,amplitude\n" (buildCsv final 0 "") ▷ "wave1d.csv"
diff --git a/examples/uncertainty/add_uncertain.goth b/examples/uncertainty/add_uncertain.goth
new file mode 100644
index 0000000..cd178d9
--- /dev/null
+++ b/examples/uncertainty/add_uncertain.goth
@@ -0,0 +1,8 @@
+# Add two uncertain values with error propagation
+# Input: a da b db (value±uncertainty for each)
+# Output: (a+b) ± √(da²+db²)
+# De Bruijn: ₃ = a, ₂ = da, ₁ = b, ₀ = db
+# Postcondition: propagated uncertainty ≥ max(da, db)
+
+╭─ main : F64 → F64 → F64 → F64 → (F64 ± F64)
+╰─ (₃ ± ₂) + (₁ ± ₀)
diff --git a/examples/uncertainty/chained_uncertain.goth b/examples/uncertainty/chained_uncertain.goth
new file mode 100644
index 0000000..a70ebe8
--- /dev/null
+++ b/examples/uncertainty/chained_uncertain.goth
@@ -0,0 +1,9 @@
+# Chained uncertainty propagation through multiple operations
+# Computes sin(√(a ± da) + (b ± db))
+# Input: a da b db
+# Output: result with fully propagated uncertainty
+# De Bruijn: ₃ = a, ₂ = da, ₁ = b, ₀ = db
+# Property: uncertainty propagates through each step automatically
+
+╭─ main : F64 → F64 → F64 → F64 → (F64 ± F64)
+╰─ sin (√(₃ ± ₂) + (₁ ± ₀))
diff --git a/examples/uncertainty/mul_uncertain.goth b/examples/uncertainty/mul_uncertain.goth
new file mode 100644
index 0000000..cdb019d
--- /dev/null
+++ b/examples/uncertainty/mul_uncertain.goth
@@ -0,0 +1,8 @@
+# Multiply two uncertain values with relative error propagation
+# Input: a da b db (value±uncertainty for each)
+# Output: (a×b) ± |a×b|×√((da/a)²+(db/b)²)
+# De Bruijn: ₃ = a, ₂ = da, ₁ = b, ₀ = db
+# Property: relative uncertainty adds in quadrature
+
+╭─ main : F64 → F64 → F64 → F64 → (F64 ± F64)
+╰─ (₃ ± ₂) × (₁ ± ₀)
diff --git a/examples/uncertainty/sin_uncertain.goth b/examples/uncertainty/sin_uncertain.goth
new file mode 100644
index 0000000..6d7fb5f
--- /dev/null
+++ b/examples/uncertainty/sin_uncertain.goth
@@ -0,0 +1,8 @@
+# Sine of an uncertain value
+# Input: x dx (value±uncertainty, x in radians)
+# Output: sin(x) ± |cos(x)|×dx
+# De Bruijn: ₁ = x, ₀ = dx
+# Property: derivative-based propagation δ(sin x) = |cos x| × δx
+
+╭─ main : F64 → F64 → (F64 ± F64)
+╰─ sin (₁ ± ₀)
diff --git a/examples/uncertainty/sqrt_uncertain.goth b/examples/uncertainty/sqrt_uncertain.goth
new file mode 100644
index 0000000..6530bac
--- /dev/null
+++ b/examples/uncertainty/sqrt_uncertain.goth
@@ -0,0 +1,9 @@
+# Square root of an uncertain value
+# Input: x dx (value±uncertainty)
+# Output: √x ± dx/(2√x)
+# De Bruijn: ₁ = x, ₀ = dx
+# Precondition: x ≥ 0
+# Property: derivative-based propagation δf = |f'(x)| × δx
+
+╭─ main : F64 → F64 → (F64 ± F64)
+╰─ √(₁ ± ₀)