feat!: Gaussian's EP operations stop wearing arithmetic's clothes

`Gaussian` publicly implemented `Mul`, `Div`, `Add` and `Sub`. They were
the EP product, cavity and variance-space convolutions, and every one of
them lies to a reader who takes the operator at face value:

    a = N(10, 2)   b = N(4, 3)   c = N(1, 1)

    a * b        N(8.15, 1.66)   not 40
    a - b        sigma GREW, 2 -> sqrt(4 + 9)
    a * N(1, 0)  mu = NaN        "multiply by one"
    a / c        pi = -0.75      mu() prints a confident 0

The last is this crate's signature defect on a public operator. `Div` is
the cavity and can legitimately leave a negative precision, which is not
a distribution — and `mu()`/`sigma()` guard `pi <= 0` and report `0.0`
and `inf`, so it comes back as a plausible number with no panic, no
`Debug` marker and nothing to test against.

The four impls are now `pub(crate)` inherent methods that say what they
do: `ep_product`, `cavity`, `convolve`, `convolve_diff`, plus `scale`
for the one operation that genuinely is arithmetic. Nothing in a user's
workflow needed operator syntax; inference did, and it still has it.

`pi()` and `tau()` follow. Storing natural parameters is a performance
decision — it makes message passing two adds — not a contract. The
public surface is now exactly: `from_ms`, `from_mv`, `mu`, `sigma`,
`variance`, `probability_below`, `probability_above`. `from_mv` and
`variance` are promoted from `pub(crate)`; they are the honest pair for
callers who already hold a variance and should not pay a round trip
through the square root.

Four integration tests asserted bit-identity on `(pi, tau)`. They assert
it on `(mu, variance)` instead — still `assert_eq!`, still exact, and
`1/pi` and `tau/pi` are deterministic, so bit-equal natural parameters
give bit-equal moments. `a_nan_sigma_passes_through_from_ms` drops its
`|| g.pi().is_nan()` half: `sigma()` substitutes for `pi <= 0` and
`pi == inf`, so NaN survives to it only from a NaN precision.

`benches/gaussian.rs` is deleted. It timed two f64 additions through the
public operators, and keeping those public solely to feed it is the same
thing #73 objected to when a benchmark was dictating five public types.
The paths it covered are exercised by `batch` and `history_converge`
through the real call chain.

Closes #71.

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

src/time_slice.rs

@@ -26,7 +26,9 @@ pub(crate) struct Skill {
 
 impl Skill {
     pub(crate) fn posterior(&self) -> Gaussian {
-        self.likelihood * self.backward * self.forward
+        self.likelihood
+            .ep_product(self.backward)
+            .ep_product(self.forward)
     }
 }
 
@@ -74,7 +76,7 @@ impl Item {
         if forward {
             Rating::new(skill.forward, r.beta, r.drift).with_drift_scale(r.drift_scale)
         } else {
-            Rating::new(skill.posterior() / self.likelihood, r.beta, r.drift)
+            Rating::new(skill.posterior().cavity(self.likelihood), r.beta, r.drift)
                 .with_drift_scale(r.drift_scale)
         }
     }
@@ -171,7 +173,7 @@ impl Event {
             for (i, item) in team.items.iter_mut().enumerate() {
                 let fresh = update.likelihoods[t][i];
                 let old_likelihood = skills.at(item.slot).likelihood;
-                let new_likelihood = (old_likelihood / item.likelihood) * fresh;
+                let new_likelihood = old_likelihood.cavity(item.likelihood).ep_product(fresh);
                 skills.at_mut(item.slot).likelihood = new_likelihood;
                 item.likelihood = fresh;
             }
@@ -435,8 +437,9 @@ impl<T: Time> TimeSlice<T> {
                 for (t, team) in event.teams.iter_mut().enumerate() {
                     for (i, item) in team.items.iter_mut().enumerate() {
                         let old_likelihood = self.skills.at(item.slot).likelihood;
-                        let new_likelihood =
-                            (old_likelihood / item.likelihood) * g.likelihoods[t][i];
+                        let new_likelihood = old_likelihood
+                            .cavity(item.likelihood)
+                            .ep_product(g.likelihoods[t][i]);
                         self.skills.at_mut(item.slot).likelihood = new_likelihood;
                         item.likelihood = g.likelihoods[t][i];
                     }
@@ -586,7 +589,7 @@ impl<T: Time> TimeSlice<T> {
 
     pub(crate) fn forward_prior_out(&self, competitor: &Index) -> Gaussian {
         let skill = self.skills.get(*competitor).unwrap();
-        skill.forward * skill.likelihood
+        skill.forward.ep_product(skill.likelihood)
     }
 
     pub(crate) fn backward_prior_out<D: Drift<T>>(
@@ -595,7 +598,7 @@ impl<T: Time> TimeSlice<T> {
         competitors: &CompetitorStore<T, D>,
     ) -> Gaussian {
         let skill = self.skills.get(*competitor).unwrap();
-        let n = skill.likelihood * skill.backward;
+        let n = skill.likelihood.ep_product(skill.backward);
         n.forget(
             competitors[*competitor]
                 .rating