perf(gaussian): drop the sqrt round-trip from variance-space operations

`Add`, `Sub`, `exclude` and `forget` combined variances by way of standard
deviations: `sigma()` takes a square root, `.powi(2)` squares it away,
`var.sqrt()` takes another, and `from_ms` squares that one back. Three roots
to compute a value that is `1/pi` all along.

They now go through `variance()` and a new `from_mv(mu, var)`, which skip
both conversions. `Sub` is the hot one — `RankDiffFactor::propagate` is
`a - b`, run for every adjacent team pair on every forward and backward
sweep of every EP iteration.

`run_chain` also stopped recomputing each team's weighted performance in the
likelihood loop; the fold is already in `arena.team_prior`, indexed by the
sorted position the loop has in hand. Each `performance()` is itself a
`forget`, so the duplicate cost scaled with players per team.

Measured on this machine, before and after, same fixtures:

    Batch::iteration          23.57us -> 19.31us   (-18%)
    scored_history_60_events   1.071ms -> 983us    (-8%)

The `Gaussian::add`/`sub` microbenchmarks cannot resolve the change: they
sit at ~234ps against a ~218ps floor that `mul`/`div` also hit, so the
harness overhead dominates a single operation.

One golden moved. Two identical competitors drawing must land on their
shared prior mean exactly, by symmetry; the root-free path now returns
25.0 where the reference transcription recorded 24.999999 — that value
rounded to six decimals. Asserting a six-decimal transcription at
epsilon 1e-6 left no headroom, so the expectation is now the exact value.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01DnsaJg74eNSva3PJjK2eej

src/gaussian.rs

@@ -35,6 +35,28 @@ impl Gaussian {
         }
     }
 
+    /// Construct from mean and *variance*, skipping the square-root round trip.
+    ///
+    /// `from_ms(mu, var.sqrt())` immediately squares the root away again to
+    /// recover `pi = 1/var`. Variance-combining operations (`Add`, `Sub`,
+    /// `exclude`, `forget`) work in variance space throughout, so they go
+    /// through here instead and never take a root.
+    #[inline]
+    pub(crate) fn from_mv(mu: f64, var: f64) -> Self {
+        if var == f64::INFINITY {
+            Self { pi: 0.0, tau: 0.0 }
+        } else if var == 0.0 {
+            // Point mass at mu; see `from_ms` for the tau convention.
+            Self {
+                pi: f64::INFINITY,
+                tau: if mu == 0.0 { 0.0 } else { f64::INFINITY },
+            }
+        } else {
+            let pi = 1.0 / var;
+            Self { pi, tau: mu * pi }
+        }
+    }
+
     /// Construct directly from natural parameters.
     #[inline]
     pub(crate) const fn from_natural(pi: f64, tau: f64) -> Self {
@@ -64,6 +86,21 @@ impl Gaussian {
         }
     }
 
+    /// Variance, `1 / pi`, without the root-and-square of `sigma().powi(2)`.
+    ///
+    /// Mirrors `sigma()`'s treatment of the improper (`pi <= 0`) and point-mass
+    /// (`pi == inf`) cases.
+    #[inline]
+    pub(crate) fn variance(&self) -> f64 {
+        if self.pi <= 0.0 {
+            f64::INFINITY
+        } else if self.pi.is_infinite() {
+            0.0
+        } else {
+            1.0 / self.pi
+        }
+    }
+
     #[inline]
     pub fn sigma(&self) -> f64 {
         // A non-positive precision is improper → infinite standard deviation. Guarding
@@ -86,22 +123,21 @@ impl Gaussian {
     }
 
     pub(crate) fn exclude(&self, other: Gaussian) -> Self {
-        let var = self.sigma().powi(2) - other.sigma().powi(2);
+        let var = self.variance() - other.variance();
         if var <= 0.0 {
             // When sigma_self ≈ sigma_other (including ULP-level rounding differences
             // from the pi→sigma accessor round-trip), the excluded contribution is N00.
             // Computing from_ms(tiny_mu, 0.0) would give {pi:inf, tau:inf}, whose
             // mu() = inf/inf = NaN.  Returning N00 is correct: when both Gaussians
             // carry the same variance, the residual is a point mass at 0.
-            return Gaussian::from_ms(0.0, 0.0);
+            return Gaussian::from_mv(0.0, 0.0);
         }
-        let mu = self.mu() - other.mu();
-        Self::from_ms(mu, var.sqrt())
+
+        Self::from_mv(self.mu() - other.mu(), var)
     }
 
     pub(crate) fn forget(&self, variance_delta: f64) -> Self {
-        let var = self.sigma().powi(2) + variance_delta;
-        Self::from_ms(self.mu(), var.sqrt())
+        Self::from_mv(self.mu(), self.variance() + variance_delta)
     }
 
     /// EP damping in natural-parameter space: `α·new + (1−α)·self`.
@@ -128,9 +164,7 @@ impl ops::Add<Gaussian> for Gaussian {
     /// Variance addition: (mu1 + mu2, sqrt(σ1² + σ2²)).
     /// Used for combining performance and noise; rare relative to mul/div.
     fn add(self, rhs: Gaussian) -> Self::Output {
-        let mu = self.mu() + rhs.mu();
-        let var = self.sigma().powi(2) + rhs.sigma().powi(2);
-        Self::from_ms(mu, var.sqrt())
+        Self::from_mv(self.mu() + rhs.mu(), self.variance() + rhs.variance())
     }
 }
 
@@ -138,9 +172,7 @@ impl ops::Sub<Gaussian> for Gaussian {
     type Output = Gaussian;
     /// (mu1 - mu2, sqrt(σ1² + σ2²)). Same sigma combination as Add.
     fn sub(self, rhs: Gaussian) -> Self::Output {
-        let mu = self.mu() - rhs.mu();
-        let var = self.sigma().powi(2) + rhs.sigma().powi(2);
-        Self::from_ms(mu, var.sqrt())
+        Self::from_mv(self.mu() - rhs.mu(), self.variance() + rhs.variance())
     }
 }
 
@@ -161,7 +193,7 @@ impl ops::Mul<f64> for Gaussian {
         if scalar == 0.0 {
             // Scaling by 0 collapses to a point mass at 0 (sigma' = 0, mu' = 0).
             // This is N00, the additive identity, NOT N_INF.
-            return Gaussian::from_ms(0.0, 0.0);
+            return Gaussian::from_mv(0.0, 0.0);
         }
         // sigma' = sigma * |scalar|  =>  pi' = pi / scalar²
         // mu'    = mu * scalar       =>  tau' = tau / scalar