diff --git a/CHANGELOG.md b/CHANGELOG.md
index 76cd143d2..07cdc910b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,19 @@
 # Changelog
 
+## Unreleased
+
+### Bugfix
+
+- Solves the Vandermonde interpolation systems of the polynomial TreeSHAP machinery (`TreeSHAPIQ`, `LinearTreeSHAP`) exactly at every depth, replacing the explicit float64 inverse:
+  - an O(n^2) Björck-Pereyra recursion in scaled-integer arithmetic with a bitwise convergence certificate returns the float64 rounding of the exact rational solution for any grid of distinct nodes,
+  - this **fixes a pre-existing silent-wrong-answer bug**: the previous inversion drifted at the ~1e-7 level from interpolation degree ~20, returned silently wrong values from ~32, and crashed with an unexplained `LinAlgError` at ~60+; the exact solves carry no conditioning error at any degree,
+  - beyond the float64 representation limit of the monomial-basis pipeline (interpolation degree ~29 for `LinearTreeSHAP` and ~25 for `TreeSHAPIQ` on the default grids, where the exact coefficients exceed ~3e10 and downstream double-precision evaluation provably cancels more than ~1% of the result) construction raises an explanatory `RepresentationLimitError` (a `ValueError`) instead of returning silently wrong values. [#545](https://github.com/mmschlk/shapiq/issues/545)
+
+### Changed
+
+- Models whose interpolation degree exceeds the float64 representation limit (~29 for `LinearTreeSHAP`, ~25 for `TreeSHAPIQ`, on the default grids) now raise `RepresentationLimitError` (a `ValueError`) at construction; previously they returned silently wrong values (degrees ~32-59) or crashed with `LinAlgError` (~60+). The degree depends on the configuration: the tree depth for `LinearTreeSHAP`, `min(depth, features in the tree)` for `TreeSHAPIQ` with SV/SII/k-SII, and the number of features in the tree for STII/FSII/BII. For order-1 SV/SII configurations `TreeExplainer` re-routes the affected trees to `TreeSHAPIQ` when its feature-bounded degree still fits; other configurations propagate the error.
+- Shapley values for interpolation degrees above ~20 may shift at the ~1e-7 level relative to previous releases: the coefficients are now the correctly rounded exact solutions rather than the output of a float64 inverse. Below the representation limit the remaining end-to-end completeness error is bounded by the downstream float64 evaluation (about `max|N| * 1e-13`, i.e. up to ~1e-4 at degree 24 and ~1e-2 at degree 28 on adversarial chain trees).
+
 ## v1.5.1 (2026-05-30)
 
 ### Bugfix
diff --git a/docs/source/_templates/autosummary/module.rst b/docs/source/_templates/autosummary/module.rst
index 40dbdeff9..ec19e6c25 100644
--- a/docs/source/_templates/autosummary/module.rst
+++ b/docs/source/_templates/autosummary/module.rst
@@ -14,6 +14,17 @@
 {%- endfor %}
 {% endif %}
 
+{% if exceptions %}
+.. rubric:: Exceptions
+
+.. autosummary::
+   :toctree:
+   :nosignatures:
+{% for item in exceptions %}
+   {{ item }}
+{%- endfor %}
+{% endif %}
+
 {% if functions %}
 .. rubric:: Functions
 
diff --git a/src/shapiq/__init__.py b/src/shapiq/__init__.py
index 140804b74..d3ecf02ed 100644
--- a/src/shapiq/__init__.py
+++ b/src/shapiq/__init__.py
@@ -79,8 +79,9 @@
 )
 from .tree import TreeExplainer
 
-# public utils functions
-from .utils import (  # sets.py  # tree.py
+# public utils functions and errors
+from .utils import (
+    RepresentationLimitError,
     get_explicit_subsets,
     powerset,
     safe_isinstance,
@@ -114,6 +115,8 @@
     "ProxySHAP",
     "ProxySPEX",
     "RegressionMSR",
+    # errors
+    "RepresentationLimitError",
     # explainers
     "Explainer",
     "TabularExplainer",
diff --git a/src/shapiq/tree/__init__.py b/src/shapiq/tree/__init__.py
index 547da8690..968fa19be 100644
--- a/src/shapiq/tree/__init__.py
+++ b/src/shapiq/tree/__init__.py
@@ -4,6 +4,9 @@
 along with supporting data structures and algorithm variants.
 """
 
+# re-exported for convenience; documented (and listed in __all__) only in shapiq.utils
+from shapiq.utils.errors import RepresentationLimitError  # noqa: F401
+
 from .base import TreeModel
 from .interventional import InterventionalGame, InterventionalTreeExplainer
 from .linear import LinearTreeSHAP
diff --git a/src/shapiq/tree/_numerics.py b/src/shapiq/tree/_numerics.py
new file mode 100644
index 000000000..924a9b3ce
--- /dev/null
+++ b/src/shapiq/tree/_numerics.py
@@ -0,0 +1,262 @@
+"""Exact solves for the Vandermonde systems of the polynomial TreeSHAP machinery.
+
+``TreeSHAPIQ`` and ``LinearTreeSHAP`` repeatedly need ``inv(vander(D[:i]).T) @ rhs``
+for prefixes ``D[:i]`` of a fixed interpolation grid. These Vandermonde systems are
+severely ill-conditioned in double precision — the conditioning is non-monotonic in
+the prefix length and peaks at interior prefixes (``i ~ n/2``), so explicit
+``float64`` inversion drifts at the ~1e-7 level from interpolation degree ~20,
+returns wrong values from ~32 on the default grids, and raises ``LinAlgError``
+for very deep trees (~60+).
+
+The conditioning is purely a floating-point artifact: the interpolation nodes are
+distinct, so the systems are exactly solvable over the rationals. This module
+therefore solves them **exactly** and returns the correctly rounded ``float64``
+result, at every depth:
+
+- **O(n^2) Björck-Pereyra recursion** instead of an O(n^3) inverse. The dual
+  Vandermonde recursion (Golub & Van Loan, Algorithm 4.6.2) factors the inverse
+  into bidiagonal steps, which keeps the exact-arithmetic operand sizes small.
+- **Scaled-integer fixed-point arithmetic** (standard-library ``int``): every node
+  and intermediate value is represented as ``round(v * 2**bits)``. With ``bits``
+  large enough the rounded result equals the exact rational solution rounded to
+  ``float64``; integer arithmetic keeps a depth-100 grid's full prefix workload in
+  the sub-second range, where naive ``fractions.Fraction`` elimination needs minutes.
+- **A convergence certificate instead of an error analysis**: each system is solved
+  at increasing precision (128, 256, ... bits) until two consecutive precision
+  levels produce bitwise-identical ``float64`` outputs. Agreement of two scaled
+  truncation grids that fine means the value has converged to the rounding of the
+  exact solution, so the returned coefficients carry no conditioning error at all.
+
+Every path either returns the exact result or refuses loudly: degenerate inputs
+(coincident, non-finite, or extreme custom nodes) raise ``ValueError``, and
+N-matrix rows whose magnitude exceeds what the downstream float64 pipeline can
+evaluate raise :class:`~shapiq.utils.errors.RepresentationLimitError`.
+"""
+
+from __future__ import annotations
+
+from fractions import Fraction
+from functools import lru_cache
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+from shapiq.utils.errors import RepresentationLimitError
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+# Precision ladder for the convergence certificate. 128 bits converges for every
+# Chebyshev-grid prefix up to at least depth 100; the 256-bit rung confirms it
+# whenever the cheap float64 cross-check cannot (ill-conditioned prefixes), and
+# the higher rungs exist for adversarial custom grids (LinearTreeSHAP(base_func=...)).
+_PRECISION_BITS = (128, 256, 512, 1024, 2048, 4096)
+
+# Largest N-matrix entry the downstream float64 pipeline can absorb. The
+# coefficients themselves are exact, but the explainers consume them in float64
+# inner products whose summands scale with the largest entry while the result is
+# O(prediction). Measured end to end on chain trees of depth 20-40, the relative
+# cancellation error tracks ``max|N| * 1e-13`` within one order of magnitude, so
+# at 3e10 the expected loss is in the 0.3%-3% range. Beyond that the explainers
+# refuse instead of returning values whose error may exceed ~1% of the
+# prediction. For the default Chebyshev grids this corresponds to an
+# interpolation degree of ~29 for LinearTreeSHAP's right-hand sides and ~25
+# for TreeSHAPIQ's (its identity N matrix, built for every index, saturates
+# first). The bound is empirical, not a worst-case proof.
+_REPRESENTATION_LIMIT = 3.0e10
+
+
+def _check_row_magnitude(row: np.ndarray) -> None:
+    """Refuse N-matrix rows whose magnitude exceeds the float64 pipeline's limit.
+
+    The Vandermonde solves themselves are exact at any depth, but the
+    monomial-basis N-matrix entries grow exponentially with the interpolation
+    degree, and the downstream float64 evaluation cancels them against each
+    other. Beyond ``_REPRESENTATION_LIMIT`` the empirically calibrated loss may
+    exceed ~1% of the result, so this raises instead of letting the explainer
+    return silently wrong values.
+
+    Raises:
+        RepresentationLimitError: If the row magnitude exceeds the limit.
+    """
+    peak = float(np.max(np.abs(row)))
+    if peak > _REPRESENTATION_LIMIT:
+        msg = (
+            "The interpolation degree of this tree is too large for the float64 "
+            f"TreeSHAP pipeline: the exact interpolation coefficients reach {peak:.1e}, "
+            "and evaluating them in double precision would lose more than ~1% of the "
+            "result to cancellation. Reduce the interpolation degree — depending on "
+            "the explainer and index this is the tree depth, min(tree depth, features "
+            "in the tree), or the number of features in the tree; the default grids "
+            "support up to ~29 with LinearTreeSHAP and ~25 with TreeSHAPIQ. This is a "
+            "structural limit of the monomial-basis representation, not of the solver."
+        )
+        raise RepresentationLimitError(msg)
+
+
+def build_n_matrix(
+    grid: np.ndarray,
+    rhs_for_row: Callable[[int], np.ndarray],
+) -> np.ndarray:
+    """Build an interpolation N matrix row by row from exact Vandermonde solves.
+
+    Row ``i`` is the solution of ``vander(grid[:i]).T @ x = rhs_for_row(i)``.
+    Every row is checked against the float64 representation limit as soon as it
+    is computed, so an over-deep tree is refused after the first offending row
+    instead of after the full (expensive) construction.
+
+    Args:
+        grid: The full interpolation grid (its length is the maximum degree).
+        rhs_for_row: Callback returning the length-``i`` right-hand side of row ``i``.
+
+    Returns:
+        The ``(len(grid) + 1, len(grid))`` N matrix (row 0 stays zero).
+
+    Raises:
+        RepresentationLimitError: If any row exceeds the representation limit.
+        ValueError: For degenerate grids (coincident, non-finite, or extreme
+            nodes), propagated from :func:`solve_vandermonde`.
+    """
+    depth = len(grid)
+    n_matrix = np.zeros((depth + 1, depth))
+    for i in range(1, depth + 1):
+        n_matrix[i, :i] = solve_vandermonde(grid[:i], rhs_for_row(i))
+        _check_row_magnitude(n_matrix[i, :i])
+    return n_matrix
+
+
+def _to_fixed_point(value: float, bits: int) -> int:
+    """Return ``value * 2**bits`` as an integer (truncated toward minus infinity).
+
+    The truncation direction is irrelevant for the convergence certificate: the
+    per-step error stays below ``2**-bits`` in either case, and the certificate
+    only accepts results that are stable under a 128-bit precision increase.
+    """
+    fraction = Fraction(value)
+    return (fraction.numerator << bits) // fraction.denominator
+
+
+def _bjorck_pereyra_fixed_point(points: np.ndarray, rhs: np.ndarray, bits: int) -> np.ndarray:
+    """Solve ``vander(points).T @ x = rhs`` in scaled-integer arithmetic.
+
+    Implements the dual Björck-Pereyra recursion on values scaled by ``2**bits``.
+    ``np.vander`` uses decreasing powers, so the right-hand side is reversed to
+    obtain the increasing-power moment problem the recursion expects.
+    """
+    n = len(points)
+    nodes = [_to_fixed_point(p, bits) for p in points.tolist()]
+    x = [_to_fixed_point(v, bits) for v in rhs.tolist()[::-1]]
+    for k in range(n - 1):
+        node_k = nodes[k]
+        for i in range(n - 1, k, -1):
+            x[i] = x[i] - (node_k * x[i - 1] >> bits)
+    for k in range(n - 2, -1, -1):
+        for i in range(k + 1, n):
+            x[i] = (x[i] << bits) // (nodes[i] - nodes[i - k - 1])
+        for i in range(k + 1, n):
+            x[i - 1] = x[i - 1] - x[i]
+    scale = 1 << bits
+    return np.array([float(Fraction(value, scale)) for value in x])
+
+
+def solve_vandermonde(points: np.ndarray, rhs: np.ndarray) -> np.ndarray:
+    """Solve ``vander(points).T @ x = rhs`` exactly, rounded to ``float64``.
+
+    The system is solved in scaled-integer arithmetic and certified by agreement
+    with an independent computation: first against a plain float64 LAPACK solve
+    (which corroborates the well-conditioned common case at negligible cost), and
+    otherwise against the next rung of the precision ladder, until two
+    computations agree bitwise. An agreeing result equals the exact rational
+    solution rounded to ``float64`` up to an overwhelmingly improbable
+    coincidence of correlated rounding (both computations would have to land in
+    the same wrong rounding bin), regardless of the conditioning of the system.
+    Works for any grid of distinct interpolation nodes.
+
+    Solves are memoized on the byte representation of ``(points, rhs)``: the
+    explainers issue identical solves for every tree of equal depth in an
+    ensemble, which therefore costs one tree's worth of work.
+
+    Args:
+        points: Interpolation nodes (a prefix of the interpolation grid). The nodes
+            must be pairwise distinct.
+        rhs: Right-hand side vector of the same length.
+
+    Returns:
+        The solution vector ``x`` (the ``float64`` rounding of the exact solution).
+
+    Raises:
+        ValueError: If the inputs are degenerate — coincident or non-finite
+            nodes, a non-finite right-hand side, or custom nodes so close
+            together (or so large) that the exact solution leaves the float64
+            range or does not stabilize within the precision ladder. The
+            default Chebyshev grids never trigger any of these.
+    """
+    points = np.asarray(points, dtype=float)
+    rhs = np.asarray(rhs, dtype=float)
+    if points.ndim != 1 or points.shape != rhs.shape:
+        msg = "The interpolation nodes and right-hand side must be 1-D arrays of equal length."
+        raise ValueError(msg)
+    if not (np.isfinite(points).all() and np.isfinite(rhs).all()):
+        msg = "The interpolation nodes and right-hand side must be finite."
+        raise ValueError(msg)
+    if len(points) == 1:
+        return rhs.copy()
+    solution = _solve_vandermonde_cached(points.tobytes(), rhs.tobytes())
+    return solution.copy()  # the cached array must stay immutable
+
+
+_DEGENERATE_NODES_MSG = (
+    "The exact solution leaves the float64 range or does not stabilize within "
+    f"{_PRECISION_BITS[-1]} fractional bits; the interpolation nodes are too close "
+    "together or too large in magnitude. Use better-separated nodes of moderate "
+    "magnitude."
+)
+
+
+def clear_solver_cache() -> None:
+    """Drop all memoized Vandermonde solutions (for cold-timing tests or memory-sensitive hosts)."""
+    _solve_vandermonde_cached.cache_clear()
+
+
+@lru_cache(maxsize=4096)
+def _solve_vandermonde_cached(points_bytes: bytes, rhs_bytes: bytes) -> np.ndarray:
+    n = len(points_bytes) // 8  # float64 buffers; lengths are equal by construction
+    points = np.frombuffer(points_bytes, dtype=float, count=n)
+    rhs = np.frombuffer(rhs_bytes, dtype=float, count=n)
+
+    previous = None
+    for bits in _PRECISION_BITS:
+        try:
+            current = _bjorck_pereyra_fixed_point(points, rhs, bits)
+        except ZeroDivisionError:
+            # Two nodes collapsed onto the same scaled integer at this precision.
+            # Genuinely distinct nodes separate at a finer rung (and stay separated
+            # at every rung above it), so climb instead of misreporting them as
+            # coincident; truly coincident nodes collapse at every rung and fall
+            # through to the error below.
+            previous = None
+            continue
+        except OverflowError:
+            # float(Fraction(...)) overflows when the exact solution leaves the
+            # float64 range (nearly-coincident or very large custom nodes).
+            raise ValueError(_DEGENERATE_NODES_MSG) from None
+        if previous is None:
+            # Cheap independent certificate: a float64 LAPACK solve corroborates
+            # the well-conditioned common case at a fraction of the cost of a
+            # second fixed-point solve. It must solve the same decreasing-power
+            # system the fixed-point recursion solves by reversing the rhs.
+            try:
+                float_solution = np.linalg.solve(np.vander(points).T, rhs)
+                if np.array_equal(current, float_solution):
+                    return current
+            except np.linalg.LinAlgError:  # float64 considers the system singular
+                pass
+        elif np.array_equal(previous, current):
+            return current
+        previous = current
+    if previous is None:
+        msg = "The interpolation nodes must be pairwise distinct."
+        raise ValueError(msg)
+    # Reachable only for custom grids with nearly-coincident nodes, where no
+    # float64 rounding of the (astronomically large) exact solution stabilizes.
+    raise ValueError(_DEGENERATE_NODES_MSG)
diff --git a/src/shapiq/tree/explainer.py b/src/shapiq/tree/explainer.py
index 6435233d3..12335cafc 100644
--- a/src/shapiq/tree/explainer.py
+++ b/src/shapiq/tree/explainer.py
@@ -11,6 +11,7 @@
 
 from shapiq.explainer.base import Explainer
 from shapiq.tree.interventional.explainer import InterventionalTreeExplainer
+from shapiq.utils.errors import RepresentationLimitError
 
 from .linear import LinearTreeSHAP
 from .treeshapiq import TreeSHAPIQ, TreeSHAPIQIndices
@@ -97,6 +98,10 @@ def __init__(
 
             **kwargs: Additional keyword arguments are ignored.
 
+        Raises:
+            RepresentationLimitError: If a tree's interpolation degree exceeds the
+                float64 representation limit and no re-route applies (see
+                :class:`~shapiq.utils.errors.RepresentationLimitError`).
         """
         super().__init__(model, index=index, max_order=max_order)
 
@@ -116,28 +121,42 @@ def __init__(
         self.mode = mode
         self._reference_dataset: np.ndarray | None = reference_dataset
 
-        # In ``"pathdependent"`` mode, build exactly one per-tree explainer list — either
-        # ``LinearTreeSHAP`` (cheap, order-1 only) or ``TreeSHAPIQ`` (any order). The dispatch
-        # decision is fixed at construction time so callers can mutate the chosen list (e.g.
-        # ``_tree.thresholds`` rounding in tests) before calling :meth:`explain`. In
+        # In ``"pathdependent"`` mode, build per-tree explainers at construction time:
+        # ``LinearTreeSHAP`` (cheap, order-1 only) where possible, ``TreeSHAPIQ`` (any
+        # order) otherwise. The two lists can BOTH be populated when individual trees
+        # are re-routed (a tree whose depth-based interpolation degree exceeds the
+        # float64 representation limit while its feature-bounded TreeSHAPIQ degree
+        # does not); :meth:`explain` then aggregates across both lists. In
         # ``"interventional"`` mode no per-tree list is created — the
-        # :class:`~shapiq.tree.interventional.explainer.InterventionalTreeExplainer` handles the
-        # full ensemble in one shot, so a per-tree list would be meaningless.
+        # :class:`~shapiq.tree.interventional.explainer.InterventionalTreeExplainer`
+        # handles the full ensemble in one shot.
         self._treeshapiq_explainers: list[TreeSHAPIQ] = []
         self._lineartreeshap_explainers: list[LinearTreeSHAP] = []
         self._interventional_explainer: InterventionalTreeExplainer | None = None
 
+        # ``index`` (the local parameter) is already narrowed to ``TreeSHAPIQIndices``;
+        # ``self.index`` is the broader ``ExplainerIndices`` and would not type-check.
+        def _build_treeshapiq(tree: TreeModel, tree_index: TreeSHAPIQIndices) -> TreeSHAPIQ:
+            return TreeSHAPIQ(model=tree, max_order=self._max_order, index=tree_index)
+
         if self.mode == "pathdependent":
             if self._can_use_lineartreeshap():
-                self._lineartreeshap_explainers = [
-                    LinearTreeSHAP(model=tree) for tree in self._trees
-                ]
+                for tree in self._trees:
+                    try:
+                        self._lineartreeshap_explainers.append(LinearTreeSHAP(model=tree))
+                    except RepresentationLimitError:
+                        # LinearTreeSHAP's interpolation degree is the full tree
+                        # depth, while TreeSHAPIQ's is min(depth, features in the
+                        # tree), so a deep tree over few features can still be
+                        # explained there — at the same order-1 Shapley values.
+                        # Only the affected tree leaves the fast path; if its
+                        # TreeSHAPIQ degree also exceeds the limit, the error
+                        # propagates. ``index="SV"`` matches LinearTreeSHAP's
+                        # output label (order-1 indices all reduce to SV).
+                        self._treeshapiq_explainers.append(_build_treeshapiq(tree, "SV"))
             else:
-                # ``index`` (the local parameter) is already narrowed to ``TreeSHAPIQIndices``;
-                # ``self.index`` is the broader ``ExplainerIndices`` and would not type-check.
                 self._treeshapiq_explainers = [
-                    TreeSHAPIQ(model=tree, max_order=self._max_order, index=index)
-                    for tree in self._trees
+                    _build_treeshapiq(tree, index) for tree in self._trees
                 ]
         elif self.mode == "interventional":
             if self._reference_dataset is None:
@@ -173,45 +192,6 @@ def _can_use_lineartreeshap(self) -> bool:
             and all(tree.n_features_in_tree >= 2 for tree in self._trees)
         )
 
-    def _explain_function_lineartreeshap(
-        self,
-        x: np.ndarray,
-        **kwargs: Any,  # noqa: ARG002
-    ) -> InteractionValues:
-        """Compute first-order Shapley values for ``x`` by aggregating the per-tree LinearTreeSHAP results.
-
-        Mirrors the per-tree aggregation done by ``_explain_function_treeshapiq``: each
-        ``LinearTreeSHAP`` in ``self._lineartreeshap_explainers`` runs against ``x``, the
-        resulting :class:`~shapiq.interaction_values.InteractionValues` are summed (which also
-        sums ``baseline_value`` and the ``()`` entry), and ``min_order`` is finally enforced via
-        :meth:`InteractionValues.get_n_order` when the user asked for a stricter minimum.
-
-        Args:
-            x: The instance to explain as a 1-dimensional array.
-            **kwargs: Additional keyword arguments are ignored.
-
-        Returns:
-            The aggregated Shapley values for the instance.
-        """
-        if len(x.shape) != 1:
-            msg = "explain expects a single instance, not a batch."
-            raise TypeError(msg)
-
-        interaction_values: list[InteractionValues] = [
-            lts.explain_function(x) for lts in self._lineartreeshap_explainers
-        ]
-        final_explanation = interaction_values[0]
-        for iv in interaction_values[1:]:
-            final_explanation += iv
-
-        if self._min_order > final_explanation.min_order:
-            final_explanation = final_explanation.get_n_order(
-                min_order=self._min_order,
-                max_order=self._max_order,
-            )
-
-        return final_explanation
-
     def _explain_function_interventionaltreeshapiq(
         self,
         x: np.ndarray,
@@ -231,7 +211,7 @@ def _explain_function_interventionaltreeshapiq(
             raise RuntimeError(msg)
         return self._interventional_explainer.explain_function(x)
 
-    def _explain_function_treeshapiq(
+    def _explain_function_pathdependent(
         self,
         x: np.ndarray,
         **kwargs: Any,  # noqa: ARG002
@@ -249,17 +229,19 @@ def _explain_function_treeshapiq(
         if len(x.shape) != 1:
             msg = "explain expects a single instance, not a batch."
             raise TypeError(msg)
-        # run treeshapiq for all trees
-        interaction_values: list[InteractionValues] = []
+        # run the per-tree explainers; both lists can be populated when trees
+        # were re-routed at construction time (the other list is then empty)
+        interaction_values: list[InteractionValues] = [
+            lts.explain_function(x) for lts in self._lineartreeshap_explainers
+        ]
         for explainer in self._treeshapiq_explainers:
             tree_explanation = explainer.explain(x)
             interaction_values.append(tree_explanation)
 
         # combine the explanations for all trees
         final_explanation = interaction_values[0]
-        if len(interaction_values) > 1:
-            for i in range(1, len(interaction_values)):
-                final_explanation += interaction_values[i]
+        for tree_explanation in interaction_values[1:]:
+            final_explanation += tree_explanation
 
         if self._min_order == 0 and final_explanation.min_order == 1:
             final_explanation.min_order = 0
@@ -295,8 +277,5 @@ def explain_function(  # type: ignore[override]
             The computed interaction index for the instance.
         """
         if self.mode == "pathdependent":
-            # Dispatch on whichever per-tree list __init__ chose to populate.
-            if self._lineartreeshap_explainers:
-                return self._explain_function_lineartreeshap(x, **kwargs)
-            return self._explain_function_treeshapiq(x, **kwargs)
+            return self._explain_function_pathdependent(x, **kwargs)
         return self._explain_function_interventionaltreeshapiq(x, **kwargs)
diff --git a/src/shapiq/tree/linear/explainer.py b/src/shapiq/tree/linear/explainer.py
index 2c565dc77..b639c6b1a 100644
--- a/src/shapiq/tree/linear/explainer.py
+++ b/src/shapiq/tree/linear/explainer.py
@@ -8,6 +8,7 @@
 import scipy.special as sp
 
 from shapiq.interaction_values import InteractionValues
+from shapiq.tree._numerics import build_n_matrix
 from shapiq.tree.conversion.edges import create_edge_tree
 from shapiq.tree.validation import validate_tree_model
 from shapiq.utils.sets import generate_interaction_lookup, powerset
@@ -36,11 +37,7 @@ def get_N_prime(max_size: int = 10) -> np.ndarray:
 
 def get_N_v2(D: np.ndarray) -> np.ndarray:
     """Get N_v2 matrix for Linear Tree Shap."""
-    depth = D.shape[0]
-    Ns = np.zeros((depth + 1, depth))
-    for i in range(1, depth + 1):
-        Ns[i, :i] = np.linalg.inv(np.vander(D[:i]).T).dot(1.0 / get_norm_weight(i - 1))
-    return Ns
+    return build_n_matrix(D, lambda i: 1.0 / get_norm_weight(i - 1))
 
 
 class LinearTreeSHAP:
@@ -76,6 +73,11 @@ def __init__(
                 :func:`~shapiq.tree.validation.validate_tree_model`.
             base_func: Callable ``int -> np.ndarray`` returning the interpolation base for the
                 given depth. Defaults to :func:`numpy.polynomial.chebyshev.chebpts2`.
+
+        Raises:
+            RepresentationLimitError: If the tree depth exceeds the float64
+                representation limit of the interpolation pipeline (see
+                :class:`~shapiq.utils.errors.RepresentationLimitError`).
         """
         self.clf = model
         self._tree = validate_tree_model(model, class_label=None)[0]
diff --git a/src/shapiq/tree/treeshapiq.py b/src/shapiq/tree/treeshapiq.py
index 450e9fe84..555fcd5a0 100644
--- a/src/shapiq/tree/treeshapiq.py
+++ b/src/shapiq/tree/treeshapiq.py
@@ -11,8 +11,10 @@
 
 from shapiq.game_theory.indices import get_computation_index
 from shapiq.interaction_values import InteractionValues
+from shapiq.utils.errors import RepresentationLimitError
 from shapiq.utils.sets import generate_interaction_lookup, powerset
 
+from ._numerics import build_n_matrix
 from .conversion.edges import create_edge_tree
 from .validation import validate_tree_model
 
@@ -153,6 +155,10 @@ def __init__(
         if self._n_features_in_tree > 0:
             try:
                 self._init_summary_polynomials()
+            except RepresentationLimitError:
+                # a structural float64 limit, never the trivial single-feature
+                # case the broader ValueError handler below absorbs
+                raise
             except ValueError:
                 if self._n_features_in_tree == 1:
                     self._trivial_computation = True  # for one feature the computation is trivial
@@ -709,23 +715,17 @@ def _get_n_matrix(interpolated_poly: np.ndarray) -> np.ndarray:
             The N matrix.
 
         """
-        depth = interpolated_poly.shape[0]
-        Ns = np.zeros((depth + 1, depth))
-        for i in range(1, depth + 1):
-            Ns[i, :i] = np.linalg.inv(np.vander(interpolated_poly[:i]).T).dot(
-                1.0 / np.array([sp.special.binom(i - 1, k) for k in range(i)])
-            )
-        return Ns
+        return build_n_matrix(
+            interpolated_poly,
+            lambda i: 1.0 / np.array([sp.special.binom(i - 1, k) for k in range(i)]),
+        )
 
     def _get_n_cii_matrix(self, interpolated_poly: np.ndarray, order: int) -> np.ndarray:
         """Computes the N matrix for the CII index."""
-        depth = interpolated_poly.shape[0]
-        Ns = np.zeros((depth + 1, depth))
-        for i in range(1, depth + 1):
-            Ns[i, :i] = np.linalg.inv(np.vander(interpolated_poly[:i]).T).dot(
-                i * np.array([self._get_subset_weight_cii(j, order) for j in range(i)]),
-            )
-        return Ns
+        return build_n_matrix(
+            interpolated_poly,
+            lambda i: i * np.array([self._get_subset_weight_cii(j, order) for j in range(i)]),
+        )
 
     def _get_subset_weight_cii(self, t: int, order: int) -> float | None:
         """Computes the weight for a given subset size and interaction order.
@@ -756,11 +756,7 @@ def _get_subset_weight_cii(self, t: int, order: int) -> float | None:
     @staticmethod
     def _get_n_id_matrix(D: np.ndarray) -> np.ndarray:
         """Computes N_id matrix."""
-        depth = D.shape[0]
-        Ns_id = np.zeros((depth + 1, depth))
-        for i in range(1, depth + 1):
-            Ns_id[i, :i] = np.linalg.inv(np.vander(D[:i]).T).dot(np.ones(i))
-        return Ns_id
+        return build_n_matrix(D, np.ones)
 
     @staticmethod
     def _cache(interpolated_poly: FloatVector) -> FloatVector:
diff --git a/src/shapiq/tree/validation.py b/src/shapiq/tree/validation.py
index c5e96456d..c38d1b0f9 100644
--- a/src/shapiq/tree/validation.py
+++ b/src/shapiq/tree/validation.py
@@ -56,7 +56,9 @@ def validate_tree_model(
 
     Raises:
         TypeError: If the model type is not supported (raised from the underlying
-            ``NotImplementedError`` of :func:`~shapiq.tree.conversion.convert_tree_model`).
+            ``NotImplementedError`` of :func:`~shapiq.tree.conversion.convert_tree_model`), or
+            if the model yields no trees (an empty list, or a list whose elements are not
+            ``TreeModel`` instances).
     """
     tree_model = []
     # direct returns for base tree models and dict as model
@@ -78,4 +80,9 @@ def validate_tree_model(
             msg = f"Model type {type(model)} is not supported."
             raise TypeError(msg) from e
         tree_model = result if isinstance(result, list) else [result]
+    if not tree_model:
+        # an empty list, or a list whose elements are not TreeModel instances,
+        # would otherwise propagate as an opaque IndexError at explain time
+        msg = "The model does not contain any trees: expected a non-empty (list of) TreeModel."
+        raise TypeError(msg)
     return tree_model
diff --git a/src/shapiq/utils/__init__.py b/src/shapiq/utils/__init__.py
index cc9747b64..6c1043305 100644
--- a/src/shapiq/utils/__init__.py
+++ b/src/shapiq/utils/__init__.py
@@ -1,7 +1,7 @@
 """Low-level utility functions for working with subsets, coalitions, and modules."""
 
 from .datasets import shuffle_data
-from .errors import raise_deprecation_warning
+from .errors import RepresentationLimitError, raise_deprecation_warning
 from .modules import check_import_module, safe_isinstance
 from .sets import (
     count_interactions,
@@ -32,5 +32,6 @@
     # datasets
     "shuffle_data",
     # errors
+    "RepresentationLimitError",
     "raise_deprecation_warning",
 ]
diff --git a/src/shapiq/utils/errors.py b/src/shapiq/utils/errors.py
index 19c8650b8..06b0caacf 100644
--- a/src/shapiq/utils/errors.py
+++ b/src/shapiq/utils/errors.py
@@ -17,3 +17,17 @@ def raise_deprecation_warning(
     )
 
     warn(message, DeprecationWarning, stacklevel=2)
+
+
+class RepresentationLimitError(ValueError):
+    """The interpolation degree exceeds what the float64 TreeSHAP pipeline can evaluate.
+
+    The exact interpolation coefficients are available at any depth, but the
+    downstream double-precision evaluation would cancel them beyond an acceptable
+    error. For order-1 SV/SII configurations, ``TreeExplainer`` catches this and
+    re-routes the affected tree to ``TreeSHAPIQ`` (whose interpolation degree is
+    bounded by the features in the tree rather than the tree depth); in every
+    other configuration — higher orders, the feature-bounded indices, or trees
+    that exceed the limit on both paths — the error propagates from the
+    constructor.
+    """
diff --git a/tests/shapiq/tests_unit/test_public_api.py b/tests/shapiq/tests_unit/test_public_api.py
index f7103e13a..8f6cd378d 100644
--- a/tests/shapiq/tests_unit/test_public_api.py
+++ b/tests/shapiq/tests_unit/test_public_api.py
@@ -78,7 +78,7 @@ def test_all_concrete_subclasses_in_all(module_path: str, base_path: str, label:
     exported = set(module.__all__)
     missing = concrete - exported
 
-    pkg_init = f"src/shapiq/{module_path.split('.')[-1]}/__init__.py"
+    pkg_init = f"src/shapiq/{module_path.rsplit('.', maxsplit=1)[-1]}/__init__.py"
     assert not missing, (
         f"Concrete {label} subclasses not listed in {module_path}.__all__: "
         f"{missing}. Add them to {pkg_init}."
diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py
index 17b318e66..58817026f 100644
--- a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py
+++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer.py
@@ -422,7 +422,7 @@ def test_xgboost_shap_error(xgb_clf_model, background_clf_data):
     # per-tree explainer list was populated so the mutation actually takes effect at explain time.
     per_tree_explainers = (
         explainer_shapiq_rounded._lineartreeshap_explainers
-        or explainer_shapiq_rounded._treeshapiq_explainers
+        + explainer_shapiq_rounded._treeshapiq_explainers
     )
     for tree_explainer in per_tree_explainers:
         tree_explainer._tree.thresholds = np.round(tree_explainer._tree.thresholds, 4)
diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py
index 89ed65829..b6c33fe77 100644
--- a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py
+++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_explainer_validate.py
@@ -33,6 +33,12 @@ def test_validate_model(dt_clf_model, dt_reg_model, rf_reg_model, rf_clf_model,
     with pytest.raises(TypeError):
         validate_tree_model("unsupported_model")
 
+    # treeless inputs raise instead of propagating an IndexError at explain time
+    with pytest.raises(TypeError, match="does not contain any trees"):
+        validate_tree_model([])
+    with pytest.raises(TypeError, match="does not contain any trees"):
+        validate_tree_model([dt_clf_model, dt_reg_model])  # a list of raw sklearn models
+
 
 @pytest.mark.external_libraries
 @pytest.mark.parametrize(("model_fixture", "model_class"), TREE_MODEL_FIXTURES)
diff --git a/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py
new file mode 100644
index 000000000..7a72ad02f
--- /dev/null
+++ b/tests/shapiq/tests_unit/tests_explainer/tests_tree_explainer/test_tree_numerics.py
@@ -0,0 +1,374 @@
+"""Tests for the exact Vandermonde solves in ``shapiq.tree._numerics``."""
+
+from __future__ import annotations
+
+from fractions import Fraction
+
+import numpy as np
+import pytest
+from numpy.polynomial.chebyshev import chebpts2
+
+from shapiq.tree import RepresentationLimitError, TreeExplainer, TreeModel
+from shapiq.tree._numerics import (
+    _solve_vandermonde_cached,
+    clear_solver_cache,
+    solve_vandermonde,
+)
+from shapiq.tree.linear.explainer import LinearTreeSHAP, get_norm_weight
+from shapiq.tree.treeshapiq import TreeSHAPIQ
+
+
+def _solve_exact_reference(points: np.ndarray, rhs: np.ndarray) -> np.ndarray:
+    """Oracle: exact rational Gaussian elimination for ``vander(points).T x = rhs``."""
+    n = len(points)
+    p = [Fraction(t) for t in points.tolist()]
+    rows = [
+        [p[j] ** (n - 1 - k) for j in range(n)] + [Fraction(v)] for k, v in enumerate(rhs.tolist())
+    ]
+    for col in range(n):
+        pivot = next(r for r in range(col, n) if rows[r][col] != 0)
+        rows[col], rows[pivot] = rows[pivot], rows[col]
+        for r in range(n):
+            if r != col and rows[r][col] != 0:
+                factor = rows[r][col] / rows[col][col]
+                rows[r] = [a - factor * b for a, b in zip(rows[r], rows[col], strict=True)]
+    return np.array([float(rows[k][n] / rows[k][k]) for k in range(n)])
+
+
+@pytest.mark.parametrize("grid_size", [2, 4, 8, 12, 16, 20])
+def test_matches_explicit_inverse_on_well_conditioned_grids(grid_size):
+    """On well-conditioned systems the result agrees with the previous inverse.
+
+    Larger grids are deliberately excluded here: above size ~20 the explicit
+    inverse starts drifting at the ~1e-7 level (at size 26 the peak prefix
+    condition number is already ~4e11), so the old formulation stops being a
+    valid reference; those sizes are covered by the exact-rational oracle test
+    below instead.
+    """
+    D = chebpts2(grid_size)
+    rng = np.random.default_rng(0)
+    for i in range(2, grid_size + 1):
+        rhs = rng.normal(size=i)
+        expected = np.linalg.inv(np.vander(D[:i]).T).dot(rhs)
+        np.testing.assert_allclose(solve_vandermonde(D[:i], rhs), expected, rtol=1e-7, atol=1e-9)
+
+
+@pytest.mark.parametrize("grid_size", [12, 28, 35, 45])
+def test_bitwise_equal_to_exact_rational_solution(grid_size):
+    """Every prefix solve equals the exact rational solution rounded to float64.
+
+    Sizes 28-45 are exactly the regime where the previous ``inv``-based code lost
+    precision (peak prefix condition numbers above 1e12) or silently returned wrong
+    values (rank-deficient to machine precision from size ~32).
+    """
+    D = chebpts2(grid_size)
+    rng = np.random.default_rng(1)
+    for i in (2, grid_size // 2, grid_size):  # interior prefixes are the worst-conditioned
+        rhs = rng.normal(size=i)
+        exact = _solve_exact_reference(D[:i], rhs)
+        result = solve_vandermonde(D[:i], rhs)
+        np.testing.assert_array_equal(result, exact)
+
+
+def test_replaced_inverse_construction_was_wrong_in_the_failure_band():
+    """Regression guard: the construction this module replaced fails where the solver is exact.
+
+    At grid size 32 with the library's own right-hand side the explicit inverse is
+    off by an absolute error of ~2e2 — entries of magnitude ~3e7 that downstream
+    inner products cancel to an O(prediction) result, so an absolute error of that
+    size destroys the result entirely. This test fails if the exact solver is ever
+    swapped back for the inverse.
+    """
+    D = chebpts2(32)
+    rhs = 1.0 / get_norm_weight(31)
+    exact = _solve_exact_reference(D, rhs)
+    np.testing.assert_array_equal(solve_vandermonde(D, rhs), exact)
+    try:
+        old = np.linalg.inv(np.vander(D).T) @ rhs
+        assert np.max(np.abs(old - exact)) > 1.0  # measured ~2e2; anything >1 is fatal downstream
+    except np.linalg.LinAlgError:
+        pass  # an outright failure is equally disqualifying for the old construction
+
+
+def test_issue_545_scenario_raises_instead_of_returning_wrong_values():
+    """The end-to-end scenario of issue #545 is refused instead of silently wrong.
+
+    A decision tree fit on one-hot-style data (every sample distinguished by its
+    own indicator feature) is forced into a depth-39 chain over 39 features, the
+    structure that previously returned silently wrong Shapley values (interpolation
+    degree in the ~32-59 failure band, no warning). Both the LinearTreeSHAP degree
+    (tree depth) and the TreeSHAPIQ degree (features in the tree) exceed the
+    representation limit, so there is no re-route and construction must raise.
+    """
+    sklearn_tree = pytest.importorskip("sklearn.tree")
+    x_data = np.eye(40)
+    y_target = np.arange(40, dtype=float)
+    model = sklearn_tree.DecisionTreeRegressor(random_state=0).fit(x_data, y_target)
+    assert model.get_depth() >= 30  # in the previously silent failure band
+    assert len({int(f) for f in model.tree_.feature if f >= 0}) >= 30  # no re-route possible
+    with pytest.raises(RepresentationLimitError):
+        TreeExplainer(model=model, index="SV", max_order=1)
+
+
+def test_single_point_system():
+    """The size-1 system is the identity."""
+    np.testing.assert_array_equal(solve_vandermonde(np.array([0.3]), np.array([2.5])), [2.5])
+
+
+def test_non_chebyshev_clustered_grid_is_solved_exactly():
+    """Custom interpolation grids (LinearTreeSHAP(base_func=...)) are handled exactly."""
+    grid = np.linspace(0.0, 1.0, 18) ** 8 + 0.01  # clustered near zero
+    rng = np.random.default_rng(2)
+    rhs = rng.normal(size=18)
+    np.testing.assert_array_equal(solve_vandermonde(grid, rhs), _solve_exact_reference(grid, rhs))
+
+
+def test_coincident_nodes_raise():
+    """Genuinely singular systems (repeated nodes) raise instead of returning garbage."""
+    with pytest.raises(ValueError, match="pairwise distinct"):
+        solve_vandermonde(np.array([0.1, 0.5, 0.5]), np.ones(3))
+
+
+def test_mismatched_input_shapes_raise():
+    """Mismatched or non-1-D inputs are rejected up front with a clear message."""
+    with pytest.raises(ValueError, match="equal length"):
+        solve_vandermonde(np.array([0.1, 0.3, 0.5]), np.ones(2))
+    with pytest.raises(ValueError, match="1-D"):
+        solve_vandermonde(np.ones((2, 2)), np.ones((2, 2)))
+
+
+def test_non_finite_inputs_raise():
+    """NaN or inf nodes/right-hand sides are rejected with a clear message."""
+    with pytest.raises(ValueError, match="must be finite"):
+        solve_vandermonde(np.array([0.1, np.nan, 0.5]), np.ones(3))
+    with pytest.raises(ValueError, match="must be finite"):
+        solve_vandermonde(np.array([0.1, 0.3, 0.5]), np.array([1.0, np.inf, 1.0]))
+
+
+def test_nearly_coincident_nodes_raise_a_documented_error():
+    """Nodes packed at consecutive float64 ulps push the exact solution beyond the
+    float64 range; the solver raises the documented ValueError instead of an opaque
+    OverflowError."""
+    nodes = [0.5]
+    for _ in range(21):
+        nodes.append(float(np.nextafter(nodes[-1], 1.0)))
+    with pytest.raises(ValueError, match="too close together"):
+        solve_vandermonde(np.array(nodes), np.ones(len(nodes)))
+
+
+def test_nodes_below_the_first_rung_resolution_climb_the_ladder():
+    """Distinct nodes that collapse at 128 fractional bits are not misreported.
+
+    Tiny-magnitude nodes can be genuinely distinct in float64 yet land on the
+    same scaled integer at the first precision rung; the solver must climb to a
+    finer rung and solve exactly rather than report them as coincident.
+    """
+    a = 2.0**-100
+    points = np.array([a, a + 2.0**-152])  # identical at 128 bits, distinct at 256
+    rhs = np.array([1.0, 3.0])
+    np.testing.assert_array_equal(
+        solve_vandermonde(points, rhs), _solve_exact_reference(points, rhs)
+    )
+
+
+def _deep_chain_tree(depth: int, n_features: int | None = None):
+    """A deterministic decision tree of exactly the requested depth.
+
+    Builds a left-leaning chain: every internal node splits on a fresh feature
+    (or cycles through ``n_features`` features when given) and has a leaf on the
+    right. The depth is exact by construction, so the deep-tree tests never
+    silently skip.
+    """
+    n_nodes = 2 * depth + 1
+    children_left = np.full(n_nodes, -1)
+    children_right = np.full(n_nodes, -1)
+    features = np.full(n_nodes, -2)
+    thresholds = np.full(n_nodes, np.nan)
+    values = np.zeros(n_nodes)
+    weights = np.ones(n_nodes)
+    rng = np.random.default_rng(0)
+    # Internal nodes 0..depth-1 form the left chain; node i's right child is the
+    # leaf depth+i, and the last internal node's left child is the leaf 2*depth.
+    for level in range(depth):
+        children_left[level] = level + 1 if level < depth - 1 else 2 * depth
+        children_right[level] = depth + level
+        features[level] = level if n_features is None else level % n_features
+        thresholds[level] = 0.5
+        weights[level] = depth - level + 1  # number of leaves below the node
+    for leaf in range(depth, 2 * depth + 1):
+        values[leaf] = float(rng.normal())
+    return TreeModel(
+        children_left=children_left,
+        children_right=children_right,
+        children_missing=children_left.copy(),
+        features=features,
+        thresholds=thresholds,
+        values=values,
+        node_sample_weight=weights,
+    )
+
+
+@pytest.mark.parametrize(("depth", "tolerance"), [(20, 1e-6), (24, 1e-4), (28, 1e-2)])
+def test_deep_tree_explanations_satisfy_completeness(depth, tolerance):
+    """Deep trees produce values satisfying the completeness axiom.
+
+    The tolerances are the documented downstream float64 bound (``max|N| * 1e-13``):
+    with exact coefficients the only remaining error is the double-precision
+    evaluation of the (exponentially large) N entries, which scales with the
+    largest entry.
+    """
+    tree = _deep_chain_tree(depth)
+    explainer = LinearTreeSHAP(tree)
+    x = np.zeros(depth)  # walks the full left chain through every split
+    explanation = explainer.explain_function(x)
+    prediction = tree.predict_one(x)
+    gap = float(np.sum(explanation.values) - prediction)
+    assert abs(gap) < tolerance * max(1.0, abs(prediction))
+
+
+@pytest.mark.parametrize(("depth", "is_supported"), [(29, True), (30, False)])
+def test_representation_limit_boundary_is_pinned(depth, is_supported):
+    """The largest supported interpolation degree on the default grids is exactly 29.
+
+    This pins the "~29" quoted in the error message, the CHANGELOG, and the
+    docstrings to the actual refuse/accept boundary, so a drift in the grids or
+    the limit cannot silently decouple the prose from the behavior.
+    """
+    tree = _deep_chain_tree(depth)
+    if is_supported:
+        explainer = LinearTreeSHAP(tree)
+        x = np.zeros(depth)
+        gap = float(np.sum(explainer.explain_function(x).values)) - tree.predict_one(x)
+        assert abs(gap) < 1e-2 * max(1.0, abs(tree.predict_one(x)))
+    else:
+        with pytest.raises(RepresentationLimitError):
+            LinearTreeSHAP(tree)
+
+
+@pytest.mark.parametrize(("depth", "is_supported"), [(25, True), (26, False)])
+def test_treeshapiq_representation_limit_boundary_is_pinned(depth, is_supported):
+    """TreeSHAPIQ's boundary is lower than LinearTreeSHAP's 29/30.
+
+    The identity N matrix (rhs of ones, built for every index) crosses the
+    magnitude limit at degree 26 already, so it is the binding constraint for
+    all TreeSHAPIQ indices; this pins the "~25" quoted in the error message.
+    """
+    tree = _deep_chain_tree(depth)
+    if is_supported:
+        TreeSHAPIQ(model=tree, max_order=2, index="SII")
+    else:
+        with pytest.raises(RepresentationLimitError):
+            TreeSHAPIQ(model=tree, max_order=2, index="SII")
+
+
+@pytest.mark.parametrize("depth", [35, 60])
+def test_too_deep_trees_are_refused_with_a_clear_error(depth):
+    """Beyond the float64 representation limit the explainer refuses loudly.
+
+    Depth 35 previously returned silently wrong values (the ~32-59 band where the
+    explicit inverse is rank-deficient to machine precision); depth 60 crashed
+    with an unexplained ``LinAlgError``. The coefficients are now exact at any
+    depth, but their magnitude exceeds what the downstream float64 pipeline can
+    evaluate without destroying the result, so construction raises an
+    explanatory error instead.
+    """
+    tree = _deep_chain_tree(depth)
+    with pytest.raises(RepresentationLimitError, match=r"interpolation degree.*too large"):
+        LinearTreeSHAP(tree)
+
+
+def test_clear_solver_cache_forces_recomputation():
+    """clear_solver_cache drops memoized solutions; a fresh solve still matches the oracle."""
+    D = chebpts2(10)
+    rhs = np.ones(10)
+    solve_vandermonde(D, rhs)
+    assert _solve_vandermonde_cached.cache_info().currsize > 0
+    clear_solver_cache()
+    assert _solve_vandermonde_cached.cache_info().currsize == 0
+    np.testing.assert_array_equal(solve_vandermonde(D, rhs), _solve_exact_reference(D, rhs))
+
+
+def test_solver_results_are_cached_and_isolated():
+    """Identical (points, rhs) solves hit the cache, and callers cannot corrupt it."""
+    D = chebpts2(14)
+    rhs = np.ones(14)
+    first = solve_vandermonde(D, rhs)
+    first[0] = 12345.0  # mutate the returned copy
+    second = solve_vandermonde(D, rhs)
+    assert second[0] != 12345.0  # the cached value must be unaffected
+    np.testing.assert_array_equal(second, _solve_exact_reference(D, rhs))
+
+
+def test_deep_narrow_tree_routes_to_treeshapiq():
+    """A deep tree over few features is re-routed instead of refused.
+
+    LinearTreeSHAP's interpolation degree is the tree depth (above the limit
+    here), but TreeSHAPIQ's is min(depth, features in the tree) = 5, so
+    ``TreeExplainer`` falls back to TreeSHAPIQ and still produces order-1
+    Shapley values satisfying completeness.
+    """
+    depth = 40
+    with pytest.raises(RepresentationLimitError):
+        LinearTreeSHAP(_deep_chain_tree(depth, n_features=5))
+
+    tree = _deep_chain_tree(depth, n_features=5)
+    explainer = TreeExplainer(model=[tree], index="SV", max_order=1)
+    x = np.zeros(depth)
+    explanation = explainer.explain(x)
+    prediction = tree.predict_one(x)
+    gap = abs(float(np.sum(explanation.values)) - prediction)
+    assert gap < 1e-9 * max(1.0, abs(prediction))
+
+
+def test_deep_narrow_tree_is_rerouted_for_sii_as_well():
+    """The re-route covers both order-1 fast-path indices, not just ``"SV"``.
+
+    ``index="SII"`` with ``max_order=1`` also takes the LinearTreeSHAP fast path;
+    the explainer validation normalizes it to ``"SV"`` (with a warning), which is
+    exactly the label the re-route assigns, so the re-routed result satisfies the
+    same completeness.
+    """
+    depth = 40
+    tree = _deep_chain_tree(depth, n_features=5)
+    with pytest.warns(UserWarning, match="SII generalizes 'SV'"):
+        explainer = TreeExplainer(model=[tree], index="SII", max_order=1)
+    assert len(explainer._treeshapiq_explainers) == 1
+    x = np.zeros(depth)
+    explanation = explainer.explain(x)
+    prediction = tree.predict_one(x)
+    gap = abs(float(np.sum(explanation.values)) - prediction)
+    assert gap < 1e-9 * max(1.0, abs(prediction))
+
+
+def test_treeshapiq_n_matrices_are_gated_too():
+    """The representation limit also fires inside TreeSHAPIQ's own N matrices.
+
+    For higher-order indices there is no re-route: a tree whose feature-bounded
+    interpolation degree exceeds the limit is refused at construction instead of
+    returning silently wrong interaction values.
+    """
+    tree = _deep_chain_tree(35)  # 35 distinct features -> TreeSHAPIQ degree 35
+    with pytest.raises(RepresentationLimitError):
+        TreeSHAPIQ(model=tree, max_order=2, index="k-SII")
+
+
+def test_mixed_ensemble_keeps_the_fast_path_for_shallow_trees():
+    """Only the over-deep tree leaves the LinearTreeSHAP fast path in an ensemble.
+
+    The shallow tree stays on LinearTreeSHAP, the deep-narrow tree is re-routed to
+    TreeSHAPIQ, and the mixed aggregation still satisfies completeness for the
+    two-tree ensemble.
+    """
+    assert issubclass(RepresentationLimitError, ValueError)
+    depth = 40
+    deep_tree = _deep_chain_tree(depth, n_features=5)
+    shallow_tree = _deep_chain_tree(8)
+    explainer = TreeExplainer(model=[shallow_tree, deep_tree], index="SV", max_order=1)
+    assert len(explainer._lineartreeshap_explainers) == 1
+    assert len(explainer._treeshapiq_explainers) == 1
+
+    x = np.zeros(depth)
+    explanation = explainer.explain(x)
+    prediction = shallow_tree.predict_one(x) + deep_tree.predict_one(x)
+    gap = abs(float(np.sum(explanation.values)) - prediction)
+    assert gap < 1e-9 * max(1.0, abs(prediction))