diff --git a/docs/mechanism_api.md b/docs/mechanism_api.md
index d7fefdc..38dfe46 100644
--- a/docs/mechanism_api.md
+++ b/docs/mechanism_api.md
@@ -280,7 +280,7 @@ composed or split as cleanly:
     bisections, whose composite `dp_event` is a composition of exponential
     mechanisms.
 
-[dp-quantiles-src]: https://github.com/google/dpsynth/blob/main/dpsynth/local_mode/_quantiles.py
+[dp-quantiles-src]: https://github.com/google/dpsynth/blob/main/dpsynth/local_mode/primitives.py
 
 Because these heterogeneous parameters cannot be directly combined,
 `configure()` accepts a single scalar `zcdp_rho` and translates it into each
diff --git a/dpsynth/local_mode/_quantiles.py b/dpsynth/local_mode/_quantiles.py
deleted file mode 100644
index 0fc6b0b..0000000
--- a/dpsynth/local_mode/_quantiles.py
+++ /dev/null
@@ -1,153 +0,0 @@
-# Copyright 2026 Google LLC
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#     http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-"""DP quantiles from dense histograms via recursive median bisection.
-
-This module computes differentially private quantile edges from a dense
-histogram of counts, using the discrete exponential mechanism.  It works purely
-in index space -- ``quantiles_from_histogram`` returns cell indices into the
-histogram, and the caller maps those indices to domain values.  The primary use
-case is a two-pass pipeline: a first pass computes a dense histogram over a
-fine-grained grid, then ``quantiles_from_histogram`` finds DP quantile indices
-from that histogram without touching individual records.
-
-Tie handling via jitter
------------------------
-Recursive median bisection needs each record assigned to one side of every
-split independently.  A "spike" of records tied on one grid cell breaks this: a
-whole-cell split sends all that mass to one side, biasing the quantiles and
-collapsing sub-ranges (dropping edges).  We fix this by breaking ties directly
-in the histogram domain rather than over the raw data values -- each cell's
-count is redistributed to nearby cells as ``Multinomial(count, kernel)`` (one
-draw per non-empty cell), which is distributionally identical to independently
-perturbing each record and so needs no extra privacy budget.  The ``refine``
-strategy uses a strictly-positive kernel over refined sub-cells (value-
-preserving); the ``symmetric`` strategy uses a symmetric kernel over neighboring
-grid cells.
-"""
-
-from __future__ import annotations
-
-from typing import Literal
-
-import numpy as np
-import scipy.special
-
-
-def _median_from_histogram(
-    rng: np.random.Generator,
-    counts: np.ndarray,
-    epsilon: float,
-) -> int:
-  """Returns the index of a DP median within a dense histogram.
-
-  Args:
-    rng: A numpy random number generator.
-    counts: Dense 1D histogram counts.
-    epsilon: Exponential mechanism privacy parameter for this level.
-
-  Returns:
-    The index of the selected median grid point within ``counts``.
-  """
-  total_points = len(counts)
-  if total_points == 0:
-    return 0
-  n = counts.sum()
-  target = n / 2.0
-  cumsum = np.cumsum(counts)
-
-  # Infinite budget = exact median, useful for testing.
-  if epsilon == np.inf:
-    return int(np.searchsorted(cumsum, target))
-
-  # Score u(v) = -dist(target, [L_v, R_v]), sensitivity 1/2.
-  left_ranks = np.r_[0, cumsum[:-1]]
-  scores = -np.maximum(0, np.maximum(left_ranks - target, target - cumsum))
-
-  probs = scipy.special.softmax(epsilon * scores)
-  return int(rng.choice(total_points, p=probs))
-
-
-def jitter_factor(num_partitions):
-  """Returns a data-independent jitter resolution m from num_partitions."""
-  # m >= num_partitions keeps each jittered cell below one partition's mass;
-  # the 4x absorbs multinomial fluctuation.
-  return max(1, 4 * num_partitions)
-
-
-def quantiles_from_histogram(
-    rng: np.random.Generator,
-    counts: np.ndarray,
-    epsilon_levels: np.ndarray,
-    jitter_strategy: Literal['symmetric', 'refine'],
-    max_records_per_user: int = 1,
-) -> list[int]:
-  """DP quantile edge indices into ``counts`` via jittered median bisection.
-
-  Operates purely in index space: it returns cell indices into ``counts`` and
-  leaves the mapping from index to domain value to the caller.
-
-  Args:
-    rng: A numpy random number generator.
-    counts: Dense 1D histogram counts.
-    epsilon_levels: Per-level exponential mechanism epsilons, ordered from the
-      deepest (finest) level to the shallowest (coarsest).
-    jitter_strategy: Specifies the pre-processing jitter strategy, -
-      'symmetric': jitter mass to +/- m//2 neighbors on the same grid. -
-      'refine': jitter mass to m equivalent sub-cells.
-    max_records_per_user: Assumed upper bound on the number of records per user.
-
-  Returns:
-    A sorted list of ``2 ** len(epsilon_levels) - 1`` cell indices.
-  """
-  if max_records_per_user != 1:
-    # The privacy analysis of this mechanism relies on parallel composition
-    # across the nodes of each level of the hierarchy. When users have
-    # multiple records, they may contirbute to multiple nodes, which would
-    # require a different privacy analysis (TBD).
-    raise NotImplementedError('max_records_per_user != 1 not yet supported.')
-  counts = np.asarray(counts)
-  m = jitter_factor(2 ** len(epsilon_levels))
-
-  if jitter_strategy == 'refine':
-    stride, offsets = m, np.arange(m)
-  else:
-    half = m // 2
-    stride, offsets = 1, np.arange(-half, half + 1)
-
-  # Scatter each cell's mass over its jittered targets: same law as perturbing
-  # each record, so it breaks ties without spending extra privacy budget.
-  num_cells = counts.size * stride
-  nz = np.flatnonzero(counts)
-  probas = np.full(offsets.size, 1.0 / offsets.size)
-  split = rng.multinomial(counts[nz].astype(np.int64), probas)
-  targets = np.clip(nz[:, None] * stride + offsets, 0, num_cells - 1)
-  jittered = np.bincount(  # pyrefly: ignore[no-matching-overload]
-      targets.flatten(), weights=split.flatten(), minlength=num_cells
-  )
-
-  def _rec(lo_idx, hi_idx, depth):
-    if depth == 0:
-      return []
-    median_idx = lo_idx + _median_from_histogram(
-        rng, jittered[lo_idx:hi_idx], epsilon_levels[depth - 1]
-    )
-    left = _rec(lo_idx, median_idx, depth - 1)
-    right = _rec(median_idx, hi_idx, depth - 1)
-    return left + [median_idx] + right
-
-  result = _rec(0, jittered.size, len(epsilon_levels))
-  if jitter_strategy == 'refine':
-    result = [idx // m for idx in result]
-  return result
diff --git a/dpsynth/local_mode/initialization.py b/dpsynth/local_mode/initialization.py
index f9443ee..9af3c8b 100644
--- a/dpsynth/local_mode/initialization.py
+++ b/dpsynth/local_mode/initialization.py
@@ -22,7 +22,6 @@
 import dp_accounting
 from dpsynth import api
 from dpsynth import domain
-from dpsynth.local_mode import _quantiles
 from dpsynth.local_mode import primitives
 from dpsynth.local_mode import vectorized_transformations as vtx
 import numpy as np
@@ -90,7 +89,7 @@ def compute_grid_spec(
   """Returns (lower, upper, grid_size) for the quantile candidate grid."""
   min_value = float(attribute.min_value)
   if attribute.dtype == 'int':
-    m = _quantiles.jitter_factor(num_partitions)
+    m = primitives.jitter_factor(num_partitions)
     budget = max(2, max_grid_size // m)
     int_range = int(attribute.max_value - attribute.min_value + 1)
     step = max(1, math.ceil(int_range / budget))
@@ -212,7 +211,7 @@ def from_summary(
   ) -> NumericalMeasurement:
     """Returns a NumericalMeasurement from pre-aggregated histogram counts."""
     jitter_strategy = 'refine' if self.attribute.dtype == 'int' else 'symmetric'
-    indices = _quantiles.quantiles_from_histogram(
+    indices = primitives.quantiles_from_histogram(
         rng,
         counts,
         epsilon_levels=np.asarray(self.epsilon_levels),
diff --git a/dpsynth/local_mode/primitives.py b/dpsynth/local_mode/primitives.py
index 5d0e4c3..b57e313 100644
--- a/dpsynth/local_mode/primitives.py
+++ b/dpsynth/local_mode/primitives.py
@@ -14,67 +14,220 @@
 
 """Differentially private primitives for local mode synthetic data generation.
 
+This module is intended to be the singular home for low-level DP building blocks
+in dpsynth, making a potential future switch to a PyDP or OpenDP backend
+simpler.
+While not all DP code in dpsynth currently goes through this module, that is the
+long-term goal (at least for tabular data).
+
 Design Decisions:
 
-1) Summary Statistics: All primitives are defined strictly in terms of summary
-   statistics (counts, histograms, etc.) rather than raw datasets. This ensures
-   computational efficiency and separation of data aggregation from DP
-   calibration.
-2) Pure Functions: Primitives are pure functions that take numpy inputs and
-return
-   numpy outputs. State and DP accounting, when needed, are managed externally.
+1) Pre-Computed Aggregates: All functions expect pre-computed aggregates
+   (counts, histograms, quality scores) rather than performing data aggregation
+   themselves. Therefore, it is the responsibility of the caller to ensure the
+   appropriate privacy assumptions are satisfied (primarily that each user
+   contributes at most one record to one bucket, or that user contributions are
+   properly bounded via max_records_per_user).
+2) Pure Functions: Primitives are pure functions that take NumPy inputs and
+   return NumPy or Python outputs. State and DP accounting, when needed, are
+   managed externally.
 
 Privacy Characterizations:
 
-- Quantiles (via `_quantiles.quantiles_from_histogram`): This mechanism is a
-composition
-  of exponential mechanisms with parameters taken from `epsilon_levels`.
-- Partition Selection (`_select_partitions_sips`): DP-SIPS mechanism for
-discovering
-  open-set vocabulary, utilizing Gaussian Thresholding combined with privacy
-  filtering.
-- Gaussian Thresholding (`select_partitions_gaussian_thresholding`): This
-mechanism adds
-  Gaussian noise to counts and tests against a threshold that provides a (0,
-  delta)
-  bound on false positives for unpopulated partitions.
-- Gaussian Noise (`add_gaussian_noise`): This is a standard Gaussian mechanism
-applied
-  to the input summary statistics (e.g. counts).
+- Exponential Mechanism (`exponential_mechanism`): Standard exponential
+  mechanism for discrete selection given candidate quality scores.
+- Quantiles (`quantiles_from_histogram`): Composition of exponential mechanisms
+  via jittered recursive median bisection over a dense histogram.
+- Gaussian Thresholding (`select_partitions_gaussian_thresholding`): Partition
+  selection mechanism that adds Gaussian noise to counts and tests against a
+  threshold bounding false positives for empty partitions at delta.
+- Gaussian Noise (`add_gaussian_noise`): Standard Gaussian mechanism applied to
+  input summary statistics (e.g., counts).
 """
 
 from __future__ import annotations
 
+from typing import Literal
 
 import numpy as np
 import scipy.stats
 
+# ---------------------------------------------------------------------------
+# Exponential Mechanism
+# ---------------------------------------------------------------------------
+
+
+def exponential_mechanism(
+    rng: np.random.Generator,
+    quality_scores: np.ndarray,
+    epsilon: float,
+    sensitivity: float = 1.0,
+    monotonic: bool = False,
+) -> int:
+  """Selects an index using the discrete exponential mechanism.
+
+  Samples a candidate index with probability proportional to
+  exp(coef * epsilon * quality_scores / sensitivity), where coef is 1.0
+  if monotonic is True and 0.5 otherwise. This is implemented via the
+  Gumbel-max trick (adding Gumbel noise with scale = sensitivity / (coef *
+  epsilon)
+  and returning argmax), which avoids computing normalizing constants and works
+  gracefully with infinite epsilon (scale = 0).
+
+  Args:
+    rng: A numpy random number generator.
+    quality_scores: 1D array of utility/quality scores for each candidate.
+    epsilon: Privacy parameter epsilon. Must be non-negative.
+    sensitivity: Upper bound on the quality score sensitivity. Must be positive.
+    monotonic: Whether the score function is monotonic with respect to dataset
+      modifications (sensitivity Delta u instead of 2 * Delta u). Defaults to
+      False.
+
+  Returns:
+    The index of the selected candidate.
+
+  Raises:
+    ValueError: If epsilon < 0, sensitivity <= 0, or quality_scores is empty.
+  """
+  if epsilon < 0:
+    raise ValueError(f'epsilon must be non-negative, got {epsilon}')
+  if sensitivity <= 0:
+    raise ValueError(f'sensitivity must be positive, got {sensitivity}')
+
+  scores = np.asarray(quality_scores, dtype=float)
+  if scores.size == 0:
+    raise ValueError('quality_scores must not be empty.')
+
+  if epsilon == 0:
+    return int(rng.choice(scores.size))
+
+  coef = 1.0 if monotonic else 0.5
+  scale = sensitivity / (coef * epsilon)
+  noise = rng.gumbel(scale=scale, size=scores.size)
+  return int(np.argmax(scores + noise))
 
-_UNCALIBRATED_MSG = (
-    '{param} has not been set. Set it directly or call calibrate().'
-)
 
+# ---------------------------------------------------------------------------
+# DP Quantiles via Recursive Median Bisection
+# ---------------------------------------------------------------------------
+
+
+def _median_from_histogram(
+    rng: np.random.Generator,
+    counts: np.ndarray,
+    epsilon: float,
+) -> int:
+  """Returns the index of a DP median within a dense histogram."""
+  total_points = len(counts)
+  if total_points == 0:
+    return 0
+  n = counts.sum()
+  target = n / 2.0
+  cumsum = np.cumsum(counts)
+
+  # Score u(v) = -dist(target, [L_v, R_v]), sensitivity 1/2.
+  left_ranks = np.r_[0, cumsum[:-1]]
+  scores = -np.maximum(0, np.maximum(left_ranks - target, target - cumsum))
+
+  return exponential_mechanism(
+      rng=rng,
+      quality_scores=scores,
+      epsilon=epsilon,
+      sensitivity=0.5,
+      monotonic=False,
+  )
+
+
+def jitter_factor(num_partitions: int) -> int:
+  """Returns a data-independent jitter resolution m from num_partitions."""
+  # m >= num_partitions keeps each jittered cell below one partition's mass;
+  # the 4x absorbs multinomial fluctuation.
+  return max(1, 4 * num_partitions)
+
+
+def quantiles_from_histogram(
+    rng: np.random.Generator,
+    counts: np.ndarray,
+    epsilon_levels: np.ndarray,
+    jitter_strategy: Literal['symmetric', 'refine'],
+    max_records_per_user: int = 1,
+) -> list[int]:
+  """DP quantile edge indices into ``counts`` via jittered median bisection.
+
+  Operates purely in index space: it returns cell indices into ``counts`` and
+  leaves the mapping from index to domain value to the caller.
+
+  Tie handling via jitter:
+  Recursive median bisection needs each record assigned to one side of every
+  split independently. A "spike" of records tied on one grid cell breaks this: a
+  whole-cell split sends all that mass to one side, biasing the quantiles and
+  collapsing sub-ranges (dropping edges). We fix this by breaking ties directly
+  in the histogram domain rather than over the raw data values -- each cell's
+  count is redistributed to nearby cells as Multinomial(count, kernel) (one
+  draw per non-empty cell), which is distributionally identical to independently
+  perturbing each record and so needs no extra privacy budget. The ``refine``
+  strategy uses a strictly-positive kernel over refined sub-cells (value-
+  preserving); the ``symmetric`` strategy uses a symmetric kernel over
+  neighboring grid cells.
+
+  Args:
+    rng: A numpy random number generator.
+    counts: Dense 1D histogram counts.
+    epsilon_levels: Per-level exponential mechanism epsilons, ordered from the
+      deepest (finest) level to the shallowest (coarsest).
+    jitter_strategy: Specifies the pre-processing jitter strategy, -
+      'symmetric': jitter mass to +/- m//2 neighbors on the same grid. -
+      'refine': jitter mass to m equivalent sub-cells.
+    max_records_per_user: Assumed upper bound on the number of records per user.
+
+  Returns:
+    A sorted list of ``2 ** len(epsilon_levels) - 1`` cell indices.
+  """
+  if max_records_per_user != 1:
+    # The privacy analysis of this mechanism relies on parallel composition
+    # across the nodes of each level of the hierarchy. When users have
+    # multiple records, they may contribute to multiple nodes, which would
+    # require a different privacy analysis (TBD).
+    raise NotImplementedError('max_records_per_user != 1 not yet supported.')
+  counts = np.asarray(counts)
+  m = jitter_factor(2 ** len(epsilon_levels))
+
+  if jitter_strategy == 'refine':
+    stride, offsets = m, np.arange(m)
+  else:
+    half = m // 2
+    stride, offsets = 1, np.arange(-half, half + 1)
+
+  # Scatter each cell's mass over its jittered targets: same law as perturbing
+  # each record, so it breaks ties without spending extra privacy budget.
+  num_cells = counts.size * stride
+  nz = np.flatnonzero(counts)
+  probas = np.full(offsets.size, 1.0 / offsets.size)
+  split = rng.multinomial(counts[nz].astype(np.int64), probas)
+  targets = np.clip(nz[:, None] * stride + offsets, 0, num_cells - 1)
+  jittered = np.bincount(  # pyrefly: ignore[no-matching-overload]
+      targets.flatten(), weights=split.flatten(), minlength=num_cells
+  )
+
+  def _rec(lo_idx, hi_idx, depth):
+    if depth == 0:
+      return []
+    median_idx = lo_idx + _median_from_histogram(
+        rng, jittered[lo_idx:hi_idx], epsilon_levels[depth - 1]
+    )
+    left = _rec(lo_idx, median_idx, depth - 1)
+    right = _rec(median_idx, hi_idx, depth - 1)
+    return left + [median_idx] + right
 
-def _contribution_bound(prng, user_ids, max_part):
-  """Return array idx where all ids appear <=max_part times in user_ids[idx]."""
-  # Sort by ID + noise to shuffle within groups. Then find where
-  # groups start/end, and select the first max_part elements of each group.
-  # Use lexsort with random keys to shuffle string/object IDs safely.
-  random_keys = prng.uniform(size=user_ids.size)
-  idx = np.lexsort((random_keys, user_ids))
-  sorted_ids = user_ids[idx]
-  diff = np.r_[True, sorted_ids[1:] != sorted_ids[:-1]]
-  kernel = np.ones(max_part, dtype=bool)
-  # This convolution determines if any of previous max_part elements are True.
-  mask = np.convolve(diff, kernel, mode='full')[: user_ids.size]
-  return idx[mask]
+  result = _rec(0, jittered.size, len(epsilon_levels))
+  if jitter_strategy == 'refine':
+    result = [idx // m for idx in result]
+  return result
 
 
-def _get_threshold(delta, sigma, max_part):
-  ks = np.arange(1, max_part + 1)
-  failure_prob = (1 - delta) ** (1 / ks)
-  thresholds = 1 / np.sqrt(ks) + sigma * scipy.stats.norm.ppf(failure_prob)
-  return thresholds.max()
+# ---------------------------------------------------------------------------
+# Partition Selection
+# ---------------------------------------------------------------------------
 
 
 def select_partitions_gaussian_thresholding(
@@ -171,7 +324,7 @@ def select_partitions_gaussian_thresholding(
   base = float(max_records_per_user + min_count - 1)
   threshold = base + stddev * scipy.stats.norm.ppf(1.0 - delta)
   passed = noisy_counts >= threshold
-  # unique_parts is sorted (see np.unique), so the output order is determinstic.
+  # unique_parts is sorted (np.unique), so the output order is deterministic.
   return unique_parts[passed], noisy_counts[passed], stddev
 
 
@@ -211,113 +364,8 @@ def ensure_public_partitions(
   return all_selected[order], all_counts[order]
 
 
-def _select_partitions_sips(
-    rng: np.random.Generator,
-    data: np.ndarray,
-    gdp_budget: float,
-    delta: float,
-    num_rounds: int | None = None,
-    user_ids: np.ndarray | None = None,
-    max_part: int = 1,
-    allocation_factor: float = 0.3,
-) -> tuple[np.ndarray, np.ndarray, float]:
-  """Implements the DP-SIPS mechanism for partition selection.
-
-  Args:
-    rng: A numpy random number generator.
-    data: 1D array of integers, where each element is a partition ID.
-    gdp_budget: Total privacy budget in terms of squared Gaussian DP mu
-      parameter (gdp_budget = mu^2 = 1 / sigma^2).
-    delta: Failure probability (false positive bound per empty partition).
-    num_rounds: Number of rounds to run the mechanism. Defaults to 1 if user_ids
-      is None, and 3 otherwise.
-    user_ids: Optional 1D array of user IDs corresponding to data. If provided,
-      user-level DP is guaranteed. If None, item-level DP is guaranteed
-      (assuming each record is a unique user).
-    max_part: Maximum number of partitions any single user can contribute to in
-      a single round.
-    allocation_factor: Factor by which to increase the budget each round.
-
-  Returns:
-    A tuple containing:
-      - selected_partitions: 1D array of unique partition IDs that passed the
-        threshold.
-      - estimated_counts: 1D array of noisy (or weighted noisy) counts for each
-        selected partition in the round it was discovered.
-      - standard_deviation: A single float representing the uniform standard
-        deviation of the noise added to the estimated counts.
-  """
-  if num_rounds is None:
-    num_rounds = 1 if user_ids is None else 3
-  if num_rounds <= 0:
-    raise ValueError(f'num_rounds ({num_rounds}) must be greater than 0.')
-  if gdp_budget <= 0 or delta <= 0 or delta > 1:
-    raise ValueError(f'{gdp_budget=} and {delta=} must be positive.')
-
-  fractions = allocation_factor ** np.arange(num_rounds)[::-1]
-  fractions /= fractions.sum()
-  gdp_rounds, delta_rounds = gdp_budget * fractions, delta * fractions
-  sigma_rounds = 1.0 / np.sqrt(gdp_rounds)
-  max_sigma = float(np.max(sigma_rounds))
-
-  if data.size == 0:
-    return np.empty(0, dtype=data.dtype), np.empty(0, dtype=float), max_sigma
-
-  if user_ids is None:
-    user_ids = np.arange(data.size)
-  if user_ids.size != data.size:
-    raise ValueError('user_ids must have the same size as data.')
-
-  combined = np.stack((user_ids, data), axis=1)
-  unique_combined = np.unique(combined, axis=0)
-  rem_user_ids = unique_combined[:, 0]
-  rem_partitions = unique_combined[:, 1]
-
-  selected_partitions = []
-  selected_counts = []
-  for i in range(num_rounds):
-    if rem_partitions.size == 0:
-      break
-
-    threshold = _get_threshold(delta_rounds[i], sigma_rounds[i], max_part)
-
-    mask = _contribution_bound(rng, rem_user_ids, max_part)
-    curr_user_ids = rem_user_ids[mask]
-    curr_partitions = rem_partitions[mask]
-
-    unique_users, user_counts = np.unique(curr_user_ids, return_counts=True)
-    user_to_count = dict(zip(unique_users, user_counts))
-    weights = np.array([1.0 / user_to_count[u] ** 0.5 for u in curr_user_ids])
-
-    unique_parts, inverse_indices = np.unique(
-        curr_partitions, return_inverse=True
-    )
-    weighted_counts = np.bincount(inverse_indices, weights=weights)
-    noised_counts = rng.normal(weighted_counts, scale=sigma_rounds[i])
-
-    passed_mask = noised_counts >= threshold
-    round_selections = unique_parts[passed_mask]
-    round_counts = noised_counts[passed_mask]
-    if round_selections.size > 0:
-      selected_partitions.append(round_selections)
-      selected_counts.append(round_counts)
-
-      mask = ~np.isin(rem_partitions, round_selections)
-      rem_user_ids = rem_user_ids[mask]
-      rem_partitions = rem_partitions[mask]
-
-  if not selected_partitions:
-    return (
-        np.empty(0, dtype=data.dtype),
-        np.empty(0, dtype=float),
-        max_sigma,
-    )
-  selected_partitions = np.concatenate(selected_partitions)
-  selected_counts = np.concatenate(selected_counts)
-  return selected_partitions, selected_counts, max_sigma
-
 # ---------------------------------------------------------------------------
-# Simple noisy counting functions
+# Gaussian Noise
 # ---------------------------------------------------------------------------
 
 
diff --git a/tests/local_mode/_quantiles_test.py b/tests/local_mode/_quantiles_test.py
deleted file mode 100644
index e5dd7ad..0000000
--- a/tests/local_mode/_quantiles_test.py
+++ /dev/null
@@ -1,136 +0,0 @@
-# Copyright 2026 Google LLC
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#     http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-
-from absl.testing import absltest
-from absl.testing import parameterized
-from dpsynth.local_mode import _quantiles
-import numpy as np
-
-
-class QuantilesFromHistogramTest(parameterized.TestCase):
-
-  def test_no_levels_returns_empty(self):
-    rng = np.random.default_rng(0)
-    counts = np.array([10])
-    for jitter_strategy in ('symmetric', 'refine'):
-      edges = _quantiles.quantiles_from_histogram(
-          rng, counts, np.array([]), jitter_strategy
-      )
-      self.assertEmpty(edges)
-
-  @parameterized.product(
-      levels=(1, 2, 3, 4),
-      jitter_strategy=('symmetric', 'refine'),
-  )
-  def test_edge_count_matches_levels(self, levels, jitter_strategy):
-    rng = np.random.default_rng(0)
-    grid_size = 10001
-    counts = rng.integers(0, 20, size=grid_size)
-    edges = _quantiles.quantiles_from_histogram(
-        rng,
-        counts,
-        epsilon_levels=np.ones(levels),
-        jitter_strategy=jitter_strategy,
-    )
-    self.assertLen(edges, 2**levels - 1)
-
-  @parameterized.parameters(1, 2, 3, 4)
-  def test_edge_count_matches_levels_with_spike(self, levels):
-    # A large tied spike must not collapse split ranges and drop edges. With
-    # whole-cell splits this dropped edges; jitter breaks up the spike so the
-    # recursion always emits the full 2**levels - 1 edges.
-    counts = np.zeros(101, dtype=np.int64)
-    counts[:40] = 1
-    counts[40] = 1000
-    counts[41:80] = 1
-    edges = _quantiles.quantiles_from_histogram(
-        np.random.default_rng(0),
-        counts,
-        epsilon_levels=np.array([np.inf] * levels),
-        jitter_strategy='refine',
-    )
-    self.assertLen(edges, 2**levels - 1)
-
-  def test_integer_edges_are_integer_indices(self):
-    # Integer attributes must return integer cell indices into ``counts``.
-    counts = np.zeros(101, dtype=np.int64)
-    counts[40] = 5000
-    counts[:40] = 50
-    counts[41:] = 50
-    edges = _quantiles.quantiles_from_histogram(
-        np.random.default_rng(0),
-        counts,
-        epsilon_levels=np.array([np.inf] * 3),
-        jitter_strategy='refine',
-    )
-    for edge in edges:
-      self.assertEqual(edge, int(edge))
-      self.assertBetween(edge, 0, counts.size - 1)
-
-  def test_exact_budget_matches_numpy_smooth(self):
-    rng = np.random.default_rng(0)
-    data = rng.integers(0, 100, size=50000)
-    counts = np.bincount(data, minlength=101)
-    edges = _quantiles.quantiles_from_histogram(
-        rng,
-        counts,
-        epsilon_levels=np.array([np.inf] * 3),
-        jitter_strategy='refine',
-    )
-    # Cell indices map 1:1 to values here (delta == 1), so compare directly.
-    expected = np.quantile(data, [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875])
-    np.testing.assert_allclose(edges, expected, atol=1.0)
-
-  def test_exact_budget_matches_numpy_with_spike(self):
-    # Reproduces the failure mode from the 'hours-per-week' column of the adult
-    # census dataset: a dominant spike at 40 carrying ~45% of the mass, with
-    # lighter integer support on either side. The pre-jitter whole-cell split
-    # lumped all tied mass into one subtree, biasing the low quantiles (e.g. the
-    # 0.25 edge came out as 18 instead of 40) and dropping upper edges. Jitter
-    # breaks up the spike so the recursive medians match numpy.
-    below = np.arange(1, 40).repeat(230)
-    spike = np.full(13500, 40)
-    above = np.arange(41, 80).repeat(190)
-    data = np.concatenate([below, spike, above])
-    counts = np.bincount(data, minlength=101)
-    edges = _quantiles.quantiles_from_histogram(
-        np.random.default_rng(0),
-        counts,
-        epsilon_levels=np.array([np.inf] * 3),
-        jitter_strategy='refine',
-    )
-    # Cell indices map 1:1 to values here (delta == 1), so compare directly.
-    expected = np.quantile(data, [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875])
-    np.testing.assert_allclose(edges, expected, atol=1.0)
-
-  def test_spike_owns_consecutive_edges(self):
-    # When a single value holds a majority of the mass, the quantiles on both
-    # sides of the median should collapse onto that value. This is the key
-    # correctness property the deterministic whole-cell split got wrong.
-    counts = np.zeros(101, dtype=np.int64)
-    counts[:40] = 20
-    counts[40] = 20000  # ~96% of the mass.
-    counts[41:80] = 20
-    edges = _quantiles.quantiles_from_histogram(
-        np.random.default_rng(0),
-        counts,
-        epsilon_levels=np.array([np.inf] * 3),
-        jitter_strategy='refine',
-    )
-    # The interior quantiles (0.25 through 0.75) must all be cell 40.
-    self.assertEqual(edges[1:6], [40, 40, 40, 40, 40])
-
-
-if __name__ == '__main__':
-  absltest.main()
diff --git a/tests/local_mode/initialization_test.py b/tests/local_mode/initialization_test.py
index 5d02c5a..0147a3f 100644
--- a/tests/local_mode/initialization_test.py
+++ b/tests/local_mode/initialization_test.py
@@ -16,8 +16,8 @@
 from absl.testing import parameterized
 import dp_accounting
 from dpsynth import domain
-from dpsynth.local_mode import _quantiles
 from dpsynth.local_mode import initialization
+from dpsynth.local_mode import primitives
 from dpsynth.local_mode import vectorized_transformations as vtx
 import numpy as np
 
@@ -134,7 +134,7 @@ def test_int_grid_reserves_budget_for_jitter_refinement(self):
     init = initialization.NumericalInitializerConfig(
         num_partitions=64, max_grid_size=max_grid_size
     )
-    m = _quantiles.jitter_factor(init.num_partitions)
+    m = primitives.jitter_factor(init.num_partitions)
     self.assertLessEqual(
         init.configure(attr, zcdp_rho=1.0).grid_size * m, max_grid_size
     )
diff --git a/tests/local_mode/primitives_test.py b/tests/local_mode/primitives_test.py
index 8fa5485..948a7fe 100644
--- a/tests/local_mode/primitives_test.py
+++ b/tests/local_mode/primitives_test.py
@@ -12,9 +12,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-"""Tests for quantiles primitives."""
-
-import unittest
+"""Tests for differentially private primitives."""
 
 from absl.testing import absltest
 from absl.testing import parameterized
@@ -22,107 +20,191 @@
 import numpy as np
 
 
-@unittest.skip(
-    "SIPS tests are broken at HEAD; will be replaced by Gaussian partition"
-    " selection."
-)
-class SelectPartitionsSipsTest(parameterized.TestCase):
+class ExponentialMechanismTest(parameterized.TestCase):
 
   def setUp(self):
     super().setUp()
     self.rng = np.random.default_rng(42)
 
-  def test_basic_operation(self):
-    data = np.array([1] * 50 + [2] * 5)
-    selected, counts, sigma = primitives._select_partitions_sips(
-        self.rng, data, gdp_budget=10.0, delta=1e-5
+  def test_basic_selection(self):
+    scores = np.array([5.0, 20.0, -10.0, 3.0])
+    idx = primitives.exponential_mechanism(
+        self.rng, scores, epsilon=1.0, sensitivity=1.0
     )
-    self.assertIn(1, selected)
-    self.assertEqual(sigma, 1.0 / np.sqrt(10.0))
-    self.assertEqual(selected.size, counts.size)
+    self.assertIn(idx, [0, 1, 2, 3])
 
-  def test_empty_data(self):
-    data = np.array([], dtype=int)
-    selected, counts, sigma = primitives._select_partitions_sips(
-        self.rng, data, gdp_budget=1.0, delta=1e-5
+  def test_infinite_budget_selects_max_score(self):
+    scores = np.array([5.0, 20.0, -10.0, 3.0])
+    idx = primitives.exponential_mechanism(
+        self.rng, scores, epsilon=np.inf, sensitivity=1.0
     )
-    self.assertEmpty(selected)
-    self.assertEmpty(counts)
-    self.assertEqual(sigma, 1.0)
+    self.assertEqual(idx, 1)
 
-  def test_infinite_budget(self):
-    data = np.array([1, 2, 3, 4, 5])
-    selected, counts, sigma = primitives._select_partitions_sips(
-        self.rng, data, gdp_budget=np.inf, delta=0.1
+  def test_high_budget_selects_max_score(self):
+    scores = np.array([5.0, 20.0, -10.0, 3.0])
+    idx = primitives.exponential_mechanism(
+        self.rng, scores, epsilon=100.0, sensitivity=1.0
     )
-    self.assertCountEqual(selected, [1, 2, 3, 4, 5])
-    self.assertEqual(sigma, 0.0)
-    np.testing.assert_array_equal(counts, np.ones(5))
+    self.assertEqual(idx, 1)
 
-  def test_zero_budget_raises(self):
-    data = np.array([1, 2, 3])
+  def test_zero_budget_uniform_distribution(self):
+    scores = np.array([100.0, 0.0, -100.0])
+    selections = [
+        primitives.exponential_mechanism(
+            self.rng, scores, epsilon=0.0, sensitivity=1.0
+        )
+        for _ in range(1500)
+    ]
+    counts = np.bincount(selections, minlength=3)
+    for c in counts:
+      self.assertBetween(c, 400, 600)
+
+  def test_monotonic_scaling(self):
+    scores = np.array([1.0, 2.0])
+    rng1 = np.random.default_rng(123)
+    rng2 = np.random.default_rng(123)
+    idx1 = primitives.exponential_mechanism(
+        rng1, scores, epsilon=1.0, sensitivity=1.0, monotonic=True
+    )
+    idx2 = primitives.exponential_mechanism(
+        rng2, scores, epsilon=2.0, sensitivity=1.0, monotonic=False
+    )
+    self.assertEqual(idx1, idx2)
+
+  def test_invalid_inputs_raise(self):
+    scores = np.array([1.0, 2.0])
+    with self.assertRaises(ValueError):
+      primitives.exponential_mechanism(
+          self.rng, scores, epsilon=-1.0, sensitivity=1.0
+      )
+    with self.assertRaises(ValueError):
+      primitives.exponential_mechanism(
+          self.rng, scores, epsilon=1.0, sensitivity=0.0
+      )
     with self.assertRaises(ValueError):
-      primitives._select_partitions_sips(
-          self.rng, data, gdp_budget=-0.1, delta=1e-5
+      primitives.exponential_mechanism(
+          self.rng, scores, epsilon=1.0, sensitivity=-0.5
       )
     with self.assertRaises(ValueError):
-      primitives._select_partitions_sips(
-          self.rng, data, gdp_budget=1.0, delta=-0.001
+      primitives.exponential_mechanism(
+          self.rng, np.array([]), epsilon=1.0, sensitivity=1.0
       )
 
-  def test_string_data_type(self):
-    data = np.array(["a", "b", "a", "c"])
-    selected, _, _ = primitives._select_partitions_sips(
-        self.rng, data, gdp_budget=10.0, delta=1e-5
-    )
-    self.assertTrue(all(isinstance(p, str) for p in selected))
 
-  def test_user_level_dp_weighting(self):
-    # Partition 1 has 10 unique users (1 to 10), each contributing 1 time.
-    # Partition 2 has 1 user (11) contributing 10 times.
-    data = np.array([1] * 10 + [2] * 10)
-    user_ids = np.array(list(range(1, 11)) + [11] * 10)
+class JitterFactorTest(absltest.TestCase):
 
-    selected, _, _ = primitives._select_partitions_sips(
-        self.rng, data, gdp_budget=10.0, delta=1e-5, user_ids=user_ids
-    )
-    self.assertIn(1, selected)
-    self.assertNotIn(2, selected)
+  def test_jitter_factor_calculation(self):
+    self.assertEqual(primitives.jitter_factor(0), 1)
+    self.assertEqual(primitives.jitter_factor(1), 4)
+    self.assertEqual(primitives.jitter_factor(16), 64)
+
+
+class QuantilesFromHistogramTest(parameterized.TestCase):
+
+  def test_no_levels_returns_empty(self):
+    rng = np.random.default_rng(0)
+    counts = np.array([10])
+    for jitter_strategy in ("symmetric", "refine"):
+      edges = primitives.quantiles_from_histogram(
+          rng, counts, np.array([]), jitter_strategy
+      )
+      self.assertEmpty(edges)
 
-  @parameterized.named_parameters(
-      ("item_level_default_rounds", None, None),
-      ("item_level_3_rounds", None, 3),
-      ("user_level_default_rounds", np.array([1, 2, 3]), None),
-      ("user_level_5_rounds", np.array([1, 2, 3]), 5),
+  @parameterized.product(
+      levels=(1, 2, 3, 4),
+      jitter_strategy=("symmetric", "refine"),
   )
-  def test_configurations(self, user_ids, num_rounds):
-    data = np.array([1, 2, 3])
-    gdp_budget = 10.0
-    _, _, sigma = primitives._select_partitions_sips(
-        self.rng,
-        data,
-        gdp_budget=gdp_budget,
-        delta=1e-5,
-        num_rounds=num_rounds,
-        user_ids=user_ids,
+  def test_edge_count_matches_levels(self, levels, jitter_strategy):
+    rng = np.random.default_rng(0)
+    grid_size = 10001
+    counts = rng.integers(0, 20, size=grid_size)
+    edges = primitives.quantiles_from_histogram(
+        rng,
+        counts,
+        epsilon_levels=np.ones(levels),
+        jitter_strategy=jitter_strategy,
     )
-    # Calculate expected max_sigma based on budget allocation
-    if num_rounds is None:
-      num_rounds = 1 if user_ids is None else 3
-    allocation_factor = 0.3  # default in primitives.py
-    fractions = allocation_factor ** np.arange(num_rounds)[::-1]
-    fractions /= fractions.sum()
-    gdp_rounds = gdp_budget * fractions
-    expected_max_sigma = float(np.max(1.0 / np.sqrt(gdp_rounds)))
-
-    self.assertAlmostEqual(sigma, expected_max_sigma)
-
-  def test_mismatched_user_ids_raises(self):
-    data = np.array([1, 2, 3])
-    user_ids = np.array([1, 2])
-    with self.assertRaises(ValueError):
-      primitives._select_partitions_sips(
-          self.rng, data, gdp_budget=10.0, delta=1e-5, user_ids=user_ids
+    self.assertLen(edges, 2**levels - 1)
+
+  @parameterized.parameters(1, 2, 3, 4)
+  def test_edge_count_matches_levels_with_spike(self, levels):
+    counts = np.zeros(101, dtype=np.int64)
+    counts[:40] = 1
+    counts[40] = 1000
+    counts[41:80] = 1
+    edges = primitives.quantiles_from_histogram(
+        np.random.default_rng(0),
+        counts,
+        epsilon_levels=np.array([np.inf] * levels),
+        jitter_strategy="refine",
+    )
+    self.assertLen(edges, 2**levels - 1)
+
+  def test_integer_edges_are_integer_indices(self):
+    counts = np.zeros(101, dtype=np.int64)
+    counts[40] = 5000
+    counts[:40] = 50
+    counts[41:] = 50
+    edges = primitives.quantiles_from_histogram(
+        np.random.default_rng(0),
+        counts,
+        epsilon_levels=np.array([np.inf] * 3),
+        jitter_strategy="refine",
+    )
+    for edge in edges:
+      self.assertEqual(edge, int(edge))
+      self.assertBetween(edge, 0, counts.size - 1)
+
+  def test_exact_budget_matches_numpy_smooth(self):
+    rng = np.random.default_rng(0)
+    data = rng.integers(0, 100, size=50000)
+    counts = np.bincount(data, minlength=101)
+    edges = primitives.quantiles_from_histogram(
+        rng,
+        counts,
+        epsilon_levels=np.array([np.inf] * 3),
+        jitter_strategy="refine",
+    )
+    expected = np.quantile(data, [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875])
+    np.testing.assert_allclose(edges, expected, atol=1.0)
+
+  def test_exact_budget_matches_numpy_with_spike(self):
+    below = np.arange(1, 40).repeat(230)
+    spike = np.full(13500, 40)
+    above = np.arange(41, 80).repeat(190)
+    data = np.concatenate([below, spike, above])
+    counts = np.bincount(data, minlength=101)
+    edges = primitives.quantiles_from_histogram(
+        np.random.default_rng(0),
+        counts,
+        epsilon_levels=np.array([np.inf] * 3),
+        jitter_strategy="refine",
+    )
+    expected = np.quantile(data, [0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875])
+    np.testing.assert_allclose(edges, expected, atol=1.0)
+
+  def test_spike_owns_consecutive_edges(self):
+    counts = np.zeros(101, dtype=np.int64)
+    counts[:40] = 20
+    counts[40] = 20000  # ~96% of the mass.
+    counts[41:80] = 20
+    edges = primitives.quantiles_from_histogram(
+        np.random.default_rng(0),
+        counts,
+        epsilon_levels=np.array([np.inf] * 3),
+        jitter_strategy="refine",
+    )
+    self.assertEqual(edges[1:6], [40, 40, 40, 40, 40])
+
+  def test_unsupported_max_records_per_user_raises(self):
+    rng = np.random.default_rng(0)
+    with self.assertRaises(NotImplementedError):
+      primitives.quantiles_from_histogram(
+          rng,
+          np.array([1, 2, 3]),
+          epsilon_levels=np.ones(2),
+          jitter_strategy="refine",
+          max_records_per_user=2,
       )
 
 
@@ -162,8 +244,6 @@ def test_high_budget_selects_all(self):
     self.assertCountEqual(selected_partitions, [1, 2, 3, 4, 5])
 
   def test_rare_items_not_selected(self):
-    # One item with many occurrences, another with just 1.
-    # With moderate budget and tight delta, the rare item should be dropped.
     data = np.array([1] * 100 + [2])
     selected_partitions, _, _ = (
         primitives.select_partitions_gaussian_thresholding(
@@ -183,7 +263,6 @@ def test_string_data_type(self):
     self.assertTrue(all(isinstance(p, str) for p in selected_partitions))
 
   def test_min_count_filters_low_count_partitions(self):
-    # Partition 1 has count 50, partition 2 has count 3.
     data = np.array([1] * 50 + [2] * 3)
     selected, _, _ = primitives.select_partitions_gaussian_thresholding(
         self.rng, data, gdp_budget=10.0, delta=1e-5, min_count=5
@@ -220,8 +299,6 @@ def test_min_count_zero_raises(self):
       )
 
   def test_min_count_increases_threshold(self):
-    # With very high budget (no noise), threshold is approximately min_count.
-    # Partitions with count exactly at min_count should pass.
     data = np.array([1] * 10 + [2] * 10)
     selected, _, _ = primitives.select_partitions_gaussian_thresholding(
         self.rng, data, gdp_budget=np.inf, delta=0.1, min_count=10
@@ -243,9 +320,9 @@ def test_missing_partitions_appended_and_sorted(self):
         self.rng, selected, counts, 0.0, public
     )
     np.testing.assert_array_equal(sel, ["a", "b", "c"])
-    self.assertEqual(cts[0], 10.0)  # count for 'a'
-    self.assertEqual(cts[1], 0.0)  # noise for 'b' (stddev=0)
-    self.assertEqual(cts[2], 20.0)  # count for 'c'
+    self.assertEqual(cts[0], 10.0)
+    self.assertEqual(cts[1], 0.0)
+    self.assertEqual(cts[2], 20.0)
 
   def test_all_present_is_noop(self):
     selected = np.array(["a", "b"])