diff --git a/AGENTS.md b/AGENTS.md
index 7001c12..602900b 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -34,6 +34,8 @@ src/morphobricks_normalization/
   bin_centres.py    float_bin_centres: bin centres that do not overflow near the float64 limits
   bin_reachability.py  float_bin_indices, reachable_float_bins, centres_map_back: which float
                     bins a value can fall in and whether centres map back (#86)
+  integer_representatives.py  integer bin representatives (nearest float64 midpoint) and
+                    their exact bound for 64-bit values float64 cannot hold (#28)
   ranges.py         HistogramSettings, RangeCheck, range_check_issues (Milestone 5)
   integrity.py      invariants every PlaneHistogram and HistogramTail checks when made
                     (counts, limits, tail sums) and the pixel type of every block
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1317811..a7163e7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -4,6 +4,59 @@ Release notes and migration notes of `morphobricks-normalization`. Numbers are p
 issues of [CellMigrationLab/morphobricks-normalization](https://github.com/CellMigrationLab/morphobricks-normalization).
 The decisions are in [docs/plan.md](docs/plan.md) §12.
 
+## Unreleased
+
+Numerical-correctness fixes. Each says which outputs change; every other output is
+bit-identical (`tools/dump_outputs.py` and `--corpus` lines, and the 400 golden fits).
+
+- **64-bit integers float64 cannot hold (#28).** A bin's representative is the float64 nearest
+  its midpoint, rounded once (it was rounded twice, `binning_min` first, which could miss a
+  representable value by 2 above `2**53`). `max_representative_error` of an integer plane
+  with a value of magnitude `2**51` or more is the exact largest distance between its
+  integers and the float64 standing in for them, rounded up, so a histogram is exact only
+  when every representative equals its value: `[2**53 + 1]` is no longer reported exact with
+  error 0 (its representative is `2**53`; error 1), and `normalization_exact` and
+  `Histogram_Exact_C<i>` are false for it. Only int64 and uint64 planes with such values
+  change; bin counts and limits never do.
+
+  **Who is affected, and by how much.** Only int64 and uint64 images with a pixel of
+  magnitude `2**51` (about `2.25e15`) or more. Every other pixel type (8-, 16- and 32-bit
+  integers, float32, float64), and 64-bit images whose values all stay below `2**51`, give
+  bit-identical numbers. For the affected images:
+  - representatives move by at most one float64 step (at most 1536 for int64 and 2048 for
+    uint64 values near `9e18` and `1.8e19`, about `1e-16` relative);
+  - `max_representative_error` is the true largest distance: for bins of one value float64
+    cannot hold it grows from 0 to at most half a float64 step (for example 1 for
+    `[2**53 + 1]`, 170 and 341 for the constant planes `3074457345618258602` and
+    `6148914691236517205`), so `is_exact`, `normalization_exact` and
+    `Histogram_Exact_C<i>` change from true to false; for wider bins it changes in its last
+    bits only;
+  - pooled means, standard deviations and percentiles, and fitted offsets and scales, move by
+    at most about `3e-14` relative in the `tools/dump_outputs.py` cases (a zscore offset
+    `9618469625711740` becomes `9618469625711942`, nearer the exact model mean
+    `9618469625711938`).
+
+  Statistics files need not be recomputed: they store counts and limits, which do not change,
+  and representatives and bounds are derived when the file is read, so an existing
+  `plane_histograms.csv` gives the new numbers. Parameters fitted with v0.4.0 from such
+  images differ from a refit by the amounts above.
+- **Standard deviations of large finite values (#33).** `PooledHistogram.std` and
+  `WeightedHistogram.std` summed squared deviations that overflow for values beyond about
+  `1e154` (`[-1e200, 1e200]` gave `inf`, so `zscore` was refused, although its standard
+  deviation is `1e200`; a tail could raise a bare `OverflowError`). When the plain sum is not
+  finite but the mean and every value are, it is redone with the values divided by a power of
+  two, which is exact; every finite result is unchanged bit for bit. A tail whose stored sum
+  of squares is infinite (schema version 4 keeps raw sums of squares) is refused by name
+  instead of giving `inf`: its spread is not known. No dump output changes. Not changed: the
+  raw sums of squares of the file lose precision by cancellation for tails of large values
+  with a small spread; a centred representation needs a new schema version.
+- **Missing conditions per time point (#41).** With `scope="biological_repeat"` and
+  `per_time_point=True`, `repeat_missing_conditions` checked each repeat group over all its
+  time points, not the groups that are fitted, so a group such as `r1::t=0` with control
+  images only gave no warning when treated images of `r1` existed at another time point.
+  The groups that are fitted are now checked, and the warning's items and message name them
+  (`r1::t=0`). Without `per_time_point` the warning is unchanged. No dump output changes.
+
 ## v0.4.0 (2026-11-10)
 
 Breaking release: `PlaneRecord.position_id` is required (NM-04) and the readers of old table layouts
diff --git a/docs/contract.md b/docs/contract.md
index 90ab53d..47b8abd 100644
--- a/docs/contract.md
+++ b/docs/contract.md
@@ -120,8 +120,14 @@ Notation in this section: `n` = `bin_count`, `lo` = `binning_min`, `hi` = `binni
 - If `W > n`: bin `i` holds the values `lo + ceil(i·W/n)` to `lo + ceil((i+1)·W/n) − 1`.
   Limits are whole numbers, consecutive bins touch without gaps or overlaps, and widths
   are `floor(W/n)` or `ceil(W/n)` values.
-- Representative value of a bin: `(first value + last value) / 2`.
+- Representative value of a bin: `(first value + last value) / 2`, the float64 nearest it.
 - Bound: every pixel is within `(widest bin − 1) / 2` of its representative.
+- Planes with a value of magnitude `2**51` or more (64-bit types, #28): float64 cannot hold
+  every value there (`2**53 + 1` is stored as `2**53`), so the bound is the exact largest
+  distance between a bin's (or tail's) integers and the float64 standing in for them, rounded
+  up. A bin is exact only when its representative equals its value, so `normalization_exact`
+  is false for a group with such a value. Counts and bin limits stay exact integers. Below
+  `2**51` representatives and bound are exact float64 arithmetic, unchanged.
 
 ##### Float planes
 
@@ -192,7 +198,7 @@ tail's mean and its count.
 | Statistic | Definition | Accuracy |
 |---|---|---|
 | min, max | Smallest `pixel_min`, largest `pixel_max` | Exact |
-| mean, std | Count-weighted mean and population std (÷N) of representatives; tails count with their exact sums and sums of squares. The weighted sums are correctly rounded sums of the products (`math.fsum`), so they never depend on the CPU thread count or BLAS library (#67) | Exact for exact histograms; otherwise very close (errors cancel) |
+| mean, std | Count-weighted mean and population std (÷N) of representatives; tails count with their exact sums and sums of squares. The weighted sums are correctly rounded sums of the products (`math.fsum`), so they never depend on the CPU thread count or BLAS library (#67). When the squares overflow float64 although every value is finite (`[-1e200, 1e200]`), the std is computed with the values divided by a power of two, exactly, so a finite spread is never reported as infinite (#33) | Exact for exact histograms; otherwise very close (errors cancel) |
 | percentile q | Smallest representative with cumulative count ≥ ceil(q/100 · N) (numpy `inverted_cdf`) | Within the group's largest per-plane bound |
 | median | 50th percentile | As percentiles |
 | MAD | Median of \|representative − median\| | Within twice the bound |
@@ -254,6 +260,10 @@ pixels keeps its extrema as limits; a 512 × 512 plane has 26 pixels in each tai
 - Percentiles, median and MAD: each non-empty tail stands as one value, its mean, with its
   count.
 - Mean and standard deviation use the tails' sums: as exact as without tails.
+- A tail whose sum of squares is not finite (its pixels' squares overflowed float64 when the
+  file was written, values beyond about `1.3e154`) has no known spread: the standard
+  deviation, and so `zscore`, is refused with a message naming the tail (#33). Mean,
+  percentiles, median and MAD do not use it.
 - Min-max parameters use the exact extrema (`pixel_min`, `pixel_max`).
 - `max_representative_error` is the larger of the bins' bound and, for each non-empty tail,
   the larger distance from the tail's mean to its two ends (`max(mean − pixel_min,
@@ -456,7 +466,9 @@ y_c    = (x_c − offset) / scale          for c = R, G, B
   forms a group with the other conditions' `single` rows, or with the repeat it is paired
   with.
 - **Incomplete groups:** a repeat group without rows from every condition of the statistics
-  gives the warning `repeat_missing_conditions` (§9.2) and is still used.
+  gives the warning `repeat_missing_conditions` (§9.2) and is still used. With
+  `per_time_point=true` each time point's group is checked on its own and named in the
+  warning as `<label>::t=<t>` (#41).
 
 ### 4.3 Pooling weights
 
diff --git a/pyproject.toml b/pyproject.toml
index eec9489..d5f89c8 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -103,6 +103,7 @@ baseline = "tools/code_rules_baseline.json"
 "morphobricks_normalization.fitting" = "pure"
 "morphobricks_normalization.histogram" = "pure"
 "morphobricks_normalization.integrity" = "pure"
+"morphobricks_normalization.integer_representatives" = "pure"
 "morphobricks_normalization.issues" = "pure"
 "morphobricks_normalization.normalization" = "pure"
 "morphobricks_normalization.parameter_checks" = "pure"
diff --git a/src/morphobricks_normalization/histogram.py b/src/morphobricks_normalization/histogram.py
index da6c977..36a3f9e 100644
--- a/src/morphobricks_normalization/histogram.py
+++ b/src/morphobricks_normalization/histogram.py
@@ -52,6 +52,7 @@
 
 from .bin_centres import float_bin_centres, float_bin_position
 from .dtypes import check_supported, default_bin_count, is_integer, validate_bin_count
+from .integer_representatives import beyond_float_exact, integer_representatives, wide_integer_error
 from .integrity import check_histogram, check_tail, checked_blocks
 from .ranges import CENTRAL_PERCENTILES, HistogramSettings, RangeCheck
 from .sums import weighted_sum
@@ -172,11 +173,7 @@ def representatives(self) -> np.ndarray:
         if self.is_empty:
             return np.full(self.bin_count, np.nan)
         if is_integer(self.dtype):
-            first, last = self._integer_bin_limits
-            return np.array(
-                [self.binning_min + (a + b) / 2 for a, b in zip(first, last, strict=True)],
-                dtype=np.float64,
-            )
+            return integer_representatives(self)
         lo, hi = float(self.binning_min), float(self.binning_max)
         return float_bin_centres(lo, hi, self.bin_count)
 
@@ -188,10 +185,12 @@ def max_representative_error(self) -> float:
         Float planes: half the bin width. Pixels of a tail are represented by
         the tail's mean, so they add the larger distance from that mean to the
         tail's ends: the plane's extremum and the binning limit. A histogram
-        with tails is therefore never exact.
+        with tails is therefore never exact. Beyond ``2**51``: see ``integer_representatives``.
         """
         if self.is_empty:
             return 0.0
+        if is_integer(self.dtype) and beyond_float_exact(self):
+            return wide_integer_error(self)
         if is_integer(self.dtype):
             first, last = self._integer_bin_limits
             widest = max(b - a + 1 for a, b in zip(first, last, strict=True))
diff --git a/src/morphobricks_normalization/integer_representatives.py b/src/morphobricks_normalization/integer_representatives.py
new file mode 100644
index 0000000..0c14f09
--- /dev/null
+++ b/src/morphobricks_normalization/integer_representatives.py
@@ -0,0 +1,65 @@
+"""Representative values of integer bins, and their distance from the bins' pixels (#28).
+
+A bin of the integers ``first .. last`` stands for its pixels by its midpoint
+``(first + last) / 2``, divided once as Python integers, so it is the float64 nearest the
+midpoint. Below ``2**51`` every integer, every half-integer midpoint and every distance
+between them is exact in float64, and the bound is ``(widest bin - 1) / 2``. From ``2**51``
+on (64-bit pixel types) float64 cannot hold every value: ``2**53 + 1`` is stored as
+``2**53``. There the bound is the exact largest distance from a bin's or tail's integers to
+the float64 standing in for them, rounded up, so a histogram is exact only when every
+representative equals the pixel values it stands for.
+"""
+
+from __future__ import annotations
+
+import math
+from fractions import Fraction
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+if TYPE_CHECKING:  # the histogram module imports this one
+    from .histogram import PlaneHistogram
+
+FLOAT_EXACT_LIMIT = 2**51
+"""Below this magnitude integer midpoints and their distances are exact float64 arithmetic."""
+
+
+def integer_representatives(histogram: PlaneHistogram) -> np.ndarray:
+    """Return the float64 nearest the midpoint of every bin of a non-empty integer plane."""
+    return np.array(
+        [(first + last) / 2 for first, last in histogram.integer_bin_limits()], dtype=np.float64
+    )
+
+
+def beyond_float_exact(histogram: PlaneHistogram) -> bool:
+    """Whether a non-empty integer plane has a pixel of magnitude ``2**51`` or more."""
+    return max(-int(histogram.pixel_min), int(histogram.pixel_max)) >= FLOAT_EXACT_LIMIT
+
+
+def wide_integer_error(histogram: PlaneHistogram) -> float:
+    """Return the exact bound of an integer plane beyond ``2**51``, rounded up to a float64.
+
+    Every bin that holds an integer, and every non-empty tail (its mean against the plane's
+    extremum and the binning limit), counts with the larger exact distance from the value
+    standing in for it to its two ends.
+    """
+    ends = [
+        (rep, first, last)
+        for rep, (first, last) in zip(
+            histogram.representatives.tolist(), histogram.integer_bin_limits(), strict=True
+        )
+        if first <= last  # a bin that holds no integer has no pixel to stand in for
+    ]
+    for tail, first, last in (
+        (histogram.underflow, histogram.pixel_min, histogram.binning_min),
+        (histogram.overflow, histogram.binning_max, histogram.pixel_max),
+    ):
+        if not tail.is_empty:
+            ends.append((tail.mean, first, last))
+    worst = max(
+        max(abs(Fraction(rep) - int(first)), abs(Fraction(rep) - int(last)))
+        for rep, first, last in ends
+    )
+    nearest = float(worst)
+    return nearest if Fraction(nearest) >= worst else math.nextafter(nearest, math.inf)
diff --git a/src/morphobricks_normalization/pooling.py b/src/morphobricks_normalization/pooling.py
index 10abd90..eebd8a2 100644
--- a/src/morphobricks_normalization/pooling.py
+++ b/src/morphobricks_normalization/pooling.py
@@ -167,9 +167,7 @@ def std(self) -> float:
             return _std_with_tails(
                 self.bin_values, self.bin_counts, self.tails, self.pixel_count, self.mean
             )
-        deviations = self.values - self.mean
-        squares = weighted_sum(self.counts, deviations * deviations)
-        return float(np.sqrt(squares / self.pixel_count))
+        return _std_with_tails(self.values, self.counts, [], self.pixel_count, self.mean)
 
     def percentile(self, q: float) -> float:
         """Inverse-CDF percentile, ``q`` in percent (0-100)."""
@@ -307,9 +305,7 @@ def std(self) -> float:
             return _std_with_tails(
                 self.bin_values, self.bin_weights, self.tails, self.total_weight, self.mean
             )
-        deviations = self.values - self.mean
-        squares = weighted_sum(self.weights, deviations * deviations)
-        return float(np.sqrt(squares / self.total_weight))
+        return _std_with_tails(self.values, self.weights, [], self.total_weight, self.mean)
 
     def percentile(self, q: float) -> float:
         """Weighted inverse-CDF percentile, ``q`` in percent (0-100)."""
@@ -393,14 +389,67 @@ def _std_with_tails(
     total: float,
     mean: float,
 ) -> float:
-    deviations = values - mean
+    """Return the population standard deviation of the bins (and tails) around ``mean`` (#33).
+
+    The squares are first summed as they are, which keeps every finite result bit for
+    bit. When that overflows although the mean and every value are finite (``[-1e200,
+    1e200]`` has standard deviation ``1e200`` but squares of ``1e400``), the sum is redone
+    with every value divided by a power of two near the largest, which is exact, so the finite
+    spread is found. A tail whose stored sum of squares is not finite is refused by name.
+    """
+    squares = [s2 for _, _, s2 in tails]
+    if not all(math.isfinite(s2) for s2 in squares):
+        raise ValueError(_non_finite_squares_message(squares))
+    try:
+        with np.errstate(over="ignore", invalid="ignore"):  # an overflow is redone below
+            std = _std_scaled(values, weights, tails, total, mean, 1.0)
+    except OverflowError:  # float ** 2 raises where NumPy gives inf
+        std = math.inf
+    if math.isfinite(std) or not (math.isfinite(mean) and np.isfinite(values).all()):
+        return std
+    largest = max(float(np.max(np.abs(values), initial=0.0)), abs(mean))
+    largest = max([largest] + [abs(s / count) for count, s, _ in tails])
+    # Values divided by 2**(exponent - 1) are below 2 in magnitude, and the scale is finite.
+    return _std_scaled(
+        values, weights, tails, total, mean, math.ldexp(1.0, math.frexp(largest)[1] - 1)
+    )
+
+
+def _std_scaled(
+    values: np.ndarray,
+    weights: np.ndarray,
+    tails: list[tuple[float, float, float]],
+    total: float,
+    mean: float,
+    scale: float,
+) -> float:
+    """Return the standard deviation computed in units of ``scale``, a power of two (1: as is)."""
+    mean = mean / scale
+    deviations = values / scale - mean
     squares = weighted_sum(weights, deviations * deviations)
     for count, s, s2 in tails:
         # Each tail adds its spread around its own mean and its mean's distance
         # from the overall mean; the spread never goes negative through rounding.
+        s, s2 = s / scale, s2 / scale / scale
         tail_mean = s / count
         squares += max(s2 - s * tail_mean, 0.0) + count * (tail_mean - mean) ** 2
-    return float(np.sqrt(squares / total))
+    return float(np.sqrt(squares / total)) * scale
+
+
+def _non_finite_squares_message(squares: list[float]) -> str:
+    """Refusal text for tails whose sum of squares overflowed float64 when it was computed."""
+    bad = ", ".join(
+        f"tail {number} of {len(squares)} has sum of squares {s2!r}"
+        for number, s2 in enumerate(squares, 1)
+        if not math.isfinite(s2)
+    )
+    return (
+        f"The standard deviation cannot be pooled: {bad}. The squares of pixels outside the "
+        "robust range overflowed float64 when the statistics file was written (pixel values "
+        "beyond about 1.3e154), so their spread is not known, although it may be finite. Use "
+        "a method that does not need the standard deviation (percentile or robust z-score), "
+        "or the min-max range."
+    )
 
 
 def _channel_histogram(record: PlaneRecord, channel: int | None) -> PlaneHistogram:
diff --git a/src/morphobricks_normalization/scopes.py b/src/morphobricks_normalization/scopes.py
index c8a201f..4253ef3 100644
--- a/src/morphobricks_normalization/scopes.py
+++ b/src/morphobricks_normalization/scopes.py
@@ -104,7 +104,8 @@ def scope_groups(
         The rows of each group, by group key in sorted order.
     issues : list of Issue
         ``repeat_missing_conditions`` when a repeat group has no rows from
-        some condition of the statistics; the group is still used.
+        some condition of the statistics; the group is still used. With
+        ``per_time_point`` each time point's group is checked (``<label>::t=<t>``).
 
     Raises
     ------
@@ -128,7 +129,7 @@ def scope_groups(
         groups.setdefault(group_key(record, settings, repeat_groups), []).append(record)
     issues = []
     if settings.scope == "biological_repeat":
-        issues.extend(_incomplete_repeat_groups(records, repeat_groups))
+        issues.extend(_incomplete_repeat_groups(groups))
     return dict(sorted(groups.items())), issues
 
 
@@ -179,17 +180,17 @@ def _missing_pairs_error(pairs: set[tuple[str, str]]) -> ValueError:
     )
 
 
-def _incomplete_repeat_groups(
-    records: list[PlaneRecord], repeat_groups: RepeatGroups | None
-) -> list[Issue]:
-    conditions = {record.condition for record in records}
-    present: dict[str, set[str]] = {}
-    for record in records:
-        present.setdefault(_repeat_group(record, repeat_groups), set()).add(record.condition)
+def _incomplete_repeat_groups(groups: dict[tuple, list[PlaneRecord]]) -> list[Issue]:
+    """Warn about every fitted group without rows from some condition of the statistics.
+
+    The groups are the ones that are fitted, so with ``per_time_point`` each time point of
+    a repeat group is checked on its own and named ``<label>::t=<t>`` (#41).
+    """
+    conditions = {record.condition for members in groups.values() for record in members}
     missing = {
-        label: sorted(conditions - found)
-        for label, found in sorted(present.items())
-        if conditions - found
+        group_name(key): sorted(conditions - {record.condition for record in members})
+        for key, members in sorted(groups.items())
+        if conditions - {record.condition for record in members}
     }
     if not missing:
         return []
diff --git a/tests/test_histogram.py b/tests/test_histogram.py
index 90be97e..27ef002 100644
--- a/tests/test_histogram.py
+++ b/tests/test_histogram.py
@@ -1,12 +1,19 @@
 import bisect
+from fractions import Fraction
 
 import numpy as np
 import pytest
 
 from morphobricks_normalization import (
+    NormalizationSettings,
+    PlaneRecord,
+    PooledHistogram,
     UnsupportedDtypeError,
     compute_block_histogram,
     compute_plane_histogram,
+    fit,
+    read_plane_histograms,
+    write_plane_histograms,
 )
 from morphobricks_normalization import histogram as histogram_module
 
@@ -344,3 +351,84 @@ def blocks():
     floats = plane.astype(np.float32)
     compute_block_histogram(lambda: (calls.append(1), iter([floats]))[1], floats.dtype, 4)
     assert len(calls) == 2  # floats are always read twice
+
+
+# 64-bit integers that float64 cannot hold (#28) ---------------------------------------
+
+UNREPRESENTABLE = [
+    ("uint64", 2**53 - 1),
+    ("uint64", 2**53),
+    ("uint64", 2**53 + 1),
+    ("uint64", 2**64 - 1),
+    ("int64", 2**63 - 1),
+    ("int64", -(2**63)),
+    ("int64", -(2**63) + 1),
+    ("int64", -(2**53) - 1),
+]
+
+
+@pytest.mark.parametrize(("dtype", "value"), UNREPRESENTABLE)
+def test_singleton_64_bit_value_is_exact_only_if_float64_holds_it(dtype, value):
+    hist = compute_plane_histogram(np.array([[value]], dtype=dtype), 256)
+    (occupied,) = np.flatnonzero(hist.counts)
+    representative = hist.representatives[occupied]
+    error = abs(Fraction(float(representative)) - value)  # exact, independent of the package
+    assert hist.counts[occupied] == 1 and hist.pixel_min == hist.pixel_max == value
+    assert hist.is_exact == (error == 0)
+    assert hist.max_representative_error == float(error)  # 0 or 1: representable
+
+
+def test_adjacent_64_bit_integers_sharing_a_float_are_not_exact():
+    plane = np.array([[2**53, 2**53 + 1, 2**53 + 1]], dtype=np.uint64)
+    hist = compute_plane_histogram(plane, 256)
+    np.testing.assert_array_equal(np.flatnonzero(hist.counts), [0, 1])  # bins stay exact
+    assert hist.representatives[0] == hist.representatives[1] == 2.0**53
+    assert not hist.is_exact and hist.max_representative_error == 1.0
+    pooled = PooledHistogram([hist])
+    assert pooled.counts.tolist() == [3] and pooled.max_representative_error == 1.0
+    record = PlaneRecord("a.tif", "c", "r", "t", [hist], position_id="a")
+    settings = NormalizationSettings(method="minmax", scope="image")
+    with pytest.raises(ValueError, match="zero"):  # 2**53 + 1 - 2**53 rounds to 0 in float64
+        fit([record], settings, group="a.tif")
+    plane = np.array([[2**53 - 2, 2**53 + 1]], dtype=np.uint64)
+    parameters = fit(
+        [
+            PlaneRecord(
+                "a.tif", "c", "r", "t", [compute_plane_histogram(plane, 256)], position_id="a"
+            )
+        ],
+        settings,
+        group="a.tif",
+    )
+    assert parameters.exact is False
+
+
+@pytest.mark.parametrize("dtype", ["int64", "uint64"])
+@pytest.mark.parametrize("bins", [2, 7, 256])
+def test_wide_64_bit_bound_covers_every_pixel_exactly(dtype, bins):
+    info = np.iinfo(dtype)
+    plane = np.random.default_rng(bins).integers(
+        info.min, info.max, (24, 31), dtype=dtype, endpoint=True
+    )
+    plane[0, :2] = (info.min, info.max)
+    hist = compute_plane_histogram(plane, bins)
+    block = compute_block_histogram(lambda: iter((plane[:5], plane[5:])), plane.dtype, bins)
+    assert block.max_representative_error == hist.max_representative_error
+    starts = [first for first, _ in hist.integer_bin_limits()]
+    representatives = [Fraction(float(r)) for r in hist.representatives]
+    worst = max(
+        abs(representatives[bisect.bisect_right(starts, int(v)) - 1] - int(v))
+        for v in plane.ravel().tolist()
+    )
+    assert Fraction(hist.max_representative_error) >= worst and not hist.is_exact
+
+
+def test_64_bit_bound_survives_the_statistics_file(tmp_path):
+    plane = np.array([[2**53 + 1, 2**53 + 3, 2**53 + 4]], dtype=np.uint64)
+    hist = compute_plane_histogram(plane, 256)
+    record = PlaneRecord("a.tif", "c", "r", "t", [hist], position_id="a")
+    (back,) = read_plane_histograms(write_plane_histograms([record], tmp_path / "p.csv"))
+    read = back.histograms[0]
+    assert read.pixel_min == 2**53 + 1 and read.pixel_max == 2**53 + 4
+    assert read.max_representative_error == hist.max_representative_error == 1.0
+    assert not read.is_exact
diff --git a/tests/test_pooling.py b/tests/test_pooling.py
index 154481d..4333506 100644
--- a/tests/test_pooling.py
+++ b/tests/test_pooling.py
@@ -2,13 +2,23 @@
 import pytest
 
 from morphobricks_normalization import (
+    NormalizationSettings,
     NoValidPixelsError,
+    PlaneRecord,
     PooledHistogram,
     compute_plane_histogram,
+    fit,
     group_records,
     pooled_statistics,
+    read_plane_histograms,
+    write_plane_histograms,
+    z_score,
+)
+from morphobricks_normalization.pooling import (
+    DEFAULT_PERCENTILES,
+    WeightedHistogram,
+    _linear_percentile,
 )
-from morphobricks_normalization.pooling import DEFAULT_PERCENTILES, _linear_percentile
 from morphobricks_normalization.ranges import HistogramSettings
 from simulation import direct_stats, make_records, simulate_experiment
 
@@ -142,3 +152,88 @@ def test_percentiles_from_pixels_near_float64_limits_are_finite_and_exact(q, lim
     assert np.isfinite(result)
     assert abs(Fraction(result) - expected) <= abs(expected) * Fraction(1, 10**14)
     assert np.isfinite(pooled.percentile(q))
+
+
+# Standard deviations of large finite values (#33) -------------------------------------
+
+
+def exact_std(values, weights, tails=()):
+    """Population standard deviation of weighted values and (count, sum, squares) tails.
+
+    Every float is taken at its exact value (Fraction) and the square root in 40 digits, so
+    nothing here can overflow or round like the package's float64 arithmetic.
+    """
+    from decimal import Decimal, localcontext
+    from fractions import Fraction
+
+    values = [Fraction(v) for v in np.asarray(values, dtype=float).tolist()]
+    weights = [Fraction(w) for w in np.asarray(weights, dtype=float).tolist()]
+    tails = [tuple(Fraction(x) for x in tail) for tail in tails]
+    total = sum(weights) + sum(c for c, _, _ in tails)
+    mean = (
+        sum(w * v for w, v in zip(weights, values, strict=True)) + sum(s for _, s, _ in tails)
+    ) / total
+    squares = sum(w * (v - mean) ** 2 for w, v in zip(weights, values, strict=True))
+    squares += sum(s2 - s * s / c + c * (s / c - mean) ** 2 for c, s, s2 in tails)
+    variance = squares / total
+    with localcontext() as context:
+        context.prec = 40
+        return float((Decimal(variance.numerator) / Decimal(variance.denominator)).sqrt())
+
+
+LARGE_PLANES = {
+    "opposite 1e200": [-1e200, 1e200],
+    "1e200 +- 1e185": [1e200 - 1e185, 1e200, 1e200 + 1e185],
+    "opposite 1e308": [-1e308, 0.0, 1e308],
+    "constant 1e200": [1e200] * 5,
+}
+
+
+@pytest.mark.parametrize("name", LARGE_PLANES)
+def test_std_of_large_finite_values_is_finite_and_exact(name):
+    hist = compute_plane_histogram(np.array(LARGE_PLANES[name]), 64)
+    pooled = PooledHistogram([hist])
+    expected = exact_std(pooled.values, pooled.counts)
+    assert np.isfinite(expected)
+    assert pooled.std == pytest.approx(expected, rel=1e-14, abs=0)
+    if name != "constant 1e200":
+        parameters = z_score(pooled)
+        assert parameters.scale == pooled.std and np.isfinite(parameters.offset)
+
+
+def test_weighted_std_of_large_finite_values_is_finite_and_exact():
+    parts = [
+        PooledHistogram([compute_plane_histogram(np.array(LARGE_PLANES[name]), 64)])
+        for name in ("opposite 1e200", "1e200 +- 1e185")
+    ]
+    weighted = WeightedHistogram(parts, [0.5, 2.0])
+    expected = exact_std(weighted.values, weighted.weights)
+    assert weighted.std == pytest.approx(expected, rel=1e-14, abs=0)
+
+
+def test_std_with_large_tails_is_found_by_scaling():
+    plane = np.random.default_rng(7).normal(0, 1, 10000) * 1e153
+    settings = HistogramSettings(histogram_range="robust", robust_quantiles=(0.05, 99.95))
+    hist = compute_plane_histogram(plane, 256, settings=settings)
+    assert hist.has_tails and np.isfinite(hist.overflow.sum_of_squares)
+    pooled = PooledHistogram([hist])
+    expected = exact_std(pooled.bin_values, pooled.bin_counts, pooled.tails)
+    assert pooled.std == pytest.approx(expected, rel=1e-12, abs=0)
+    assert pooled.std == pytest.approx(np.std(plane / 1e153) * 1e153, rel=1e-4)
+
+
+def test_tail_whose_squares_overflowed_is_a_named_refusal(tmp_path):
+    plane = np.random.default_rng(1).normal(0, 1, 10000) * 1e150
+    plane[:5] = (1e200, 2e200, -1e200, 3e200, -2e200)
+    with np.errstate(over="ignore"):  # the tail's squares overflow as the file is written
+        hist = compute_plane_histogram(plane, 256, settings=ROBUST)
+    assert hist.overflow.sum_of_squares == np.inf
+    record = PlaneRecord("a.tif", "c", "r", "t", [hist], position_id="a")
+    (back,) = read_plane_histograms(write_plane_histograms([record], tmp_path / "p.csv"))
+    pooled = PooledHistogram([back.histograms[0]])
+    assert np.isfinite(pooled.mean) and np.isfinite(pooled.percentile(50))
+    message = "pooled: tail 1 of 2 has sum of squares inf, tail 2 of 2 has sum of squares inf"
+    with pytest.raises(ValueError, match=message):
+        pooled.summary()  # the standard deviation; mean and percentiles stay available
+    with pytest.raises(ValueError, match="percentile or robust z-score"):
+        fit([back], NormalizationSettings(method="zscore", scope="image"), group="a.tif")
diff --git a/tests/test_scopes.py b/tests/test_scopes.py
index 23e49db..fa76728 100644
--- a/tests/test_scopes.py
+++ b/tests/test_scopes.py
@@ -184,3 +184,42 @@ def test_existing_label_groups_keep_their_meaning():
         ("control", "r1"),
         ("treated", "r1"),
     ]
+
+
+REPEAT_PER_T = NormalizationSettings(scope="biological_repeat", per_time_point=True)
+
+
+def missing_items(records, settings, pairing=None):
+    _, issues = scope_groups(records, settings, pairing)
+    return [item for issue in issues for item in issue.items]
+
+
+def test_per_time_point_repeat_groups_are_checked_at_each_time_point():
+    # Issue #41: each condition only at one time point of r1, so both fitted groups lack one.
+    records = [record("control", "r1", t=0), record("treated", "r1", t=1)]
+    groups, issues = scope_groups(records, REPEAT_PER_T)
+    assert list(groups) == [("r1", 0), ("r1", 1)]
+    (issue,) = issues
+    assert (issue.code, issue.items) == ("repeat_missing_conditions", ("r1::t=0", "r1::t=1"))
+    assert "r1::t=0: no images from treated" in issue.message
+    assert "r1::t=1: no images from control" in issue.message
+    # Pooled over time the same rows form one complete group: no warning, as before.
+    assert missing_items(records, REPEAT) == []
+
+
+def test_per_time_point_completeness_controls():
+    complete = [record(c, "r1", t=t) for c in ("control", "treated") for t in (0, 1)]
+    assert missing_items(complete, REPEAT_PER_T) == []
+    # A condition with no images at all at one time point, in every repeat.
+    staggered = [r for r in complete if (r.condition, r.t) != ("treated", 1)]
+    staggered += [record("control", "r2", t=0), record("treated", "r2", t=0)]
+    staggered.append(record("control", "r2", t=1))
+    assert missing_items(staggered, REPEAT_PER_T) == ["r1::t=1", "r2::t=1"]
+    assert missing_items(staggered, REPEAT) == []
+    # Renamed pairs are checked under the label of the pairing.
+    pairing = {("control", "r1"): "A", ("treated", "r1"): "A"}
+    pairing |= {("control", "r2"): "B", ("treated", "r2"): "B"}
+    assert missing_items(staggered, REPEAT_PER_T, pairing) == ["A::t=1", "B::t=1"]
+    # The image scope has no repeat check, with or without time points.
+    per_t_image = NormalizationSettings(scope="image", per_time_point=True)
+    assert missing_items(staggered, per_t_image) == []
diff --git a/tools/code_rules_baseline.json b/tools/code_rules_baseline.json
index 4848f84..41c1cdf 100644
--- a/tools/code_rules_baseline.json
+++ b/tools/code_rules_baseline.json
@@ -32,10 +32,9 @@
     "src/morphobricks_normalization/histogram.py _integer_bins: # Ruff findings too-many-arguments": 1,
     "src/morphobricks_normalization/histogram.py _integer_bins: # parameters (max #)": 9,
     "src/morphobricks_normalization/histogram.py _integer_bins: private function of # lines without a docstring": 25,
-    "src/morphobricks_normalization/histogram.py: # lines (max #)": 862,
+    "src/morphobricks_normalization/histogram.py: # lines (max #)": 861,
     "src/morphobricks_normalization/pooling.py __init__: # Ruff findings missing-return-type-special-method": 1,
     "src/morphobricks_normalization/records.py _check_channels: private function of # lines without a docstring": 19,
-    "src/morphobricks_normalization/scopes.py _incomplete_repeat_groups: private function of # lines without a docstring": 26,
     "src/morphobricks_normalization/settings.py _percentile_pair: # Ruff findings any-type": 1,
     "src/morphobricks_normalization/settings.py _positive_number: # Ruff findings any-type": 1,
     "src/morphobricks_normalization/summaries.py _saturated_fraction: # Ruff findings missing-type-function-argument": 1,