fix(stats): keep r² within [0, 1] when y sits far from zero

The centred sums lose their last digits when every y is close to one
large value, and a two-point fit came out at r² 1.00000004. Clamp the
ratio to 1 and the y sum of squares at 0, and pin the two shrunk
counterexamples the nightly property run found.

src/lib/stats/__tests__/stats.property.test.ts

@@ -244,11 +244,30 @@ describe("least-squares line", () => {
           return;
         }
         expect(fit.r2).toBeGreaterThanOrEqual(0);
-        expect(fit.r2).toBeLessThanOrEqual(1 + 1e-9);
+        expect(fit.r2).toBeLessThanOrEqual(1);
       }),
     );
   });
 
+  // Shrunk counterexamples of the property above from two nightly seeds
+  // (802398975 and 1070990592): two points on a line far from zero, where
+  // the centred sums rounded r² above 1.
+  it.each([
+    [
+      [0, 1],
+      [-999.87, -999.97],
+    ],
+    [
+      [0, -1],
+      [9996.28, 9993.55],
+    ],
+  ])("keeps r² at most 1 for x %j, y %j", (xs, ys) => {
+    const fit = linearFit(xs, ys)!;
+    expect(fit.r2).not.toBeNull();
+    expect(fit.r2!).toBeLessThanOrEqual(1);
+    expect(fit.r2!).toBeGreaterThan(0.999);
+  });
+
   it("refuses a vertical line instead of inventing a slope", () => {
     assertProperty(
       fc.property(

src/lib/stats/linear-fit.ts

@@ -48,7 +48,9 @@ export function fitFromSums(sums: FitSums): LinearFit | null {
   if (n < 2) return null;
   const sxx = sumXx - (sumX * sumX) / n;
   const sxy = sumXy - (sumX * sumY) / n;
-  const syy = sumYy - (sumY * sumY) / n;
+  // A sum of squares is never negative; the subtraction can still land a
+  // hair below zero when every y sits close to one large value.
+  const syy = Math.max(0, sumYy - (sumY * sumY) / n);
   if (sxx === 0) return null;
   const slope = sxy / sxx;
   const meanX = sumX / n;
@@ -56,7 +58,10 @@ export function fitFromSums(sums: FitSums): LinearFit | null {
   return {
     slope,
     intercept: meanY - slope * meanX,
-    r2: syy === 0 ? null : (sxy * sxy) / (sxx * syy),
+    // r² is at most 1 by Cauchy–Schwarz. With y far from zero and little
+    // spread (two weighings of about 1000, a tenth apart) the centred sums
+    // lose their last digits and the ratio came out at 1.00000004.
+    r2: syy === 0 ? null : Math.min(1, (sxy * sxy) / (sxx * syy)),
     n,
     meanX,
     meanY,