diff --git a/docs/inference/model_selection.rst b/docs/inference/model_selection.rst
index 6e66dd0..6fefa30 100644
--- a/docs/inference/model_selection.rst
+++ b/docs/inference/model_selection.rst
@@ -8,9 +8,11 @@ Model Selection
 Classical model-selection criteria for choosing the order (state count) of a
 fitted generator when a fully-Bayesian evidence is unavailable or undesirable.
 The module scores any fitted :class:`~sofic.generators.base.HiddenMarkovModel`
-using the natural-log likelihood from
+using the log-likelihood (in bits) from
 :func:`sofic.generators.hmm_inference.log_likelihood` and a free-parameter count
-read off the transition graph:
+read off the transition graph. Log-likelihoods, log scores, and the MDL code
+length are in bits; AIC, AICc, BIC, and WAIC are computed from the natural
+log-likelihood so they keep their standard deviance scale:
 
 * **AIC** :cite:`Akaike1974` and the small-sample-corrected **AICc**
   :cite:`HurvichTsai1989`,
diff --git a/docs/quickstart.rst.txt b/docs/quickstart.rst.txt
index 9b3cf8f..61eed6d 100644
--- a/docs/quickstart.rst.txt
+++ b/docs/quickstart.rst.txt
@@ -77,7 +77,7 @@ topological entropy (the log of the golden ratio) :cite:`LindMarcus1995`:
 
    @doctest float
    In [4]: tmc.topological_entropy()
-   Out[4]: 0.4812118250596035
+   Out[4]: 0.6942419136306174
 
 The Parry measure turns the shift into the maximum-entropy ``MealyHMM`` on it
 (``tmc.parry_measure()``), the measure that attains this topological entropy.
diff --git a/docs/shifts/covers.rst b/docs/shifts/covers.rst
index 5ce62bd..a935dc8 100644
--- a/docs/shifts/covers.rst
+++ b/docs/shifts/covers.rst
@@ -5,18 +5,36 @@
 Covers
 ******
 
-Fischer and Krieger covers convert a Sofic shift into unifilar presentations
-:cite:`Fischer1975,Krieger1984,LindMarcus1995`.
-
-* :class:`LeftFischerCover`, :class:`RightFischerCover` — implemented.
-* :class:`LeftKriegerCover`, :class:`RightKriegerCover` — construction raises
-  :exc:`NotImplementedError`.
+Fischer and Krieger covers are canonical resolving presentations of a sofic
+shift :cite:`Fischer1975,Krieger1984,LindMarcus1995`. Following Lind and Marcus,
+*right* means right-resolving (deterministic reading forward):
+
+* :class:`RightFischerCover` -- the unique minimal right-resolving presentation
+  of an irreducible sofic shift. Reducible shifts raise
+  :exc:`~sofic.exceptions.SoficValidationError`, since their minimal
+  right-resolving presentation need not be unique.
+* :class:`RightKriegerCover` -- the future cover, with one vertex per follower
+  set of a left-infinite ray. Defined for every sofic shift; for an irreducible
+  shift the Fischer cover is its unique terminal component.
+* :class:`LeftFischerCover`, :class:`LeftKriegerCover` -- the left-resolving
+  mirror images, built from the reversed shift.
+
+Both constructions are exact: the Fischer cover is the terminal component of the
+follower-merged subset construction, and the Krieger cover's vertices are the
+images of path relations lying on cycles of the finite relation monoid.
 
 .. ipython::
 
-   In [1]: from sofic.examples import golden_mean_shift_parry; parry = golden_mean_shift_parry()
+   In [1]: from sofic.graph import ATTR_SYMBOL
+
+   In [2]: from sofic.shifts import RightKriegerCover, SoficShift
+
+   In [3]: even = SoficShift(symbol_alphabet=frozenset("01"))
 
-   In [2]: parry.validate()
+   In [4]: for s, t, a in [("A", "A", "0"), ("A", "B", "1"), ("B", "A", "1")]:
+      ...:     even.graph.add_transition(s, t, **{ATTR_SYMBOL: a})
+
+   In [5]: len(list(RightKriegerCover.from_sofic(even).states()))
 
 API
 ===
@@ -26,7 +44,11 @@ API
 .. autoclass:: RightFischerCover
    :members: from_sofic
 .. autoclass:: LeftKriegerCover
+   :members: from_sofic
 .. autoclass:: RightKriegerCover
+   :members: from_sofic
 
 .. autofunction:: sofic.shifts.cover_construction.left_fischer_from_sofic
 .. autofunction:: sofic.shifts.cover_construction.right_fischer_from_sofic
+.. autofunction:: sofic.shifts.cover_construction.left_krieger_from_sofic
+.. autofunction:: sofic.shifts.cover_construction.right_krieger_from_sofic
diff --git a/docs/shifts/topological_anatomy.rst b/docs/shifts/topological_anatomy.rst
index 4d6359e..c309518 100644
--- a/docs/shifts/topological_anatomy.rst
+++ b/docs/shifts/topological_anatomy.rst
@@ -67,8 +67,8 @@ entropy across both parts, while the sofic even shift is purely bound
    Out[10]: 0.5527864045001022
 
 The parts add up to :math:`h_\mathrm{top}`, which equals
-:meth:`~sofic.shifts.sofic.SoficShift.topological_entropy` divided by
-:math:`\ln 2` on a right-resolving presentation.
+:meth:`~sofic.shifts.sofic.SoficShift.topological_entropy` (in bits) on a
+right-resolving presentation.
 
 API
 ===
diff --git a/docs/shifts/topological_markov_chain.rst b/docs/shifts/topological_markov_chain.rst
index bba680d..e50ffe3 100644
--- a/docs/shifts/topological_markov_chain.rst
+++ b/docs/shifts/topological_markov_chain.rst
@@ -15,7 +15,7 @@ adjacency matrix. Its Parry measure is the maximum-entropy stochastic generator
 
    @doctest float
    In [2]: tmc.topological_entropy()
-   Out[2]: 0.48121182505960347
+   Out[2]: 0.6942419136306174
 
 API
 ===
diff --git a/docs/shifts/wheeler.rst b/docs/shifts/wheeler.rst
index 37f3905..0110c37 100644
--- a/docs/shifts/wheeler.rst
+++ b/docs/shifts/wheeler.rst
@@ -43,10 +43,9 @@ excluded outright, and :func:`wheeler_cover` raises
 
 The cover is a :class:`~sofic.shifts.covers.WheelerCover`, a
 :class:`~sofic.shifts.sofic.SoficShift` subclass that sits beside the Fischer
-and Krieger covers. It does not fill the Krieger stubs in
-:mod:`sofic.shifts.cover_construction`: a Krieger cover's states are *all*
-sets of pasts closed under the follower relation, whereas a Wheeler cover
-carries only those that happen to be recency intervals.
+and Krieger covers. It is not a Krieger cover: a Krieger cover's states are
+*all* follower sets of left-infinite pasts, whereas a Wheeler cover carries
+only those that happen to be recency intervals.
 
 Once a shift has a Wheeler cover, :func:`wheeler_index_of_shift` gives
 ``O(|w| log |A|)`` factor-language membership in place of scanning
diff --git a/sofic/automata/algorithms.py b/sofic/automata/algorithms.py
index bb58d58..b7b41a6 100644
--- a/sofic/automata/algorithms.py
+++ b/sofic/automata/algorithms.py
@@ -244,9 +244,15 @@ def minimize_hopcroft(dfa: DFA, *, alphabet: frozenset[Any] | None = None) -> DF
 def equivalent(
     aut1: LabeledAutomaton,
     aut2: LabeledAutomaton,
-    alphabet: frozenset[Any],
+    alphabet: frozenset[Any] | None = None,
 ) -> bool:
-    """Return whether two automata recognize the same language over ``alphabet``."""
+    """Return whether two automata recognize the same language.
+
+    The comparison runs over ``alphabet`` together with every symbol either
+    automaton declares or uses, so a too-small ``alphabet`` cannot hide a
+    difference on the omitted symbols.
+    """
+    alphabet = frozenset(alphabet or ()) | _transition_alphabet(aut1) | _transition_alphabet(aut2)
     d1 = minimize(_to_nfa(aut1), alphabet=alphabet, algorithm="hopcroft")
     d2 = minimize(_to_nfa(aut2), alphabet=alphabet, algorithm="hopcroft")
     return _isomorphic_minimal_dfa(d1, d2, alphabet)
@@ -281,6 +287,15 @@ def _effective_alphabet(aut: LabeledAutomaton) -> frozenset[Any]:
     return frozenset(symbols)
 
 
+def _transition_alphabet(aut: LabeledAutomaton) -> frozenset[Any]:
+    symbols = {symbol for symbol in aut.input_alphabet if symbol is not EPSILON}
+    for transition in aut.transitions():
+        symbol = transition.data.get(ATTR_SYMBOL)
+        if symbol is not None and symbol is not EPSILON:
+            symbols.add(symbol)
+    return frozenset(symbols)
+
+
 def _forward_reachable(aut: LabeledAutomaton) -> set[Hashable]:
     if not aut.initial_states:
         return set()
diff --git a/sofic/automata/buchi_simulation.py b/sofic/automata/buchi_simulation.py
index 26da207..f6f9973 100644
--- a/sofic/automata/buchi_simulation.py
+++ b/sofic/automata/buchi_simulation.py
@@ -9,7 +9,12 @@
 
 
 def accepts_lasso_buchi(ba: BuchiAutomaton, prefix: Sequence[Any], loop: Sequence[Any]) -> bool:
-    """Accept if repeating ``loop`` after ``prefix`` visits an accepting state infinitely often."""
+    """Accept if repeating ``loop`` after ``prefix`` visits an accepting state infinitely often.
+
+    ``loop`` must be non-empty: ``prefix loop^omega`` is an infinite word only then.
+    """
+    if not loop:
+        raise ValueError("an ultimately periodic omega-word needs a non-empty loop")
     post = ba._run_nfa(prefix)
     if not post:
         return False
diff --git a/sofic/examples/epsilon_machines.py b/sofic/examples/epsilon_machines.py
index b30e777..6673ea6 100644
--- a/sofic/examples/epsilon_machines.py
+++ b/sofic/examples/epsilon_machines.py
@@ -27,6 +27,7 @@
 
 import numpy as np
 
+from sofic.exceptions import SoficError
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.mealy import MealyHMM
 from sofic.graph import ATTR_EMISSION, ATTR_FUTURE_SYMBOL, ATTR_PROB, TransitionGraph
@@ -1004,7 +1005,7 @@ def tent_map_misiurewicz_bidirectional_fig8(a: Any | None = None):
     else:
         try:
             reverse = EpsilonMachine.from_hmm(reverse_raw)
-        except Exception:
+        except SoficError:
             from sofic.generators.epsilon_machine import _row_normalized_presentation
 
             reverse = _row_normalized_presentation(reverse_raw)
diff --git a/sofic/examples/processes.py b/sofic/examples/processes.py
index 928ba77..b0b11e3 100644
--- a/sofic/examples/processes.py
+++ b/sofic/examples/processes.py
@@ -64,7 +64,7 @@ def _stationary_initial(
     transition = transition / row_sums[:, None]
     try:
         pi = stationary_distribution_from_transition(transition)
-    except Exception:
+    except (ValueError, np.linalg.LinAlgError):
         return _uniform_initial(states)
     return {state: float(pi[i]) for i, state in enumerate(states)}
 
diff --git a/sofic/generators/alternative_complexity.py b/sofic/generators/alternative_complexity.py
index c65d118..68e5ad0 100644
--- a/sofic/generators/alternative_complexity.py
+++ b/sofic/generators/alternative_complexity.py
@@ -100,5 +100,5 @@ def _stationary_vector(transition: np.ndarray) -> np.ndarray | None:
 
     try:
         return stationary_distribution_from_transition(transition)
-    except Exception:
+    except (ValueError, np.linalg.LinAlgError):
         return None
diff --git a/sofic/generators/block_convergence.py b/sofic/generators/block_convergence.py
index 7800981..a3b592b 100644
--- a/sofic/generators/block_convergence.py
+++ b/sofic/generators/block_convergence.py
@@ -8,6 +8,7 @@
 
 import numpy as np
 
+from sofic.exceptions import SoficError
 from sofic.generators._word_measures import (
     _block_caekl,
     _block_coinformation,
@@ -203,7 +204,7 @@ def _anatomy_curves(
 def _exact_anatomy_scalars(machine: EpsilonMachine) -> dict[str, float] | None:
     try:
         bidir = machine.to_bidirectional()
-    except Exception:
+    except (SoficError, ValueError, NotImplementedError, np.linalg.LinAlgError):
         return None
     h_mu = float(bidir.entropy_rate())
     rho_mu = float(bidir.predicted_information())
diff --git a/sofic/generators/block_entropy.py b/sofic/generators/block_entropy.py
index 64a8ca0..318fee3 100644
--- a/sofic/generators/block_entropy.py
+++ b/sofic/generators/block_entropy.py
@@ -9,6 +9,8 @@
 
 import numpy as np
 
+from sofic.exceptions import SoficError
+
 if TYPE_CHECKING:
     from sofic.generators.epsilon_machine import EpsilonMachine
 
@@ -327,7 +329,10 @@ def block_entropy_estimates(
 
     entropy_asymptote = excess_entropy + h_mu_l
     crypticity_estimate = state_block_entropy - block_state_entropy
-    crypticity = float(crypticity_estimate[-1]) if crypticity_estimate.size else 0.0
+    if use_exact:
+        crypticity = statistical_complexity - excess_entropy
+    else:
+        crypticity = float(crypticity_estimate[-1]) if crypticity_estimate.size else 0.0
 
     cm = _cm_extension_curves(
         lengths,
@@ -453,7 +458,7 @@ def _excess_entropy(
 ) -> float:
     try:
         return float(machine.to_bidirectional().excess_entropy())
-    except Exception:
+    except (SoficError, ValueError, NotImplementedError, np.linalg.LinAlgError):
         return _block_entropy_excess_entropy(machine, entropy_rate=entropy_rate, block_entropy=block_entropy)
 
 
@@ -498,7 +503,7 @@ def _estimated_excess_entropy(machine: EpsilonMachine, estimate: np.ndarray, *,
     if use_exact:
         try:
             return float(machine.excess_entropy())
-        except Exception:
+        except (SoficError, ValueError, NotImplementedError, np.linalg.LinAlgError):
             pass
     return float(estimate[-1]) if estimate.size else 0.0
 
diff --git a/sofic/generators/epsilon_inference.py b/sofic/generators/epsilon_inference.py
index 488be70..ad3ab3c 100644
--- a/sofic/generators/epsilon_inference.py
+++ b/sofic/generators/epsilon_inference.py
@@ -832,7 +832,8 @@ def suggest_lmax(
     select_markov_order = getattr(dit.inference, "select_markov_order", None)
     if select_markov_order is None:  # pragma: no cover - depends on the installed dit
         raise ImportError("suggest_lmax requires a dit release with dit.inference.select_markov_order")
-    seq = [repr(symbol) for symbol in sequence]
+    codes: dict[Any, str] = {}
+    seq = [codes.setdefault(symbol, str(len(codes))) for symbol in sequence]
     if max_order is None:
         k = max(2, len(set(seq)))
         max_order = max(1, min(10, int(np.log(max(len(seq), 1) / 5) / np.log(k)) - 1))
diff --git a/sofic/generators/hmm_inference.py b/sofic/generators/hmm_inference.py
index f83728a..48f575b 100644
--- a/sofic/generators/hmm_inference.py
+++ b/sofic/generators/hmm_inference.py
@@ -8,6 +8,7 @@
 
 from __future__ import annotations
 
+import warnings
 from collections import defaultdict
 from collections.abc import Hashable, Iterable, Sequence
 from typing import Any
@@ -131,7 +132,7 @@ def _stationary_emission_tensors(
             return limited, joint
     try:
         pi = stationary_distribution_from_transition(transition)
-    except Exception:
+    except (np.linalg.LinAlgError, ValueError):
         pi = pi_initial
     return pi, joint
 
@@ -143,7 +144,7 @@ def _forward_scaled(
 ) -> tuple[np.ndarray, np.ndarray]:
     """Return per-step-normalized forward messages and log scaling factors.
 
-    ``alpha_hat[t]`` sums to one; ``log P(obs) = log_scales.sum()``. A ``-inf``
+    ``alpha_hat[t]`` sums to one; ``log2 P(obs) = log_scales.sum()``. A ``-inf``
     entry in ``log_scales`` marks an impossible step. Normalizing each step avoids
     the underflow that makes the raw forward product vanish for long sequences.
     """
@@ -155,7 +156,7 @@ def _forward_scaled(
         log_scales[0] = -np.inf
         return alpha_hat, log_scales
     alpha_hat[0] = pi / total0
-    log_scales[0] = float(np.log(total0))
+    log_scales[0] = float(np.log2(total0))
     for t, symbol in enumerate(obs):
         matrix = joint.get(symbol)
         if matrix is None:
@@ -167,7 +168,7 @@ def _forward_scaled(
             log_scales[t + 1] = -np.inf
             continue
         alpha_hat[t + 1] = row / scale
-        log_scales[t + 1] = float(np.log(scale))
+        log_scales[t + 1] = float(np.log2(scale))
     return alpha_hat, log_scales
 
 
@@ -239,7 +240,7 @@ def _backward_scaled(joint: dict[Any, np.ndarray], obs: list[Any], n_states: int
 
 
 def log_likelihood(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> float:
-    """Natural-log likelihood ``log P(observations)``.
+    """Log-likelihood ``log2 P(observations)`` in bits.
 
     Uses the per-step-scaled forward recursion so the result stays finite for long
     sequences instead of underflowing to ``-inf``.
@@ -316,7 +317,7 @@ def _expected_edge_counts(
     - ``source_totals[i] = \sum_{t=0}^{n-1} P(X_t = i \mid Y)`` is the expected
       number of transitions out of state ``i`` (the Baum-Welch denominator);
     - ``gamma0`` is the smoothed marginal of the initial state ``X_0``;
-    - ``loglik`` is the natural-log likelihood of the sequence.
+    - ``loglik`` is the log-likelihood of the sequence in bits.
 
     Only edges present in ``joint`` (structural support) receive mass, so the
     statistics preserve the model topology.
@@ -391,7 +392,7 @@ def baum_welch(
     so the fit is returned as a plain :class:`~sofic.generators.mealy.MealyHMM`.
 
     Returns ``(fitted_model, loglik_trace)`` where ``loglik_trace`` is the
-    non-decreasing sequence of total natural-log likelihoods observed before each
+    non-decreasing sequence of total log-likelihoods (bits) observed before each
     parameter update.
 
     EM converges to a local maximum of the likelihood. With ``n_restarts > 1`` the
@@ -483,15 +484,25 @@ def _baum_welch_run(
         total_source = np.zeros(n_states, dtype=float)
         gamma0_sum = np.zeros(n_states, dtype=float)
         total_ll = 0.0
+        skipped = 0
         for obs in seqs:
             edge_counts, source_totals, gamma0, loglik = _expected_edge_counts(pi, joint, obs)
             if not np.isfinite(loglik):
+                skipped += 1
                 continue
             for key, value in edge_counts.items():
                 total_edge_counts[key] += value
             total_source += source_totals
             gamma0_sum += gamma0
             total_ll += loglik
+        if seqs and skipped == len(seqs):
+            raise ValueError("every observation sequence has zero probability under the model")
+        if skipped and _iteration == 0:
+            warnings.warn(
+                f"{skipped} of {len(seqs)} sequences have zero probability under the model and are ignored",
+                RuntimeWarning,
+                stacklevel=3,
+            )
         loglik_trace.append(total_ll)
         if prev_ll is not None and abs(total_ll - prev_ll) < tol:
             break
@@ -522,6 +533,10 @@ def score(hmm: HiddenMarkovModel, observations: Sequence[Any]) -> dict[tuple[Has
     ``\nabla \log L(\theta) = E[\nabla \log f(X, Y; \theta) \mid Y]`` (Cappe,
     Moulines & Ryden, 2005, Section 10.2.3), evaluated in the raw (unconstrained)
     joint-edge parameters. Keys are ``(source, symbol, target)`` state labels.
+
+    The score is the gradient of the *natural* log-likelihood
+    :math:`\ln P(Y) = \ln 2 \cdot` :func:`log_likelihood`, the convention under
+    which the observed information and standard errors are standard.
     """
     mealy = hmm.to_mealy()
     idx = mealy.reindex()
@@ -761,7 +776,10 @@ def sample(
     mealy = _as_mealy_hmm(hmm)
     idx = mealy.reindex()
     pi, joint = _emission_transition_tensors_from_mealy(mealy)
-    state = int(generator.choice(len(idx), p=pi / pi.sum()))
+    total = float(np.sum(pi))
+    if not total > 0.0:
+        raise ValueError("cannot sample: the initial state distribution has no mass")
+    state = int(generator.choice(len(idx), p=pi / total))
 
     observations: list[Any] = []
     states: list[Hashable] = []
diff --git a/sofic/generators/markov.py b/sofic/generators/markov.py
index e2d61dd..36d6882 100644
--- a/sofic/generators/markov.py
+++ b/sofic/generators/markov.py
@@ -82,10 +82,23 @@ def words_of_length(self, length: int) -> dict[tuple[Hashable, ...], float]:
         return markov_words_of_length(self, length)
 
     def sample_path(self, n: int, rng: np.random.Generator | None = None) -> list[Hashable]:
+        """Sample a state path of length ``n``.
+
+        Starts from :attr:`initial_distribution`, or from the stationary
+        distribution when no initial distribution is given.
+        """
         generator = rng if rng is not None else np.random.default_rng()
         idx = self.reindex()
-        pi = self.stationary_distribution()
-        state = int(generator.choice(len(idx), p=pi))
+        if self.initial_distribution:
+            pi = np.zeros(len(idx), dtype=float)
+            for state, mass in self.initial_distribution.items():
+                pi[idx.index(state)] = float(mass)
+        else:
+            pi = np.asarray(self.stationary_distribution(), dtype=float)
+        total = float(pi.sum())
+        if not total > 0.0:
+            raise ValueError("cannot sample: the initial state distribution has no mass")
+        state = int(generator.choice(len(idx), p=pi / total))
 
         path: list[Hashable] = []
         for _ in range(n):
diff --git a/sofic/generators/measures.py b/sofic/generators/measures.py
index 62bb5f0..d6deee6 100644
--- a/sofic/generators/measures.py
+++ b/sofic/generators/measures.py
@@ -239,7 +239,7 @@ def entropy_rate_markov(chain: MarkovChain) -> Any:
 
 
 def collision_entropy(quasi_model: QuasiStochasticModel) -> float:
-    """Second Renyi entropy rate from quasi transition matrices."""
+    """Second Renyi entropy rate (bits) from quasi transition matrices."""
     matrices = quasi_model.transition_matrices()
     pi = quasi_model.stationary_quasidistribution()
     total = 0.0
@@ -247,7 +247,7 @@ def collision_entropy(quasi_model: QuasiStochasticModel) -> float:
         total += float(pi @ (matrix @ matrix) @ np.ones(len(pi)))
     if total <= 0.0:
         return 0.0
-    return float(-np.log(total))
+    return float(-np.log2(total))
 
 
 def process_negativity(quasi_model: QuasiStochasticModel) -> float:
diff --git a/sofic/generators/words.py b/sofic/generators/words.py
index 171a999..6fcffec 100644
--- a/sofic/generators/words.py
+++ b/sofic/generators/words.py
@@ -164,11 +164,13 @@ def markov_words_of_length(chain: MarkovChain, length: int) -> dict[tuple[Hashab
     if length < 0:
         raise ValueError("length must be nonnegative")
     states = tuple(chain.states())
+    start = _markov_start(chain)
     if length == 0:
-        return {(): 1.0}
+        total = float(sum(start.values()))
+        return {(): total} if total > _TOL else {}
     distribution: dict[tuple[Hashable, ...], float] = {}
     for word in product(states, repeat=length):
-        probability = _markov_path_probability(chain, word)
+        probability = _markov_path_probability(chain, word, start)
         if probability > _TOL:
             distribution[word] = probability
     return distribution
@@ -211,10 +213,19 @@ def _start_vector(
     return vector.copy()
 
 
-def _markov_path_probability(chain: MarkovChain, path: tuple[Hashable, ...]) -> float:
+def _markov_start(chain: MarkovChain) -> dict[Hashable, float]:
+    """Initial law of ``chain``, or its stationary law when none is given."""
+    if chain.initial_distribution:
+        return {state: float(mass) for state, mass in chain.initial_distribution.items()}
+    idx = chain.reindex()
+    pi = chain.stationary_distribution()
+    return {idx.state(i): float(mass) for i, mass in enumerate(pi)}
+
+
+def _markov_path_probability(chain: MarkovChain, path: tuple[Hashable, ...], start: Mapping[Hashable, float]) -> float:
     if not path:
-        return 1.0
-    probability = float(chain.initial_distribution.get(path[0], 0.0))
+        return float(sum(start.values()))
+    probability = float(start.get(path[0], 0.0))
     for source, target in zip(path, path[1:], strict=False):
         edge_probability = 0.0
         for transition in chain.graph.out_transitions(source):
diff --git a/sofic/inference/bayesian/hdp_hmm.py b/sofic/inference/bayesian/hdp_hmm.py
index 300e722..64b3a6b 100644
--- a/sofic/inference/bayesian/hdp_hmm.py
+++ b/sofic/inference/bayesian/hdp_hmm.py
@@ -47,7 +47,7 @@ class HDPHMMPosterior:
     state_counts
         Number of occupied states in each retained draw (``len == len(samples)``).
     log_likelihoods
-        Data log-likelihood (natural log) of each retained draw.
+        Data log-likelihood (bits) of each retained draw.
     alphabet
         Sorted observation alphabet used by the sampler.
     """
@@ -110,7 +110,7 @@ def _ffbs(
     emit: np.ndarray,
     rng: np.random.Generator,
 ) -> tuple[np.ndarray, float]:
-    """Forward-filter backward-sample one sequence; return states and log-likelihood."""
+    """Forward-filter backward-sample one sequence; return states and log-likelihood (bits)."""
     n_states = trans.shape[0]
     length = obs_idx.shape[0]
     alpha = np.empty((length, n_states))
@@ -122,7 +122,7 @@ def _ffbs(
         weights = emit[:, obs_idx[0]].copy()
         scale = weights.sum()
     alpha[0] = weights / scale
-    loglik += np.log(scale)
+    loglik += np.log2(scale)
 
     for t in range(1, length):
         predicted = alpha[t - 1] @ trans
@@ -132,7 +132,7 @@ def _ffbs(
             weights = emit[:, obs_idx[t]].copy()
             scale = weights.sum()
         alpha[t] = weights / scale
-        loglik += np.log(scale)
+        loglik += np.log2(scale)
 
     states = np.empty(length, dtype=int)
     states[length - 1] = rng.choice(n_states, p=alpha[length - 1])
diff --git a/sofic/inference/model_selection.py b/sofic/inference/model_selection.py
index c108140..912c177 100644
--- a/sofic/inference/model_selection.py
+++ b/sofic/inference/model_selection.py
@@ -1,4 +1,4 @@
-"""Classical model-selection criteria for stochastic generators.
+r"""Classical model-selection criteria for stochastic generators.
 
 Point-estimate information criteria -- AIC :cite:`Akaike1974`, the
 small-sample-corrected AICc :cite:`HurvichTsai1989`, BIC :cite:`Schwarz1978`,
@@ -7,11 +7,15 @@
 criterion (WAIC) :cite:`Watanabe2010`. These complement the exact Bayesian
 evidences of :mod:`sofic.inference.bayesian`: they score any fitted
 :class:`~sofic.generators.base.HiddenMarkovModel` (ε-machine, Mealy HMM, Markov
-chain) using the natural-log likelihood from
+chain) using the log-likelihood (in bits) from
 :func:`sofic.generators.hmm_inference.log_likelihood` and a free-parameter count
 read off the transition graph, so they are likelihood-agnostic and apply
 directly to discrete-emission models.
 
+Log-likelihoods, log scores, and the MDL code length are reported in **bits**.
+AIC, AICc, BIC, and WAIC keep their standard deviance scale: they are computed
+from the natural log-likelihood :math:`\ln L = \ln 2 \cdot \log_2 L`.
+
 All information criteria follow the convention **lower is better**;
 cross-validated and WAIC log scores follow **higher is better** for the raw log
 score (WAIC itself is reported on the deviance scale, lower is better).
@@ -123,7 +127,7 @@ def _smoothed_log_likelihood(
     smoothing: float,
     alphabet_size: int,
 ) -> float:
-    """Natural-log likelihood with each one-step prediction mixed with the uniform law.
+    """Log-likelihood (bits) with each one-step prediction mixed with the uniform law.
 
     ``P'(x_t | x_{0:t}) = (1 - smoothing) P(x_t | x_{0:t}) + smoothing / alphabet_size``,
     computed by forward filtering. After a symbol the model forbids, the belief is
@@ -140,7 +144,7 @@ def _smoothed_log_likelihood(
         matrix = joint.get(symbol)
         unnormalized = belief @ matrix if matrix is not None else np.zeros_like(belief)
         predicted = float(unnormalized.sum())
-        total += float(np.log((1.0 - smoothing) * predicted + smoothing / alphabet_size))
+        total += float(np.log2((1.0 - smoothing) * predicted + smoothing / alphabet_size))
         if predicted > 0.0:
             belief = unnormalized / predicted
         else:
@@ -158,8 +162,9 @@ def score_model(
     """Score ``model`` on ``data`` with AIC, AICc, BIC, and MDL.
 
     ``data`` may be a single observation sequence or an iterable of sequences.
-    The scores use natural-log likelihoods; the number of observations is the
-    total symbol count. When the data has zero probability under the model the
+    ``log_likelihood`` and ``mdl`` are in bits; AIC, AICc, and BIC use the natural
+    log-likelihood so they keep their usual scale. The number of observations is
+    the total symbol count. When the data has zero probability under the model the
     likelihood is ``-inf`` and every criterion is ``+inf``.
     """
     sequences = _normalize_sequences(data)
@@ -171,12 +176,12 @@ def score_model(
         inf = float("inf")
         return ModelScores(float("-inf"), k, n, inf, inf, inf, inf)
 
-    aic = 2.0 * k - 2.0 * ll
+    ln_l = ll * np.log(2.0)
+    aic = 2.0 * k - 2.0 * ln_l
     denom = n - k - 1
     aicc = aic + (2.0 * k * (k + 1)) / denom if denom > 0 else float("inf")
-    log_n = np.log(n) if n > 0 else 0.0
-    bic = k * log_n - 2.0 * ll
-    mdl = 0.5 * k * log_n - ll
+    bic = k * (np.log(n) if n > 0 else 0.0) - 2.0 * ln_l
+    mdl = 0.5 * k * (np.log2(n) if n > 0 else 0.0) - ll
     return ModelScores(float(ll), k, n, float(aic), float(aicc), float(bic), float(mdl))
 
 
@@ -229,7 +234,7 @@ def cross_validated_log_likelihood(
     gap: int = 0,
     smoothing: float = 0.0,
 ) -> float:
-    """Return the total held-out natural-log likelihood under ``folds``-fold CV.
+    """Return the total held-out log-likelihood (bits) under ``folds``-fold CV.
 
     ``fit(train_sequences)`` must fit and return a model from a list of training
     sequences. When ``data`` is a collection of sequences the folds partition the
@@ -249,7 +254,7 @@ def cross_validated_log_likelihood(
     smoothing
         Mix each one-step held-out prediction with the uniform distribution over
         the observed alphabet, with this weight. Then a single transition that the
-        fitted model forbids costs ``log(smoothing / |A|)`` instead of making the
+        fitted model forbids costs ``log2(smoothing / |A|)`` instead of making the
         whole fold ``-inf``, so models can still be compared.
     """
     if gap < 0:
@@ -307,17 +312,20 @@ def waic(pointwise_log_likelihoods: np.ndarray) -> WAICResult:
     r"""Widely applicable information criterion from posterior samples.
 
     ``pointwise_log_likelihoods`` has shape ``(n_samples, n_points)`` with entry
-    ``[s, i] = log p(y_i | theta_s)`` for posterior draw ``theta_s``. Returns the
-    WAIC on the deviance scale (lower is better),
+    ``[s, i] = log2 p(y_i | theta_s)`` (bits) for posterior draw ``theta_s``, as
+    returned by :func:`posterior_pointwise_log_likelihoods`. Returns the WAIC on
+    the standard (natural-log) deviance scale (lower is better),
     ``WAIC = -2 (lppd - p_waic)`` with the log pointwise predictive density
     ``lppd = sum_i log mean_s p(y_i | theta_s)`` and effective parameter count
-    ``p_waic = sum_i Var_s log p(y_i | theta_s)`` :cite:`Watanabe2010`.
+    ``p_waic = sum_i Var_s log p(y_i | theta_s)`` :cite:`Watanabe2010`. ``waic``,
+    ``p_waic``, and ``standard_error`` use natural logs; ``lppd`` is in bits.
     """
     from scipy.special import logsumexp
 
     matrix = np.asarray(pointwise_log_likelihoods, dtype=float)
     if matrix.ndim != 2 or matrix.size == 0:
         raise ValueError("pointwise_log_likelihoods must be a non-empty (n_samples, n_points) array")
+    matrix = matrix * np.log(2.0)
     n_samples = matrix.shape[0]
     lppd_pointwise = logsumexp(matrix, axis=0) - np.log(n_samples)
     p_waic_pointwise = matrix.var(axis=0, ddof=1) if n_samples > 1 else np.zeros(matrix.shape[1])
@@ -327,7 +335,7 @@ def waic(pointwise_log_likelihoods: np.ndarray) -> WAICResult:
     standard_error = float(np.sqrt(n_points * np.var(-2.0 * elpd_pointwise, ddof=0))) if n_points > 1 else 0.0
     return WAICResult(
         waic=waic_value,
-        lppd=float(lppd_pointwise.sum()),
+        lppd=float(lppd_pointwise.sum() / np.log(2.0)),
         p_waic=float(p_waic_pointwise.sum()),
         standard_error=standard_error,
     )
diff --git a/sofic/properties.py b/sofic/properties.py
index ed67507..4edf4c0 100644
--- a/sofic/properties.py
+++ b/sofic/properties.py
@@ -214,7 +214,7 @@ def is_detailed_balance(model: StateMachine, *, rtol: float = 1e-8, atol: float
     """Return whether stationary labeled flows satisfy detailed balance."""
     try:
         pi = np.asarray(model.stationary_distribution(), dtype=float)
-    except Exception:
+    except (AttributeError, ValueError, np.linalg.LinAlgError):
         pi, _transition = _initial_vector_and_transition(model)
 
     matrices = _labeled_or_internal_matrices(model)
diff --git a/sofic/shifts/algorithms.py b/sofic/shifts/algorithms.py
index 2214cf6..514855a 100644
--- a/sofic/shifts/algorithms.py
+++ b/sofic/shifts/algorithms.py
@@ -8,7 +8,7 @@
 
 import numpy as np
 
-from sofic.graph import ATTR_SYMBOL
+from sofic.graph import ATTR_MULTIPLICITY, ATTR_SYMBOL
 from sofic.shifts.base import SymbolicModel
 
 
@@ -54,16 +54,19 @@ def adjacency_matrix(model: SymbolicModel) -> tuple[np.ndarray, tuple[Hashable,
     for transition in model.transitions():
         i = index[transition.source]
         j = index[transition.target]
-        matrix[i, j] += 1.0
+        matrix[i, j] += float(transition.data.get(ATTR_MULTIPLICITY, 1))
     return matrix, states
 
 
 def topological_entropy_from_matrix(matrix: np.ndarray) -> float:
+    """Return ``log2`` of the spectral radius of ``matrix`` (bits per symbol)."""
     if matrix.size == 0:
         return 0.0
     eigenvalues = np.linalg.eigvals(matrix)
     spectral_radius = float(np.max(np.abs(eigenvalues)))
-    return float(np.log(max(spectral_radius, 0.0)))
+    if spectral_radius <= 0.0:
+        return 0.0
+    return float(np.log2(spectral_radius))
 
 
 def _forward_reachable(model: SymbolicModel) -> set[Hashable]:
diff --git a/sofic/shifts/cover_construction.py b/sofic/shifts/cover_construction.py
index cb4daf1..e1a0b5f 100644
--- a/sofic/shifts/cover_construction.py
+++ b/sofic/shifts/cover_construction.py
@@ -1,76 +1,204 @@
-"""Fischer and Krieger cover constructions."""
+"""Fischer and Krieger cover constructions.
+
+Both covers are computed exactly from a presentation ``G`` with vertex set
+``Q``. Write ``S . w`` for the set of vertices reached from ``S`` along paths
+labeled ``w``, and ``F(S)`` for the follower set (future language) of ``S``.
+
+* The **right Fischer cover** of an irreducible sofic shift is its unique
+  minimal right-resolving presentation :cite:`Fischer1975`
+  :cite:`LindMarcus1995`. It is the unique terminal strongly connected
+  component of the subset construction from ``Q`` after merging subsets with
+  equal follower sets: an intrinsically synchronizing word ``m`` sends every
+  subset to the follower class ``F(m)``, so that class is reachable from all
+  others.
+* The **right Krieger cover** (future cover) has one vertex per follower set
+  ``F(x^-)`` of a left-infinite ray, with ``F(x^-) --a--> F(x^- a)``
+  :cite:`Krieger1984` :cite:`LindMarcus1995`. ``F(x^-) = F(T(x^-))`` where
+  ``T(x^-) = Q . s`` for every long enough suffix ``s`` of ``x^-``. Reading a
+  ray right to left composes path relations ``rho_{cs} = rho_c o rho_s`` in a
+  finite monoid, so the sets ``T(x^-)`` are exactly the images ``Q . rho`` of
+  relations ``rho`` that lie on a cycle reachable from the identity.
+
+The left covers are the mirror images: the right cover of the reversed shift,
+reversed back.
+"""
 
 from __future__ import annotations
 
-from collections import defaultdict
+from collections import defaultdict, deque
+from collections.abc import Hashable, Iterable
 from typing import Any
 
+import networkx as nx
+
+from sofic.exceptions import SoficValidationError
 from sofic.graph import ATTR_SYMBOL, TransitionGraph
 from sofic.shifts.covers import LeftFischerCover, LeftKriegerCover, RightFischerCover, RightKriegerCover
 from sofic.shifts.sofic import SoficShift
 from sofic.states import sequential_labels
 
+Subset = frozenset[Hashable]
+Relation = frozenset[tuple[Hashable, Hashable]]
 
-def _follower_language(shift: SoficShift, vertex: Any, max_len: int = 8) -> frozenset[tuple[Any, ...]]:
-    from collections import deque
 
-    seen: set[tuple[Any, ...]] = set()
-    queue: deque[tuple[Any, tuple[Any, ...]]] = deque([(vertex, ())])
+def _labeled_successors(shift: SoficShift) -> dict[Hashable, dict[Any, set[Hashable]]]:
+    successors: dict[Hashable, dict[Any, set[Hashable]]] = {state: defaultdict(set) for state in shift.states()}
+    for transition in shift.transitions():
+        symbol = transition.data.get(ATTR_SYMBOL)
+        if symbol is not None:
+            successors[transition.source][symbol].add(transition.target)
+    return successors
+
+
+def _subset_automaton(
+    successors: dict[Hashable, dict[Any, set[Hashable]]],
+) -> dict[Subset, dict[Any, Subset]]:
+    """Deterministic subset automaton reachable from the full vertex set."""
+    start: Subset = frozenset(successors)
+    delta: dict[Subset, dict[Any, Subset]] = {}
+    queue = deque([start])
     while queue:
-        state, prefix = queue.popleft()
-        if len(prefix) > max_len:
+        subset = queue.popleft()
+        if subset in delta:
             continue
-        if prefix:
-            seen.add(prefix)
-        for transition in shift.graph.out_transitions(state):
-            symbol = transition.data.get(ATTR_SYMBOL)
-            if symbol is None:
-                continue
-            queue.append((transition.target, prefix + (symbol,)))
-    return frozenset(seen)
-
-
-def _build_left_fischer(shift: SoficShift) -> LeftFischerCover:
-    followers: dict[Any, frozenset[tuple[Any, ...]]] = {
-        vertex: _follower_language(shift, vertex) for vertex in shift.states()
+        moves: dict[Any, set[Hashable]] = defaultdict(set)
+        for state in subset:
+            for symbol, targets in successors[state].items():
+                moves[symbol] |= targets
+        delta[subset] = {symbol: frozenset(targets) for symbol, targets in moves.items() if targets}
+        queue.extend(target for target in delta[subset].values() if target not in delta)
+    return delta
+
+
+def _follower_classes(delta: dict[Subset, dict[Any, Subset]]) -> dict[Subset, int]:
+    """Moore refinement: subsets with equal follower sets share a class."""
+    block = dict.fromkeys(delta, 0)
+    while True:
+        signatures = {
+            subset: (block[subset], tuple(sorted(((repr(a), block[t]) for a, t in moves.items()))))
+            for subset, moves in delta.items()
+        }
+        ids: dict[Any, int] = {}
+        refined = {subset: ids.setdefault(signature, len(ids)) for subset, signature in signatures.items()}
+        if len(ids) == len(set(block.values())):
+            return refined
+        block = refined
+
+
+def _labels(count: int) -> tuple[Hashable, ...]:
+    return sequential_labels(count) if count <= 26 else tuple(range(count))
+
+
+def _quotient_shift(
+    cls: type[SoficShift],
+    shift: SoficShift,
+    vertices: Iterable[Subset],
+    delta: dict[Subset, dict[Any, Subset]],
+    classes: dict[Subset, int],
+) -> SoficShift:
+    keep = set(vertices)
+    used = sorted({classes[subset] for subset in keep})
+    name = dict(zip(used, _labels(len(used)), strict=True))
+    graph = TransitionGraph()
+    for class_id in used:
+        graph.add_state(name[class_id])
+    edges = {
+        (classes[subset], symbol, classes[target])
+        for subset in keep
+        for symbol, target in delta[subset].items()
+        if target in keep
     }
-    classes: dict[frozenset[tuple[Any, ...]], list[Any]] = defaultdict(list)
-    for vertex, language in followers.items():
-        classes[language].append(vertex)
+    for source, symbol, target in sorted(edges, key=repr):
+        graph.add_transition(name[source], name[target], **{ATTR_SYMBOL: symbol})
+    return cls(graph=graph, symbol_alphabet=shift.symbol_alphabet)
 
-    graph = TransitionGraph()
-    class_for_vertex = {vertex: language for vertex, language in followers.items()}
-    state_ids = {language: sequential_labels(len(classes))[index] for index, language in enumerate(classes)}
-    for _language, state_id in state_ids.items():
-        graph.add_state(state_id)
 
-    for transition in shift.transitions():
-        source_lang = class_for_vertex[transition.source]
-        target_lang = class_for_vertex[transition.target]
-        symbol = transition.data.get(ATTR_SYMBOL)
-        graph.add_transition(state_ids[source_lang], state_ids[target_lang], **{ATTR_SYMBOL: symbol})
+def _mirror(cls: type[SoficShift], cover: SoficShift) -> SoficShift:
+    return cls(graph=cover.graph.reverse(), symbol_alphabet=cover.symbol_alphabet)
 
-    return LeftFischerCover(
-        graph=graph,
-        symbol_alphabet=shift.symbol_alphabet,
-    )
+
+def right_fischer_from_sofic(shift: SoficShift) -> RightFischerCover:
+    """Return the minimal right-resolving presentation of an irreducible sofic shift.
+
+    Raises :class:`~sofic.exceptions.SoficValidationError` when the shift is
+    reducible (more than one terminal component), in which case a minimal
+    right-resolving presentation need not be unique :cite:`LindMarcus1995`.
+    """
+    trimmed = shift.trim_transient()
+    delta = _subset_automaton(_labeled_successors(trimmed))
+    if not delta or not any(delta.values()):
+        return RightFischerCover(symbol_alphabet=shift.symbol_alphabet)
+    classes = _follower_classes(delta)
+
+    quotient = nx.DiGraph()
+    quotient.add_nodes_from(set(classes.values()))
+    for subset, moves in delta.items():
+        for target in moves.values():
+            quotient.add_edge(classes[subset], classes[target])
+    condensation = nx.condensation(quotient)
+    terminal = [node for node in condensation.nodes if condensation.out_degree(node) == 0]
+    if len(terminal) != 1:
+        raise SoficValidationError("the Fischer cover is defined for irreducible sofic shifts; this shift is reducible")
+    members = set(condensation.nodes[terminal[0]]["members"])
+    if not any(classes[s] in members and delta[s] for s in delta):
+        return RightFischerCover(symbol_alphabet=shift.symbol_alphabet)
+    vertices = [subset for subset in delta if classes[subset] in members]
+    return _quotient_shift(RightFischerCover, shift, vertices, delta, classes)
 
 
 def left_fischer_from_sofic(shift: SoficShift) -> LeftFischerCover:
-    return _build_left_fischer(shift.trim_transient())
+    """Return the minimal left-resolving presentation (mirror of the right Fischer cover)."""
+    return _mirror(LeftFischerCover, right_fischer_from_sofic(shift.reverse()))
 
 
-def right_fischer_from_sofic(shift: SoficShift) -> RightFischerCover:
-    left = left_fischer_from_sofic(shift.reverse())
-    return RightFischerCover(
-        graph=left.graph.copy(),
-        symbol_alphabet=left.symbol_alphabet,
-    )
+def _ray_terminal_sets(successors: dict[Hashable, dict[Any, set[Hashable]]]) -> set[Subset]:
+    """Return ``{T(x^-)}``: images of path relations lying on reachable cycles."""
+    symbols = {symbol for moves in successors.values() for symbol in moves}
+    letter: dict[Any, Relation] = {
+        symbol: frozenset((p, q) for p, moves in successors.items() for q in moves.get(symbol, ()))
+        for symbol in symbols
+    }
+    identity: Relation = frozenset((q, q) for q in successors)
 
+    def prepend(symbol: Any, relation: Relation) -> Relation:
+        after: dict[Hashable, set[Hashable]] = defaultdict(set)
+        for q, r in relation:
+            after[q].add(r)
+        return frozenset((p, r) for p, q in letter[symbol] for r in after.get(q, ()))
 
-def left_krieger_from_sofic(shift: SoficShift) -> LeftKriegerCover:
-    raise NotImplementedError("Left Krieger cover construction is not yet implemented; use left_fischer_from_sofic")
+    graph = nx.DiGraph()
+    graph.add_node(identity)
+    queue = deque([identity])
+    while queue:
+        relation = queue.popleft()
+        for symbol in symbols:
+            extended = prepend(symbol, relation)
+            if not extended:
+                continue
+            if extended not in graph:
+                queue.append(extended)
+            graph.add_edge(relation, extended)
+
+    recurrent: set[Relation] = set()
+    for component in nx.strongly_connected_components(graph):
+        node = next(iter(component))
+        if len(component) > 1 or graph.has_edge(node, node):
+            recurrent |= component
+    return {frozenset(r for _q, r in relation) for relation in recurrent}
 
 
 def right_krieger_from_sofic(shift: SoficShift) -> RightKriegerCover:
-    raise NotImplementedError("Right Krieger cover construction is not yet implemented; use right_fischer_from_sofic")
+    """Return the right Krieger (future) cover of ``shift``."""
+    trimmed = shift.trim_transient()
+    successors = _labeled_successors(trimmed)
+    if not successors:
+        return RightKriegerCover(symbol_alphabet=shift.symbol_alphabet)
+    delta = _subset_automaton(successors)
+    classes = _follower_classes(delta)
+    vertices = [subset for subset in _ray_terminal_sets(successors) if subset in delta]
+    return _quotient_shift(RightKriegerCover, shift, vertices, delta, classes)
+
+
+def left_krieger_from_sofic(shift: SoficShift) -> LeftKriegerCover:
+    """Return the left Krieger (past) cover: the mirror of the right Krieger cover."""
+    return _mirror(LeftKriegerCover, right_krieger_from_sofic(shift.reverse()))
diff --git a/sofic/shifts/sft_construction.py b/sofic/shifts/sft_construction.py
index 05fe2a9..559d433 100644
--- a/sofic/shifts/sft_construction.py
+++ b/sofic/shifts/sft_construction.py
@@ -4,6 +4,7 @@
 
 from typing import Any
 
+from sofic.exceptions import SoficValidationError
 from sofic.graph import ATTR_SYMBOL, TransitionGraph
 from sofic.shifts.sft import ShiftOfFiniteType
 
@@ -14,7 +15,12 @@ def from_forbidden_words(
     *,
     max_states: int = 256,
 ) -> ShiftOfFiniteType:
-    """Build an SFT presentation via a follower automaton on allowed prefixes."""
+    """Build an SFT presentation via a follower automaton on allowed prefixes.
+
+    Raises :class:`~sofic.exceptions.SoficValidationError` when the presentation
+    would need more than ``max_states`` states rather than returning a truncated
+    automaton.
+    """
     alphabet = tuple(symbol_alphabet)
     forbidden_set = set(forbidden)
     max_len = max((len(word) for word in forbidden_set), default=0)
@@ -30,16 +36,16 @@ def is_allowed(prefix: tuple[Any, ...]) -> bool:
 
     while queue:
         prefix = queue.pop(0)
-        if len(seen) >= max_states:
-            break
         for symbol in alphabet:
             extended = prefix + (symbol,)
             if not is_allowed(extended):
                 continue
-            trimmed = extended
-            if max_len > 0:
-                trimmed = extended[-max_len:]
+            trimmed = extended[-max_len:] if max_len > 0 else ()
             if trimmed not in seen:
+                if len(seen) >= max_states:
+                    raise SoficValidationError(
+                        f"SFT presentation needs more than max_states={max_states} states; raise max_states"
+                    )
                 seen.add(trimmed)
                 graph.add_state(trimmed)
                 queue.append(trimmed)
diff --git a/sofic/shifts/sliding_block_code.py b/sofic/shifts/sliding_block_code.py
index 3f249b1..577e91c 100644
--- a/sofic/shifts/sliding_block_code.py
+++ b/sofic/shifts/sliding_block_code.py
@@ -73,22 +73,43 @@ def apply_word(self, word: Sequence[Any]) -> tuple[Any, ...]:
     def apply(self, shift: Any) -> SoficShift:
         """Return the image subshift ``Phi(shift)`` as a sofic presentation.
 
-        Uses the higher-block construction: vertices are allowed
-        ``(window - 1)``-blocks of ``shift`` and each allowed ``window``-block
-        contributes an edge labeled by its image symbol.
+        Vertices pair a vertex ``q`` of ``shift``'s presentation with the last
+        ``window - 1`` symbols read along a path into ``q``; each edge of
+        ``shift`` out of ``q`` reading ``a`` emits ``Phi(context + a)``. Tracking
+        the presentation vertex keeps every constraint of ``shift``, not just
+        those visible in ``window``-blocks. Raises :class:`ValueError` when a
+        ``window``-block of ``shift`` is missing from ``block_map``.
         """
+        outgoing: dict[Any, list[tuple[Any, Any]]] = {state: [] for state in shift.states()}
+        for transition in shift.transitions():
+            symbol = transition.data.get(ATTR_SYMBOL)
+            if symbol is not None:
+                outgoing[transition.source].append((symbol, transition.target))
+
+        context_length = self.window - 1
+        frontier = {(state, ()) for state in outgoing}
+        for _ in range(context_length):
+            frontier = {
+                (target, (*context, symbol)) for state, context in frontier for symbol, target in outgoing[state]
+            }
+
         image = SoficShift(symbol_alphabet=frozenset(self.output_alphabet))
-        blocks = list(shift.factor_language(self.window))
-        vertices = {block[:-1] for block in blocks} | {block[1:] for block in blocks}
-        for vertex in vertices:
-            image.graph.add_state(vertex)
+        seen = set(frontier)
+        queue = list(frontier)
         used_outputs: set[Any] = set()
-        for block in blocks:
-            output = self.block_map.get(block)
-            if output is None:
-                continue
-            image.add_transition(block[:-1], block[1:], output)
-            used_outputs.add(output)
+        while queue:
+            state, context = queue.pop()
+            image.graph.add_state((state, context))
+            for symbol, target in outgoing[state]:
+                block = (*context, symbol)
+                if block not in self.block_map:
+                    raise ValueError(f"block {block!r} of the shift is not in block_map")
+                successor = (target, block[1:])
+                image.add_transition((state, context), successor, self.block_map[block])
+                used_outputs.add(self.block_map[block])
+                if successor not in seen:
+                    seen.add(successor)
+                    queue.append(successor)
         image.symbol_alphabet = frozenset(used_outputs)
         return image.trim_transient()
 
diff --git a/sofic/shifts/topological_anatomy.py b/sofic/shifts/topological_anatomy.py
index 4615b62..3909ec2 100644
--- a/sofic/shifts/topological_anatomy.py
+++ b/sofic/shifts/topological_anatomy.py
@@ -87,10 +87,9 @@ def _right_resolving(shift: SoficShift) -> SoficShift:
     If ``shift`` is already unifilar it is returned unchanged. Otherwise the right
     Fischer cover (:meth:`~sofic.shifts.covers.RightFischerCover.from_sofic`) is
     built and its duplicate labeled edges merged (:func:`_dedup_symbol_edges`).
-    The cover construction uses a bounded follower language, so it is not
-    guaranteed to determinize every presentation; if the result is still not
-    unifilar a :class:`~sofic.exceptions.UnifilarityError` is raised asking for a
-    right-resolving input.
+    The cover is exact for irreducible shifts; a reducible non-unifilar
+    presentation raises :class:`~sofic.exceptions.SoficValidationError` from the
+    cover construction.
     """
     if shift.is_unifilar():
         return shift
diff --git a/tests/test_block_entropy.py b/tests/test_block_entropy.py
index 4849f05..b449687 100644
--- a/tests/test_block_entropy.py
+++ b/tests/test_block_entropy.py
@@ -10,6 +10,7 @@
 from hypothesis import given, settings
 
 from sofic.examples import fair_coin, golden_mean
+from sofic.exceptions import SoficValidationError
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.generators.topological_epsilon_enumeration import idfa_string_to_epsilon_machine
 from sofic.testing.strategies import epsilon_machines
@@ -144,7 +145,7 @@ def test_golden_mean_block_entropy_estimates_match_finite_order_values():
 
 
 def test_block_entropy_estimates_fallback_when_exact_excess_entropy_fails():
-    with patch.object(EpsilonMachine, "excess_entropy", side_effect=RuntimeError("no bidirectional")):
+    with patch.object(EpsilonMachine, "excess_entropy", side_effect=SoficValidationError("no bidirectional")):
         estimates = golden_mean(0.5).block_entropy_estimates(3, use_exact=True)
 
     assert np.isfinite(estimates.excess_entropy)
@@ -163,7 +164,7 @@ def test_block_entropy_diagram_uses_exact_excess_entropy_when_available():
 
 
 def test_block_entropy_diagram_fallback_when_exact_excess_entropy_fails():
-    with patch.object(EpsilonMachine, "to_bidirectional", side_effect=RuntimeError("no bidirectional")):
+    with patch.object(EpsilonMachine, "to_bidirectional", side_effect=SoficValidationError("no bidirectional")):
         diagram = golden_mean(0.5).block_entropy_diagram(1)
 
     assert diagram.excess_entropy == pytest.approx(0.25162916738782304, abs=1e-12)
diff --git a/tests/test_buchi.py b/tests/test_buchi.py
index 16041d2..0c80703 100644
--- a/tests/test_buchi.py
+++ b/tests/test_buchi.py
@@ -53,3 +53,8 @@ def test_accepts_omega_non_periodic_raises():
     ba = _accepting_loop_ba()
     with pytest.raises(NotImplementedError):
         ba.accepts_omega(("a", "b", "a"))
+
+
+def test_empty_loop_is_not_an_omega_word():
+    with pytest.raises(ValueError, match="non-empty loop"):
+        _accepting_loop_ba().accepts_lasso(("a",), ())
diff --git a/tests/test_correctness_regressions.py b/tests/test_correctness_regressions.py
new file mode 100644
index 0000000..4b3e040
--- /dev/null
+++ b/tests/test_correctness_regressions.py
@@ -0,0 +1,108 @@
+"""Regression tests for correctness fixes found in the package review."""
+
+import numpy as np
+import pytest
+
+from sofic.automata.algorithms import equivalent
+from sofic.automata.dfa import DFA
+from sofic.examples import even_process, golden_mean
+from sofic.exceptions import SoficValidationError
+from sofic.generators.markov import MarkovChain
+from sofic.generators.mealy import MealyHMM
+from sofic.graph import ATTR_MULTIPLICITY, ATTR_SYMBOL
+from sofic.shifts.sft import ShiftOfFiniteType
+from sofic.shifts.tmc import TopologicalMarkovChain
+
+
+def test_tmc_multiplicity_counts_in_entropy_and_parry_measure():
+    tmc = TopologicalMarkovChain(symbol_alphabet=frozenset({"a"}))
+    tmc.graph.add_state("s")
+    tmc.graph.add_transition("s", "s", **{ATTR_SYMBOL: "a", ATTR_MULTIPLICITY: 2})
+    assert tmc.topological_entropy() == pytest.approx(1.0)
+
+    matrix_tmc = TopologicalMarkovChain.from_adjacency(np.array([[2, 1], [1, 0]]), symbol_alphabet=frozenset("ab"))
+    expected = np.log2(max(abs(np.linalg.eigvals(np.array([[2, 1], [1, 0]])))))
+    assert matrix_tmc.topological_entropy() == pytest.approx(expected)
+    parry = matrix_tmc.parry_measure()
+    parry.validate()
+    weights = {}
+    for t in parry.transitions():
+        weights[(t.source, t.target)] = weights.get((t.source, t.target), 0.0) + t.data["prob"]
+    for state in parry.states():
+        assert sum(w for (s, _), w in weights.items() if s == state) == pytest.approx(1.0)
+
+
+def test_golden_mean_topological_entropy_is_log2_golden_ratio():
+    tmc = TopologicalMarkovChain.from_adjacency(np.array([[1, 1], [1, 0]]), symbol_alphabet=frozenset("01"))
+    assert tmc.topological_entropy() == pytest.approx(np.log2((1 + np.sqrt(5)) / 2))
+
+
+def test_sft_construction_raises_instead_of_truncating():
+    forbidden = {tuple("0" * 9)}
+    with pytest.raises(SoficValidationError, match="max_states"):
+        ShiftOfFiniteType.from_forbidden_words(forbidden, frozenset("01"), max_states=8)
+
+
+def test_markov_words_of_length_zero_respects_initial_distribution():
+    chain = MarkovChain(initial_distribution={"a": 1.0})
+    chain.graph.add_state("a")
+    chain.graph.add_state("b")
+    chain.add_transition("a", "b", 1.0)
+    chain.add_transition("b", "a", 1.0)
+    assert chain.words_of_length(0) == {(): pytest.approx(1.0)}
+    assert chain.words_of_length(2) == {("a", "b"): pytest.approx(1.0)}
+
+    stationary = MarkovChain()
+    stationary.graph.add_state("a")
+    stationary.graph.add_state("b")
+    stationary.add_transition("a", "b", 1.0)
+    stationary.add_transition("b", "a", 1.0)
+    words = stationary.words_of_length(1)
+    assert words == {("a",): pytest.approx(0.5), ("b",): pytest.approx(0.5)}
+
+
+def test_markov_sample_path_starts_from_initial_distribution():
+    chain = MarkovChain(initial_distribution={"b": 1.0})
+    chain.graph.add_state("a")
+    chain.graph.add_state("b")
+    chain.add_transition("a", "b", 1.0)
+    chain.add_transition("b", "a", 1.0)
+    for seed in range(5):
+        assert chain.sample_path(3, rng=np.random.default_rng(seed)) == ["b", "a", "b"]
+
+
+def test_sample_without_initial_mass_raises():
+    hmm = MealyHMM(initial_distribution={"A": 0.0}, observation_alphabet=frozenset({0}))
+    hmm.graph.add_state("A")
+    hmm.add_transition("A", "A", 0, 1.0)
+    with pytest.raises(ValueError, match="no mass"):
+        hmm.sample(3, rng=np.random.default_rng(0))
+
+
+def test_block_entropy_estimates_exact_crypticity_is_cmu_minus_excess_entropy():
+    machine = even_process()
+    estimates = machine.block_entropy_estimates(3, use_exact=True)
+    assert estimates.crypticity == pytest.approx(estimates.statistical_complexity - estimates.excess_entropy)
+    assert estimates.crypticity == pytest.approx(machine.crypticity(), abs=1e-9)
+
+
+def test_log_likelihood_is_in_bits():
+    machine = golden_mean(0.5)
+    observations = [0, 0, 0, 0]
+    assert 2.0 ** machine.log_likelihood(observations) == pytest.approx(machine.word_probability(observations))
+
+
+def _single_symbol_dfa(symbols: set[str]) -> DFA:
+    dfa = DFA(input_alphabet=frozenset(symbols), initial_states=frozenset({0}), accepting_states=frozenset({1}))
+    dfa.graph.add_state(0)
+    dfa.graph.add_state(1)
+    for symbol in symbols:
+        dfa.add_transition(0, 1, symbol)
+    return dfa
+
+
+def test_equivalent_does_not_hide_differences_outside_the_given_alphabet():
+    only_a = _single_symbol_dfa({"a"})
+    a_or_b = _single_symbol_dfa({"a", "b"})
+    assert not equivalent(only_a, a_or_b, frozenset({"a"}))
+    assert equivalent(only_a, _single_symbol_dfa({"a"}))
diff --git a/tests/test_covers.py b/tests/test_covers.py
index 9172871..b1129dc 100644
--- a/tests/test_covers.py
+++ b/tests/test_covers.py
@@ -1,30 +1,104 @@
 """Tests for Fischer and Krieger covers."""
 
+import networkx as nx
 import pytest
 
+from sofic.exceptions import SoficValidationError
 from sofic.graph import ATTR_SYMBOL
 from sofic.shifts.covers import LeftFischerCover, LeftKriegerCover, RightFischerCover, RightKriegerCover
 from sofic.shifts.sofic import SoficShift
 
 
-def _golden_mean() -> SoficShift:
-    shift = SoficShift(symbol_alphabet=frozenset({"0", "1"}))
-    shift.graph.add_state("A")
-    shift.graph.add_state("B")
-    shift.graph.add_transition("A", "B", **{ATTR_SYMBOL: "1"})
-    shift.graph.add_transition("B", "A", **{ATTR_SYMBOL: "0"})
-    shift.graph.add_transition("B", "B", **{ATTR_SYMBOL: "1"})
+def _shift(edges, alphabet=("0", "1")) -> SoficShift:
+    shift = SoficShift(symbol_alphabet=frozenset(alphabet))
+    for source, target, symbol in edges:
+        shift.graph.add_state(source)
+        shift.graph.add_state(target)
+        shift.graph.add_transition(source, target, **{ATTR_SYMBOL: symbol})
     return shift
 
 
-@pytest.mark.parametrize("cls", [LeftFischerCover, RightFischerCover])
-def test_fischer_cover_from_sofic(cls):
-    cover = cls.from_sofic(_golden_mean())
+def _golden_mean() -> SoficShift:
+    return _shift([("A", "B", "1"), ("B", "A", "0"), ("B", "B", "0"), ("A", "A", "0")])
+
+
+def _even_shift() -> SoficShift:
+    """Runs of 1s between 0s have even length; Krieger cover has three vertices."""
+    return _shift([("A", "A", "0"), ("A", "B", "1"), ("B", "A", "1")])
+
+
+def _nondeterministic_even_shift() -> SoficShift:
+    # Two copies of the even-shift presentation glued nondeterministically.
+    return _shift(
+        [
+            ("A", "A", "0"),
+            ("A", "B", "1"),
+            ("B", "A", "1"),
+            ("A", "C", "0"),
+            ("C", "D", "1"),
+            ("D", "C", "1"),
+            ("C", "A", "0"),
+        ]
+    )
+
+
+def _language(shift: SoficShift, max_length: int = 8) -> set[tuple]:
+    return {word for n in range(max_length + 1) for word in shift.factor_language(n)}
+
+
+def _is_left_resolving(shift: SoficShift) -> bool:
+    return SoficShift(graph=shift.graph.reverse(), symbol_alphabet=shift.symbol_alphabet).is_unifilar()
+
+
+@pytest.mark.parametrize("builder", [_golden_mean, _even_shift, _nondeterministic_even_shift])
+def test_right_fischer_cover_is_minimal_right_resolving_and_presents_shift(builder):
+    shift = builder()
+    cover = RightFischerCover.from_sofic(shift)
     cover.validate()
-    assert len(list(cover.states())) >= 1
+    assert cover.is_unifilar()
+    assert nx.is_strongly_connected(cover.graph.nx)
+    assert _language(cover) == _language(shift)
+    assert len(list(cover.states())) == 2
+
+
+@pytest.mark.parametrize("builder", [_golden_mean, _even_shift, _nondeterministic_even_shift])
+def test_left_fischer_cover_is_left_resolving_and_presents_shift(builder):
+    shift = builder()
+    cover = LeftFischerCover.from_sofic(shift)
+    assert _is_left_resolving(cover)
+    assert _language(cover) == _language(shift)
+
+
+def test_even_shift_krieger_cover_has_three_vertices_and_contains_fischer_cover():
+    shift = _even_shift()
+    krieger = RightKriegerCover.from_sofic(shift)
+    fischer = RightFischerCover.from_sofic(shift)
+    assert krieger.is_unifilar()
+    assert len(list(krieger.states())) == 3
+    assert _language(krieger) == _language(shift)
+    condensation = nx.condensation(krieger.graph.nx)
+    terminal = [n for n in condensation.nodes if condensation.out_degree(n) == 0]
+    assert len(terminal) == 1
+    assert len(condensation.nodes[terminal[0]]["members"]) == len(list(fischer.states()))
+
+
+def test_golden_mean_krieger_cover_equals_fischer_cover_size():
+    shift = _golden_mean()
+    assert len(list(RightKriegerCover.from_sofic(shift).states())) == 2
+    assert len(list(LeftKriegerCover.from_sofic(shift).states())) == 2
+
+
+@pytest.mark.parametrize("builder", [_golden_mean, _even_shift, _nondeterministic_even_shift])
+def test_left_krieger_cover_is_left_resolving_and_presents_shift(builder):
+    shift = builder()
+    cover = LeftKriegerCover.from_sofic(shift)
+    assert _is_left_resolving(cover)
+    assert _language(cover) == _language(shift)
 
 
-@pytest.mark.parametrize("cls", [LeftKriegerCover, RightKriegerCover])
-def test_krieger_cover_not_implemented(cls):
-    with pytest.raises(NotImplementedError, match="Krieger"):
-        cls.from_sofic(_golden_mean())
+def test_fischer_cover_rejects_reducible_shift():
+    reducible = _shift([("A", "A", "0"), ("B", "B", "1")])
+    with pytest.raises(SoficValidationError, match="irreducible"):
+        RightFischerCover.from_sofic(reducible)
+    krieger = RightKriegerCover.from_sofic(reducible)
+    assert _language(krieger) == _language(reducible)
diff --git a/tests/test_examples.py b/tests/test_examples.py
index 60b2f2f..30b5232 100644
--- a/tests/test_examples.py
+++ b/tests/test_examples.py
@@ -81,7 +81,7 @@ def test_golden_mean_shift_parry_entropy():
         np.array([[1, 1], [1, 0]], dtype=float),
         symbol_alphabet=frozenset({0, 1}),
     )
-    assert parry.entropy_rate() == pytest.approx(tmc.topological_entropy() / np.log(2), rel=0.05)
+    assert parry.entropy_rate() == pytest.approx(tmc.topological_entropy(), rel=0.05)
 
 
 def test_butterfly_statistical_complexity():
diff --git a/tests/test_hmm_inference.py b/tests/test_hmm_inference.py
index db6668c..390a3a3 100644
--- a/tests/test_hmm_inference.py
+++ b/tests/test_hmm_inference.py
@@ -32,7 +32,7 @@ def test_forward_coin_initial_and_likelihood():
     alpha = forward(coin, observations)
     assert alpha.shape == (4, 1)
     assert alpha[0].sum() == pytest.approx(1.0, abs=1e-9)
-    assert alpha[-1].sum() == pytest.approx(np.exp(log_likelihood(coin, observations)), abs=1e-9)
+    assert alpha[-1].sum() == pytest.approx(2.0 ** log_likelihood(coin, observations), abs=1e-9)
 
 
 def test_backward_coin():
@@ -92,7 +92,7 @@ def test_log_likelihood_long_sequence_stays_finite():
     observations = ["0", "1"] * 1500
     ll = log_likelihood(coin, observations)
     assert np.isfinite(ll)
-    assert ll == pytest.approx(-3000 * np.log(2), rel=1e-9)
+    assert ll == pytest.approx(-3000.0, rel=1e-9)
 
 
 def test_forward_scaled_rows_are_normalized():
@@ -204,7 +204,7 @@ def loglik_entry(symbol: int, i: int, j: int, value: float) -> float:
         perturbed = {sym: matrix.copy() for sym, matrix in joint.items()}
         perturbed[symbol][i, j] = value
         _alpha, log_scales = _forward_scaled(pi, perturbed, list(obs))
-        return float(log_scales.sum())
+        return float(log_scales.sum()) * np.log(2)
 
     analytic = score(gm, obs)
     h = 1e-6
@@ -226,7 +226,7 @@ def loglik_theta(theta: float) -> float:
         perturbed[0][a, a] = theta
         perturbed[1][a, b] = 1.0 - theta
         _alpha, log_scales = _forward_scaled(pi, perturbed, list(obs))
-        return float(log_scales.sum())
+        return float(log_scales.sum()) * np.log(2)
 
     assert free_parameter_labels(gm) == [("A", 0, "A")]
     theta0 = 0.5
@@ -254,7 +254,7 @@ def loglik_free(theta: np.ndarray) -> float:
         perturbed[1][0, 0] = theta[1]
         perturbed[2][0, 0] = 1.0 - theta[0] - theta[1]
         _alpha, log_scales = _forward_scaled(pi, perturbed, list(obs))
-        return float(log_scales.sum())
+        return float(log_scales.sum()) * np.log(2)
 
     base = np.array([0.2, 0.3])
     h = 1e-5
@@ -380,3 +380,22 @@ def test_baum_welch_restarts_reproducible_and_validated():
     assert log_likelihood(first, data) == pytest.approx(log_likelihood(second, data))
     with pytest.raises(ValueError):
         baum_welch(_symmetric_two_state(), data, n_restarts=0)
+
+
+def test_viterbi_impossible_after_first_step_has_no_path():
+    gm = golden_mean(0.5)
+    observations = [0, 1, 1, 0]
+    assert log_likelihood(gm, observations) == float("-inf")
+    assert viterbi(gm, observations) == []
+
+
+def test_baum_welch_rejects_data_with_zero_probability():
+    gm = golden_mean(0.5)
+    with pytest.raises(ValueError, match="zero probability"):
+        baum_welch(gm, [[1, 1], [0, 1, 1]])
+
+
+def test_baum_welch_warns_when_some_sequences_are_impossible():
+    gm = golden_mean(0.5)
+    with pytest.warns(RuntimeWarning, match="zero probability"):
+        baum_welch(gm, [[0, 1, 0, 0], [1, 1]], max_iter=3)
diff --git a/tests/test_measures.py b/tests/test_measures.py
index 8fb4c43..43600f5 100644
--- a/tests/test_measures.py
+++ b/tests/test_measures.py
@@ -136,7 +136,7 @@ def test_golden_mean_shift_parry_entropy_rate():
         np.array([[1, 1], [1, 0]], dtype=float),
         symbol_alphabet=frozenset({0, 1}),
     )
-    assert parry.entropy_rate() == pytest.approx(tmc.topological_entropy() / np.log(2), rel=0.05)
+    assert parry.entropy_rate() == pytest.approx(tmc.topological_entropy(), rel=0.05)
 
 
 def test_collision_entropy_nmachine():
diff --git a/tests/test_model_selection.py b/tests/test_model_selection.py
index 39d661b..e05f1be 100644
--- a/tests/test_model_selection.py
+++ b/tests/test_model_selection.py
@@ -63,9 +63,9 @@ def test_score_model_relationships():
     scores = score_model(golden_mean(0.3), data)
     k, n, ll = scores.num_parameters, scores.num_observations, scores.log_likelihood
     assert n == 2000
-    assert scores.aic == pytest.approx(2 * k - 2 * ll)
-    assert scores.bic == pytest.approx(k * np.log(n) - 2 * ll)
-    assert scores.mdl == pytest.approx(0.5 * k * np.log(n) - ll)
+    assert scores.aic == pytest.approx(2 * k - 2 * ll * np.log(2))
+    assert scores.bic == pytest.approx(k * np.log(n) - 2 * ll * np.log(2))
+    assert scores.mdl == pytest.approx(0.5 * k * np.log2(n) - ll)
     assert scores.aicc == pytest.approx(scores.aic + 2 * k * (k + 1) / (n - k - 1))
 
 
@@ -139,7 +139,7 @@ def test_waic_zero_variance_matches_deviance():
     result = waic(matrix)
     assert result.p_waic == pytest.approx(0.0, abs=1e-12)
     assert result.lppd == pytest.approx(-3.5)
-    assert result.waic == pytest.approx(7.0)
+    assert result.waic == pytest.approx(7.0 * np.log(2))
 
 
 def test_waic_positive_effective_parameters():
@@ -189,6 +189,7 @@ def test_rank_topological_epsilon_machines_prefers_two_states():
     assert best.criterion_value < ranked[-1].criterion_value
 
 
+@pytest.mark.filterwarnings("ignore:.*zero probability under the model:RuntimeWarning")
 def test_cross_validation_smoothing_keeps_forbidden_folds_finite():
     """The golden mean forbids 11; a held-out '11' makes an unsmoothed fold -inf."""
     rng = np.random.default_rng(6)
diff --git a/tests/test_sliding_block_code.py b/tests/test_sliding_block_code.py
index 3d3ac25..25ca3f4 100644
--- a/tests/test_sliding_block_code.py
+++ b/tests/test_sliding_block_code.py
@@ -66,3 +66,27 @@ def test_memoryless_transducer_round_trip():
     code = BitFlip().to_sliding_block_code()
     assert code.memory == 0
     assert code.apply_word(["0"]) == ("1",)
+
+
+def _golden_mean_shift() -> SoficShift:
+    shift = SoficShift(symbol_alphabet=frozenset({"0", "1"}))
+    shift.graph.add_state("A")
+    shift.graph.add_state("B")
+    shift.add_transition("A", "A", "0")
+    shift.add_transition("A", "B", "1")
+    shift.add_transition("B", "A", "0")
+    return shift
+
+
+def test_apply_keeps_constraints_longer_than_the_window():
+    identity = SlidingBlockCode({("0",): "0", ("1",): "1"})
+    shift = _golden_mean_shift()
+    image = identity.apply(shift)
+    for length in range(1, 7):
+        assert set(image.factor_language(length)) == set(shift.factor_language(length))
+
+
+def test_apply_rejects_partial_block_map():
+    partial = SlidingBlockCode({("0",): "0"}, input_alphabet={"0", "1"})
+    with pytest.raises(ValueError, match="not in block_map"):
+        partial.apply(_golden_mean_shift())
diff --git a/tests/test_tmc.py b/tests/test_tmc.py
index 69d3b8a..0daa712 100644
--- a/tests/test_tmc.py
+++ b/tests/test_tmc.py
@@ -31,4 +31,4 @@ def test_parry_measure():
     tmc = TopologicalMarkovChain.from_adjacency(np.array([[0, 1], [1, 1]]), symbol_alphabet=frozenset({"0", "1"}))
     parry = tmc.parry_measure()
     parry.validate()
-    assert parry.entropy_rate() == pytest.approx(tmc.topological_entropy() / np.log(2), rel=0.1)
+    assert parry.entropy_rate() == pytest.approx(tmc.topological_entropy(), rel=0.1)
diff --git a/tests/test_topological_anatomy.py b/tests/test_topological_anatomy.py
index 56855da..0a7a2b3 100644
--- a/tests/test_topological_anatomy.py
+++ b/tests/test_topological_anatomy.py
@@ -7,7 +7,6 @@
 import numpy as np
 import pytest
 
-from sofic.exceptions import UnifilarityError
 from sofic.generators.epsilon_machine import EpsilonMachine
 from sofic.shifts.sofic import SoficShift
 from sofic.shifts.tmc import TopologicalMarkovChain
@@ -88,7 +87,7 @@ def test_split_refines_h_top():
 def test_h_top_matches_topological_entropy_and_is_additive(builder):
     shift = builder()
     anatomy = shift.topological_anatomy()
-    assert anatomy["h_top"] == pytest.approx(shift.topological_entropy() / np.log(2), abs=1e-9)
+    assert anatomy["h_top"] == pytest.approx(shift.topological_entropy(), abs=1e-9)
     assert anatomy["h_top"] == pytest.approx(anatomy["b_top"] + anatomy["r_top"], abs=1e-9)
 
 
@@ -142,8 +141,10 @@ def test_determinizes_nondeterministic_full_shift():
     assert anatomy["b_top"] == pytest.approx(0.0, abs=1e-9)
 
 
-def test_non_auto_determinizable_presentation_raises():
-    """A nondeterministic presentation the bounded Fischer cover cannot resolve raises."""
+def test_nondeterministic_presentation_is_resolved_by_exact_fischer_cover():
+    """A presentation the old bounded follower-set cover could not resolve now works."""
+    from sofic.shifts.covers import RightFischerCover
+
     shift = SoficShift(symbol_alphabet=frozenset({0, 1}))
     for state in ("u", "v"):
         shift.graph.add_state(state)
@@ -152,5 +153,8 @@ def test_non_auto_determinizable_presentation_raises():
     shift.add_transition("u", "v", 1)
     shift.add_transition("v", "u", 0)
     assert not shift.is_unifilar()
-    with pytest.raises(UnifilarityError):
-        shift.topological_anatomy()
+    cover = RightFischerCover.from_sofic(shift)
+    assert cover.is_unifilar()
+    anatomy = shift.topological_anatomy()
+    assert anatomy["h_top"] == pytest.approx(cover.topological_entropy(), abs=1e-9)
+    assert anatomy["h_top"] == pytest.approx(anatomy["b_top"] + anatomy["r_top"], abs=1e-9)
diff --git a/tests/test_yaml.py b/tests/test_yaml.py
index da8bc67..b5d8237 100644
--- a/tests/test_yaml.py
+++ b/tests/test_yaml.py
@@ -342,3 +342,29 @@ def test_read_write_yaml_file(tmp_path):
 
     assert type(restored) is type(eps)
     assert model_to_dict(restored) == model_to_dict(eps)
+
+
+def test_cover_and_symbolic_models_round_trip():
+    from sofic.shifts.base import SymbolicModel
+    from sofic.shifts.covers import (
+        LeftFischerCover,
+        LeftKriegerCover,
+        RightFischerCover,
+        RightKriegerCover,
+        WheelerCover,
+    )
+
+    shift = SoficShift(symbol_alphabet=frozenset({"0", "1"}))
+    for source, target, symbol in (("A", "A", "0"), ("A", "B", "1"), ("B", "A", "0")):
+        shift.graph.add_state(source)
+        shift.graph.add_transition(source, target, **{ATTR_SYMBOL: symbol})
+    for cls in (LeftFischerCover, RightFischerCover, LeftKriegerCover, RightKriegerCover, WheelerCover):
+        cover = cls.from_sofic(shift)
+        restored = _round_trip(cover)
+        assert type(restored) is cls
+        assert sorted(map(repr, restored.states())) == sorted(map(repr, cover.states()))
+
+    bare = SymbolicModel(symbol_alphabet=frozenset({"x"}))
+    bare.graph.add_state(0)
+    bare.graph.add_transition(0, 0, **{ATTR_SYMBOL: "x"})
+    assert type(_round_trip(bare)) is SymbolicModel