\section{Exact optimum of the nonscalar four-level model}
\label{app:nonscalaroptimal}

We now adapt the scalar construction to the nonscalar average used in the
main text. Throughout, \(k=1000001/1000000\) is fixed as in
Eq.~\eqref{eq:nonscalar}; no claim is made for an arbitrary average state.
Write \(s_1=s_2=p,s_3=q,s_4=t\), and let the three upper entries of
\(B_k\) be \(w_0=p,m_0=2p,n_0=3(q+t)\).
The finite and square-prior models have the same \(S_k,B_k,T_{0,k}\);
only \(\lambda\) changes. We retain the dimensionless identity loss.

\subsection{A lower comparison at every estimate}

The factorization must now account for the unequal diagonal entries of
\(S_k\). We choose the following matrix so that its weighted row products
vanish off the diagonal.
For real \(u,v\), put
\begin{align}
 d_k&=s_1s_2v(n_0+v)-s_3s_4u(w_0+u),\nonumber\\
 a&=\frac{u m_0s_1s_4(n_0+v)}{d_k},\qquad
 b=\frac{u m_0s_3s_4(w_0+u)}{d_k},\nonumber\\
 R(u,v)&=\begin{pmatrix}
 0&u&0&a\\
 -w_0-u&0&b&0\\
 0&-m_0-b&0&v\\
 -a&0&-n_0-v&0
 \end{pmatrix}.
 \label{eq:weightedgram}
\end{align}
Whenever \(d_k\ne0\), one has \(R+R^T=-(B_k+B_k^\dagger)\).
Only the \(13\) and \(24\) entries of \(RS_k^{-1}R^T\) can be
off-diagonal. Substitution gives
\[
 -\frac{u(m_0+b)}{s_2}+\frac{av}{s_4}=0,\qquad
 \frac{a(w_0+u)}{s_1}-\frac{b(n_0+v)}{s_3}=0.
\]
Define \(c_i=(RS_k^{-1}R^T)_{ii}\) and
\(A_k(u,v)=T_{0,k}-\operatorname{diag}(c_i)\).
Expanding the weighted Gram product proves
Eq.~\eqref{eq:nonscalaroptimalfactor} with \(R(u,v)\) in place of
\(R_\diamond\), for all \(z=rw\), including \(z=0\).
The phase identity is unchanged because \(S_k,T_{0,k}\) and the diagonal
Gram matrix commute with \(D(w)\).
Consequently every normalized POVM and measurable estimator has risk at least
\(\lambda-\sum_i c_i\). The Radon--Nikodym and full-rank integrability
argument above applies verbatim with \(S_k\) in place of \(I/4\).

\subsection{The allowed reports and effect normalization}

Set \(H=RS_k^{-1}\), and reorder its indices as \(1,3,2,4\), so that
\(H=\left(\begin{smallmatrix}0&K\\L&0\end{smallmatrix}\right)\). Define
\[
 \sigma=\Tr(KL),\qquad \kappa=\det(KL),\qquad
 h(x)=x^2-\sigma x+\kappa .
\]
Block elimination, or the characteristic polynomial of \(H\), gives
\[
 \det(rS_k+R)=\det(S_k)h(r^2).
\]
If \(\sigma>0,\kappa>0,\sigma^2>4\kappa\), the two roots
\(x_-<x_+\) of \(h\) are positive. Put \(r_s=\sqrt{x_s}\).
The real tridiagonal matrix \(N(r)=F_k(r)-A_k(u,v)\) has diagonal
\(s_i r^2+c_i\) and neighboring entries \(-r(w_0,m_0,n_0)\).
At \(x=x_s\) define
\begin{align*}
 \nu_1&=1,&
 \nu_2&=\frac{s_1x+c_1}{r_s w_0},\\
 \nu_3&=\frac{(s_2x+c_2)(s_1x+c_1)-xw_0^2}{xw_0m_0},&
 \nu_4&=\frac{(s_3x+c_3)\nu_3-m_0(s_1x+c_1)/w_0}{r_s n_0}.
\end{align*}
The first three row equations fix these components.
The determinant identity
\(\det N(r)=\det(S_k)h(r^2)^2\), together with the nonzero neighboring
entries, then forces the last row equation. Equivalently,
\((r_sS_k+R)^T\nu_s=0\): these are left-null vectors of the Gram factor,
not right-null vectors.

Reduce each square modulo \(h\), writing
\(\nu_i(x)^2=\eta_i+\xi_i x\) at its roots, with
\(x^{-1}=(\sigma-x)/\kappa\). In particular,
\[
 \eta_2=\frac{2s_1c_1+c_1^2\sigma/\kappa}{w_0^2},
 \qquad
 \xi_2=\frac{s_1^2-c_1^2/\kappa}{w_0^2}.
\]
For completeness, the other rational coefficients can be evaluated without
square roots. Put
\begin{align*}
 A_3&=\frac{s_1c_2+s_2c_1-w_0^2+c_1c_2\sigma/\kappa}{w_0m_0},&
 B_3&=\frac{s_1s_2-c_1c_2/\kappa}{w_0m_0},\\
 A_4&=c_3A_3-s_3B_3\kappa-m_0c_1/w_0,&
 B_4&=s_3A_3+c_3B_3+s_3B_3\sigma-m_0s_1/w_0 .
\end{align*}
Then
\begin{align*}
 \eta_3&=A_3^2-B_3^2\kappa,&
 \xi_3&=2A_3B_3+B_3^2\sigma,\\
 \eta_4&=\frac{(A_4^2-B_4^2\kappa)\sigma/\kappa
                    +2A_4B_4+B_4^2\sigma}{n_0^2},&
 \xi_4&=-\frac{A_4^2-B_4^2\kappa}{\kappa n_0^2}.
\end{align*}
Where \(\xi_2\ne0\), let \(y=(1-\eta_2)/\xi_2\).
The required branch is a root of the two rational equations
\begin{equation}
 g_3(u,v)=\eta_3+\xi_3y-1=0,\qquad
 g_4(u,v)=\eta_4+\xi_4y-1=0.
 \label{eq:nonscalarnormalization}
\end{equation}
If also \(h(y)<0\), the weights
\[
 \alpha_+=\frac{y-x_-}{x_+-x_-},\qquad \alpha_-=1-\alpha_+
\]
are strictly positive. They satisfy
\(\alpha_-\nu_{-,i}^2+\alpha_+\nu_{+,i}^2=1\) for every \(i\);
the first component uses \(\alpha_-+\alpha_+=1\), the second the definition
of \(y\), and the remaining two Eq.~\eqref{eq:nonscalarnormalization}.
Thus the eight rank-one effects
\[
 M_{s,w}=\frac{\alpha_s}{4}D(w)\nu_s\nu_s^TD(w)^\dagger,\qquad
 z_{s,w}=r_s w,\qquad w\in\{1,i,-1,-i\},
\]
are positive and sum to \(I_4\). Each annihilates its comparison residual.
Their risk is therefore exactly \(\lambda-\sum_i c_i\), matching the
unrestricted lower bound. Covariance constructs one attaining measurement;
it does not restrict the lower-bound quantifier.

\subsection{Specifying the algebraic solution and bounding its risk}

Let the rational center be the terminating decimals
\[
 \begin{split}
 b_u&=-0.06554647804483970189007943245440558,\\
 b_v&=-0.75260434123165148908923759763752985 .
 \end{split}
\]
In Eq.~\eqref{eq:weightedgram} select the unique root
\((u_\diamond,v_\diamond)\) of Eq.~\eqref{eq:nonscalarnormalization}
in the closed box \(\mathcal B=b+[-10^{-30},10^{-30}]^2\).
This box, not the rounded approximations to a measurement, defines the branch.
We now justify existence and uniqueness.

With \(g=(g_3,g_4)^T\), calculate the exact rational matrix
\(C_0=Jg(b)^{-1}\). Rational interval evaluation on the whole box certifies
that all evaluated denominators exclude zero and that
\[
 \|I-C_0Jg(\mathcal B)\|_\infty<2\,10^{-23}<1,\qquad
 b-C_0g(b)+(I-C_0Jg(\mathcal B))(\mathcal B-b)
       \subset\operatorname{int}\mathcal B .
\]
The norm bound is the maximum absolute row-sum enclosure.
The map \(T(x)=x-C_0g(x)\) is consequently a contraction of the closed
box into itself. The mean-value formula proves the displayed image inclusion
for every point, so the contraction theorem gives exactly one fixed point.
Since \(C_0\) is invertible, that point solves \(g=0\).
The same bound makes \(Jg\) nonsingular on the box.
Clearing its nonzero rational-function denominators produces polynomial
equations over \(\mathbb Q\) with a nonsingular isolated solution.
Thus the selected coordinates, radii, weights and risk are algebraic.
No assertion of minimal polynomial degree is needed.

The whole-box certificate also verifies
\(\sigma>0,\kappa>0,\sigma^2-4\kappa>0,\xi_2\ne0,h(y)<0\), and the strict
risk bounds. With \(R_\diamond=R(u_\diamond,v_\diamond)\), they give
\[
 \begin{split}
 94.28592480080068802537609221
  &<\mathcal L_{\min,k}<94.28592480080068802537609222,\\
 0.03592633552227579562708351
  &<\Gamma_k<0.03592633552227579562708352 .
 \end{split}
\]
The first line uses \(\lambda=100\). For the square prior, \(\lambda=600\)
and the minimum risk is larger by \(500\); \(\Gamma_k\) is unchanged.
In both cases the exact definition is
\(\mathcal L_{\min,k}=\lambda-\sum_i c_i\), with the branch just specified.

The ancillary certificate \texttt{nonscalar-root-certificate.json} contains
only \(k\), the rational center and box radius.
\texttt{weighted\_exact\_checker.py} uses exact fractions, forward
differentiation and outward-rounded \(180\)-bit dyadic intervals.
\texttt{independent\_symbolic\_checker.py} instead uses the explicit
component formulas above, symbolic differentiation and outward-rounded
decimal intervals. Each independently verifies contraction, strict
self-inclusion, positive radii and mixture weights, and the risk enclosure.
Neither treats a numerical residual as proof or imports the numerical
proposal generator. A wrong model parameter is rejected by both.