def budgeted_torque_extrema(Q: np.ndarray, D: np.ndarray, budget: float) -> dict[str, object]:
    """Largest and smallest signed torque over free complex excitations with ``c^H D c = budget``.

    Parameters
    ----------
    Q : numpy.ndarray
        Hermitian torque matrix of a qualified model, possibly indefinite [N m / phasor^2].
    D : numpy.ndarray
        Hermitian positive-definite excitation metric [-].
    budget : float
        The excitation budget ``P_0 > 0`` [phasor^2].

    Returns
    -------
    extrema : dict
        ``T_max``/``T_min`` [N m], ``c_max``/``c_min`` meeting the budget
        [phasor], and the generalized ``eigenvalues`` ascending [N m / phasor^2].

    Raises
    ------
    ValueError
        ``D`` is not positive definite, shapes disagree, or the budget is not
        positive.

    Processing steps
    ----------------
    1. Solve ``Q v = lambda D v`` (generalized Hermitian eigenproblem).
    2. Scale the extreme eigenvectors to ``v^H D v = P_0``; ``T = lambda P_0``.

    Convention
    ----------
    Signed: with an indefinite ``Q`` the minimum is a negative torque, not a
    small magnitude.  The returned phasors are defined up to a common phase.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    A mathematical optimum of the phasor space: actual sector or circuit
    current limits, winding signs and fixed group amplitudes are not imposed
    and must be checked separately.

    Provenance
    ----------
    .. [1887] Issue #1887, Sec. 8 (free complex excitation with a budget).
    """
    from scipy.linalg import eigh

    q = _hermitian(Q)
    d = _hermitian(D)
    if d.shape != q.shape:
        raise ValueError("Q and D differ in shape")
    if not budget > 0:
        raise ValueError("the budget must be positive")
    if np.linalg.eigvalsh(d).min() <= 0:
        raise ValueError("D must be positive definite")
    values, vectors = eigh(q, d)

    def scaled(v):
        return v * np.sqrt(budget / float(np.real(v.conj() @ d @ v)))

    return {
        "T_max": float(values[-1] * budget),
        "T_min": float(values[0] * budget),
        "c_max": scaled(vectors[:, -1]),
        "c_min": scaled(vectors[:, 0]),
        "eigenvalues": values,
    }

def coil_phasors(amplitudes: Sequence[float], phases: Sequence[float]) -> np.ndarray:
    """Complex excitation phasors ``A_k exp(+i phi_k)`` of the coil groups.

    Parameters
    ----------
    amplitudes : sequence of float
        Declared amplitude normalization of each group, not negative [phasor].
    phases : sequence of float
        Phase of each group [rad].

    Returns
    -------
    c : numpy.ndarray
        Complex phasor per group [phasor].

    Raises
    ------
    ValueError
        The two sequences differ in length or an amplitude is negative or not
        finite.

    Convention
    ----------
    ``c_k = A_k exp(+i phi_k)``; the sign of the exponent is the one
    :func:`hermitian_torque_matrix` recovers cross terms in.

    Applicability
    -------------
    Machine-independent.

    Provenance
    ----------
    .. [1887] Issue #1887, phasor convention.
    """
    a = np.asarray(amplitudes, dtype=float).ravel()
    phi = np.asarray(phases, dtype=float).ravel()
    if a.shape != phi.shape:
        raise ValueError(f"{a.size} amplitudes for {phi.size} phases")
    if not np.all(np.isfinite(a)) or np.any(a < 0):
        raise ValueError("amplitudes must be finite and not negative")
    return a * np.exp(1j * phi)

def hermitian_torque_matrix(torques: Mapping[str, float], groups: int) -> np.ndarray:
    """The Hermitian ``Q`` of ``T = c^H Q c`` from the probe torques.

    Parameters
    ----------
    torques : mapping of str to float
        Signed torque of each probe of :func:`quadratic_probe_design`, keyed
        by its label [N m].
    groups : int
        Number of coil groups ``K`` [-].

    Returns
    -------
    Q : numpy.ndarray
        ``K x K`` Hermitian matrix [N m / phasor^2].

    Raises
    ------
    KeyError
        A probe label is missing.
    ValueError
        A probe torque is not finite (a failed evaluation cannot identify ``Q``).

    Processing steps
    ----------------
    1. ``Q_ii = T(e_i)``.
    2. ``Re Q_ij = [T(e_i + e_j) - T(e_i) - T(e_j)] / 2``.
    3. ``Im Q_ij = [T(e_i) + T(e_j) - T(e_i + i e_j)] / 2``.
    4. ``Q_ji = conj(Q_ij)``.

    Convention
    ----------
    With ``c = e_i + i e_j``, ``c^H Q c = Q_ii + Q_jj - 2 Im Q_ij``, which is
    step 3; it holds for the ``exp(+i phi)`` phasor of :func:`coil_phasors`.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    A ``Q`` that reproduces its own probes exactly is not thereby validated:
    the probes determine it.  Only :func:`qualify_torque_matrix` on held-out
    evaluations can qualify it.

    Provenance
    ----------
    .. [1887] Issue #1887, Sec. 6.
    """
    k = int(groups)
    labels = [label for label, _ in quadratic_probe_design(k)]
    values = {label: float(torques[label]) for label in labels}
    bad = sorted(label for label, value in values.items() if not np.isfinite(value))
    if bad:
        raise ValueError(f"probe torques are not finite: {bad}")
    q = np.zeros((k, k), dtype=complex)
    for i in range(k):
        q[i, i] = values[f"e{i}"]
    for i, j in combinations(range(k), 2):
        real = 0.5 * (values[f"e{i}+e{j}"] - values[f"e{i}"] - values[f"e{j}"])
        imag = 0.5 * (values[f"e{i}"] + values[f"e{j}"] - values[f"e{i}+ie{j}"])
        q[i, j] = real + 1j * imag
        q[j, i] = real - 1j * imag
    return q

def quadratic_invariance_checks(
    base_torque: float,
    *,
    scaled: Sequence[tuple[float, float]] = (),
    reversed_torque: float | None = None,
    rotated: Sequence[tuple[float, float]] = (),
    atol: float,
    rtol: float,
) -> tuple[QuadraticTorqueCheck, ...]:
    """Tests of a pure quadratic response that need no ``Q``: scaling, sign, common phase.

    Parameters
    ----------
    base_torque : float
        Signed torque of the reference excitation ``c`` [N m].
    scaled : sequence of (float, float)
        ``(a, T(a c))`` for positive real factors ``a`` [-, N m].
    reversed_torque : float, optional
        ``T(-c)`` [N m].
    rotated : sequence of (float, float)
        ``(chi, T(exp(i chi) c))`` for common phase rotations [rad, N m].
    atol : float
        Absolute tolerance, the scale below which a torque is near zero [N m].
    rtol : float
        Relative tolerance away from zero [-].

    Returns
    -------
    checks : tuple of QuadraticTorqueCheck
        ``amplitude_scaling`` (expects ``a^2 T(c)``), ``sign_reversal``
        (expects ``T(c)``) and ``common_phase`` (expects ``T(c)``), each only
        when its samples were given [-].

    Raises
    ------
    ValueError
        A tolerance is negative, a factor ``a`` is not positive, or a torque
        is not finite.

    Processing steps
    ----------------
    1. Expected torque per sample from the homogeneous-quadratic hypothesis.
    2. A sample passes when ``|T - T_expected| <= atol + rtol |T_expected|``.
    3. One check per kind, passed only if every sample of that kind passes.

    Convention
    ----------
    The common-phase test assumes a single toroidal mode ``n``; a mixed-``n``
    excitation is not phase-invariant and should not be given here.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    A failed check rejects the *pure* quadratic form; a constant or linear
    term (a background field) is examined by :func:`fit_inhomogeneous_torque`.
    The tolerances are the caller's: no universal NTV tolerance exists.

    Provenance
    ----------
    .. [1887] Issue #1887, Sec. 5.
    """
    _tolerances(atol, rtol)
    base = _finite(base_torque, "base_torque")
    checks: list[QuadraticTorqueCheck] = []
    if scaled:
        factors = np.array([float(a) for a, _ in scaled])
        if np.any(factors <= 0):
            raise ValueError("scaling factors must be positive")
        values = np.array([_finite(t, "scaled torque") for _, t in scaled])
        checks.append(_check("amplitude_scaling", values, factors**2 * base, atol, rtol))
    if reversed_torque is not None:
        checks.append(_check("sign_reversal", np.array([_finite(reversed_torque, "reversed_torque")]),
                             np.array([base]), atol, rtol))
    if rotated:
        values = np.array([_finite(t, "rotated torque") for _, t in rotated])
        checks.append(_check("common_phase", values, np.full(values.shape, base), atol, rtol))
    return tuple(checks)

def quadratic_probe_design(groups: int) -> tuple[tuple[str, np.ndarray], ...]:
    """The ``K^2`` unit excitations that determine a Hermitian ``Q``.

    Parameters
    ----------
    groups : int
        Number of independently driven coil groups ``K``, at least 1 [-].

    Returns
    -------
    probes : tuple of (str, numpy.ndarray)
        ``("e0", e_0)``, ..., then for every ``i < j`` ``("e{i}+e{j}", e_i + e_j)``
        and ``("e{i}+ie{j}", e_i + i e_j)``: ``K + 2 C(K, 2) = K^2`` phasors [phasor].

    Raises
    ------
    ValueError
        ``groups`` is below 1.

    Convention
    ----------
    The mixed probes put ``+i`` on the higher-index group, matching the
    ``exp(+i phi)`` phasor of :func:`coil_phasors`.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    ``K^2`` evaluations identify ``Q`` only; validating it needs further
    held-out evaluations (:func:`qualify_torque_matrix`).

    Provenance
    ----------
    .. [1887] Issue #1887, Sec. 6 (probe identities and their count).
    """
    k = int(groups)
    if k < 1:
        raise ValueError("at least one coil group is needed")
    unit = np.eye(k, dtype=complex)
    probes = [(f"e{i}", unit[i]) for i in range(k)]
    for i, j in combinations(range(k), 2):
        probes.append((f"e{i}+e{j}", unit[i] + unit[j]))
        probes.append((f"e{i}+ie{j}", unit[i] + 1j * unit[j]))
    return tuple(probes)

def quadratic_torque(Q: np.ndarray, phasors: np.ndarray) -> np.ndarray:
    """Signed torque ``c^H Q c`` of one excitation or a stack of them.

    Parameters
    ----------
    Q : numpy.ndarray
        ``K x K`` Hermitian torque matrix [N m / phasor^2].
    phasors : numpy.ndarray
        One phasor vector of length ``K``, or an ``(N, K)`` stack [phasor].

    Returns
    -------
    torque : numpy.ndarray
        Real signed torque, scalar or length ``N`` [N m].

    Raises
    ------
    ValueError
        ``Q`` is not square and Hermitian, or the phasors do not match it.

    Convention
    ----------
    ``T = c^H Q c`` with ``c_k = A_k exp(+i phi_k)``; the imaginary part,
    zero for a Hermitian ``Q``, is dropped after the check.

    Applicability
    -------------
    Machine-independent.
    """
    q = _hermitian(Q)
    c = np.asarray(phasors, dtype=complex)
    if c.shape[-1] != q.shape[0]:
        raise ValueError(f"phasors of length {c.shape[-1]} for a {q.shape[0]}-group Q")
    return np.real(np.einsum("...i,ij,...j->...", c.conj(), q, c))

def two_group_phase_extrema(Q: np.ndarray, amplitudes: Sequence[float]) -> dict[str, object]:
    """Relative phases of the largest and smallest signed torque of two groups.

    Parameters
    ----------
    Q : numpy.ndarray
        ``2 x 2`` Hermitian torque matrix of a qualified model [N m / phasor^2].
    amplitudes : sequence of float
        Fixed amplitudes ``(A_1, A_2)``, not negative [phasor].

    Returns
    -------
    extrema : dict
        ``delta_phi_max`` and ``delta_phi_min`` in ``[0, 2 pi)`` [rad], or
        ``None`` when the torque does not depend on the phase; ``T_max`` and
        ``T_min`` [N m]; ``phase_independent`` [-].

    Raises
    ------
    ValueError
        ``Q`` is not ``2 x 2`` Hermitian or an amplitude is negative.

    Processing steps
    ----------------
    1. ``T(dphi) = Q_11 A_1^2 + Q_22 A_2^2 + 2 A_1 A_2 |Q_12| cos(dphi + arg Q_12)``.
    2. Maximum at ``dphi = -arg Q_12``, minimum at ``pi - arg Q_12``.
    3. With ``Q_12 = 0`` or a zero amplitude no phase is preferred.

    Convention
    ----------
    ``c = (A_1, A_2 exp(+i dphi))``; the extrema are of the *signed* torque,
    not of ``|T|`` or of a braking criterion.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    Valid only where :func:`qualify_torque_matrix` qualified ``Q``.

    Provenance
    ----------
    .. [1887] Issue #1887, Sec. 8.
    """
    q = _hermitian(Q)
    if q.shape != (2, 2):
        raise ValueError("two_group_phase_extrema needs a 2 x 2 Q")
    a1, a2 = (float(a) for a in amplitudes)
    if a1 < 0 or a2 < 0:
        raise ValueError("amplitudes must not be negative")
    base = float(np.real(q[0, 0])) * a1**2 + float(np.real(q[1, 1])) * a2**2
    swing = 2 * a1 * a2 * abs(q[0, 1])
    scale = max(abs(base), abs(q[0, 0]) * a1**2, abs(q[1, 1]) * a2**2, np.finfo(float).tiny)
    if swing <= 64 * np.finfo(float).eps * scale:  # below the roundoff of the diagonal terms
        return {"delta_phi_max": None, "delta_phi_min": None, "T_max": base, "T_min": base,
                "phase_independent": True}
    arg = float(np.angle(q[0, 1]))
    return {
        "delta_phi_max": float(np.mod(-arg, 2 * np.pi)),
        "delta_phi_min": float(np.mod(np.pi - arg, 2 * np.pi)),
        "T_max": base + swing,
        "T_min": base - swing,
        "phase_independent": False,
    }