diff --git a/benchmarks/benchmark_training_throughput.py b/benchmarks/benchmark_training_throughput.py
index 46b28100a2..a2d3e87dc2 100644
--- a/benchmarks/benchmark_training_throughput.py
+++ b/benchmarks/benchmark_training_throughput.py
@@ -20,6 +20,7 @@
 from transformers.optimization import get_cosine_schedule_with_warmup
 
 import fla
+from benchmarks.distributions import sample_lognormal_packed_lengths
 
 classes = [getattr(fla.models, i) for i in fla.models.__all__]
 configs = {i.model_type: i() for i in classes if issubclass(i, PretrainedConfig)}
@@ -66,19 +67,46 @@ def prepare_inputs(
     varlen: bool,
     vocab_size: int,
     device: torch.device,
+    length_distribution: str = 'random',
+    num_sequences: int | None = None,
+    length_sigma: float = 1.0,
+    generator: torch.Generator | None = None,
 ):
     if varlen:
-        tokens = torch.randint(high=vocab_size, size=(1, batch_size * seq_len), device=device)
-        cu_seqlens = torch.cat([
-            torch.tensor([0]),
-            torch.randperm(batch_size * seq_len - 16)[:torch.randint(8, 64, size=(1,))] + 16,
-            torch.tensor([batch_size * seq_len]),
-        ], 0).sort()[0].to(dtype=torch.int32, device=device)
-        if context_len is not None:
-            cu_seqlens = torch.cat(
-                [torch.arange(i, j, context_len) for i, j in zip(cu_seqlens[:-1].tolist(), cu_seqlens[1:].tolist())] +
-                [torch.tensor([len(tokens[0])])],
-            ).to(dtype=torch.int32, device=device)
+        total_tokens = batch_size * seq_len
+        tokens = torch.randint(high=vocab_size, size=(1, total_tokens), device=device)
+        if length_distribution == 'random':
+            num_cuts = int(torch.randint(8, 64, size=(1,), generator=generator))
+            cut_points = torch.randperm(total_tokens - 16, generator=generator)[:num_cuts] + 16
+            cu_seqlens = torch.cat([
+                torch.tensor([0]),
+                cut_points,
+                torch.tensor([total_tokens]),
+            ], 0).sort()[0]
+            if context_len is not None:
+                cu_seqlens = torch.cat(
+                    [torch.arange(i, j, context_len) for i, j in zip(cu_seqlens[:-1].tolist(), cu_seqlens[1:].tolist())] +
+                    [torch.tensor([total_tokens])],
+                )
+        elif length_distribution == 'lognormal':
+            max_length = context_len or seq_len
+            min_sequences = (total_tokens + max_length - 1) // max_length
+            max_sequences = total_tokens // 16
+            if num_sequences is None:
+                default_sequences = min(64, max(8, batch_size * 4))
+                num_sequences = min(max(default_sequences, min_sequences), max_sequences)
+            lengths = sample_lognormal_packed_lengths(
+                total_tokens=total_tokens,
+                num_sequences=num_sequences,
+                max_length=max_length,
+                min_length=16,
+                sigma=length_sigma,
+                generator=generator,
+            )
+            cu_seqlens = torch.cat([torch.zeros(1, dtype=torch.long), lengths.cumsum(0)])
+        else:
+            raise ValueError(f"unsupported length_distribution: {length_distribution!r}")
+        cu_seqlens = cu_seqlens.to(dtype=torch.int32, device=device)
     else:
         tokens = torch.randint(high=vocab_size, size=(batch_size, seq_len), device=device)
         cu_seqlens = None
@@ -106,8 +134,13 @@ def profile(
     enable_profile: bool = False,
     profile_steps: int = 64,
     profile_trace: str | None = None,
+    length_distribution: str = 'random',
+    num_sequences: int | None = None,
+    length_sigma: float = 1.0,
+    seed: int = 42,
 ):
     device = torch.device('cuda')
+    torch.manual_seed(seed)
     config = configs[name] if name in configs else AutoConfig.from_pretrained(name)
     if num_heads is not None:
         if not hasattr(config, 'num_heads'):
@@ -149,7 +182,8 @@ def _mark(override):
     _print_run_header({
         'model':    name,
         'arch':     ' '.join(arch_parts),
-        'data':     f"B={batch_size} T={seq_len} ctx={context_len} varlen={varlen}",
+        'data':     f"B={batch_size} T={seq_len} ctx={context_len} varlen={varlen} "
+                    f"lengths={length_distribution} n={num_sequences} sigma={length_sigma} seed={seed}",
         'training': f"{dtype} (mixed={mixed_precision}) compile={compile} "
         f"warmup={warmup_steps} steps={steps}",
         'profile':  profile_str,
@@ -178,6 +212,7 @@ def _mark(override):
     bar = trange(warmup_steps)
 
     model, optimizer, scheduler = accelerator.prepare(model, optimizer, scheduler)
+    length_generator = torch.Generator().manual_seed(seed)
     torch.cuda.synchronize(device)
     for _ in bar:
         # forward pass
@@ -188,6 +223,10 @@ def _mark(override):
             varlen=varlen,
             vocab_size=config.vocab_size,
             device=device,
+            length_distribution=length_distribution,
+            num_sequences=num_sequences,
+            length_sigma=length_sigma,
+            generator=length_generator,
         )
         outputs = model(tokens, labels=tokens, cu_seqlens=cu_seqlens)
         # backward pass
@@ -209,6 +248,10 @@ def _mark(override):
             varlen=varlen,
             vocab_size=config.vocab_size,
             device=device,
+            length_distribution=length_distribution,
+            num_sequences=num_sequences,
+            length_sigma=length_sigma,
+            generator=length_generator,
         )
         outputs = model(tokens, labels=tokens, cu_seqlens=cu_seqlens)
         # backward pass
@@ -246,6 +289,10 @@ def _mark(override):
                         varlen=varlen,
                         vocab_size=config.vocab_size,
                         device=device,
+                        length_distribution=length_distribution,
+                        num_sequences=num_sequences,
+                        length_sigma=length_sigma,
+                        generator=length_generator,
                     )
                     outputs = model(tokens, labels=tokens, cu_seqlens=cu_seqlens)
                     accelerator.backward(outputs.loss)
@@ -276,6 +323,10 @@ def _mark(override):
     parser.add_argument("--seq_len", default=4096, type=int)
     parser.add_argument("--context_len", default=None, type=int)
     parser.add_argument("--varlen", action='store_true')
+    parser.add_argument("--length_distribution", choices=['random', 'lognormal'], default='random')
+    parser.add_argument("--num_sequences", default=None, type=int)
+    parser.add_argument("--length_sigma", default=1.0, type=float)
+    parser.add_argument("--seed", default=42, type=int)
     parser.add_argument("--num_heads", default=None, type=int)
     parser.add_argument("--head_dim", default=None, type=int)
     parser.add_argument("--num_hidden_layers", default=None, type=int)
@@ -291,6 +342,10 @@ def _mark(override):
         seq_len=args.seq_len,
         context_len=args.context_len,
         varlen=args.varlen,
+        length_distribution=args.length_distribution,
+        num_sequences=args.num_sequences,
+        length_sigma=args.length_sigma,
+        seed=args.seed,
         num_heads=args.num_heads,
         head_dim=args.head_dim,
         num_hidden_layers=args.num_hidden_layers,
diff --git a/benchmarks/distributions.py b/benchmarks/distributions.py
new file mode 100644
index 0000000000..9edf769a09
--- /dev/null
+++ b/benchmarks/distributions.py
@@ -0,0 +1,61 @@
+# Copyright (c) 2023-2026, Songlin Yang, Yu Zhang, Zhiyuan Li
+#
+# This source code is licensed under the MIT license found in the
+# LICENSE file in the root directory of this source tree.
+# For a list of all contributors, visit:
+#   https://github.com/fla-org/flash-linear-attention/graphs/contributors
+
+from __future__ import annotations
+
+import torch
+
+
+def sample_lognormal_packed_lengths(
+    total_tokens: int,
+    num_sequences: int,
+    max_length: int,
+    min_length: int = 16,
+    sigma: float = 1.0,
+    generator: torch.Generator | None = None,
+) -> torch.Tensor:
+    """Sample right-skewed sequence lengths with an exact packed-token budget.
+
+    The log-normal weights are a configurable workload proxy, not a claim about any particular dataset.
+    Pass an observed sequence count, length bounds, and sigma from the workload being measured.
+    """
+    if total_tokens <= 0:
+        raise ValueError(f"total_tokens must be positive, got {total_tokens}")
+    if num_sequences <= 0:
+        raise ValueError(f"num_sequences must be positive, got {num_sequences}")
+    if not 0 < min_length <= max_length:
+        raise ValueError(f"expected 0 < min_length <= max_length, got {min_length=} and {max_length=}")
+    if sigma <= 0:
+        raise ValueError(f"sigma must be positive, got {sigma}")
+
+    min_tokens = num_sequences * min_length
+    max_tokens = num_sequences * max_length
+    if not min_tokens <= total_tokens <= max_tokens:
+        raise ValueError(
+            f"cannot pack {total_tokens} tokens into {num_sequences} sequences with "
+            f"lengths in [{min_length}, {max_length}]"
+        )
+
+    weights = torch.empty(num_sequences, dtype=torch.float64).log_normal_(
+        mean=0.0,
+        std=sigma,
+        generator=generator,
+    )
+    lengths = torch.full((num_sequences,), min_length, dtype=torch.long)
+    capacity = torch.full_like(lengths, max_length - min_length)
+    remaining = total_tokens - min_tokens
+
+    while remaining:
+        active_weights = weights.masked_fill(capacity == 0, 0)
+        picks = torch.multinomial(active_weights, remaining, replacement=True, generator=generator)
+        requested = torch.bincount(picks, minlength=num_sequences)
+        added = torch.minimum(requested, capacity)
+        lengths += added
+        capacity -= added
+        remaining -= int(added.sum())
+
+    return lengths
diff --git a/benchmarks/ops/registry.py b/benchmarks/ops/registry.py
index 34e267aa16..c8133a466b 100644
--- a/benchmarks/ops/registry.py
+++ b/benchmarks/ops/registry.py
@@ -14,6 +14,7 @@
 from __future__ import annotations
 
 import logging
+import math
 from collections.abc import Callable
 from dataclasses import dataclass, field
 from typing import Any
@@ -81,6 +82,42 @@ def logsigmoid_clamp(t):
     return F.logsigmoid(t).clamp_min(-5)
 
 
+def silu_transform(t):
+    return F.silu(t)
+
+
+def normalized_logsigmoid(t, normalizer: float = 16.0):
+    return F.logsigmoid(t) / normalizer
+
+
+def _inv_softplus(t):
+    return t + torch.log(-torch.expm1(-t))
+
+
+def learned_decay_transform(t):
+    """Match the activated decay range used by GDN and Comba layers."""
+    feature_shape = t.shape[2:]
+    log_dt = torch.empty(feature_shape, dtype=torch.float32, device=t.device).uniform_(math.log(0.001), math.log(0.1))
+    dt = log_dt.exp()
+    dt_bias = _inv_softplus(dt)
+    scale_shape = (t.shape[2],) + (1,) * (t.ndim - 3)
+    scale = torch.empty(scale_shape, dtype=torch.float32, device=t.device).uniform_(1, 16)
+    return (-scale * F.softplus(t.float() + dt_bias)).to(t.dtype)
+
+
+def safe_decay_transform(t, lower_bound: float = -5.0):
+    """Match KDA's bounded gate with its log-uniform step-size initialization."""
+    feature_shape = t.shape[2:]
+    log_dt = torch.empty(feature_shape, dtype=torch.float32, device=t.device).uniform_(
+        math.log(0.001),
+        math.log(0.1),
+    )
+    dt_bias = _inv_softplus(log_dt.exp())
+    scale_shape = (t.shape[2],) + (1,) * (t.ndim - 3)
+    scale = torch.empty(scale_shape, dtype=torch.float32, device=t.device).uniform_(1, 16)
+    return (lower_bound * torch.sigmoid(scale * (t.float() + dt_bias))).to(t.dtype)
+
+
 RWKV7_W_MIN = -0.6065306597126334
 
 
@@ -103,11 +140,13 @@ class TensorSpec:
         requires_grad:  whether the tensor needs gradients
         dtype:          'default' inherits from the benchmark, or 'float32'/'long'
         transform:      applied after randn, e.g. F.logsigmoid
+        realistic_transform: layer-shaped transform selected by the realistic input profile
     """
     shape_fn: Callable
     requires_grad: bool = True
     dtype: str = 'default'
     transform: Callable | None = None
+    realistic_transform: Callable | None = None
 
 
 # ---------------------------------------------------------------------------
@@ -216,6 +255,7 @@ def generate_inputs(
     B: int, T: int, H: int, D: int,
     dtype: torch.dtype = torch.bfloat16,
     device: str | torch.device = 'cuda',
+    input_profile: str = 'synthetic',
     **extra_shape_kw,
 ) -> dict[str, torch.Tensor]:
     """Create input tensors for *config* at the given shape.
@@ -223,6 +263,9 @@ def generate_inputs(
     Returns a dict mapping parameter names to tensors.
     Raises ValueError if dim_constraints are not satisfied (caller should skip).
     """
+    if input_profile not in ('synthetic', 'realistic'):
+        raise ValueError(f"input_profile must be 'synthetic' or 'realistic', got {input_profile!r}")
+
     # Check dim constraints
     if config.dim_constraints:
         shape_vals = {'B': B, 'T': T, 'H': H, 'D': D, **extra_shape_kw}
@@ -252,8 +295,10 @@ def generate_inputs(
         else:
             tensor = torch.randn(shape, dtype=tensor_dtype, device=device)
 
-        if spec.transform is not None:
-            tensor = spec.transform(tensor)
+        transform = spec.realistic_transform if input_profile == 'realistic' else None
+        transform = transform or spec.transform
+        if transform is not None:
+            tensor = transform(tensor)
 
         if spec.requires_grad and tensor.is_floating_point():
             tensor = tensor.requires_grad_(True)
@@ -262,7 +307,7 @@ def generate_inputs(
 
     # Custom post-init mutation
     if config.post_init is not None:
-        config.post_init(inputs, B=B, T=T, H=H, D=D, **extra_shape_kw)
+        config.post_init(inputs, B=B, T=T, H=H, D=D, input_profile=input_profile, **extra_shape_kw)
 
     return inputs
 
@@ -279,6 +324,12 @@ def generate_inputs(
     'v': TensorSpec(shape_BTHD),
 }
 
+_silu_qkv = {
+    'q': TensorSpec(shape_BTHD, realistic_transform=silu_transform),
+    'k': TensorSpec(shape_BTHD, realistic_transform=silu_transform),
+    'v': TensorSpec(shape_BTHD, realistic_transform=silu_transform),
+}
+
 register_op(OpConfig(
     name='chunk_retention',
     import_path='fla.ops.retention',
@@ -300,7 +351,11 @@ def generate_inputs(
     import_path='fla.ops.gla',
     inputs={
         **_simple_qkv,
-        'g': TensorSpec(shape_BTHD, transform=logsigmoid_clamp),
+        'g': TensorSpec(
+            shape_BTHD,
+            transform=logsigmoid_clamp,
+            realistic_transform=normalized_logsigmoid,
+        ),
     },
     category='elem_gate',
 ))
@@ -311,7 +366,7 @@ def generate_inputs(
     name='chunk_delta_rule',
     import_path='fla.ops.delta_rule',
     inputs={
-        **_simple_qkv,
+        **_silu_qkv,
         'beta': TensorSpec(shape_BTH, transform=sigmoid_transform),
     },
     category='beta',
@@ -324,8 +379,8 @@ def generate_inputs(
     name='chunk_gdn',
     import_path='fla.ops.gated_delta_rule',
     inputs={
-        **_simple_qkv,
-        'g': TensorSpec(shape_BTH, transform=logsigmoid),
+        **_silu_qkv,
+        'g': TensorSpec(shape_BTH, transform=logsigmoid, realistic_transform=learned_decay_transform),
         'beta': TensorSpec(shape_BTH, transform=sigmoid_transform),
     },
     func_name='chunk_gated_delta_rule',
@@ -334,12 +389,20 @@ def generate_inputs(
     test_file='tests/ops/test_gdn.py',
 ))
 
+
+def _comba_post_init(inputs, B, T, H, D, input_profile='synthetic', **kw):
+    """Comba derives p from k, optionally scaled by a learned per-head decay."""
+    if input_profile == 'realistic':
+        decay = torch.ones(H, dtype=inputs['k'].dtype, device=inputs['k'].device).sigmoid()
+        inputs['p'] = (inputs['k'].detach() * decay[None, None, :, None]).requires_grad_(True)
+
+
 register_op(OpConfig(
     name='chunk_kda',
     import_path='fla.ops.kda',
     inputs={
-        **_simple_qkv,
-        'g': TensorSpec(shape_BTHD, transform=logsigmoid),
+        **_silu_qkv,
+        'g': TensorSpec(shape_BTHD, transform=logsigmoid, realistic_transform=safe_decay_transform),
         'beta': TensorSpec(shape_BTH, transform=sigmoid_transform),
     },
     extra_kwargs={'use_qk_l2norm_in_kernel': True, 'safe_gate': True, 'lower_bound': -5},
@@ -386,7 +449,7 @@ def generate_inputs(
     import_path='fla.ops.simple_gla',
     inputs={
         **_simple_qkv,
-        'g': TensorSpec(shape_BTH, transform=logsigmoid),
+        'g': TensorSpec(shape_BTH, transform=logsigmoid, realistic_transform=normalized_logsigmoid),
     },
     category='head_gate',
 ))
@@ -436,12 +499,13 @@ def _rwkv7_post_init(inputs, B, T, H, D, **kw):
     name='chunk_comba',
     import_path='fla.ops.comba',
     inputs={
-        **_simple_qkv,
+        **_silu_qkv,
         'p': TensorSpec(shape_BTHD),
-        'g': TensorSpec(shape_BTH, transform=logsigmoid),
+        'g': TensorSpec(shape_BTH, transform=logsigmoid, realistic_transform=learned_decay_transform),
         'beta': TensorSpec(shape_BTH, transform=sigmoid_transform),
     },
     extra_kwargs={'use_qk_l2norm_in_kernel': True},
+    post_init=_comba_post_init,
     category='comba',
 ))
 
diff --git a/benchmarks/ops/run.py b/benchmarks/ops/run.py
index be7a93da32..a780082246 100644
--- a/benchmarks/ops/run.py
+++ b/benchmarks/ops/run.py
@@ -46,6 +46,10 @@
     python -m benchmarks.ops.run --op chunk_gla \\
         --custom-shapes '{"test": {"B": 2, "T": 4096, "H": 32, "D": 128}}'
 
+    # Layer-shaped q/k/v and gate distributions with a reproducible seed
+    python -m benchmarks.ops.run --op chunk_gla chunk_gdn chunk_kda \\
+        --input-profile realistic --seed 42
+
     # List all registered ops
     python -m benchmarks.ops.run --list
 
@@ -131,6 +135,8 @@
 4. **Isolation**: HEAD and ``--base`` each run in a subprocess. The parent
    never holds accelerator tensors, so a large HEAD shape cannot starve the
    baseline process of HBM (and vice versa).
+5. ``--input-profile realistic`` follows the transforms and parameter ranges used by FLA layers.
+   ``--seed`` is reapplied per op/shape so input samples are stable across runs and baseline refs.
 """
 
 from __future__ import annotations
@@ -279,6 +285,8 @@ def benchmark_op(
     shapes: dict[str, dict[str, int]],
     modes: list[str] | None = None,
     backend: str | None = None,
+    input_profile: str = 'synthetic',
+    seed: int = 42,
 ) -> list[dict]:
     """Benchmark a single op across all *shapes* and *modes*.
 
@@ -309,6 +317,8 @@ def benchmark_op(
     elif backend == 'triton':
         for env in backend_env.values():
             os.environ[env] = '0'
+    if input_profile != 'synthetic':
+        op_label = f"{op_label}[{input_profile}]"
 
     if config.skip_backward and 'fwdbwd' in modes:
         modes = [m for m in modes if m != 'fwdbwd']
@@ -343,11 +353,22 @@ def benchmark_op(
     # Phase 1: warmup ALL shapes before timing ANY
     print(f"\n  [{op_name}] Warming up {len(valid_shapes)} shape(s)...")
     failed_shapes = set()
-    for shape_name, shape_dict in valid_shapes.items():
+    for shape_index, (shape_name, shape_dict) in enumerate(valid_shapes.items()):
         B, T, H, D = shape_dict['B'], shape_dict['T'], shape_dict['H'], shape_dict['D']
         extra_shape_kw = {k: v for k, v in shape_dict.items() if k not in ('B', 'T', 'H', 'D')}
         try:
-            inputs = generate_inputs(config, B, T, H, D, dtype=dtype, device=device_name, **extra_shape_kw)
+            torch.manual_seed(seed + shape_index)
+            inputs = generate_inputs(
+                config,
+                B,
+                T,
+                H,
+                D,
+                dtype=dtype,
+                device=device_name,
+                input_profile=input_profile,
+                **extra_shape_kw,
+            )
             out = op_fn(**inputs, **call_kwargs)
             out_tensor = out[0] if config.output_is_tuple else out
             do = torch.randn_like(out_tensor)
@@ -368,11 +389,22 @@ def _fwdbwd_fn(inputs=inputs, do=do):
 
     # Phase 2: timing
     results = []
-    for shape_name, shape_dict in list(valid_shapes.items()):
+    for shape_index, (shape_name, shape_dict) in enumerate(valid_shapes.items()):
         B, T, H, D = shape_dict['B'], shape_dict['T'], shape_dict['H'], shape_dict['D']
         extra_shape_kw = {k: v for k, v in shape_dict.items() if k not in ('B', 'T', 'H', 'D')}
         try:
-            inputs = generate_inputs(config, B, T, H, D, dtype=dtype, device=device_name, **extra_shape_kw)
+            torch.manual_seed(seed + shape_index)
+            inputs = generate_inputs(
+                config,
+                B,
+                T,
+                H,
+                D,
+                dtype=dtype,
+                device=device_name,
+                input_profile=input_profile,
+                **extra_shape_kw,
+            )
         except Exception as e:
             logger.warning(f"Input generation failed for {op_name} @ {shape_name}: {e}")
             continue
@@ -563,11 +595,21 @@ def _find_project_root() -> str:
 _WORKER_ENV = 'FLA_BENCH_WORKER'
 
 
-def _isolated_bench_cmd(runner, op_names, shape_configs, modes, backend, out_json):
+def _isolated_bench_cmd(
+    runner,
+    op_names,
+    shape_configs,
+    modes,
+    backend,
+    out_json,
+    input_profile='synthetic',
+    seed=42,
+):
     """Worker argv. ``--no-base`` plus ``FLA_BENCH_WORKER`` block nested compares."""
     cmd = [sys.executable, runner, '--op', *op_names,
            '--custom-shapes', json.dumps(shape_configs),
-           '--modes', *modes, '--json', out_json, '--no-base']
+           '--modes', *modes, '--input-profile', input_profile,
+           '--seed', str(seed), '--json', out_json, '--no-base']
     if backend is not None:
         cmd += ['--backend', backend]
     return cmd
@@ -587,7 +629,14 @@ def _read_bench_json(out_json):
     return None, None
 
 
-def _bench_current(op_names, shape_configs, modes, backend=None):
+def _bench_current(
+    op_names,
+    shape_configs,
+    modes,
+    backend=None,
+    input_profile='synthetic',
+    seed=42,
+):
     """Run the current working tree in a subprocess, then exit to free HBM.
 
     Returns (results_list, machine_info_dict) or (None, None) on failure.
@@ -600,7 +649,16 @@ def _bench_current(op_names, shape_configs, modes, backend=None):
         out_json = os.path.join(tmpdir, 'results.json')
         _run_isolated_bench(
             project_root,
-            _isolated_bench_cmd(runner, op_names, shape_configs, modes, backend, out_json),
+            _isolated_bench_cmd(
+                runner,
+                op_names,
+                shape_configs,
+                modes,
+                backend,
+                out_json,
+                input_profile,
+                seed,
+            ),
         )
         return _read_bench_json(out_json)
     except Exception as e:
@@ -610,7 +668,7 @@ def _bench_current(op_names, shape_configs, modes, backend=None):
         shutil.rmtree(tmpdir, ignore_errors=True)
 
 
-def _bench_at_ref(ref, op_names, shape_configs, modes, backend=None):
+def _bench_at_ref(ref, op_names, shape_configs, modes, backend=None, input_profile='synthetic', seed=42):
     """Run benchmarks at a git ref using a temporary worktree.
 
     Returns (results_list, machine_info_dict) or (None, None) on failure.
@@ -643,7 +701,16 @@ def _bench_at_ref(ref, op_names, shape_configs, modes, backend=None):
         out_json = os.path.join(tmpdir, 'results.json')
         _run_isolated_bench(
             worktree_dir,
-            _isolated_bench_cmd(runner, op_names, shape_configs, modes, backend, out_json),
+            _isolated_bench_cmd(
+                runner,
+                op_names,
+                shape_configs,
+                modes,
+                backend,
+                out_json,
+                input_profile,
+                seed,
+            ),
         )
         return _read_bench_json(out_json)
     except Exception as e:
@@ -685,6 +752,14 @@ def main():
         choices=['fwd', 'fwdbwd'],
         help='Benchmark modes (default: fwd fwdbwd)',
     )
+    parser.add_argument(
+        '--input-profile', choices=['synthetic', 'realistic'], default='synthetic',
+        help='Input distribution profile (default: synthetic)',
+    )
+    parser.add_argument(
+        '--seed', type=int, default=42,
+        help='Input seed reapplied for each op/shape (default: 42)',
+    )
     parser.add_argument(
         '--json', dest='json_file', default=None,
         help='Output file path for JSON results',
@@ -750,6 +825,7 @@ def main():
     print(f"Git: {_get_git_label()}")
     print(f"Shapes: {len(shape_configs)} configs")
     print(f"Ops: {op_names}")
+    print(f"Input profile: {args.input_profile} (seed={args.seed})")
 
     if is_worker:
         machine_info = _get_machine_info()
@@ -758,14 +834,26 @@ def main():
         for op_name in op_names:
             try:
                 all_results.extend(
-                    benchmark_op(op_name, shape_configs, modes=args.modes, backend=args.backend),
+                    benchmark_op(
+                        op_name,
+                        shape_configs,
+                        modes=args.modes,
+                        backend=args.backend,
+                        input_profile=args.input_profile,
+                        seed=args.seed,
+                    ),
                 )
             except Exception as e:
                 logger.error(f"Failed to benchmark {op_name}: {e}")
         baseline, baseline_info = None, None
     else:
         all_results, machine_info = _bench_current(
-            op_names, shape_configs, args.modes, backend=args.backend,
+            op_names,
+            shape_configs,
+            args.modes,
+            backend=args.backend,
+            input_profile=args.input_profile,
+            seed=args.seed,
         )
         all_results = all_results or []
         machine_info = machine_info or {'git_label': _get_git_label()}
@@ -781,7 +869,14 @@ def main():
         baseline, baseline_info = None, None
         if base_ref:
             baseline, baseline_info = _bench_at_ref(
-                base_ref, op_names, shape_configs, args.modes, backend=args.backend)
+                base_ref,
+                op_names,
+                shape_configs,
+                args.modes,
+                backend=args.backend,
+                input_profile=args.input_profile,
+                seed=args.seed,
+            )
 
     # Sort by (mode, L, B, T, H, D, op) so the table groups by mode first
     # and (when present) by L so different residual-source counts cluster.
diff --git a/tests/test_benchmark_distributions.py b/tests/test_benchmark_distributions.py
new file mode 100644
index 0000000000..0c33bd8443
--- /dev/null
+++ b/tests/test_benchmark_distributions.py
@@ -0,0 +1,113 @@
+# Copyright (c) 2023-2026, Songlin Yang, Yu Zhang, Zhiyuan Li
+#
+# This source code is licensed under the MIT license found in the
+# LICENSE file in the root directory of this source tree.
+# For a list of all contributors, visit:
+#   https://github.com/fla-org/flash-linear-attention/graphs/contributors
+
+import pytest
+import torch
+
+from benchmarks.distributions import sample_lognormal_packed_lengths
+from benchmarks.ops.registry import generate_inputs, get_op
+
+
+def test_lognormal_lengths_are_seeded_bounded_and_exact():
+    first = sample_lognormal_packed_lengths(
+        total_tokens=8192,
+        num_sequences=32,
+        min_length=16,
+        max_length=1024,
+        sigma=1.0,
+        generator=torch.Generator().manual_seed(42),
+    )
+    second = sample_lognormal_packed_lengths(
+        total_tokens=8192,
+        num_sequences=32,
+        min_length=16,
+        max_length=1024,
+        sigma=1.0,
+        generator=torch.Generator().manual_seed(42),
+    )
+
+    assert torch.equal(first, second)
+    assert first.sum() == 8192
+    assert first.min() >= 16
+    assert first.max() <= 1024
+    assert first.float().quantile(0.9) > 2 * first.float().median()
+
+
+@pytest.mark.parametrize(
+    'kwargs',
+    [
+        {'total_tokens': 0, 'num_sequences': 1, 'max_length': 16},
+        {'total_tokens': 32, 'num_sequences': 0, 'max_length': 16},
+        {'total_tokens': 32, 'num_sequences': 2, 'min_length': 17, 'max_length': 16},
+        {'total_tokens': 33, 'num_sequences': 2, 'min_length': 16, 'max_length': 16},
+    ],
+)
+def test_lognormal_lengths_reject_invalid_packing(kwargs):
+    with pytest.raises(ValueError):
+        sample_lognormal_packed_lengths(**kwargs)
+
+
+def test_realistic_gla_gate_matches_layer_normalizer():
+    config = get_op('chunk_gla')
+    torch.manual_seed(42)
+    synthetic = generate_inputs(config, 1, 64, 2, 16, dtype=torch.float32, device='cpu')
+    torch.manual_seed(42)
+    realistic = generate_inputs(
+        config,
+        1,
+        64,
+        2,
+        16,
+        dtype=torch.float32,
+        device='cpu',
+        input_profile='realistic',
+    )
+
+    assert realistic['g'].max() <= 0
+    assert torch.allclose(realistic['g'], synthetic['g'] / 16)
+
+
+@pytest.mark.parametrize('op_name', ['chunk_gdn', 'chunk_comba'])
+def test_realistic_learned_decay_is_nonsymmetric(op_name):
+    torch.manual_seed(42)
+    inputs = generate_inputs(
+        get_op(op_name),
+        1,
+        256,
+        4,
+        16,
+        dtype=torch.float32,
+        device='cpu',
+        input_profile='realistic',
+    )
+
+    assert inputs['q'].mean() > 0
+    assert inputs['g'].max() < 0
+    assert inputs['g'].std() > 0
+
+    if op_name == 'chunk_comba':
+        expected_p = inputs['k'].detach() * torch.ones(4).sigmoid()[None, None, :, None]
+        assert torch.allclose(inputs['p'], expected_p)
+
+
+def test_realistic_kda_gate_respects_safe_bound():
+    torch.manual_seed(42)
+    inputs = generate_inputs(
+        get_op('chunk_kda'),
+        1,
+        256,
+        4,
+        16,
+        dtype=torch.float32,
+        device='cpu',
+        input_profile='realistic',
+    )
+
+    assert inputs['g'].min() >= -5
+    assert inputs['g'].max() <= 0
+    assert inputs['g'].min() < 0
+    assert inputs['g'].std() > 0