diff --git a/benchmarks/pandas/bench_case_when.py b/benchmarks/pandas/bench_case_when.py
new file mode 100644
index 00000000..6dfb8761
--- /dev/null
+++ b/benchmarks/pandas/bench_case_when.py
@@ -0,0 +1,35 @@
+"""Benchmark: case_when — conditional value selection on 100k-element Series (pandas 2.2+)"""
+import json, time
+import numpy as np
+import pandas as pd
+
+ROWS = 100_000
+WARMUP = 5
+ITERATIONS = 20
+
+data = np.arange(ROWS, dtype=float) % 100
+s = pd.Series(data)
+cond1 = s < 25
+cond2 = s < 50
+cond3 = s < 75
+
+caselist = [
+    (cond1, "low"),
+    (cond2, "medium-low"),
+    (cond3, "medium-high"),
+]
+
+for _ in range(WARMUP):
+    s.case_when(caselist)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    s.case_when(caselist)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "case_when",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_contingency.py b/benchmarks/pandas/bench_contingency.py
new file mode 100644
index 00000000..431f58cb
--- /dev/null
+++ b/benchmarks/pandas/bench_contingency.py
@@ -0,0 +1,78 @@
+"""
+Benchmark: contingency — expected_freq, relative_risk, odds_ratio, association
+Pure-numpy equivalents (no scipy) matching the TypeScript benchmark.
+Dataset: same 4x4 table as the TypeScript benchmark.
+"""
+import json
+import time
+import numpy as np
+
+WARMUP = 10
+ITERS = 50
+
+observed = np.array([
+    [120, 80, 40, 60],
+    [90, 110, 70, 30],
+    [50, 60, 100, 90],
+    [40, 50, 90, 120],
+], dtype=float)
+
+two_by_two = np.array([
+    [60.0, 40.0],
+    [30.0, 70.0],
+])
+
+
+def expected_freq(obs):
+    row_sums = obs.sum(axis=1, keepdims=True)
+    col_sums = obs.sum(axis=0, keepdims=True)
+    total = obs.sum()
+    return row_sums * col_sums / total
+
+
+def relative_risk(obs):
+    a, b = obs[0, 0], obs[0, 1]
+    c, d = obs[1, 0], obs[1, 1]
+    risk0 = a / (a + b)
+    risk1 = c / (c + d)
+    return risk0 / risk1
+
+
+def odds_ratio(obs):
+    a, b = obs[0, 0], obs[0, 1]
+    c, d = obs[1, 0], obs[1, 1]
+    return (a * d) / (b * c)
+
+
+def association_cramer(obs):
+    exp = expected_freq(obs)
+    chi2 = np.sum((obs - exp) ** 2 / exp)
+    n = obs.sum()
+    r, c = obs.shape
+    return np.sqrt(chi2 / n / min(r - 1, c - 1))
+
+
+# Warm up
+for _ in range(WARMUP):
+    expected_freq(observed)
+    relative_risk(two_by_two)
+    odds_ratio(two_by_two)
+    association_cramer(observed)
+
+# Measure
+start = time.perf_counter()
+for _ in range(ITERS):
+    expected_freq(observed)
+    relative_risk(two_by_two)
+    odds_ratio(two_by_two)
+    association_cramer(observed)
+total_s = time.perf_counter() - start
+total_ms = total_s * 1000
+mean_ms = total_ms / ITERS
+
+print(json.dumps({
+    "function": "contingency",
+    "mean_ms": round(mean_ms, 4),
+    "iterations": ITERS,
+    "total_ms": round(total_ms, 4),
+}))
diff --git a/benchmarks/pandas/bench_entropy.py b/benchmarks/pandas/bench_entropy.py
new file mode 100644
index 00000000..5c3f3f05
--- /dev/null
+++ b/benchmarks/pandas/bench_entropy.py
@@ -0,0 +1,44 @@
+import numpy as np
+import json
+import time
+
+N = 100
+WARMUP = 5
+ITERS = 50
+
+p = np.arange(1, N + 1, dtype=float)
+q = np.arange(N, 0, -1, dtype=float)
+
+# Normalise
+p_norm = p / p.sum()
+q_norm = q / q.sum()
+
+
+def entropy_fn(pk):
+    pk = pk / pk.sum()
+    return -np.sum(pk * np.log(pk + 1e-300))
+
+
+def kl_divergence(pk, qk):
+    pk = pk / pk.sum()
+    qk = qk / qk.sum()
+    mask = pk > 0
+    return np.sum(pk[mask] * np.log(pk[mask] / (qk[mask] + 1e-300)))
+
+
+for _ in range(WARMUP):
+    entropy_fn(p)
+    kl_divergence(p, q)
+
+t0 = time.perf_counter()
+for _ in range(ITERS):
+    entropy_fn(p)
+    kl_divergence(p, q)
+total_ms = (time.perf_counter() - t0) * 1000
+
+print(json.dumps({
+    "function": "entropy_klDivergence",
+    "mean_ms": total_ms / ITERS,
+    "iterations": ITERS,
+    "total_ms": total_ms,
+}))
diff --git a/benchmarks/pandas/bench_feather.py b/benchmarks/pandas/bench_feather.py
new file mode 100644
index 00000000..9a6424e8
--- /dev/null
+++ b/benchmarks/pandas/bench_feather.py
@@ -0,0 +1,42 @@
+"""
+Benchmark: read_feather / to_feather — Arrow Feather v2 round-trip on 10k rows
+"""
+import json
+import time
+import io
+import pandas as pd
+import numpy as np
+
+ROWS = 10_000
+WARMUP = 3
+ITERATIONS = 20
+
+rng = np.random.default_rng(42)
+df = pd.DataFrame({
+    "id": np.arange(ROWS, dtype=np.int64),
+    "value": np.arange(ROWS, dtype=np.float64) * 1.1,
+    "label": [f"cat_{i % 50}" for i in range(ROWS)],
+})
+
+# Warm up
+for _ in range(WARMUP):
+    buf = io.BytesIO()
+    df.to_feather(buf)
+    buf.seek(0)
+    pd.read_feather(buf)
+
+# Measure round-trip
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    buf = io.BytesIO()
+    df.to_feather(buf)
+    buf.seek(0)
+    pd.read_feather(buf)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "feather",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_flags_options.py b/benchmarks/pandas/bench_flags_options.py
new file mode 100644
index 00000000..3b105d7d
--- /dev/null
+++ b/benchmarks/pandas/bench_flags_options.py
@@ -0,0 +1,73 @@
+"""
+Benchmark: flags and options (pandas equivalent)
+
+Measures:
+ - DataFrame.flags.allows_duplicate_labels get+set
+ - pd.get_option / pd.set_option / pd.reset_option for multiple keys
+ - pd.options proxy read
+
+Dataset: 10,000-row Series and DataFrame; 20 measured iterations.
+"""
+
+import json
+import time
+
+import numpy as np
+import pandas as pd
+
+N = 10_000
+WARMUP = 5
+ITERS = 20
+
+data = np.arange(N, dtype=np.float64)
+s = pd.Series(data)
+df = pd.DataFrame({"a": data, "b": data})
+
+
+def bench_flags_options():
+    sink = 0
+    for _ in range(ITERS + WARMUP):
+        # flags on Series
+        prev = s.flags.allows_duplicate_labels
+        s.flags.allows_duplicate_labels = not prev
+        s.flags.allows_duplicate_labels = prev
+        sink ^= int(s.flags.allows_duplicate_labels)
+
+        # flags on DataFrame
+        prev_df = df.flags.allows_duplicate_labels
+        df.flags.allows_duplicate_labels = not prev_df
+        df.flags.allows_duplicate_labels = prev_df
+        sink ^= int(df.flags.allows_duplicate_labels)
+
+        # options get/set/reset
+        v = pd.get_option("display.max_rows")
+        pd.set_option("display.max_rows", v + 1)
+        pd.reset_option("display.max_rows")
+        sink ^= int(pd.options.display.max_rows > 0)
+
+        pd.set_option("display.max_columns", 20)
+        pd.reset_option("display.max_columns")
+        sink ^= int(pd.options.display.max_columns > 0)
+
+    return sink
+
+
+# Warm-up
+bench_flags_options()
+
+# Measure
+t0 = time.perf_counter()
+for _ in range(ITERS):
+    bench_flags_options()
+total = (time.perf_counter() - t0) * 1000  # ms
+
+print(
+    json.dumps(
+        {
+            "function": "flags_options",
+            "mean_ms": total / ITERS,
+            "iterations": ITERS,
+            "total_ms": total,
+        }
+    )
+)
diff --git a/benchmarks/pandas/bench_floating_array.py b/benchmarks/pandas/bench_floating_array.py
new file mode 100644
index 00000000..297062e2
--- /dev/null
+++ b/benchmarks/pandas/bench_floating_array.py
@@ -0,0 +1,42 @@
+"""
+Benchmark: FloatingArray (pandas Float64 nullable float array).
+N=100_000 elements with ~10% nulls. Tests from/sum/mean/min/max/add/fillna.
+"""
+import json
+import time
+
+import pandas as pd
+import numpy as np
+
+N = 100_000
+WARMUP = 3
+ITERATIONS = 20
+
+raw = [None if i % 10 == 0 else (i % 1000) * 0.001 - 0.5 for i in range(N)]
+
+for _ in range(WARMUP):
+    a = pd.array(raw, dtype="Float64")
+    float(a.sum())
+    float(a.mean())
+    float(a.min())
+    float(a.max())
+    a + 1.0
+    a.fillna(0.0)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    a = pd.array(raw, dtype="Float64")
+    float(a.sum())
+    float(a.mean())
+    float(a.min())
+    float(a.max())
+    a + 1.0
+    a.fillna(0.0)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "floating_array",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_gaussianKDE.py b/benchmarks/pandas/bench_gaussianKDE.py
new file mode 100644
index 00000000..353cb996
--- /dev/null
+++ b/benchmarks/pandas/bench_gaussianKDE.py
@@ -0,0 +1,45 @@
+"""Benchmark: Gaussian KDE evaluate on 1000-point dataset at 200 grid points.
+Uses Scott bandwidth (numpy equivalent of Silverman for univariate data).
+Pure numpy implementation to match pages workflow constraints.
+"""
+import json
+import time
+import numpy as np
+
+N = 1_000
+GRID = 200
+WARMUP = 3
+ITERATIONS = 20
+
+# Create dataset: mix of two gaussians (same as TS benchmark)
+i_vals = np.arange(N)
+data = np.sin(i_vals * 0.01) * 2 + np.where(i_vals % 2 == 0, 0.0, 5.0)
+
+xmin, xmax = -5.0, 10.0
+grid = np.linspace(xmin, xmax, GRID)
+
+
+def gaussian_kde_evaluate(data: np.ndarray, points: np.ndarray) -> np.ndarray:
+    """Gaussian KDE with Silverman bandwidth, pure numpy."""
+    n = len(data)
+    std = np.std(data, ddof=1)
+    bw = (4.0 / (3.0 * n)) ** 0.2 * std  # Silverman rule
+    diff = points[:, np.newaxis] - data[np.newaxis, :]  # (GRID, N)
+    kernels = np.exp(-0.5 * (diff / bw) ** 2) / (bw * np.sqrt(2 * np.pi))
+    return kernels.mean(axis=1)
+
+
+for _ in range(WARMUP):
+    gaussian_kde_evaluate(data, grid)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    gaussian_kde_evaluate(data, grid)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "gaussianKDE",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_hypothesis_tests.py b/benchmarks/pandas/bench_hypothesis_tests.py
new file mode 100644
index 00000000..ee6863b6
--- /dev/null
+++ b/benchmarks/pandas/bench_hypothesis_tests.py
@@ -0,0 +1,103 @@
+"""Benchmark: hypothesis_tests — scipy-style hypothesis tests vs pandas/scipy equivalents."""
+import json
+import math
+import time
+
+WARMUP = 3
+ITERS = 20
+N = 1000
+
+
+def make_data(n: int, seed: int) -> list:
+    arr = []
+    x = seed
+    for _ in range(n):
+        x = (x * 1664525 + 1013904223) & 0xFFFFFFFF
+        arr.append((x & 0xFFFFFFFF) / 0x100000000)
+    return [v * 4 + 2 for v in arr]
+
+
+a = make_data(N, 42)
+b = make_data(N, 99)
+
+# Pure-numpy/stdlib implementations matching what tsb does
+try:
+    import numpy as np
+
+    def ttest1samp_np(x, popmean):
+        x = np.asarray(x)
+        n = len(x)
+        mean = x.mean()
+        se = x.std(ddof=1) / math.sqrt(n)
+        t = (mean - popmean) / se
+        return t
+
+    def ttestind_np(x, y):
+        x, y = np.asarray(x), np.asarray(y)
+        nx, ny = len(x), len(y)
+        vx, vy = x.var(ddof=1), y.var(ddof=1)
+        se = math.sqrt(vx / nx + vy / ny)
+        t = (x.mean() - y.mean()) / se
+        return t
+
+    def ttestrel_np(x, y):
+        d = np.asarray(x) - np.asarray(y)
+        n = len(d)
+        t = d.mean() / (d.std(ddof=1) / math.sqrt(n))
+        return t
+
+    def foneway_np(*groups):
+        grand = np.concatenate(groups)
+        grand_mean = grand.mean()
+        ss_between = sum(len(g) * (np.mean(g) - grand_mean) ** 2 for g in groups)
+        ss_within = sum(((np.asarray(g) - np.mean(g)) ** 2).sum() for g in groups)
+        df_between = len(groups) - 1
+        df_within = len(grand) - len(groups)
+        F = (ss_between / df_between) / (ss_within / df_within)
+        return F
+
+    def pearsonr_np(x, y):
+        x, y = np.asarray(x), np.asarray(y)
+        r = np.corrcoef(x, y)[0, 1]
+        return r
+
+    def spearmanr_np(x, y):
+        x, y = np.asarray(x), np.asarray(y)
+        rx = np.argsort(np.argsort(x)).astype(float)
+        ry = np.argsort(np.argsort(y)).astype(float)
+        r = np.corrcoef(rx, ry)[0, 1]
+        return r
+
+    def mannwhitneyu_np(x, y):
+        nx, ny = len(x), len(y)
+        all_vals = sorted([(v, 0) for v in x] + [(v, 1) for v in y], key=lambda t: t[0])
+        ranks = list(range(1, nx + ny + 1))
+        u1 = sum(ranks[i] for i, (_, g) in enumerate(all_vals) if g == 0)
+        u1 -= nx * (nx + 1) / 2
+        return u1
+
+    HAS_NP = True
+except ImportError:
+    HAS_NP = False
+
+
+def bench():
+    total = 0.0
+    for i in range(WARMUP + ITERS):
+        t0 = time.perf_counter()
+        if HAS_NP:
+            ttest1samp_np(a, 2.5)
+            ttestind_np(a, b)
+            ttestrel_np(a, b)
+            foneway_np(a, b)
+            pearsonr_np(a, b)
+            spearmanr_np(a, b)
+            mannwhitneyu_np(a, b)
+        elapsed = (time.perf_counter() - t0) * 1000
+        if i >= WARMUP:
+            total += elapsed
+    mean_ms = total / ITERS
+    print(json.dumps({"function": "hypothesis_tests", "mean_ms": mean_ms, "iterations": ITERS, "total_ms": total}))
+
+
+bench()
diff --git a/benchmarks/pandas/bench_integer_array.py b/benchmarks/pandas/bench_integer_array.py
new file mode 100644
index 00000000..fb481445
--- /dev/null
+++ b/benchmarks/pandas/bench_integer_array.py
@@ -0,0 +1,41 @@
+"""Benchmark: IntegerArray — nullable integer extension array operations.
+N=100_000 elements with ~10% nulls using pandas IntegerArray.
+Tests: from_sequence, sum, mean, min, max, add scalar, fillna.
+"""
+import json
+import time
+import numpy as np
+import pandas as pd
+
+N = 100_000
+WARMUP = 3
+ITERATIONS = 20
+
+# Build input with ~10% nulls (same pattern as TS version)
+raw = [(None if i % 10 == 0 else int((i % 1000) - 500)) for i in range(N)]
+
+
+def run():
+    a = pd.array(raw, dtype="Int32")
+    _ = a.sum(skipna=True)
+    _ = a.mean(skipna=True)
+    _ = a.min(skipna=True)
+    _ = a.max(skipna=True)
+    _ = a + 1
+    _ = a.fillna(0)
+
+
+for _ in range(WARMUP):
+    run()
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    run()
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "integer_array",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_linregress_polyfit.py b/benchmarks/pandas/bench_linregress_polyfit.py
new file mode 100644
index 00000000..c4c08a97
--- /dev/null
+++ b/benchmarks/pandas/bench_linregress_polyfit.py
@@ -0,0 +1,43 @@
+"""Benchmark: linregress and polyfit — simple linear regression and polynomial fit.
+Dataset: 10,000 points, 20 iterations.
+"""
+import json
+import time
+import numpy as np
+
+N = 10_000
+WARMUP = 3
+ITERATIONS = 20
+
+x = np.arange(N, dtype=float) / N
+y = 2.5 * x + 1.0 + np.sin(np.arange(N) * 0.01) * 0.1
+
+
+def linregress_numpy(xs, ys):
+    """Simple OLS linear regression matching scipy.stats.linregress."""
+    n = len(xs)
+    sx = xs.sum()
+    sy = ys.sum()
+    sxx = (xs * xs).sum()
+    sxy = (xs * ys).sum()
+    slope = (n * sxy - sx * sy) / (n * sxx - sx * sx)
+    intercept = (sy - slope * sx) / n
+    return slope, intercept
+
+
+for _ in range(WARMUP):
+    linregress_numpy(x, y)
+    np.polyfit(x, y, 2)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    linregress_numpy(x, y)
+    np.polyfit(x, y, 2)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "linregress_polyfit",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_lreshape.py b/benchmarks/pandas/bench_lreshape.py
new file mode 100644
index 00000000..96b02a3e
--- /dev/null
+++ b/benchmarks/pandas/bench_lreshape.py
@@ -0,0 +1,37 @@
+"""Benchmark: lreshape — wide-to-long reshape using named column groups.
+Dataset: 10,000 rows with 4 value columns (v1..v4), 50 iterations.
+"""
+import json
+import time
+import numpy as np
+import pandas as pd
+
+N = 10_000
+WARMUP = 3
+ITERATIONS = 50
+
+ids = np.arange(N)
+data = {
+    "id": ids,
+    "v1": ids * 1.0,
+    "v2": ids * 2.0,
+    "v3": ids * 3.0,
+    "v4": ids * 4.0,
+}
+df = pd.DataFrame(data)
+groups = {"value": ["v1", "v2", "v3", "v4"]}
+
+for _ in range(WARMUP):
+    pd.lreshape(df, groups)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    pd.lreshape(df, groups)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "lreshape",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_multivariate.py b/benchmarks/pandas/bench_multivariate.py
new file mode 100644
index 00000000..4922a1d5
--- /dev/null
+++ b/benchmarks/pandas/bench_multivariate.py
@@ -0,0 +1,53 @@
+"""Benchmark: multivariate statistics — Mahalanobis distance, covariance matrix, PCA
+Dataset: 500 observations x 5 features (matching the TypeScript benchmark)
+"""
+import json
+import time
+import numpy as np
+
+N = 500
+P = 5
+WARMUP = 3
+ITERATIONS = 20
+
+# Generate deterministic dataset matching TS version
+i_idx = np.arange(N)
+j_idx = np.arange(P)
+X = np.sin(i_idx[:, None] * 0.1 + j_idx[None, :]) * 10 + j_idx[None, :] * 2  # (N, P)
+
+u = X[0]
+v = X[1]
+
+# Pre-compute covariance and diagonal inverse covariance
+cov = np.cov(X, rowvar=False)  # (P, P) sample covariance
+diag_inv_cov = np.diag(1.0 / np.maximum(np.diag(cov), 1e-10))  # diagonal approx of VI
+
+
+def run_iteration(X, u, v, diag_inv_cov):
+    # Mahalanobis distance with pre-computed VI
+    diff = u - v
+    dist = np.sqrt(diff @ diag_inv_cov @ diff)
+    # Covariance matrix
+    cov_m = np.cov(X, rowvar=False)
+    # PCA via SVD (3 components)
+    X_centered = X - X.mean(axis=0)
+    _, _, Vt = np.linalg.svd(X_centered, full_matrices=False)
+    components = Vt[:3]
+    scores = X_centered @ components.T
+    return dist, cov_m, scores
+
+
+for _ in range(WARMUP):
+    run_iteration(X, u, v, diag_inv_cov)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    run_iteration(X, u, v, diag_inv_cov)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "multivariate",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_mutual_information.py b/benchmarks/pandas/bench_mutual_information.py
new file mode 100644
index 00000000..2d0bafe1
--- /dev/null
+++ b/benchmarks/pandas/bench_mutual_information.py
@@ -0,0 +1,68 @@
+import numpy as np
+import json
+import time
+
+N = 1000
+WARMUP = 5
+ITERS = 50
+CATS = 10
+
+# Same paired observations as the TS benchmark
+xs = np.array([i % CATS for i in range(N)])
+ys = np.array([(i % CATS) + (i // CATS) % 3 for i in range(N)])
+
+
+def mutual_information(xs, ys):
+    """Compute mutual information I(X;Y) from paired observations."""
+    n = len(xs)
+    ux, cx = np.unique(xs, return_counts=True)
+    uy, cy = np.unique(ys, return_counts=True)
+    px = cx / n
+    py = cy / n
+
+    # Joint counts
+    joint_counts = {}
+    for x, y in zip(xs, ys):
+        key = (int(x), int(y))
+        joint_counts[key] = joint_counts.get(key, 0) + 1
+
+    mi = 0.0
+    for (xi, yi), cnt in joint_counts.items():
+        pxy = cnt / n
+        pxi = px[np.searchsorted(ux, xi)]
+        pyi = py[np.searchsorted(uy, yi)]
+        if pxy > 0:
+            mi += pxy * np.log(pxy / (pxi * pyi + 1e-300))
+    return mi
+
+
+def normalized_mi(xs, ys):
+    """Normalized mutual information (arithmetic normalization)."""
+    mi = mutual_information(xs, ys)
+    n = len(xs)
+    _, cx = np.unique(xs, return_counts=True)
+    _, cy = np.unique(ys, return_counts=True)
+    px = cx / n
+    py = cy / n
+    hx = -np.sum(px * np.log(px + 1e-300))
+    hy = -np.sum(py * np.log(py + 1e-300))
+    denom = (hx + hy) / 2
+    return mi / denom if denom > 0 else 0.0
+
+
+for _ in range(WARMUP):
+    mutual_information(xs, ys)
+    normalized_mi(xs, ys)
+
+t0 = time.perf_counter()
+for _ in range(ITERS):
+    mutual_information(xs, ys)
+    normalized_mi(xs, ys)
+total_ms = (time.perf_counter() - t0) * 1000
+
+print(json.dumps({
+    "function": "mutual_information",
+    "mean_ms": total_ms / ITERS,
+    "iterations": ITERS,
+    "total_ms": total_ms,
+}))
diff --git a/benchmarks/pandas/bench_pipe_apply.py b/benchmarks/pandas/bench_pipe_apply.py
new file mode 100644
index 00000000..0025b60d
--- /dev/null
+++ b/benchmarks/pandas/bench_pipe_apply.py
@@ -0,0 +1,56 @@
+"""Benchmark: pipe / apply / applymap on 10,000-row datasets.
+
+Exercises three functional-pipeline operations:
+  - pipe: chain 4 transforms on a Series
+  - apply: element-wise function on 10k-element Series
+  - applymap: element-wise function on 10k x 3 DataFrame
+"""
+import json
+import time
+import numpy as np
+import pandas as pd
+
+N = 10_000
+WARMUP = 5
+ITERATIONS = 20
+
+raw = np.array([(i % 100) + 1 for i in range(N)], dtype=float)
+series = pd.Series(raw, name="x")
+df = pd.DataFrame({
+    "a": [(i % 50) + 1 for i in range(N)],
+    "b": [(i % 30) + 1 for i in range(N)],
+    "c": [(i % 20) + 1 for i in range(N)],
+}, dtype=float)
+
+
+def run_pipe(s):
+    return s.pipe(lambda x: x + 1).pipe(lambda x: x * 2).pipe(lambda x: x - 1).pipe(lambda x: x / 2)
+
+
+def apply_fn(v):
+    return v * 2 + 1
+
+
+def applymap_fn(v):
+    return v * 2
+
+
+# Warm-up
+for _ in range(WARMUP):
+    run_pipe(series)
+    series.apply(apply_fn)
+    df.map(applymap_fn)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    run_pipe(series)
+    series.apply(apply_fn)
+    df.map(applymap_fn)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "pipe_apply",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_readStata.py b/benchmarks/pandas/bench_readStata.py
new file mode 100644
index 00000000..997a8762
--- /dev/null
+++ b/benchmarks/pandas/bench_readStata.py
@@ -0,0 +1,37 @@
+"""Benchmark: readStata / toStata round-trip on a 10k-row DataFrame"""
+import json, time, tempfile, os
+import numpy as np
+import pandas as pd
+
+ROWS = 10_000
+WARMUP = 3
+ITERATIONS = 20
+
+ids = np.arange(ROWS, dtype=np.int32)
+values = np.sin(np.arange(ROWS) * 0.01) * 1000
+categories = np.array([f"cat_{i % 5}" for i in range(ROWS)])
+
+df = pd.DataFrame({"id": ids, "value": values, "category": categories})
+
+# Write to a temp Stata file so read_stata benchmarks read from disk
+tmp = tempfile.NamedTemporaryFile(suffix=".dta", delete=False)
+tmp.close()
+df.to_stata(tmp.name, write_index=False)
+
+# Warm up
+for _ in range(WARMUP):
+    pd.read_stata(tmp.name)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    pd.read_stata(tmp.name)
+total = (time.perf_counter() - start) * 1000
+
+os.unlink(tmp.name)
+
+print(json.dumps({
+    "function": "readStata",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_read_fwf.py b/benchmarks/pandas/bench_read_fwf.py
new file mode 100644
index 00000000..b13760f6
--- /dev/null
+++ b/benchmarks/pandas/bench_read_fwf.py
@@ -0,0 +1,37 @@
+"""Benchmark: read_fwf — parse a fixed-width formatted text file into a DataFrame.
+Dataset: 10,000 rows x 4 columns (id, name, value, flag).
+"""
+import json, time, io
+import numpy as np
+import pandas as pd
+
+ROWS = 10_000
+WARMUP = 3
+ITERATIONS = 10
+
+# Build fixed-width text matching the TypeScript benchmark
+lines = ["id    name      value     flag"]
+for i in range(ROWS):
+    id_col = str(i).rjust(6)
+    name_col = ("item" + str(i % 500)).ljust(10)
+    value_col = f"{np.sin(i * 0.01) * 1000:.3f}".rjust(10)
+    flag_col = ("Y" if i % 2 == 0 else "N").ljust(4)
+    lines.append(id_col + name_col + value_col + flag_col)
+text = "\n".join(lines)
+
+colspecs = [(0, 6), (6, 16), (16, 26), (26, 30)]
+
+for _ in range(WARMUP):
+    pd.read_fwf(io.StringIO(text), colspecs=colspecs, header=0)
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    pd.read_fwf(io.StringIO(text), colspecs=colspecs, header=0)
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "readFwf",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_read_sas.py b/benchmarks/pandas/bench_read_sas.py
new file mode 100644
index 00000000..9af9c561
--- /dev/null
+++ b/benchmarks/pandas/bench_read_sas.py
@@ -0,0 +1,148 @@
+"""
+Benchmark: pd.read_sas — parse a 1,000-row SAS XPORT (XPT) file.
+Outputs JSON: {"function": "read_sas", "mean_ms": ..., "iterations": ..., "total_ms": ...}
+"""
+import json
+import time
+import io
+import struct
+import math
+import numpy as np
+import pandas as pd
+
+ROWS = 1_000
+WARMUP = 3
+ITERATIONS = 20
+
+# ── IBM 370 double encoder ────────────────────────────────────────────────────
+
+def ibm_encode(val: float) -> bytes:
+    if val == 0.0:
+        return b"\x00" * 8
+    if not math.isfinite(val):
+        return b"\x2e" + b"\x00" * 7
+    sign = 1 if val < 0 else 0
+    abs_val = abs(val)
+    exp = 0
+    mant = abs_val
+    while mant >= 1.0:
+        mant /= 16
+        exp += 1
+    while mant < 1.0 / 16 and mant > 0:
+        mant *= 16
+        exp -= 1
+    mant_int = round(mant * 2**56)
+    out = bytearray(8)
+    out[0] = (sign << 7) | ((exp + 64) & 0x7f)
+    for i in range(1, 8):
+        out[i] = (mant_int >> ((7 - i) * 8)) & 0xff
+    return bytes(out)
+
+# ── Minimal XPORT v5 builder ─────────────────────────────────────────────────
+
+def build_xpt(num_vars, char_vars, rows_data):
+    RECORD = 80
+
+    def pad80(s):
+        return s.ljust(RECORD).encode("ascii")[:RECORD]
+
+    def write_u16(val):
+        return struct.pack(">H", val)
+
+    def write_u32(val):
+        return struct.pack(">I", val)
+
+    # Compute variable metadata
+    metas = []
+    pos = 0
+    for name in num_vars:
+        metas.append({"type": 1, "name": name, "len": 8, "pos": pos})
+        pos += 8
+    for name, length in char_vars:
+        metas.append({"type": 2, "name": name, "len": length, "pos": pos})
+        pos += length
+    row_len = pos
+
+    chunks = bytearray()
+
+    # Library header (5 × 80 bytes)
+    chunks += pad80("HEADER RECORD*******LIBRARY HEADER RECORD!!!!!!!000000000000000000000000000000  ")
+    chunks += pad80("SAS     SAS     SASLIB  6.06    ASCII")
+    chunks += pad80("20240101")
+    chunks += pad80("")
+    chunks += pad80("")
+
+    # Member header (3 × 80 bytes)
+    chunks += pad80("HEADER RECORD*******MEMBER  HEADER RECORD!!!!!!!000000000000000000000000000001600000000140  ")
+    chunks += pad80("SAS     BENCH   SASDATA 6.06    ASCII")
+    chunks += pad80("")
+
+    # Namestr header
+    nvar = len(metas)
+    chunks += pad80(f"HEADER RECORD*******NAMESTR HEADER RECORD!!!!!!!{nvar:06d}00000000000000000000  ")
+
+    # Namestr records (140 bytes each)
+    ns_buf = bytearray(nvar * 140)
+    for i, m in enumerate(metas):
+        off = i * 140
+        ns_buf[off:off+2] = write_u16(m["type"])
+        ns_buf[off+2:off+4] = write_u16(140)
+        name_bytes = m["name"].encode("ascii").ljust(8)[:8]
+        ns_buf[off+4:off+12] = name_bytes
+        ns_buf[off+52:off+54] = write_u16(m["len"])
+        ns_buf[off+84:off+88] = write_u32(m["pos"])
+    padded_ns = math.ceil(len(ns_buf) / RECORD) * RECORD
+    ns_buf_padded = ns_buf.ljust(padded_ns, b"\x00")
+    chunks += ns_buf_padded
+
+    # Obs header
+    chunks += pad80("HEADER RECORD*******OBS     HEADER RECORD!!!!!!!000000000000000000000000000000  ")
+
+    # Observations
+    padded_row_len = math.ceil(row_len / RECORD) * RECORD
+    obs_buf = bytearray(len(rows_data) * padded_row_len)
+    for r, row in enumerate(rows_data):
+        base = r * padded_row_len
+        for m in metas:
+            val = row.get(m["name"])
+            if m["type"] == 1:
+                encoded = ibm_encode(float(val) if val is not None else 0.0)
+                obs_buf[base + m["pos"]:base + m["pos"] + 8] = encoded
+            else:
+                s = str(val) if val is not None else ""
+                b = s.encode("ascii")[:m["len"]].ljust(m["len"], b" ")
+                obs_buf[base + m["pos"]:base + m["pos"] + m["len"]] = b
+
+    chunks += obs_buf
+    return bytes(chunks)
+
+
+# ── Build dataset ─────────────────────────────────────────────────────────────
+
+rows_data = [
+    {"id": float(i), "value": i * 1.5, "score": math.sin(i * 0.01), "label": f"item_{i % 100}"}
+    for i in range(ROWS)
+]
+
+xpt_bytes = build_xpt(
+    ["id", "value", "score"],
+    [("label", 12)],
+    rows_data,
+)
+
+# ── Benchmark ─────────────────────────────────────────────────────────────────
+
+for _ in range(WARMUP):
+    pd.read_sas(io.BytesIO(xpt_bytes), format="xport")
+
+start = time.perf_counter()
+for _ in range(ITERATIONS):
+    pd.read_sas(io.BytesIO(xpt_bytes), format="xport")
+total = (time.perf_counter() - start) * 1000
+
+print(json.dumps({
+    "function": "read_sas",
+    "mean_ms": total / ITERATIONS,
+    "iterations": ITERATIONS,
+    "total_ms": total,
+}))
diff --git a/benchmarks/pandas/bench_sql.py b/benchmarks/pandas/bench_sql.py
new file mode 100644
index 00000000..1b419767
--- /dev/null
+++ b/benchmarks/pandas/bench_sql.py
@@ -0,0 +1,53 @@
+"""Benchmark: read_sql / to_sql on 10k-row DataFrames using SQLite in-memory"""
+import json
+import time
+import math
+import sqlite3
+import pandas as pd
+
+ROWS = 10_000
+WARMUP = 3
+ITERATIONS = 20
+
+# ── Build a matching dataset ──────────────────────────────────────────────────
+data = {
+    "id": list(range(ROWS)),
+    "value": [math.sin(i * 0.01) * 1000 for i in range(ROWS)],
+    "label": [f"item_{i % 100}" for i in range(ROWS)],
+}
+df = pd.DataFrame(data)
+
+# ── SQLite in-memory database ─────────────────────────────────────────────────
+con = sqlite3.connect(":memory:")
+df.to_sql("mock_table", con, index=False, if_exists="replace")
+
+# ── Warm-up reads ─────────────────────────────────────────────────────────────
+for _ in range(WARMUP):
+    pd.read_sql_query("SELECT * FROM mock_table", con)
+
+# ── read_sql_query benchmark ──────────────────────────────────────────────────
+start_read = time.perf_counter()
+for _ in range(ITERATIONS):
+    pd.read_sql_query("SELECT * FROM mock_table", con)
+total_read = (time.perf_counter() - start_read) * 1000
+
+# ── Warm-up writes ────────────────────────────────────────────────────────────
+for _ in range(WARMUP):
+    df.to_sql("bench_table", con, index=False, if_exists="replace")
+
+# ── to_sql benchmark ──────────────────────────────────────────────────────────
+start_write = time.perf_counter()
+for _ in range(ITERATIONS):
+    df.to_sql("bench_table", con, index=False, if_exists="replace")
+total_write = (time.perf_counter() - start_write) * 1000
+
+con.close()
+
+print(json.dumps({
+    "function": "sql",
+    "mean_ms": (total_read + total_write) / (2 * ITERATIONS),
+    "iterations": ITERATIONS,
+    "total_ms": total_read + total_write,
+    "read_mean_ms": total_read / ITERATIONS,
+    "write_mean_ms": total_write / ITERATIONS,
+}))
diff --git a/benchmarks/pandas/bench_to_excel.py b/benchmarks/pandas/bench_to_excel.py
new file mode 100644
index 00000000..05b04472
--- /dev/null
+++ b/benchmarks/pandas/bench_to_excel.py
@@ -0,0 +1,27 @@
+"""Benchmark: to_excel — write a DataFrame to an XLSX buffer (BytesIO)."""
+import json, time, io
+import pandas as pd
+
+ROWS = 5_000
+WARMUP = 3
+ITERATIONS = 20
+
+df = pd.DataFrame({
+    "name": [f"name_{i % 1000}" for i in range(ROWS)],
+    "value": [i * 1.5 for i in range(ROWS)],
+    "flag": [i % 2 == 0 for i in range(ROWS)],
+})
+
+for _ in range(WARMUP):
+    buf = io.BytesIO()
+    df.to_excel(buf, index=True)
+
+times = []
+for _ in range(ITERATIONS):
+    buf = io.BytesIO()
+    t0 = time.perf_counter()
+    df.to_excel(buf, index=True)
+    times.append((time.perf_counter() - t0) * 1000)
+
+total_ms = sum(times)
+print(json.dumps({"function": "to_excel", "mean_ms": round(total_ms / ITERATIONS, 3), "iterations": ITERATIONS, "total_ms": round(total_ms, 3)}))
diff --git a/benchmarks/pandas/bench_xml.py b/benchmarks/pandas/bench_xml.py
new file mode 100644
index 00000000..0c515f7b
--- /dev/null
+++ b/benchmarks/pandas/bench_xml.py
@@ -0,0 +1,65 @@
+"""Benchmark: read_xml / to_xml — parse and serialize XML
+
+Creates a 1,000-row XML document, then benchmarks:
+  - pd.read_xml (parse XML string → DataFrame)
+  - df.to_xml  (DataFrame → XML string)
+"""
+import json
+import time
+import io
+import numpy as np
+import pandas as pd
+
+ROWS = 1_000
+WARMUP = 3
+ITERATIONS = 20
+
+# Build XML string with ROWS row elements (matching TS benchmark)
+lines = ['<?xml version="1.0"?>', "<data>"]
+for i in range(ROWS):
+    lines.append(
+        f'  <row id="{i}" value="{i * 1.1:.4f}" label="cat_{i % 50}" />'
+    )
+lines.append("</data>")
+xml_string = "\n".join(lines)
+
+# Build a DataFrame for to_xml benchmarks
+df = pd.DataFrame(
+    {
+        "id": np.arange(ROWS, dtype=np.int64),
+        "value": np.arange(ROWS, dtype=np.float64) * 1.1,
+        "label": [f"cat_{i % 50}" for i in range(ROWS)],
+    }
+)
+
+# Warm up
+for _ in range(WARMUP):
+    pd.read_xml(io.StringIO(xml_string))
+    df.to_xml()
+
+# Benchmark read_xml
+t0 = time.perf_counter()
+for _ in range(ITERATIONS):
+    pd.read_xml(io.StringIO(xml_string))
+read_total = (time.perf_counter() - t0) * 1000
+
+# Benchmark to_xml
+t1 = time.perf_counter()
+for _ in range(ITERATIONS):
+    df.to_xml()
+write_total = (time.perf_counter() - t1) * 1000
+
+total = read_total + write_total
+
+print(
+    json.dumps(
+        {
+            "function": "xml",
+            "mean_ms": total / (ITERATIONS * 2),
+            "iterations": ITERATIONS * 2,
+            "total_ms": total,
+            "read_mean_ms": read_total / ITERATIONS,
+            "write_mean_ms": write_total / ITERATIONS,
+        }
+    )
+)
diff --git a/benchmarks/tsb/bench_case_when.ts b/benchmarks/tsb/bench_case_when.ts
new file mode 100644
index 00000000..3435c0d2
--- /dev/null
+++ b/benchmarks/tsb/bench_case_when.ts
@@ -0,0 +1,39 @@
+/**
+ * Benchmark: caseWhen — conditional value selection on 100k-element Series
+ */
+import { Series, caseWhen } from "../../src/index.js";
+
+const ROWS = 100_000;
+const WARMUP = 5;
+const ITERATIONS = 20;
+
+const data = Float64Array.from({ length: ROWS }, (_, i) => i % 100);
+const s = new Series(data);
+const cond1 = s.map((v) => (v as number) < 25);
+const cond2 = s.map((v) => (v as number) < 50);
+const cond3 = s.map((v) => (v as number) < 75);
+
+const caselist: [Series, string][] = [
+  [cond1 as Series, "low"],
+  [cond2 as Series, "medium-low"],
+  [cond3 as Series, "medium-high"],
+];
+
+for (let i = 0; i < WARMUP; i++) {
+  caseWhen(s, caselist);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  caseWhen(s, caselist);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "case_when",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_contingency.ts b/benchmarks/tsb/bench_contingency.ts
new file mode 100644
index 00000000..609776c4
--- /dev/null
+++ b/benchmarks/tsb/bench_contingency.ts
@@ -0,0 +1,50 @@
+/**
+ * Benchmark: contingency — expectedFreq, relativeRisk, oddsRatio, association
+ * Dataset: 4×4 contingency table built from 100,000 categorised observations.
+ */
+import { expectedFreq, relativeRisk, oddsRatio, association } from "../../src/index.js";
+
+const WARMUP = 10;
+const ITERS = 50;
+
+// Build a 4×4 observed count table
+const observed: readonly (readonly number[])[] = [
+  [120, 80, 40, 60],
+  [90, 110, 70, 30],
+  [50, 60, 100, 90],
+  [40, 50, 90, 120],
+];
+
+// 2×2 table for relativeRisk / oddsRatio (requires exactly 2 rows × 2 cols)
+const twoByTwo: readonly (readonly number[])[] = [
+  [60, 40],
+  [30, 70],
+];
+
+// Warm up
+for (let i = 0; i < WARMUP; i++) {
+  expectedFreq(observed);
+  relativeRisk(twoByTwo);
+  oddsRatio(twoByTwo);
+  association(observed, "cramer");
+}
+
+// Measure
+const start = performance.now();
+for (let i = 0; i < ITERS; i++) {
+  expectedFreq(observed);
+  relativeRisk(twoByTwo);
+  oddsRatio(twoByTwo);
+  association(observed, "cramer");
+}
+const total_ms = performance.now() - start;
+const mean_ms = total_ms / ITERS;
+
+console.log(
+  JSON.stringify({
+    function: "contingency",
+    mean_ms: parseFloat(mean_ms.toFixed(4)),
+    iterations: ITERS,
+    total_ms: parseFloat(total_ms.toFixed(4)),
+  }),
+);
diff --git a/benchmarks/tsb/bench_entropy.ts b/benchmarks/tsb/bench_entropy.ts
new file mode 100644
index 00000000..f4783389
--- /dev/null
+++ b/benchmarks/tsb/bench_entropy.ts
@@ -0,0 +1,31 @@
+import { entropy, klDivergence } from "../../src/index.js";
+
+const N = 100;
+const WARMUP = 5;
+const ITERS = 50;
+
+// Build two probability distributions of length N
+const p: number[] = Array.from({ length: N }, (_, i) => i + 1);
+const q: number[] = Array.from({ length: N }, (_, i) => N - i);
+
+let t0 = performance.now();
+for (let i = 0; i < WARMUP; i++) {
+  entropy(p);
+  klDivergence(p, q);
+}
+t0 = performance.now();
+
+for (let i = 0; i < ITERS; i++) {
+  entropy(p);
+  klDivergence(p, q);
+}
+const total_ms = performance.now() - t0;
+
+console.log(
+  JSON.stringify({
+    function: "entropy_klDivergence",
+    mean_ms: total_ms / ITERS,
+    iterations: ITERS,
+    total_ms,
+  }),
+);
diff --git a/benchmarks/tsb/bench_feather.ts b/benchmarks/tsb/bench_feather.ts
new file mode 100644
index 00000000..b701fa1b
--- /dev/null
+++ b/benchmarks/tsb/bench_feather.ts
@@ -0,0 +1,38 @@
+/**
+ * Benchmark: readFeather / toFeather — Arrow IPC Feather v2 round-trip on 10k rows
+ */
+import { DataFrame, toFeather, readFeather } from "../../src/index.js";
+
+const ROWS = 10_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Build a DataFrame with int, float, and string columns
+const ids = Array.from({ length: ROWS }, (_, i) => i);
+const values = Array.from({ length: ROWS }, (_, i) => i * 1.1);
+const labels = Array.from({ length: ROWS }, (_, i) => `cat_${i % 50}`);
+
+const df = new DataFrame({ id: ids, value: values, label: labels });
+
+// Warm up
+for (let i = 0; i < WARMUP; i++) {
+  const buf = toFeather(df);
+  readFeather(buf);
+}
+
+// Measure round-trip
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  const buf = toFeather(df);
+  readFeather(buf);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "feather",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_flags_options.ts b/benchmarks/tsb/bench_flags_options.ts
new file mode 100644
index 00000000..d17204ce
--- /dev/null
+++ b/benchmarks/tsb/bench_flags_options.ts
@@ -0,0 +1,76 @@
+/**
+ * Benchmark: flags and options
+ *
+ * Measures:
+ *  - getFlags / allowsDuplicateLabels get+set (Series & DataFrame)
+ *  - getOption / setOption / resetOption for multiple keys
+ *  - options proxy read
+ *
+ * Dataset: 10,000-row Series and DataFrame; 20 measured iterations.
+ */
+
+import {
+  Series,
+  DataFrame,
+  getFlags,
+  getOption,
+  setOption,
+  resetOption,
+  options,
+} from "../../src/index.js";
+
+const N = 10_000;
+const WARMUP = 5;
+const ITERS = 20;
+
+// Build test data once
+const data = Float64Array.from({ length: N }, (_, i) => i);
+const s = new Series(data);
+const df = DataFrame.fromColumns({ a: Array.from(data), b: Array.from(data) });
+
+function benchFlagsOptions(): number {
+  let sink = 0;
+  for (let i = 0; i < ITERS + WARMUP; i++) {
+    // flags on Series
+    const sf = getFlags(s);
+    const prev = sf.allowsDuplicateLabels;
+    sf.allowsDuplicateLabels = !prev;
+    sf.allowsDuplicateLabels = prev;
+    sink ^= sf.allowsDuplicateLabels ? 1 : 0;
+
+    // flags on DataFrame
+    const dff = getFlags(df);
+    const prevDf = dff.allowsDuplicateLabels;
+    dff.allowsDuplicateLabels = !prevDf;
+    dff.allowsDuplicateLabels = prevDf;
+    sink ^= dff.allowsDuplicateLabels ? 1 : 0;
+
+    // options get/set/reset
+    const v = getOption("display.max_rows") as number;
+    setOption("display.max_rows", v + 1);
+    resetOption("display.max_rows");
+    sink ^= (options.display as Record<string, unknown>).max_rows ? 1 : 0;
+
+    setOption("display.max_columns", 20);
+    resetOption("display.max_columns");
+    sink ^= (options.display as Record<string, unknown>).max_columns ? 1 : 0;
+  }
+  return sink;
+}
+
+// Warm-up
+benchFlagsOptions();
+
+// Measure
+const t0 = performance.now();
+for (let i = 0; i < ITERS; i++) benchFlagsOptions();
+const total = performance.now() - t0;
+
+console.log(
+  JSON.stringify({
+    function: "flags_options",
+    mean_ms: total / ITERS,
+    iterations: ITERS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_floating_array.ts b/benchmarks/tsb/bench_floating_array.ts
new file mode 100644
index 00000000..dedc4829
--- /dev/null
+++ b/benchmarks/tsb/bench_floating_array.ts
@@ -0,0 +1,45 @@
+/**
+ * Benchmark: FloatingArray — nullable float64 extension array operations.
+ * N=100_000 elements with ~10% nulls. Tests from/sum/mean/min/max/add/fillna.
+ */
+import { arrays } from "../../src/index.js";
+
+const N = 100_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Build input with ~10% nulls
+const raw: (number | null)[] = Array.from({ length: N }, (_, i) =>
+  i % 10 === 0 ? null : (i % 1000) * 0.001 - 0.5,
+);
+
+for (let i = 0; i < WARMUP; i++) {
+  const a = arrays.FloatingArray.from(raw, "Float64");
+  a.sum();
+  a.mean();
+  a.min();
+  a.max();
+  a.add(1.0);
+  a.fillna(0.0);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  const a = arrays.FloatingArray.from(raw, "Float64");
+  a.sum();
+  a.mean();
+  a.min();
+  a.max();
+  a.add(1.0);
+  a.fillna(0.0);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "floating_array",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_gaussianKDE.ts b/benchmarks/tsb/bench_gaussianKDE.ts
new file mode 100644
index 00000000..9e3f4de8
--- /dev/null
+++ b/benchmarks/tsb/bench_gaussianKDE.ts
@@ -0,0 +1,44 @@
+/**
+ * Benchmark: Gaussian KDE evaluate on 1000-point dataset at 200 grid points.
+ * Uses Silverman bandwidth (default).
+ */
+import { gaussianKDE } from "../../src/index.js";
+
+const N = 1_000;
+const GRID = 200;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Create dataset: mix of two gaussians
+const data: number[] = [];
+for (let i = 0; i < N; i++) {
+  const x = Math.sin(i * 0.01) * 2 + (i % 2 === 0 ? 0 : 5);
+  data.push(x);
+}
+
+// Grid points to evaluate KDE at
+const xmin = -5;
+const xmax = 10;
+const grid: number[] = Array.from({ length: GRID }, (_, i) => xmin + (i / (GRID - 1)) * (xmax - xmin));
+
+// Warm-up
+for (let i = 0; i < WARMUP; i++) {
+  const kde = gaussianKDE(data);
+  kde.evaluate(grid);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  const kde = gaussianKDE(data);
+  kde.evaluate(grid);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "gaussianKDE",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_hypothesis_tests.ts b/benchmarks/tsb/bench_hypothesis_tests.ts
new file mode 100644
index 00000000..05d36cb7
--- /dev/null
+++ b/benchmarks/tsb/bench_hypothesis_tests.ts
@@ -0,0 +1,41 @@
+import { ttest1samp, ttestInd, ttestRel, fOneway, pearsonr, spearmanr, mannWhitneyU } from "../../src/index.js";
+
+const WARMUP = 3;
+const ITERS = 20;
+const N = 1000;
+
+// Seeded deterministic data
+function makeData(n: number, seed: number): number[] {
+  const arr: number[] = [];
+  let x = seed;
+  for (let i = 0; i < n; i++) {
+    x = (x * 1664525 + 1013904223) & 0xffffffff;
+    arr.push((x >>> 0) / 0x100000000);
+  }
+  return arr;
+}
+
+const a = makeData(N, 42).map((v) => v * 4 + 2);
+const b = makeData(N, 99).map((v) => v * 4 + 2.5);
+
+function bench(): void {
+  let total = 0;
+  for (let i = 0; i < WARMUP + ITERS; i++) {
+    const t0 = performance.now();
+    ttest1samp(a, 2.5);
+    ttestInd(a, b);
+    ttestRel(a, b);
+    fOneway([a, b]);
+    pearsonr(a, b);
+    spearmanr(a, b);
+    mannWhitneyU(a, b);
+    const elapsed = performance.now() - t0;
+    if (i >= WARMUP) total += elapsed;
+  }
+  const mean_ms = total / ITERS;
+  console.log(
+    JSON.stringify({ function: "hypothesis_tests", mean_ms, iterations: ITERS, total_ms: total }),
+  );
+}
+
+bench();
diff --git a/benchmarks/tsb/bench_integer_array.ts b/benchmarks/tsb/bench_integer_array.ts
new file mode 100644
index 00000000..a2eaea3d
--- /dev/null
+++ b/benchmarks/tsb/bench_integer_array.ts
@@ -0,0 +1,45 @@
+/**
+ * Benchmark: IntegerArray — nullable integer extension array operations.
+ * N=100_000 elements with ~10% nulls. Tests from/sum/mean/min/max/add/fillna.
+ */
+import { arrays } from "../../src/index.js";
+
+const N = 100_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Build input with ~10% nulls
+const raw: (number | null)[] = Array.from({ length: N }, (_, i) =>
+  i % 10 === 0 ? null : (i % 1000) - 500,
+);
+
+for (let i = 0; i < WARMUP; i++) {
+  const a = arrays.IntegerArray.from(raw, "Int32");
+  a.sum();
+  a.mean();
+  a.min();
+  a.max();
+  a.add(1);
+  a.fillna(0);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  const a = arrays.IntegerArray.from(raw, "Int32");
+  a.sum();
+  a.mean();
+  a.min();
+  a.max();
+  a.add(1);
+  a.fillna(0);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "integer_array",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_linregress_polyfit.ts b/benchmarks/tsb/bench_linregress_polyfit.ts
new file mode 100644
index 00000000..b1f100be
--- /dev/null
+++ b/benchmarks/tsb/bench_linregress_polyfit.ts
@@ -0,0 +1,33 @@
+/**
+ * Benchmark: linregress and polyfit — simple linear regression and polynomial fit.
+ * Dataset: 10,000 points, 20 iterations.
+ */
+import { linregress, polyfit, polyval } from "../../src/index.js";
+
+const N = 10_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+const x = Array.from({ length: N }, (_, i) => i / N);
+const y = Array.from({ length: N }, (_, i) => 2.5 * (i / N) + 1.0 + Math.sin(i * 0.01) * 0.1);
+
+for (let i = 0; i < WARMUP; i++) {
+  linregress(x, y);
+  polyfit(x, y, 2);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  linregress(x, y);
+  polyfit(x, y, 2);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "linregress_polyfit",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_lreshape.ts b/benchmarks/tsb/bench_lreshape.ts
new file mode 100644
index 00000000..58f94459
--- /dev/null
+++ b/benchmarks/tsb/bench_lreshape.ts
@@ -0,0 +1,37 @@
+/**
+ * Benchmark: lreshape — wide-to-long reshape using named column groups.
+ * Dataset: 10,000 rows with 4 value columns (v1..v4), 50 iterations.
+ */
+import { DataFrame, lreshape } from "../../src/index.js";
+
+const N = 10_000;
+const WARMUP = 3;
+const ITERATIONS = 50;
+
+const ids = Array.from({ length: N }, (_, i) => i);
+const v1 = Array.from({ length: N }, (_, i) => i * 1.0);
+const v2 = Array.from({ length: N }, (_, i) => i * 2.0);
+const v3 = Array.from({ length: N }, (_, i) => i * 3.0);
+const v4 = Array.from({ length: N }, (_, i) => i * 4.0);
+
+const df = DataFrame.fromColumns({ id: ids, v1, v2, v3, v4 });
+const groups = { value: ["v1", "v2", "v3", "v4"] };
+
+for (let i = 0; i < WARMUP; i++) {
+  lreshape(df, groups);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  lreshape(df, groups);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "lreshape",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_multivariate.ts b/benchmarks/tsb/bench_multivariate.ts
new file mode 100644
index 00000000..199ec9e7
--- /dev/null
+++ b/benchmarks/tsb/bench_multivariate.ts
@@ -0,0 +1,49 @@
+/**
+ * Benchmark: multivariate statistics — mahalanobis distance, covMatrix, PCA
+ * Dataset: 500 observations × 5 features (realistic small-to-medium size)
+ */
+import { mahalanobis, covMatrix, PCA } from "../../src/index.js";
+
+const N = 500;
+const P = 5;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Generate a deterministic dataset
+const X: number[][] = Array.from({ length: N }, (_, i) =>
+  Array.from({ length: P }, (_, j) => Math.sin(i * 0.1 + j) * 10 + j * 2),
+);
+
+const u = X[0]!;
+const v = X[1]!;
+
+// Pre-compute inverse covariance for mahalanobis
+const cov = covMatrix(X);
+// Simple diagonal approximation for VI (invertMatrix is tested via mahalanobis internals)
+const VI: number[][] = Array.from({ length: P }, (_, i) =>
+  Array.from({ length: P }, (_, j) => (i === j ? 1 / Math.max(cov[i]![i]!, 1e-10) : 0)),
+);
+
+// Warm up
+for (let i = 0; i < WARMUP; i++) {
+  mahalanobis(u, v, VI);
+  covMatrix(X);
+  new PCA({ n_components: 3 }).fit(X);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  mahalanobis(u, v, VI);
+  covMatrix(X);
+  new PCA({ n_components: 3 }).fit(X);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "multivariate",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_mutual_information.ts b/benchmarks/tsb/bench_mutual_information.ts
new file mode 100644
index 00000000..0f8df611
--- /dev/null
+++ b/benchmarks/tsb/bench_mutual_information.ts
@@ -0,0 +1,34 @@
+import { mutualInformation, normalizedMI } from "../../src/index.js";
+
+const N = 1000;
+const WARMUP = 5;
+const ITERS = 50;
+
+// Build paired observations: two correlated categorical variables (10 categories each)
+const CATS = 10;
+const pairs: [number, number][] = Array.from({ length: N }, (_, i) => [
+  i % CATS,
+  (i % CATS) + Math.floor(i / CATS) % 3,
+]);
+
+let t0 = performance.now();
+for (let i = 0; i < WARMUP; i++) {
+  mutualInformation(pairs);
+  normalizedMI(pairs, "arithmetic");
+}
+t0 = performance.now();
+
+for (let i = 0; i < ITERS; i++) {
+  mutualInformation(pairs);
+  normalizedMI(pairs, "arithmetic");
+}
+const total_ms = performance.now() - t0;
+
+console.log(
+  JSON.stringify({
+    function: "mutual_information",
+    mean_ms: total_ms / ITERS,
+    iterations: ITERS,
+    total_ms,
+  }),
+);
diff --git a/benchmarks/tsb/bench_pipe_apply.ts b/benchmarks/tsb/bench_pipe_apply.ts
new file mode 100644
index 00000000..29213f06
--- /dev/null
+++ b/benchmarks/tsb/bench_pipe_apply.ts
@@ -0,0 +1,45 @@
+/**
+ * Benchmark: pipe / seriesApply / dataFrameApplyMap on 10,000-row datasets.
+ *
+ * Exercises three functional-pipeline operations:
+ *   - pipe: chain 4 transforms on a Series
+ *   - seriesApply: element-wise function on 10k-element Series
+ *   - dataFrameApplyMap: element-wise function on 10k × 3 DataFrame
+ */
+import { Series, DataFrame, pipe, seriesApply, dataFrameApplyMap } from "../../src/index.js";
+
+const N = 10_000;
+const WARMUP = 5;
+const ITERATIONS = 20;
+
+const raw = Float64Array.from({ length: N }, (_, i) => (i % 100) + 1);
+const series = new Series(raw, { name: "x" });
+const df = DataFrame.fromColumns({
+  a: Array.from({ length: N }, (_, i) => (i % 50) + 1),
+  b: Array.from({ length: N }, (_, i) => (i % 30) + 1),
+  c: Array.from({ length: N }, (_, i) => (i % 20) + 1),
+});
+
+// Warm-up
+for (let i = 0; i < WARMUP; i++) {
+  pipe(series, (s) => s.add(1), (s) => s.mul(2), (s) => s.sub(1), (s) => s.div(2));
+  seriesApply(series, (v) => (v as number) * 2 + 1);
+  dataFrameApplyMap(df, (v) => (v as number) * 2);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  pipe(series, (s) => s.add(1), (s) => s.mul(2), (s) => s.sub(1), (s) => s.div(2));
+  seriesApply(series, (v) => (v as number) * 2 + 1);
+  dataFrameApplyMap(df, (v) => (v as number) * 2);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "pipe_apply",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_readStata.ts b/benchmarks/tsb/bench_readStata.ts
new file mode 100644
index 00000000..a6703324
--- /dev/null
+++ b/benchmarks/tsb/bench_readStata.ts
@@ -0,0 +1,42 @@
+/**
+ * Benchmark: readStata / toStata round-trip on a 10k-row DataFrame
+ */
+import { DataFrame, Series, readStata, toStata } from "../../src/index.js";
+
+const ROWS = 10_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Build a DataFrame with numeric and string columns
+const ids = Int32Array.from({ length: ROWS }, (_, i) => i);
+const values = Float64Array.from({ length: ROWS }, (_, i) => Math.sin(i * 0.01) * 1000);
+const categories = Array.from({ length: ROWS }, (_, i) => `cat_${i % 5}`);
+
+const df = new DataFrame({
+  id: new Series(ids),
+  value: new Series(values),
+  category: new Series(categories),
+});
+
+// Serialize once so readStata benchmarks read from a pre-built buffer
+const buf = toStata(df);
+
+// Warm up
+for (let i = 0; i < WARMUP; i++) {
+  readStata(buf);
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  readStata(buf);
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "readStata",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_read_fwf.ts b/benchmarks/tsb/bench_read_fwf.ts
new file mode 100644
index 00000000..f6d7eace
--- /dev/null
+++ b/benchmarks/tsb/bench_read_fwf.ts
@@ -0,0 +1,39 @@
+/**
+ * Benchmark: readFwf — parse a fixed-width formatted text file into a DataFrame.
+ * Dataset: 10,000 rows × 4 columns (id, name, value, flag).
+ */
+import { readFwf } from "../../src/index.js";
+
+const ROWS = 10_000;
+const WARMUP = 3;
+const ITERATIONS = 10;
+
+// Build a fixed-width text: id(6), name(10), value(10), flag(4)
+const lines: string[] = ["id    name      value     flag"];
+for (let i = 0; i < ROWS; i++) {
+  const id = String(i).padStart(6);
+  const name = ("item" + (i % 500)).padEnd(10);
+  const value = (Math.sin(i * 0.01) * 1000).toFixed(3).padStart(10);
+  const flag = (i % 2 === 0 ? "Y" : "N").padEnd(4);
+  lines.push(id + name + value + flag);
+}
+const text = lines.join("\n");
+
+for (let i = 0; i < WARMUP; i++) {
+  readFwf(text, { colspecs: [[0, 6], [6, 16], [16, 26], [26, 30]] });
+}
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  readFwf(text, { colspecs: [[0, 6], [6, 16], [16, 26], [26, 30]] });
+}
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "readFwf",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_read_sas.ts b/benchmarks/tsb/bench_read_sas.ts
new file mode 100644
index 00000000..dd873ef7
--- /dev/null
+++ b/benchmarks/tsb/bench_read_sas.ts
@@ -0,0 +1,148 @@
+/**
+ * Benchmark: readSas — parse a 1,000-row SAS XPORT (XPT) file.
+ * Outputs JSON: {"function": "read_sas", "mean_ms": ..., "iterations": ..., "total_ms": ...}
+ */
+import { readSas } from "../../src/index.ts";
+
+// ─── IBM 370 floating-point encoder ──────────────────────────────────────────
+
+function ibmEncode(val: number): Uint8Array {
+  const out = new Uint8Array(8);
+  if (val === 0) return out;
+  if (!Number.isFinite(val)) { out[0] = 0x2e; return out; }
+  const sign = val < 0 ? 1 : 0;
+  const abs = Math.abs(val);
+  let exp = 0;
+  let mant = abs;
+  while (mant >= 1) { mant /= 16; exp++; }
+  while (mant < 1 / 16 && mant > 0) { mant *= 16; exp--; }
+  const mantInt = BigInt(Math.round(mant * 2 ** 56));
+  out[0] = (sign << 7) | ((exp + 64) & 0x7f);
+  for (let i = 1; i <= 7; i++) {
+    out[i] = Number((mantInt >> BigInt((7 - i) * 8)) & 0xffn);
+  }
+  return out;
+}
+
+// ─── Minimal XPORT v5 builder ────────────────────────────────────────────────
+
+function buildXpt(
+  numVars: readonly string[],
+  charVars: readonly { name: string; len: number }[],
+  rows: readonly Readonly<Record<string, number | string>>[],
+): Uint8Array {
+  const RECORD = 80;
+
+  function encodeAscii(s: string, maxLen: number): Uint8Array {
+    const buf = new Uint8Array(maxLen);
+    for (let i = 0; i < Math.min(s.length, maxLen); i++) buf[i] = s.charCodeAt(i) & 0x7f;
+    return buf;
+  }
+  function padTo80(s: string): Uint8Array { return encodeAscii(s.padEnd(RECORD, " "), RECORD); }
+  function writeU16(b: Uint8Array, o: number, v: number) { b[o] = (v >> 8) & 0xff; b[o + 1] = v & 0xff; }
+  function writeU32(b: Uint8Array, o: number, v: number) {
+    b[o] = (v >> 24) & 0xff; b[o + 1] = (v >> 16) & 0xff; b[o + 2] = (v >> 8) & 0xff; b[o + 3] = v & 0xff;
+  }
+
+  type Meta = { type: 1 | 2; name: string; len: number; pos: number };
+  const metas: Meta[] = [];
+  let pos = 0;
+  for (const name of numVars) { metas.push({ type: 1, name, len: 8, pos }); pos += 8; }
+  for (const { name, len } of charVars) { metas.push({ type: 2, name, len, pos }); pos += len; }
+  const rowLen = pos;
+
+  const chunks: Uint8Array[] = [];
+
+  // Library header (5 × 80 bytes)
+  chunks.push(padTo80("HEADER RECORD*******LIBRARY HEADER RECORD!!!!!!!000000000000000000000000000000  "));
+  chunks.push(padTo80("SAS     SAS     SASLIB  6.06    ASCII"));
+  chunks.push(padTo80("20240101"));
+  chunks.push(padTo80(""));
+  chunks.push(padTo80(""));
+
+  // Member header (3 × 80 bytes)
+  chunks.push(padTo80("HEADER RECORD*******MEMBER  HEADER RECORD!!!!!!!000000000000000000000000000001600000000140  "));
+  chunks.push(padTo80("SAS     BENCH   SASDATA 6.06    ASCII"));
+  chunks.push(padTo80(""));
+
+  // Namestr header
+  const nvar = metas.length;
+  chunks.push(padTo80(`HEADER RECORD*******NAMESTR HEADER RECORD!!!!!!!${String(nvar).padStart(6, "0")}00000000000000000000  `));
+
+  // Namestr records (140 bytes each)
+  const nsBuf = new Uint8Array(nvar * 140);
+  for (let i = 0; i < metas.length; i++) {
+    const m = metas[i]!;
+    const off = i * 140;
+    writeU16(nsBuf, off, m.type);
+    writeU16(nsBuf, off + 2, 140);
+    nsBuf.set(encodeAscii(m.name, 8), off + 4);
+    writeU16(nsBuf, off + 52, m.len);
+    writeU32(nsBuf, off + 84, m.pos);
+  }
+  const nsPadded = Math.ceil(nsBuf.length / RECORD) * RECORD;
+  const nsPaddedBuf = new Uint8Array(nsPadded);
+  nsPaddedBuf.set(nsBuf);
+  chunks.push(nsPaddedBuf);
+
+  // Obs header
+  chunks.push(padTo80("HEADER RECORD*******OBS     HEADER RECORD!!!!!!!000000000000000000000000000000  "));
+
+  // Observations
+  const paddedRowLen = Math.ceil(rowLen / RECORD) * RECORD;
+  const obsBuf = new Uint8Array(rows.length * paddedRowLen);
+  for (let r = 0; r < rows.length; r++) {
+    const base = r * paddedRowLen;
+    const row = rows[r]!;
+    for (const m of metas) {
+      const val = row[m.name];
+      if (m.type === 1) {
+        const encoded = ibmEncode(typeof val === "number" ? val : 0);
+        obsBuf.set(encoded, base + m.pos);
+      } else {
+        const s = typeof val === "string" ? val : "";
+        obsBuf.set(encodeAscii(s, m.len), base + m.pos);
+      }
+    }
+  }
+  chunks.push(obsBuf);
+
+  let total = 0;
+  for (const c of chunks) total += c.length;
+  const out = new Uint8Array(total);
+  let off = 0;
+  for (const c of chunks) { out.set(c, off); off += c.length; }
+  return out;
+}
+
+// ─── Build dataset ────────────────────────────────────────────────────────────
+
+const ROWS = 1_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+const rows: Readonly<Record<string, number | string>>[] = Array.from({ length: ROWS }, (_, i) => ({
+  id: i,
+  value: i * 1.5,
+  score: Math.sin(i * 0.01),
+  label: `item_${i % 100}`,
+}));
+
+const xpt = buildXpt(["id", "value", "score"], [{ name: "label", len: 12 }], rows);
+
+// ─── Benchmark ────────────────────────────────────────────────────────────────
+
+for (let i = 0; i < WARMUP; i++) readSas(xpt);
+
+const start = performance.now();
+for (let i = 0; i < ITERATIONS; i++) readSas(xpt);
+const total = performance.now() - start;
+
+console.log(
+  JSON.stringify({
+    function: "read_sas",
+    mean_ms: total / ITERATIONS,
+    iterations: ITERATIONS,
+    total_ms: total,
+  }),
+);
diff --git a/benchmarks/tsb/bench_sql.ts b/benchmarks/tsb/bench_sql.ts
new file mode 100644
index 00000000..5889d38c
--- /dev/null
+++ b/benchmarks/tsb/bench_sql.ts
@@ -0,0 +1,85 @@
+/**
+ * Benchmark: readSqlQuery / toSql on 10k-row result sets
+ */
+import { DataFrame, readSqlQuery, toSql } from "../../src/index.js";
+import type { SqlConnection, SqlResult, SqlRow, SqlValue, IfExistsStrategy } from "../../src/index.js";
+
+const ROWS = 10_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// ── Shared result set ─────────────────────────────────────────────────────────
+const columns: string[] = ["id", "value", "label"];
+const rows: SqlRow[] = Array.from({ length: ROWS }, (_, i) => ({
+  id: i,
+  value: Math.sin(i * 0.01) * 1000,
+  label: `item_${i % 100}`,
+}));
+
+// ── Mock adapter for reads ────────────────────────────────────────────────────
+class ReadAdapter implements SqlConnection {
+  query(_sql: string, _params?: readonly SqlValue[]): SqlResult {
+    return { columns, rows };
+  }
+  listTables(): readonly string[] {
+    return ["mock_table"];
+  }
+}
+
+const readConn = new ReadAdapter();
+
+// ── Warm-up reads ─────────────────────────────────────────────────────────────
+for (let i = 0; i < WARMUP; i++) {
+  readSqlQuery("SELECT * FROM mock_table", readConn);
+}
+
+// ── readSqlQuery benchmark ────────────────────────────────────────────────────
+const startRead = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  readSqlQuery("SELECT * FROM mock_table", readConn);
+}
+const totalRead = performance.now() - startRead;
+
+// ── Mock adapter for writes ───────────────────────────────────────────────────
+class WriteAdapter implements SqlConnection {
+  query(_sql: string, _params?: readonly SqlValue[]): SqlResult {
+    return { columns: [], rows: [] };
+  }
+  listTables(): readonly string[] {
+    return [];
+  }
+  insert(
+    _tableName: string,
+    _rows: readonly SqlRow[],
+    _columns: readonly string[],
+    _ifExists: IfExistsStrategy,
+  ): number {
+    return _rows.length;
+  }
+}
+
+const writeConn = new WriteAdapter();
+const df = readSqlQuery("SELECT * FROM mock_table", readConn);
+
+// ── Warm-up writes ────────────────────────────────────────────────────────────
+for (let i = 0; i < WARMUP; i++) {
+  toSql(df, "bench_table", writeConn, { ifExists: "replace" });
+}
+
+// ── toSql benchmark ───────────────────────────────────────────────────────────
+const startWrite = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  toSql(df, "bench_table", writeConn, { ifExists: "replace" });
+}
+const totalWrite = performance.now() - startWrite;
+
+console.log(
+  JSON.stringify({
+    function: "sql",
+    mean_ms: (totalRead + totalWrite) / (2 * ITERATIONS),
+    iterations: ITERATIONS,
+    total_ms: totalRead + totalWrite,
+    read_mean_ms: totalRead / ITERATIONS,
+    write_mean_ms: totalWrite / ITERATIONS,
+  }),
+);
diff --git a/benchmarks/tsb/bench_to_excel.ts b/benchmarks/tsb/bench_to_excel.ts
new file mode 100644
index 00000000..941535c2
--- /dev/null
+++ b/benchmarks/tsb/bench_to_excel.ts
@@ -0,0 +1,38 @@
+/**
+ * Benchmark: toExcel — serialize a DataFrame to an XLSX binary buffer.
+ * Outputs JSON: {"function": "to_excel", "mean_ms": ..., "iterations": ..., "total_ms": ...}
+ */
+import { DataFrame, toExcel } from "../../src/index.ts";
+
+const ROWS = 5_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+const index = Array.from({ length: ROWS }, (_, i) => i);
+const colA = Array.from({ length: ROWS }, (_, i) => `name_${i % 1000}`);
+const colB = Array.from({ length: ROWS }, (_, i) => i * 1.5);
+const colC = Array.from({ length: ROWS }, (_, i) => i % 2 === 0);
+
+const df = new DataFrame({ name: colA, value: colB, flag: colC }, { index });
+
+for (let i = 0; i < WARMUP; i++) {
+  toExcel(df);
+}
+
+const times: number[] = [];
+for (let i = 0; i < ITERATIONS; i++) {
+  const start = performance.now();
+  toExcel(df);
+  times.push(performance.now() - start);
+}
+
+const totalMs = times.reduce((a, b) => a + b, 0);
+const meanMs = totalMs / ITERATIONS;
+console.log(
+  JSON.stringify({
+    function: "to_excel",
+    mean_ms: Math.round(meanMs * 1000) / 1000,
+    iterations: ITERATIONS,
+    total_ms: Math.round(totalMs * 1000) / 1000,
+  }),
+);
diff --git a/benchmarks/tsb/bench_xml.ts b/benchmarks/tsb/bench_xml.ts
new file mode 100644
index 00000000..a10f1eb6
--- /dev/null
+++ b/benchmarks/tsb/bench_xml.ts
@@ -0,0 +1,65 @@
+/**
+ * Benchmark: readXml / toXml — parse and serialize XML
+ *
+ * Creates a 1,000-row XML document, then benchmarks:
+ *   - readXml (parse XML string → DataFrame)
+ *   - toXml (DataFrame → XML string)
+ */
+import { readXml, toXml, DataFrame, Series } from "../../src/index.js";
+
+const ROWS = 1_000;
+const WARMUP = 3;
+const ITERATIONS = 20;
+
+// Build XML string with ROWS row elements
+const lines: string[] = ['<?xml version="1.0"?>', "<data>"];
+for (let i = 0; i < ROWS; i++) {
+  lines.push(
+    `  <row id="${i}" value="${(i * 1.1).toFixed(4)}" label="cat_${i % 50}" />`,
+  );
+}
+lines.push("</data>");
+const xmlString = lines.join("\n");
+
+// Build a DataFrame for toXml benchmarks
+const ids = Array.from({ length: ROWS }, (_, i) => i);
+const values = Array.from({ length: ROWS }, (_, i) => i * 1.1);
+const labels = Array.from({ length: ROWS }, (_, i) => `cat_${i % 50}`);
+const df = new DataFrame({
+  id: new Series(ids),
+  value: new Series(values),
+  label: new Series(labels),
+});
+
+// Warm up
+for (let i = 0; i < WARMUP; i++) {
+  readXml(xmlString);
+  toXml(df);
+}
+
+// Benchmark readXml
+const t0 = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  readXml(xmlString);
+}
+const readTotal = performance.now() - t0;
+
+// Benchmark toXml
+const t1 = performance.now();
+for (let i = 0; i < ITERATIONS; i++) {
+  toXml(df);
+}
+const writeTotal = performance.now() - t1;
+
+const total = readTotal + writeTotal;
+
+console.log(
+  JSON.stringify({
+    function: "xml",
+    mean_ms: total / (ITERATIONS * 2),
+    iterations: ITERATIONS * 2,
+    total_ms: total,
+    read_mean_ms: readTotal / ITERATIONS,
+    write_mean_ms: writeTotal / ITERATIONS,
+  }),
+);