\section{Reference implementations}\label{app:algorithms}

The replication package (\cref{app:replication}) contains the package \texttt{pamle} with the following components, all exercised by the test suite.
\begin{itemize}
\item \texttt{jackknife.py}: exact rational Vandermonde weights for any pattern, the constants $\betar$ and $\Amp$, condition numbers, and the combination of level averages.
\item \texttt{likelihoods.py}: exponential, logistic and Poisson models with batched Newton fitting over blocks and sub-blocks (step halving, convergence on the score, status codes for non-convergence, singular curvature and divergence, and a separation safeguard for logistic regression that refuses a stationary point with fitted $\abs{\eta} > 30$ or a numerically singular information matrix). Derivatives are checked against finite differences.
\item \texttt{algorithms.py}: local split fits, corrected worker estimates, size-weighted pilots, global Newton corrections with exact, stale, single-machine and quantized curvature (linear solves, condition numbers and step norms as diagnostics; no explicit inverse in the Newton step), the aggregated estimating-equation estimator of \citet{lin2011aggregated}, the surrogate-likelihood estimator of \citet{jordan2019communication} with its full local maximization, and the analytic first-order bias correction of canonical GLMs \citep{cox1968general,cordeiro1991bias}.
\item \texttt{communication.py}: rounds, messages, numbers and bytes for every variant of \cref{tab:protocols}.
\item \texttt{protocol.py}, \texttt{worker.py}: the coordinator--worker protocol over TCP with length-prefixed messages, per-stage wall-clock, per-worker compute times, latency percentiles and byte counts (\cref{alg:protocol}).
\end{itemize}

\paragraph{Analytic bias correction used as a baseline.} For a canonical-link GLM with per-observation score $\psi_i(\beta) = x_i(y_i - \mu_i)$, the $O(1/n)$ bias of the MLE conditional on the design is $b_1/n$ with
\[
b_1 = \tfrac12\, I^{-1}\Big[n^{-1}\textstyle\sum_i \nabla^2\psi_i\,[I^{-1}]\Big] = -\tfrac12\, I^{-1}\, n^{-1}\textstyle\sum_i x_i\,\mu''(\eta_i)\,x_i^\top I^{-1}x_i,
\]
where $I = n^{-1}\sum_i x_ix_i^\top\mu'(\eta_i)$; this is the form of \citet{cordeiro1991bias} and reduces to $b_1 = \theta_0$ in the exponential model. The corrected local estimate $\thetahat_{D_k} - \hat b_1/m$ needs one local fit and no sub-blocks; it is the local analogue of the debiased weighted estimator of \citet{gu2023heterogeneity} and provides, in \cref{sec:mc}, a baseline that attains the order-$2$ exponent at one third of the local fits of the half-block jackknife whenever the model's bias formula is available and correct.