def dimension_matrix(variables: Sequence[str], dimensions: Mapping[str, Sequence] = PHYSICAL_DIMENSIONS) -> np.ndarray:
    r"""Dimension matrix $D$ of a list of variables, exactly.

    $$D_{kj} = \text{exponent of base dimension } k \text{ in variable } j$$

    Parameters
    ----------
    variables : sequence of str
        Variable names, keys of ``dimensions`` [-].
    dimensions : Mapping, optional
        Base-dimension exponents $(L, M, T, I)$ per name; default
        :data:`PHYSICAL_DIMENSIONS` [-].

    Returns
    -------
    numpy.ndarray
        Object array of :class:`fractions.Fraction`, one row per base
        dimension of :data:`BASE_DIMENSIONS` and one column per variable [-].

    Raises
    ------
    KeyError
        A variable with no declared dimensions.
    ValueError
        No variable, or a dimension vector of the wrong length.

    Convention
    ----------
    Rows are $(L, M, T, I)$; temperature is an energy, so there is no
    temperature row.

    References
    ----------
    .. [1] E. Buckingham, Phys. Rev. 4 (1914) 345.
    """
    names = list(variables)
    if not names:
        raise ValueError("at least one variable is needed")
    cols = []
    for name in names:
        dim = tuple(dimensions[name])
        if len(dim) != len(BASE_DIMENSIONS):
            raise ValueError(f"{name}: {len(dim)} base-dimension exponents, expected {len(BASE_DIMENSIONS)}")
        cols.append([_frac(v) for v in dim])
    return np.array([[col[k] for col in cols] for k in range(len(BASE_DIMENSIONS))], dtype=object)

def dimensionless_groups(variables: Sequence[str], dimensions: Mapping[str, Sequence] = PHYSICAL_DIMENSIONS) -> list[np.ndarray]:
    r"""A basis of the Buckingham-Pi groups of a set of variables: $\ker D$.

    $$\Pi = \prod_j x_j^{a_j}\ \text{dimensionless} \iff D\mathbf a = 0$$

    Parameters
    ----------
    variables : sequence of str
        Variable names, keys of ``dimensions`` [-].
    dimensions : Mapping, optional
        Base-dimension exponents per name; default :data:`PHYSICAL_DIMENSIONS` [-].

    Returns
    -------
    list of numpy.ndarray
        Integer exponent vectors over ``variables``, one per group; there are
        ``len(variables) - rank(D)`` of them [-].

    Raises
    ------
    KeyError
        A variable with no declared dimensions.

    Convention
    ----------
    The basis is not unique; :func:`express_in_basis` changes it to a
    conventional one.

    References
    ----------
    .. [1] E. Buckingham, Phys. Rev. 4 (1914) 345.
    """
    return rational_null_space(dimension_matrix(variables, dimensions))

def express_in_basis(target: Sequence, basis: Sequence[Sequence]) -> np.ndarray:
    r"""Exact coefficients $\mathbf c$ with $\mathbf t = \sum_k c_k \mathbf b_k$.

    $$\mathbf t = \sum_k c_k\,\mathbf b_k$$

    Parameters
    ----------
    target : sequence
        Exponent vector to express, e.g. a conventional group [-].
    basis : sequence of sequences
        Basis vectors of the same length, e.g. from :func:`dimensionless_groups` [-].

    Returns
    -------
    numpy.ndarray
        Object array of :class:`fractions.Fraction`, one coefficient per basis
        vector [-].

    Raises
    ------
    ValueError
        Lengths that disagree, linearly dependent basis vectors, or a target
        outside their span.

    Convention
    ----------
    In multiplicative form $\Pi_t = \prod_k \Pi_k^{c_k}$: the coefficients
    are the powers of the basis groups that build the target group.

    References
    ----------
    .. [1] G. Strang, *Linear Algebra and Its Applications*, 4th ed.,
           Thomson (2006), Sec. 2.3.
    """
    b = _matrix([list(v) for v in basis])
    t = [_frac(v) for v in target]
    if len(t) != len(b[0]):
        raise ValueError(f"target of length {len(t)} for basis vectors of length {len(b[0])}")
    k = len(b)
    # Solve B^T c = t as an augmented system over the rationals.
    aug = [[b[j][i] for j in range(k)] + [t[i]] for i in range(len(t))]
    reduced, pivots = _rref(aug)
    if len(pivots) and pivots[-1] == k:
        raise ValueError("target is not in the span of the basis")
    if [p for p in pivots if p < k] != list(range(k)):
        raise ValueError("basis vectors are linearly dependent")
    return np.array([reduced[i][k] for i in range(k)], dtype=object)

def monomial_dimension(exponents: Sequence, variables: Sequence[str],
                       dimensions: Mapping[str, Sequence] = PHYSICAL_DIMENSIONS) -> np.ndarray:
    r"""Base dimensions of the monomial $\prod_j x_j^{a_j}$: $D\mathbf a$.

    $$\dim\Big(\prod_j x_j^{a_j}\Big) = D\mathbf a$$

    Parameters
    ----------
    exponents : sequence
        Exponent $a_j$ per variable [-].
    variables : sequence of str
        Variable names, keys of ``dimensions`` [-].
    dimensions : Mapping, optional
        Base-dimension exponents per name; default :data:`PHYSICAL_DIMENSIONS` [-].

    Returns
    -------
    numpy.ndarray
        Object array of :class:`fractions.Fraction`, one entry per base
        dimension; all zero for a dimensionless monomial [-].

    Raises
    ------
    ValueError
        ``exponents`` and ``variables`` of different lengths.

    Convention
    ----------
    Base dimensions in the order of :data:`BASE_DIMENSIONS`.

    References
    ----------
    .. [1] E. Buckingham, Phys. Rev. 4 (1914) 345.
    """
    a = [_frac(v) for v in exponents]
    names = list(variables)
    if len(a) != len(names):
        raise ValueError(f"{len(a)} exponents for {len(names)} variables")
    d = dimension_matrix(names, dimensions)
    return np.array([sum((d[k][j] * a[j] for j in range(len(a))), Fraction(0)) for k in range(d.shape[0])],
                    dtype=object)

def rational_null_space(matrix) -> list[np.ndarray]:
    r"""Exact basis of the null space $\{\mathbf a : M\mathbf a = 0\}$ of a rational matrix.

    $$M\mathbf a = 0$$

    Parameters
    ----------
    matrix : array-like
        Rows of integers, fractions or floats; a float is read as the decimal
        it prints as [-].

    Returns
    -------
    list of numpy.ndarray
        One object array of integer-valued :class:`fractions.Fraction` per
        basis vector, coprime and with a positive first non-zero entry; empty
        when $M$ has full column rank [-].

    Raises
    ------
    ValueError
        An empty or ragged matrix, or a non-finite entry.

    Convention
    ----------
    The basis is the one Gauss-Jordan elimination gives (one vector per free
    column, that column set to one); any invertible combination of it spans
    the same space.

    References
    ----------
    .. [1] G. Strang, *Linear Algebra and Its Applications*, 4th ed.,
           Thomson (2006), Sec. 2.2.
    """
    m = _matrix(np.asarray(matrix, dtype=object).tolist())
    ncol = len(m[0])
    reduced, pivots = _rref(m)
    basis = []
    for free in (c for c in range(ncol) if c not in pivots):
        v = [Fraction(0)] * ncol
        v[free] = Fraction(1)
        for row, pc in zip(reduced, pivots):
            v[pc] = -row[free]
        basis.append(_integer_scaled(v))
    return basis

def similarity_constraint_vectors(kept_groups: Sequence[Sequence]) -> list[np.ndarray]:
    r"""Linear constraints a law depending only on ``kept_groups`` imposes on its exponents.

    $$\mathbf e \in \mathrm{span}\{\mathbf g_k\} \iff \mathbf c\cdot\mathbf e = 0\ \ \forall\,\mathbf c \perp \mathbf g_k$$

    Parameters
    ----------
    kept_groups : sequence of sequences
        Exponent vectors, over the varied variables, of the dimensionless
        groups the similarity assumption keeps [-].

    Returns
    -------
    list of numpy.ndarray
        Integer constraint vectors $\mathbf c$, a basis of the orthogonal
        complement of the groups' span; empty when the groups span the whole
        space (no constraint) [-].

    Raises
    ------
    ValueError
        An empty or ragged group list.

    Convention
    ----------
    $\mathbf e$ is the exponent vector of the dimensionless response (e.g.
    $\Omega_i\tau_E$) over the *varied* variables only; the constants and the
    variables held fixed (shape, $q$, ion mass) are not coordinates.

    References
    ----------
    .. [1] J. W. Connor and J. B. Taylor, Nucl. Fusion 17 (1977) 1047.
    .. [2] B. B. Kadomtsev, Sov. J. Plasma Phys. 1 (1975) 295.
    """
    return rational_null_space([list(g) for g in kept_groups])

def state_exponents_from_engineering_exponents(a_I, a_B, a_P, a_n=0, a_R=0) -> np.ndarray:
    r"""Exponents of $\tau_E$ over $(n, T, B, R)$ implied by an engineering scaling, exactly.

    $$\tau_E^{1+\alpha_P} \propto n^{\alpha_n+\alpha_P}\,T^{\alpha_P}\,
      B^{\alpha_I+\alpha_B}\,R^{\alpha_I+3\alpha_P+\alpha_R}$$

    Parameters
    ----------
    a_I : float or Fraction
        Plasma-current exponent $\alpha_I$ [-].
    a_B : float or Fraction
        Toroidal-field exponent $\alpha_B$ [-].
    a_P : float or Fraction
        Loss-power exponent $\alpha_P$ [-].
    a_n : float or Fraction, optional
        Density exponent $\alpha_n$, default 0 [-].
    a_R : float or Fraction, optional
        Size exponent $\alpha_R$, default 0 [-].

    Returns
    -------
    numpy.ndarray
        The exponents of $n$, $T$, $B$ and $R$ in $\tau_E$: exact
        :class:`fractions.Fraction` when every input is an int, Fraction or
        rational string, float otherwise [-].

    Raises
    ------
    ValueError
        $\alpha_P = -1$, where $\tau_E$ cancels and the map is singular.

    Convention
    ----------
    At fixed $\epsilon$, $\kappa$, $q$ and ion mass: $I_p \propto B R$
    (cylindrical $q$) and $P = W/\tau_E$ with $W \propto n T R^3$.  Every
    exponent carries the factor $1/(1+\alpha_P)$, the conditioning of the
    whole engineering-to-dimensionless map.

    References
    ----------
    .. [1] J. W. Connor and J. B. Taylor, Nucl. Fusion 17 (1977) 1047.
    .. [2] ITER Physics Expert Groups, Nucl. Fusion 39 (1999) 2175, Ch. 2, Sec. 6.
    """
    values = (a_I, a_B, a_P, a_n, a_R)
    exact = all(isinstance(v, (int, np.integer, Fraction, str)) for v in values)
    if exact:
        i, b, p, n, r = (_frac(v) for v in values)
    else:
        i, b, p, n, r = (float(v) for v in values)
        if not np.all(np.isfinite([i, b, p, n, r])):
            raise ValueError("non-finite exponent")
    d = 1 + p
    if d == 0:
        raise ValueError("alpha_P = -1: tau_E cancels from P = W/tau_E and the map is singular")
    out = [(n + p) / d, p / d, (i + b) / d, (i + 3 * p + r) / d]
    return np.array(out, dtype=object if exact else float)