zorch.sumcheck.domain

Evaluation domains for sumcheck round polynomials.

A round polynomial goes on the wire as ascending coefficients (the form verifier.CoeffsSumcheckRound checks). EvalDomain names the points a prover samples it at and owns the map from those samples to coefficients: a finite node set — the Gruen set {0, 1, *extra, eq_root(z)}, or the naturals — optionally led by the value at infinity (the leading coefficient, cheap for a product since it is the product of the factor slopes). extend_to_round_domain lifts a linear pair onto the Û_d sample domain; product_round_poly / product_round_coeffs build a product round message over it (the prover Rounds that emit them live in sqrt_space).

EvalDomain dataclass

The sample points of a round polynomial and the map from those samples to ascending coefficients.

nodes are the finite sample points — the Gruen set {0, 1, *extra, eq_root(z)}, or (when None) the naturals {0..}. inf_index is the index at which the value at infinity — the leading coefficient (cheap for a product: the product of the factor slopes) — sits among the samples: 0 first (Û), -1 last (the compressed [s(node), s(∞)] wire), None for no ∞ sample. Degree and field come from the samples, so only the node shape is fixed here.

Source code in zorch/sumcheck/domain.py

@dataclass(frozen=True)
class EvalDomain:
    """The sample points of a round polynomial and the map from those samples to
    ascending coefficients.

    nodes are the finite sample points — the Gruen set {0, 1, *extra, eq_root(z)},
    or (when None) the naturals {0..}. inf_index is the index at which the value
    at infinity — the leading coefficient (cheap for a product: the product of the
    factor slopes) — sits among the samples: 0 first (Û), -1 last (the compressed
    `[s(node), s(∞)]` wire), None for no ∞ sample. Degree and field come from the
    samples, so only the node shape is fixed here."""

    nodes: Array | None = None
    inf_index: int | None = None

    def coeff_matrix(self) -> Array:
        """Value → coefficient matrix for the explicit finite nodes — what a driver
        precomputes once per round. The leading / naturals domain has no fixed size,
        so it reads its degree off the values instead: use to_coeffs."""
        assert self.nodes is not None, "leading / naturals domain has no fixed matrix"
        return _finite_coeff_matrix(self.nodes)

    def to_coeffs(self, values: Array) -> Array:
        """Ascending coefficients from this domain's samples of a round polynomial;
        the degree (len−1) and field come from values."""
        if self.inf_index is None:
            return fnp.dot(self.coeff_matrix(), values)
        # The ∞ sample (at `inf_index`) is the leading coefficient c_d; the finite
        # samples (naturals {0..d−1} unless given) interpolate the residual p − c_d·xᵈ.
        pos = self.inf_index % values.shape[0]
        v_inf = values[pos]
        finite = fnp.concatenate([values[:pos], values[pos + 1 :]])
        d = finite.shape[0]
        if self.nodes is None:
            # Naturals domain: _finite_coeff_matrix(naturals(d)) is exactly the
            # inverse Vandermonde — its Lagrange build over the naturals collapses to
            # the identity — so skip that identity vmap + matmul and use inv_vand.
            nodes, cmat = _interp_constants(d - 1, values.dtype)
        else:
            nodes, cmat = self.nodes, _finite_coeff_matrix(self.nodes)
        low = fnp.dot(cmat, finite - v_inf * nodes**d)
        return fnp.concatenate([low, fnp.atleast_1d(v_inf)])

    def sample(self, p0: Array, p1: Array) -> Array:
        """Sample the linear factor p(x) = p0 + x·(p1−p0) at this domain's points
        [∞ if leading, *nodes]: p(∞) = slope p1−p0, p(node) = p0 + node·slope. The
        leading axis indexes the domain. Requires explicit nodes — a sampling domain
        must be concrete (the naturals-sized domain is an output-only coeff map)."""
        assert self.nodes is not None, "a sampling domain needs an explicit node set"
        diff = p1 - p0
        finite = p0[None] + self.nodes.reshape((-1,) + (1,) * p0.ndim) * diff[None]
        if self.inf_index is None:
            return finite
        pos = self.inf_index % (finite.shape[0] + 1)
        return fnp.concatenate([finite[:pos], diff[None], finite[pos:]])

coeff_matrix

coeff_matrix() -> Array

Value → coefficient matrix for the explicit finite nodes — what a driver precomputes once per round. The leading / naturals domain has no fixed size, so it reads its degree off the values instead: use to_coeffs.

Source code in zorch/sumcheck/domain.py

def coeff_matrix(self) -> Array:
    """Value → coefficient matrix for the explicit finite nodes — what a driver
    precomputes once per round. The leading / naturals domain has no fixed size,
    so it reads its degree off the values instead: use to_coeffs."""
    assert self.nodes is not None, "leading / naturals domain has no fixed matrix"
    return _finite_coeff_matrix(self.nodes)

to_coeffs

to_coeffs(values: Array) -> Array

Ascending coefficients from this domain's samples of a round polynomial; the degree (len−1) and field come from values.

Source code in zorch/sumcheck/domain.py

def to_coeffs(self, values: Array) -> Array:
    """Ascending coefficients from this domain's samples of a round polynomial;
    the degree (len−1) and field come from values."""
    if self.inf_index is None:
        return fnp.dot(self.coeff_matrix(), values)
    # The ∞ sample (at `inf_index`) is the leading coefficient c_d; the finite
    # samples (naturals {0..d−1} unless given) interpolate the residual p − c_d·xᵈ.
    pos = self.inf_index % values.shape[0]
    v_inf = values[pos]
    finite = fnp.concatenate([values[:pos], values[pos + 1 :]])
    d = finite.shape[0]
    if self.nodes is None:
        # Naturals domain: _finite_coeff_matrix(naturals(d)) is exactly the
        # inverse Vandermonde — its Lagrange build over the naturals collapses to
        # the identity — so skip that identity vmap + matmul and use inv_vand.
        nodes, cmat = _interp_constants(d - 1, values.dtype)
    else:
        nodes, cmat = self.nodes, _finite_coeff_matrix(self.nodes)
    low = fnp.dot(cmat, finite - v_inf * nodes**d)
    return fnp.concatenate([low, fnp.atleast_1d(v_inf)])

sample

sample(p0: Array, p1: Array) -> Array

Sample the linear factor p(x) = p0 + x·(p1−p0) at this domain's points [∞ if leading, *nodes]: p(∞) = slope p1−p0, p(node) = p0 + node·slope. The leading axis indexes the domain. Requires explicit nodes — a sampling domain must be concrete (the naturals-sized domain is an output-only coeff map).

Source code in zorch/sumcheck/domain.py

def sample(self, p0: Array, p1: Array) -> Array:
    """Sample the linear factor p(x) = p0 + x·(p1−p0) at this domain's points
    [∞ if leading, *nodes]: p(∞) = slope p1−p0, p(node) = p0 + node·slope. The
    leading axis indexes the domain. Requires explicit nodes — a sampling domain
    must be concrete (the naturals-sized domain is an output-only coeff map)."""
    assert self.nodes is not None, "a sampling domain needs an explicit node set"
    diff = p1 - p0
    finite = p0[None] + self.nodes.reshape((-1,) + (1,) * p0.ndim) * diff[None]
    if self.inf_index is None:
        return finite
    pos = self.inf_index % (finite.shape[0] + 1)
    return fnp.concatenate([finite[:pos], diff[None], finite[pos:]])

UnivariateSkipDomain dataclass

The univariate skip's round-0 evaluation domain: the order-2^skip_rounds two-adic subgroup D, with the value↔coeff map carried by the NTT (subgroup_to_coeffs / subgroup_evals) rather than the inverse Vandermonde a finite EvalDomain uses. Holds the knob (skip_rounds, the number of collapsed leading rounds) and the field; the round-0 arithmetic runs in the base field (nodes and the transforms are base-field), extension arithmetic entering only once r₀ is bound at round 1.

Source code in zorch/sumcheck/domain.py

@dataclass(frozen=True)
class UnivariateSkipDomain:
    """The univariate skip's round-0 evaluation domain: the order-2^skip_rounds two-adic
    subgroup D, with the value↔coeff map carried by the NTT (`subgroup_to_coeffs` /
    `subgroup_evals`) rather than the inverse Vandermonde a finite `EvalDomain` uses.
    Holds the knob (`skip_rounds`, the number of collapsed leading rounds) and the
    field; the round-0 arithmetic runs in the base field (`nodes` and the transforms are
    base-field), extension arithmetic entering only once r₀ is bound at round 1."""

    skip_rounds: int
    dtype: Any

    @property
    def size(self) -> int:
        """|D| = 2^skip_rounds."""
        return 1 << self.skip_rounds

    def nodes(self) -> Array:
        """The |D| subgroup points in `lax.ntt` order — `ntt(e₁)` reads them off the
        same transform the map uses (mirrors `coding.reed_solomon.eval_domain`)."""
        base = base_field(self.dtype)
        if self.size == 1:
            return fnp.ones((1,), base)
        e1 = fnp.zeros(self.size, base).at[1].set(fnp.ones((), base))
        return lax.ntt(e1, ntt_type="NTT", ntt_length=self.size)

    def to_coeffs(self, values: Array) -> Array:
        """Value→coeff over D via the iNTT (`subgroup_to_coeffs`)."""
        return subgroup_to_coeffs(values)

    def sum_over_subgroup(self, coeffs: Array) -> Array:
        """Σ_{z∈D} p(z) from ascending coeffs (`subgroup_sum`)."""
        return subgroup_sum(coeffs, self.skip_rounds)

size property

size: int

|D| = 2^skip_rounds.

nodes

nodes() -> Array

The |D| subgroup points in lax.ntt order — ntt(e₁) reads them off the same transform the map uses (mirrors coding.reed_solomon.eval_domain).

Source code in zorch/sumcheck/domain.py

def nodes(self) -> Array:
    """The |D| subgroup points in `lax.ntt` order — `ntt(e₁)` reads them off the
    same transform the map uses (mirrors `coding.reed_solomon.eval_domain`)."""
    base = base_field(self.dtype)
    if self.size == 1:
        return fnp.ones((1,), base)
    e1 = fnp.zeros(self.size, base).at[1].set(fnp.ones((), base))
    return lax.ntt(e1, ntt_type="NTT", ntt_length=self.size)

to_coeffs

to_coeffs(values: Array) -> Array

Value→coeff over D via the iNTT (subgroup_to_coeffs).

Source code in zorch/sumcheck/domain.py

def to_coeffs(self, values: Array) -> Array:
    """Value→coeff over D via the iNTT (`subgroup_to_coeffs`)."""
    return subgroup_to_coeffs(values)

sum_over_subgroup

sum_over_subgroup(coeffs: Array) -> Array

Σ_{z∈D} p(z) from ascending coeffs (subgroup_sum).

Source code in zorch/sumcheck/domain.py

def sum_over_subgroup(self, coeffs: Array) -> Array:
    """Σ_{z∈D} p(z) from ascending coeffs (`subgroup_sum`)."""
    return subgroup_sum(coeffs, self.skip_rounds)

extend_to_round_domain

extend_to_round_domain(
    p0: Array, p1: Array, d: int, *, skip_one: bool = False
) -> Array

Lift a linear pair (p(0), p(1)) onto U_d = [∞, 0, 1, …, d−1] (or Û_d, u=1 dropped, when skip_one). p(∞) is the slope p(1)−p(0); p(u) = p(0) + u·(p(1)−p(0)). The leading axis indexes the domain.

Source code in zorch/sumcheck/domain.py

def extend_to_round_domain(
    p0: Array, p1: Array, d: int, *, skip_one: bool = False
) -> Array:
    """Lift a linear pair (p(0), p(1)) onto U_d = [∞, 0, 1, …, d−1] (or Û_d, u=1
    dropped, when skip_one). p(∞) is the slope p(1)−p(0);
    p(u) = p(0) + u·(p(1)−p(0)). The leading axis indexes the domain."""
    diff = p1 - p0
    if d == 1:  # U_1 = [∞, 0]: u=1 is out of range {0..d−1}, so p1 is not a node
        return fnp.stack([diff, p0])
    base = fnp.stack([diff, p0]) if skip_one else fnp.stack([diff, p0, p1])
    if d == 2:
        return base
    # Python-int multiplier avoids a field-dtype iota (unsupported in frx).
    rest = fnp.stack([p0 + diff * u for u in range(2, d)], axis=0)
    return fnp.concatenate([base, rest], axis=0)

natural_domain

natural_domain(degree: int, dtype: Any) -> EvalDomain

The natural evaluation domain {0, 1, …, degree}: the round poly sent as its plain values [s(0), …, s(degree)] — the wire form verifier.SumcheckRound checks and the default domain of the generic StandardRound. Nodes live in the base field (an integer node is a base-field element; extension factors promote at multiply), reproducing the list prover's per-point lift byte-for-byte.

Source code in zorch/sumcheck/domain.py

def natural_domain(degree: int, dtype: Any) -> EvalDomain:
    """The natural evaluation domain {0, 1, …, degree}: the round poly sent as its
    plain values [s(0), …, s(degree)] — the wire form verifier.SumcheckRound checks
    and the default domain of the generic StandardRound. Nodes live in the base
    field (an integer node is a base-field element; extension factors promote at
    multiply), reproducing the list prover's per-point lift byte-for-byte."""
    return EvalDomain(naturals(degree + 1, dtype))

uhat_domain

uhat_domain(degree: int, dtype: Any) -> EvalDomain

The compressed product round domain Û_degree = {∞, 0, 2, …, degree−1}: ∞-leading, u=1 dropped (the verifier recovers s(1) from s(0)+s(1)=claim). The default sampling domain for the eq-poly / sqrt-space engines.

Source code in zorch/sumcheck/domain.py

def uhat_domain(degree: int, dtype: Any) -> EvalDomain:
    """The compressed product round domain Û_degree = {∞, 0, 2, …, degree−1}:
    ∞-leading, u=1 dropped (the verifier recovers s(1) from s(0)+s(1)=claim). The
    default sampling domain for the eq-poly / sqrt-space engines."""
    nat = naturals(degree, dtype)
    return EvalDomain(fnp.concatenate([nat[:1], nat[2:]]), inf_index=0)

compressed_domain

compressed_domain(node: int, dtype: Any) -> EvalDomain

The two-point compressed product-round domain [s(node), s(∞)]: one finite node (0 → s(0)=c_0, 1 → s(1)) with s(∞) trailing. The third value of the degree-2 round poly stays off the wire — the verifier recovers it from s(0)+s(1)=claim.

Source code in zorch/sumcheck/domain.py

def compressed_domain(node: int, dtype: Any) -> EvalDomain:
    """The two-point compressed product-round domain `[s(node), s(∞)]`: one finite
    node (`0` → `s(0)=c_0`, `1` → `s(1)`) with `s(∞)` trailing. The third value of
    the degree-2 round poly stays off the wire — the verifier recovers it from
    `s(0)+s(1)=claim`."""
    if node not in (0, 1):
        raise ValueError(f"compressed node must be 0 or 1, got {node}")
    return EvalDomain(naturals(node + 1, dtype)[node:], inf_index=-1)

split_halves

split_halves(arr: Array) -> tuple[Array, Array]

Split the last variable MSB-first into contiguous halves (arr[..., :N/2], arr[..., N/2:]) — the dense bind. ndim-agnostic; the MSB dual of split_pairs.

Source code in zorch/sumcheck/domain.py

def split_halves(arr: Array) -> tuple[Array, Array]:
    """Split the last variable MSB-first into contiguous halves
    `(arr[..., :N/2], arr[..., N/2:])` — the dense bind. ndim-agnostic; the MSB
    dual of split_pairs."""
    half = arr.shape[-1] // 2
    return arr[..., :half], arr[..., half:]

split_pairs

split_pairs(arr: Array) -> tuple[Array, Array]

Split the last variable LSB-first into stride-2 consecutive pairs (arr[..., 0::2], arr[..., 1::2]). The jagged engines bind LSB-first — a batch-major layout makes the pair the in-segment dimension, so a fold never crosses a segment boundary — while the dense drivers split MSB-first (split_halves). Its split-only form also serves the LogUp paired_evals, which needs both halves without folding.

Source code in zorch/sumcheck/domain.py

def split_pairs(arr: Array) -> tuple[Array, Array]:
    """Split the last variable LSB-first into stride-2 consecutive pairs
    `(arr[..., 0::2], arr[..., 1::2])`. The jagged engines bind LSB-first — a
    batch-major layout makes the pair the in-segment dimension, so a fold never
    crosses a segment boundary — while the dense drivers split MSB-first
    (split_halves). Its split-only form also serves the LogUp `paired_evals`,
    which needs both halves without folding."""
    return arr[..., 0::2], arr[..., 1::2]

summand_evals

summand_evals(
    stacked: Array,
    combine: Callable[..., Array],
    domain: EvalDomain,
    *,
    weight: Array | None = None,
    msb: bool = True
) -> Array

Σ_x' weight(x')·combine(f₁, …, f_m)(x') per point of domain: the m stacked factors sampled at the domain's points (domain.sample), combined, optionally weighted per hypercube point, then summed. The one reduction body the round-message builders share — generic over the summand's combine (Πₖ, LogUp, …), the evaluation domain (∞-leading Û, Gruen, {0,½}, {0,2,4}, …), and the bind order.

weight (a length-N/2 vector over the folded hypercube) is the eq-weight of an eq-weighted sumcheck — Σ eq·Π in one pass, no eq factor stacked into the state. msb=False binds the LOW variable (split_pairs) instead of the high (split_halves) — the LSB order the jagged / GHASH engines fold in.

A leading (∞) node encodes s(∞) = combine(*slopes), the true leading coefficient only for a HOMOGENEOUS combine — every monomial a product of exactly degree factors (a plain product, or the LogUp combine); a mixed-degree combine like eq·(â◦b̂ − ĉ) is not. A finite domain carries no such restriction — any summand samples cleanly on it.

Source code in zorch/sumcheck/domain.py

def summand_evals(
    stacked: Array,
    combine: Callable[..., Array],
    domain: EvalDomain,
    *,
    weight: Array | None = None,
    msb: bool = True,
) -> Array:
    """Σ_x' weight(x')·combine(f₁, …, f_m)(x') per point of `domain`: the m stacked
    factors sampled at the domain's points (domain.sample), combined, optionally
    weighted per hypercube point, then summed. The one reduction body the
    round-message builders share — generic over the summand's `combine` (Πₖ, LogUp,
    …), the evaluation domain (∞-leading Û, Gruen, {0,½}, {0,2,4}, …), and the bind
    order.

    `weight` (a length-N/2 vector over the folded hypercube) is the eq-weight of an
    eq-weighted sumcheck — `Σ eq·Π` in one pass, no eq factor stacked into the state.
    `msb=False` binds the LOW variable (split_pairs) instead of the high
    (split_halves) — the LSB order the jagged / GHASH engines fold in.

    A leading (∞) node encodes s(∞) = combine(*slopes), the true leading coefficient
    only for a HOMOGENEOUS combine — every monomial a product of exactly `degree`
    factors (a plain product, or the LogUp combine); a mixed-degree combine like
    eq·(â◦b̂ − ĉ) is not. A finite domain carries no such restriction — any summand
    samples cleanly on it."""
    p0, p1 = split_halves(stacked) if msb else split_pairs(stacked)
    combined = combine(*frx.vmap(domain.sample)(p0, p1))
    if weight is not None:
        combined = combined * weight
    return fnp.sum(combined, axis=1)

product_round_poly

product_round_poly(stacked: Array) -> Array

Round message s = Σₓ Πₖ fₖ over Û_m for the m stacked multilinears, shape (m,) — the product summand on the compressed domain.

Source code in zorch/sumcheck/domain.py

def product_round_poly(stacked: Array) -> Array:
    """Round message s = Σₓ Πₖ fₖ over Û_m for the m stacked multilinears, shape
    (m,) — the product summand on the compressed domain."""
    return summand_evals(
        stacked, _product, uhat_domain(stacked.shape[0], stacked.dtype)
    )

product_round_coeffs

product_round_coeffs(stacked: Array) -> Array

Ascending coefficients of the degree-m product round polynomial for m stacked factors: the same Σ_x' Πₖ fₖ as product_round_poly but sampled over the full round domain [∞, 0, 1, …, m−1] so it is fully determined, then mapped to coefficients — the wire form verifier.CoeffsSumcheckRound checks.

Source code in zorch/sumcheck/domain.py

def product_round_coeffs(stacked: Array) -> Array:
    """Ascending coefficients of the degree-m product round polynomial for m stacked
    factors: the same Σ_x' Πₖ fₖ as product_round_poly but sampled over the full round
    domain [∞, 0, 1, …, m−1] so it is fully determined, then mapped to coefficients —
    the wire form verifier.CoeffsSumcheckRound checks."""
    m = stacked.shape[0]
    full = EvalDomain(naturals(m, stacked.dtype), inf_index=0)  # [∞, 0, 1, …, m−1]
    return EvalDomain(inf_index=0).to_coeffs(summand_evals(stacked, _product, full))

fold

fold(arr: Array, r: Array, *, msb: bool = True) -> Array

Fold the last variable at challenge r: P0 + r*(P1 - P0), halving the last axis. msb (the dense default) splits contiguous halves [low | high]; msb=False splits stride-2 consecutive pairs — the jagged bind, where a batch-major layout makes the pair the in-segment dimension, so the fold never crosses a segment boundary. ndim-agnostic: the leading factor/batch axes broadcast.

Source code in zorch/sumcheck/domain.py

def fold(arr: Array, r: Array, *, msb: bool = True) -> Array:
    """Fold the last variable at challenge `r`: P0 + r*(P1 - P0), halving the last
    axis. `msb` (the dense default) splits contiguous halves [low | high]; msb=False
    splits stride-2 consecutive pairs — the jagged bind, where a batch-major layout
    makes the pair the in-segment dimension, so the fold never crosses a segment
    boundary. ndim-agnostic: the leading factor/batch axes broadcast."""
    p0, p1 = split_halves(arr) if msb else split_pairs(arr)
    return p0 + r * (p1 - p0)

subgroup_to_coeffs

subgroup_to_coeffs(values: Array) -> Array

Ascending coefficients of the degree-<|D| univariate whose values on the order-|D| two-adic subgroup D are values (last axis = D, |D| a power of two): the iNTT over D, batched over the leading axes (the native op transforms the last axis, like coding.reed_solomon.encode). The subgroup analogue of EvalDomain.to_coeffs — the inverse Vandermonde does not scale to |D| = 2^skip_rounds.

Source code in zorch/sumcheck/domain.py

def subgroup_to_coeffs(values: Array) -> Array:
    """Ascending coefficients of the degree-<|D| univariate whose values on the
    order-|D| two-adic subgroup D are `values` (last axis = D, |D| a power of two):
    the iNTT over D, batched over the leading axes (the native op transforms the last
    axis, like `coding.reed_solomon.encode`). The subgroup analogue of
    `EvalDomain.to_coeffs` — the inverse Vandermonde does not scale to
    |D| = 2^skip_rounds."""
    return lax.ntt(values, ntt_type="INTT", ntt_length=values.shape[-1])

subgroup_evals

subgroup_evals(coeffs: Array, size: int) -> Array

Evaluate the univariate coeffs (ascending, last axis) on the order-size two-adic subgroup, batched over the leading axes — zero-pad to size then NTT. With size > len(coeffs) this is the low-degree extension onto a superset D' ⊇ D the round-0 message needs, since deg s₀ = degree·(|D|−1) outgrows |D|; size must be a power of two ≥ len(coeffs).

Source code in zorch/sumcheck/domain.py

def subgroup_evals(coeffs: Array, size: int) -> Array:
    """Evaluate the univariate `coeffs` (ascending, last axis) on the order-`size`
    two-adic subgroup, batched over the leading axes — zero-pad to `size` then NTT.
    With `size > len(coeffs)` this is the low-degree extension onto a superset D' ⊇ D
    the round-0 message needs, since `deg s₀ = degree·(|D|−1)` outgrows |D|; `size` must
    be a power of two ≥ len(coeffs)."""
    pad = coeffs.shape[-1]
    if size < pad:
        raise ValueError(f"subgroup_evals size {size} < coeff count {pad}")
    tail = coeffs.shape[:-1] + (size - pad,)
    padded = fnp.concatenate([coeffs, fnp.zeros(tail, coeffs.dtype)], axis=-1)
    return lax.ntt(padded, ntt_type="NTT", ntt_length=size)

subgroup_sum

subgroup_sum(coeffs: Array, skip_rounds: int) -> Array

Σ_{z∈D} p(z) for the univariate coeffs (ascending, last axis) over the order-2^skip_rounds subgroup D. Σ_{z∈D} zᵏ = |D| when |D| divides k, else 0 — so the sum reads off the coefficients at multiples of |D|, scaled by |D|. |D| is built in the field (a bare Python int would not Montgomery-encode).

Source code in zorch/sumcheck/domain.py

def subgroup_sum(coeffs: Array, skip_rounds: int) -> Array:
    """Σ_{z∈D} p(z) for the univariate `coeffs` (ascending, last axis) over the
    order-2^skip_rounds subgroup D. Σ_{z∈D} zᵏ = |D| when |D| divides k, else 0 — so the
    sum reads off the coefficients at multiples of |D|, scaled by |D|. |D| is built in
    the field (a bare Python int would not Montgomery-encode)."""
    size = 1 << skip_rounds
    d_field = fnp.asarray(size, coeffs.dtype)
    return d_field * fnp.sum(coeffs[..., ::size], axis=-1)