diff --git a/.Rbuildignore b/.Rbuildignore
index aa9007f..b2c7ee0 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -6,3 +6,6 @@
 ^README\.Rmd$
 ^.*\.Rcheck$
 ^.*\.tar\.gz$
+^doc$
+^Meta$
+^docs$
diff --git a/.gitignore b/.gitignore
index a2ccfd5..740d4af 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,6 @@
 .Ruserdata
 *.Rcheck/
 *.tar.gz
+/doc/
+/Meta/
+/docs/
diff --git a/DESCRIPTION b/DESCRIPTION
index fa786ca..599ea04 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,25 +1,29 @@
 Package: numops
-Type: Package
 Title: Lightweight Numerical Operations
-Version: 1.0.1
-Authors@R: c(
-    person(
-      given = c("Joao", "Claudio"),
-      family = "Macosso",
-      email = "joaoclaudiomacosso@gmail.com",
-      role = c("aut", "cre"),
-      comment = c(ORCID = "0009-0006-5051-9312")
-    )
-  )
+Version: 1.0.0
+Authors@R: 
+    person(given = c("Joao", "Claudio"), 
+          family = "Macosso",
+          email = "joaoclaudiomacosso@gmail.com",
+          role = c("aut", "cre"),
+          comment = c(ORCID = "0009-0006-5051-9312"))
 URL: https://github.com/Macosso/numops
 BugReports: https://github.com/Macosso/numops/issues
-Description: Provides dependency-free helpers for common numerical operations
-    on vectors, matrices, and arrays. Includes bounds, interpolation,
-    remapping, division, norms, normalization, and adjacent differences.
+Description: Provides dependency-free helpers for recurring numerical tasks on
+    vectors, matrices, and arrays. Operations cover bounds, interpolation,
+    remapping, division, Euclidean norms, normalization, and adjacent
+    differences. Multi-input operations use strict scalar recycling, reject
+    incompatible lengths, and preserve names, dimensions, and dimension names
+    where applicable. Explicit handling of invalid intervals, zero
+    denominators, and zero norms gives consistent behavior for common edge
+    cases.
 License: GPL-3
 Encoding: UTF-8
 Suggests:
+    knitr,
+    rmarkdown,
     testthat (>= 3.0.0)
+VignetteBuilder: knitr
 Roxygen: list(markdown = TRUE)
 Config/roxygen2/version: 8.1.0
 Config/testthat/edition: 3
diff --git a/vignettes/clamp-test-wrap-numeric-values.Rmd b/vignettes/clamp-test-wrap-numeric-values.Rmd
new file mode 100644
index 0000000..ed8c93b
--- /dev/null
+++ b/vignettes/clamp-test-wrap-numeric-values.Rmd
@@ -0,0 +1,306 @@
+---
+title: "Clamp, Test, and Wrap Numeric Values in R"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Clamp, Test, and Wrap Numeric Values in R}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+Learn how to clamp numeric values to limits, test whether values fall within
+an inclusive range, and wrap angles or other periodic measurements in R with
+the dependency-free numops package.
+
+```{r}
+library(numops)
+```
+
+## Choose the right operation
+
+The bounds functions are related, but they answer different questions.
+
+| Goal | Function | Interval |
+|---|---|---|
+| Replace values outside fixed limits | `clamp()` | `[lower, upper]` |
+| Restrict probabilities or proportions | `clamp01()` | `[0, 1]` |
+| Identify values inside fixed limits | `in_range()` | `[lower, upper]` |
+| Map periodic values to one cycle | `wrap()` | `[lower, upper)` |
+
+`clamp()` and `clamp01()` modify values. `in_range()` returns logical results
+without changing its input. `wrap()` uses modular arithmetic rather than
+replacing values at the nearest boundary.
+
+## Clamp values in R
+
+### Clamp values to a numeric interval
+
+`clamp()` restricts each value to the closed interval `[lower, upper]`. Its
+element-wise formula is
+
+```text
+min(max(x, lower), upper)
+```
+
+Values inside the interval remain unchanged. Values below or above it are
+replaced by the nearest boundary.
+
+```{r}
+x <- c(-3, -1, 0.5, 4, 8)
+
+clamp(x, lower = -1, upper = 4)
+```
+
+The equivalent base R expression requires nested parallel extrema.
+
+```{r}
+pmin(pmax(x, -1), 4)
+```
+
+Both boundaries are included, so values equal to `lower` or `upper` remain
+unchanged.
+
+```{r}
+clamp(c(-1, 4), lower = -1, upper = 4)
+```
+
+### Preserve matrix and array shape
+
+Clamping a matrix preserves its dimensions and dimnames.
+
+```{r}
+x_matrix <- matrix(
+  c(-2, 0, 3, 8),
+  nrow = 2,
+  dimnames = list(c("a", "b"), c("x", "y"))
+)
+
+clamp(x_matrix, lower = 0, upper = 5)
+```
+
+The same rule applies to higher-dimensional arrays.
+
+## Clamp probabilities to zero and one
+
+`clamp01()` is equivalent to `clamp(x, 0, 1)`. It is convenient when small
+numerical errors produce probabilities or proportions just outside their
+valid interval.
+
+```{r}
+probabilities <- c(-0.02, 0.25, 0.8, 1.03)
+
+clamp01(probabilities)
+```
+
+Clamping changes invalid values, while range testing only identifies them.
+
+```{r}
+in_range(probabilities, lower = 0, upper = 1)
+clamp01(probabilities)
+```
+
+When unexpected values may indicate a data problem, test and investigate them
+before deciding whether clamping is appropriate.
+
+## Test whether values are within a range
+
+`in_range()` evaluates the inclusive condition
+
+```text
+x >= lower & x <= upper
+```
+
+```{r}
+x <- c(-2, 0, 3, 5, 8)
+
+in_range(x, lower = 0, upper = 5)
+```
+
+The logical result can be used directly for filtering.
+
+```{r}
+x[in_range(x, lower = 0, upper = 5)]
+```
+
+Both endpoints belong to the interval.
+
+```{r}
+in_range(c(0, 5), lower = 0, upper = 5)
+```
+
+Missing inputs produce missing logical results rather than `TRUE` or `FALSE`.
+
+```{r}
+in_range(c(1, NA, 5), lower = 0, upper = 4)
+```
+
+## Wrap angles and periodic values in R
+
+`wrap()` maps values to the half-open interval `[lower, upper)`. Its conceptual
+formula is
+
+```text
+lower + (x - lower) %% (upper - lower)
+```
+
+The implementation uses an equivalent calculation that avoids unnecessary
+overflow.
+
+### Wrap angles to one rotation
+
+Angles outside a canonical rotation can be wrapped to `[0, 360)`.
+
+```{r}
+angles <- c(-370, -10, 0, 360, 370, 725)
+
+wrap(angles, lower = 0, upper = 360)
+```
+
+The lower boundary is included and the upper boundary is excluded. Therefore,
+an angle of 360 degrees maps to zero rather than remaining 360.
+
+This differs from `in_range()`, which always includes both boundaries. To test
+whether an angle is already in the canonical half-open interval, use an
+explicit upper comparison.
+
+```{r}
+angles >= 0 & angles < 360
+```
+
+### Wrap clock times and phases
+
+The same operation applies to clock times.
+
+```{r}
+hours <- c(-2, 0, 12, 24, 27, 49)
+
+wrap(hours, lower = 0, upper = 24)
+```
+
+A symmetric interval is often useful for phase angles.
+
+```{r}
+phases <- c(-2 * pi, -pi, 0, pi, 2 * pi)
+
+wrap(phases, lower = -pi, upper = pi)
+```
+
+### Visualize periodic wrapping
+
+Wrapping produces a repeating sawtooth pattern. Each complete cycle returns
+the result to the lower boundary.
+
+```{r wrap-plot, fig.alt = "Wrapped angle against original angle"}
+angle_sequence <- seq(-720, 720, length.out = 500)
+
+plot(
+  angle_sequence,
+  wrap(angle_sequence, lower = 0, upper = 360),
+  type = "l",
+  xlab = "Original angle",
+  ylab = "Wrapped angle",
+  main = "Wrapping angles to [0, 360)"
+)
+```
+
+## Scalar recycling and vectorized bounds
+
+Bounds may have length one or the same length as the values being processed.
+Scalar bounds are recycled to the shared length.
+
+```{r}
+clamp(c(-2, 5, 20), lower = 0, upper = 10)
+```
+
+Vectorized bounds allow each position to use a different interval.
+
+```{r}
+x <- c(-2, 5, 20)
+lower <- c(0, 0, 10)
+upper <- c(1, 10, 15)
+
+clamp(x, lower, upper)
+in_range(x, lower, upper)
+```
+
+Every argument must have length one or a shared length. Other combinations are
+errors rather than partial recycling. Names, dimensions, and dimnames come
+from the first input already having the shared length.
+
+## Missing and infinite values
+
+The bounds functions use consistent rules for missing and non-finite values.
+
+| Condition | Behavior |
+|---|---|
+| Missing value in `x` | Produces a missing result |
+| Missing bound | Produces an error |
+| `lower > upper` | Produces an error |
+| Infinite bound in `clamp()` or `in_range()` | Allowed |
+| Infinite bound in `wrap()` | Produces an error |
+| Infinite value passed to `wrap()` | Produces `NaN` |
+| Equal lower and upper bounds in `wrap()` | Produces an error |
+
+An empty numeric input is returned with length zero, provided its bounds are
+valid.
+
+## A validation-and-correction workflow
+
+Consider measurements containing probabilities and angles. First record which
+probabilities are valid before applying any correction.
+
+```{r}
+measurements <- data.frame(
+  probability = c(-0.02, 0.35, 1.04, NA),
+  angle = c(-10, 45, 360, 725)
+)
+
+measurements$probability_valid <- in_range(
+  measurements$probability,
+  lower = 0,
+  upper = 1
+)
+
+measurements
+```
+
+If the out-of-range probabilities are known numerical artifacts, clamp them to
+the valid interval. Wrap the angles to a canonical rotation at the same time.
+
+```{r}
+measurements$probability <- clamp01(
+  measurements$probability
+)
+
+measurements$angle <- wrap(
+  measurements$angle,
+  lower = 0,
+  upper = 360
+)
+
+measurements
+```
+
+The validation column preserves which probabilities required attention, while
+the transformed columns are ready for downstream calculations.
+
+## Interval semantics at a glance
+
+The most important distinction among these functions is whether the upper
+boundary is included.
+
+```text
+clamp():    [lower, upper]
+in_range(): [lower, upper]
+wrap():     [lower, upper)
+```
+
+Use `clamp()` to enforce limits, `in_range()` to validate or filter values, and
+`wrap()` to represent periodic values in a single cycle. Use `clamp01()` when
+the required limits are specifically zero and one.
diff --git a/vignettes/numerical-operations-in-r.Rmd b/vignettes/numerical-operations-in-r.Rmd
new file mode 100644
index 0000000..8f4a675
--- /dev/null
+++ b/vignettes/numerical-operations-in-r.Rmd
@@ -0,0 +1,363 @@
+---
+title: "Common Numerical Operations in R with numops"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Common Numerical Operations in R with numops}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r setup, include = FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>"
+)
+```
+
+numops provides dependency-free helpers for common numerical operations on
+vectors, matrices, and arrays. This vignette shows how to clamp, wrap,
+interpolate, remap, divide, and normalize numeric data with consistent input
+validation and shape preservation.
+
+```{r}
+library(numops)
+```
+
+## Why use numops?
+
+Many numerical tasks in R require short but easily repeated expressions. For
+example, clamping a value to an interval requires nested `pmin()` and `pmax()`
+calls, while row normalization requires combining `rowSums()` and `sweep()`.
+numops gives these operations concise names and consistent behavior.
+
+The package has no runtime dependencies. Functions that combine inputs use
+strict scalar recycling, and operations on matrices and arrays preserve their
+dimensions and dimnames.
+
+## Quick reference
+
+| Problem | Function |
+|---|---|
+| Restrict values to an interval | `clamp()` |
+| Restrict probabilities to `[0, 1]` | `clamp01()` |
+| Test whether values are within bounds | `in_range()` |
+| Wrap periodic values | `wrap()` |
+| Interpolate between endpoints | `lerp()` |
+| Find a relative position between endpoints | `inv_lerp()` |
+| Map values between intervals | `remap()` |
+| Calculate a midpoint | `midpoint()` |
+| Divide with a fallback for zero | `divide_or()` |
+| Calculate Euclidean length | `l2_norm()` |
+| Normalize vectors or array slices | `normalize_l2()` |
+| Calculate length-preserving differences | `adjacent_difference()` |
+
+## Clamp and wrap values in R
+
+### Clamp values to an interval
+
+`clamp()` restricts each value to the closed interval `[lower, upper]`. Its
+element-wise formula is
+
+```text
+min(max(x, lower), upper)
+```
+
+```{r}
+x <- c(-3, -1, 0.5, 4, 8)
+
+clamp(x, lower = -1, upper = 4)
+```
+
+Use `clamp01()` when the interval is `[0, 1]`, as is common for probabilities
+and proportions.
+
+```{r}
+clamp01(c(-0.2, 0.35, 1.4, NA))
+```
+
+Missing values in `x` remain missing. Bounds may be infinite, but the lower
+bound cannot exceed the upper bound.
+
+### Test whether values are within a range
+
+`in_range()` performs the inclusive test
+`x >= lower & x <= upper`.
+
+```{r}
+in_range(1:7, lower = 3, upper = 5)
+```
+
+### Wrap angles and periodic values
+
+`wrap()` maps values to the half-open interval `[lower, upper)`. A value at the
+upper boundary wraps to the lower boundary.
+
+```{r}
+wrap(c(-10, 0, 360, 370), lower = 0, upper = 360)
+```
+
+This is useful for angles, clock times, phases, and other periodic quantities.
+
+## Linear interpolation and range remapping in R
+
+### Interpolate between endpoints
+
+`lerp()` implements linear interpolation:
+
+```text
+a + t * (b - a)
+```
+
+At `t = 0`, the result is `a`; at `t = 1`, it is `b`. Values of `t` outside
+`[0, 1]` extrapolate beyond the endpoints.
+
+```{r}
+lerp(a = 10, b = 20, t = c(0, 0.25, 0.5, 1, 1.5))
+```
+
+`inv_lerp()` performs the inverse calculation. It returns the position of `x`
+relative to `a` and `b`:
+
+```text
+(x - a) / (b - a)
+```
+
+```{r}
+inv_lerp(a = 10, b = 20, x = c(5, 10, 15, 20, 25))
+```
+
+Results below zero or above one indicate that `x` lies outside the endpoints.
+
+### Remap one interval to another
+
+`remap()` combines inverse interpolation and interpolation. It maps `x` from
+the interval `from` to the interval `to`:
+
+```text
+to[1] + (x - from[1]) / (from[2] - from[1]) * (to[2] - to[1])
+```
+
+```{r}
+remap(
+  c(0, 25, 50, 75, 100),
+  from = c(0, 100),
+  to = c(-1, 1)
+)
+```
+
+Values outside `from` are extrapolated. Call `clamp()` separately when the
+output must remain inside the destination interval.
+
+### Calculate midpoints safely
+
+`midpoint()` calculates the value halfway between corresponding endpoints. It
+uses equivalent formulas chosen to avoid unnecessary overflow for large finite
+values.
+
+```{r}
+midpoint(c(0, 10), c(10, 30))
+midpoint(-.Machine$double.xmax, .Machine$double.xmax)
+```
+
+## Handle division by zero in R
+
+`divide_or()` evaluates `x / y` except where `y` is zero. At those positions,
+it returns `default`.
+
+```{r}
+divide_or(
+  x = c(12, 8, 5),
+  y = c(3, 0, 2),
+  default = NA_real_
+)
+```
+
+Only a zero denominator triggers the fallback. A missing denominator produces
+a missing result according to ordinary R division.
+
+## Normalize vectors, matrices, and arrays
+
+### Calculate the L2 norm of a vector
+
+For values $x_1, \ldots, x_n$, the Euclidean or L2 norm is
+
+$$
+\lVert x \rVert_2 = \sqrt{\sum_{i = 1}^{n} x_i^2}.
+$$
+
+```{r}
+l2_norm(c(3, 4))
+```
+
+The calculation is scaled internally to avoid unnecessary overflow and
+underflow.
+
+```{r}
+l2_norm(c(1e308, 1e308))
+```
+
+### Normalize a vector to unit length
+
+`normalize_l2()` divides a vector by its L2 norm.
+
+```{r}
+unit_vector <- normalize_l2(c(3, 4))
+
+unit_vector
+l2_norm(unit_vector)
+```
+
+### Normalize matrix rows and columns
+
+The `margin` argument identifies the dimensions that index separate slices.
+For a matrix, `margin = 1` operates on rows and `margin = 2` operates on
+columns.
+
+```{r}
+x <- matrix(
+  c(3, 4, 0, 1, 2, 2),
+  nrow = 2,
+  byrow = TRUE,
+  dimnames = list(c("a", "b"), c("x", "y", "z"))
+)
+
+l2_norm(x, margin = 1)
+normalize_l2(x, margin = 1)
+```
+
+Column normalization uses the same interface.
+
+```{r}
+normalize_l2(x, margin = 2)
+```
+
+The approach extends to higher-dimensional arrays by supplying one or more
+dimensions in `margin`.
+
+### Handle zero-length slices
+
+A zero vector cannot be scaled to unit length. The `zero` argument controls
+the result:
+
+- `"keep"` leaves the zero slice unchanged.
+- `"na"` replaces the slice with missing values.
+- `"error"` stops the calculation.
+
+```{r}
+normalize_l2(c(0, 0), zero = "keep")
+normalize_l2(c(0, 0), zero = "na")
+```
+
+## Calculate adjacent differences
+
+`adjacent_difference()` keeps the first value and then records the change from
+each preceding value:
+
+```text
+result[1] = x[1]
+result[i] = x[i] - x[i - 1]
+```
+
+```{r}
+values <- c(10, 13, 12, 18)
+
+changes <- adjacent_difference(values)
+changes
+cumsum(changes)
+```
+
+Unlike `diff()`, the result has the same length as the input. Applying
+`cumsum()` reconstructs the original values when ordinary arithmetic is
+reversible.
+
+## Recycling and output shape
+
+Functions with multiple numeric inputs use strict scalar recycling. Every
+input must have length one or a shared length. Length-one inputs are recycled;
+other length combinations produce an error.
+
+```{r}
+lerp(a = 0, b = c(10, 20, 30), t = 0.5)
+```
+
+Names, dimensions, and dimnames come from the first input already having the
+shared length.
+
+```{r}
+x <- matrix(
+  1:4,
+  nrow = 2,
+  dimnames = list(c("first", "second"), c("a", "b"))
+)
+
+clamp(x, lower = 2, upper = 3)
+```
+
+This rule avoids the partial recycling that base R permits for some length
+combinations.
+
+## A complete data-preparation example
+
+Consider three sensor readings containing temperature, relative humidity, and
+wind direction. The variables use different units and need separate numerical
+transformations before they can be combined as features.
+
+```{r}
+sensors <- data.frame(
+  temperature = c(18, 22, 30),
+  humidity = c(45, 0, 75),
+  wind_direction = c(-10, 360, 450)
+)
+
+sensors
+```
+
+Map temperatures from 0--40 degrees to `[0, 1]`, convert percentage humidity
+to a proportion, and wrap wind directions to `[0, 360)`.
+
+```{r}
+temperature_score <- remap(
+  sensors$temperature,
+  from = c(0, 40),
+  to = c(0, 1)
+)
+
+humidity_score <- clamp01(
+  divide_or(sensors$humidity, 100, default = NA_real_)
+)
+
+wind_direction <- wrap(
+  sensors$wind_direction,
+  lower = 0,
+  upper = 360
+)
+
+transformed <- data.frame(
+  temperature_score,
+  humidity_score,
+  wind_direction
+)
+
+transformed
+```
+
+The two unitless scores can then be combined and normalized by row.
+
+```{r}
+features <- as.matrix(
+  transformed[c("temperature_score", "humidity_score")]
+)
+
+normalize_l2(features, margin = 1)
+```
+
+This workflow makes each numerical step explicit: remapping changes units,
+clamping enforces valid proportions, wrapping handles periodic values, and L2
+normalization scales feature vectors to unit length.
+
+## Summary
+
+numops supplies a compact vocabulary for recurring numerical tasks in R. The
+functions are dependency-free, work across vectors, matrices, and arrays, and
+share consistent rules for validation, recycling, missing values, and output
+shape. See the individual function help pages for complete edge-case behavior.