diff --git a/src/history.rs b/src/history.rs
index fb56932..aea04e2 100644
--- a/src/history.rs
+++ b/src/history.rs
@@ -389,7 +389,15 @@ impl<T: Time, D: Drift<T>, O: Observer<T>, K: Eq + Hash + Clone> History<T, D, O
 
     /// Draw-probability quality metric for the given teams (key slices).
     ///
-    /// Values range roughly [0, 1]; 1 == perfectly matched.
+    /// Values range roughly [0, 1]; 1 == perfectly matched. Supports any
+    /// number of teams.
+    ///
+    /// # Panics
+    ///
+    /// Panics if fewer than two teams are supplied, or if a team resolves to
+    /// no known competitors — keys absent from the history, or competitors
+    /// with no recorded skill, are dropped, so a team of entirely-unknown
+    /// keys becomes empty. Use `lookup` to check keys first.
     pub fn predict_quality(&self, teams: &[&[&K]]) -> f64 {
         let groups: Vec<Vec<Gaussian>> = teams
             .iter()
diff --git a/src/lib.rs b/src/lib.rs
index cdacbfa..227c4f1 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -247,7 +247,26 @@ pub(crate) fn sort_time<T: Copy + Ord>(xs: &[T], reverse: bool) -> Vec<usize> {
 }
 
 /// Calculates the match quality of the given rating groups. A result is the draw probability in the association
+///
+/// Supports any number of groups. Values range roughly `[0, 1]`; 1 means a
+/// perfectly balanced match.
+///
+/// # Panics
+///
+/// Panics if fewer than two rating groups are supplied, or if any group is
+/// empty — match quality is a property of a contest between at least two
+/// non-empty sides.
 pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
+    assert!(
+        rating_groups.len() >= 2,
+        "quality() requires at least 2 rating groups, got {}",
+        rating_groups.len()
+    );
+    assert!(
+        rating_groups.iter().all(|group| !group.is_empty()),
+        "quality() requires every rating group to be non-empty"
+    );
+
     let flatten_ratings = rating_groups
         .iter()
         .flat_map(|group| group.iter())
@@ -271,8 +290,10 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
 
     let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length);
 
+    // Row `row` contrasts group `row` (+weight) against group `row + 1`
+    // (-weight). `t` is the column where the current group's players start;
+    // the negative block begins immediately after it.
     let mut t = 0;
-    let mut x = 0;
 
     for (row, group) in rating_groups.windows(2).enumerate() {
         let current = group[0];
@@ -280,17 +301,13 @@ pub fn quality(rating_groups: &[&[Gaussian]], beta: f64) -> f64 {
 
         for n in t..t + current.len() {
             rotated_a_matrix[(row, n)] = flatten_weights[n];
-
-            x += 1;
         }
 
         t += current.len();
 
-        for n in x..x + next.len() {
+        for n in t..t + next.len() {
             rotated_a_matrix[(row, n)] = -flatten_weights[n];
         }
-
-        x += next.len();
     }
 
     let a_matrix = rotated_a_matrix.transpose();
diff --git a/src/matrix.rs b/src/matrix.rs
index cf74061..cf29777 100644
--- a/src/matrix.rs
+++ b/src/matrix.rs
@@ -1,29 +1,13 @@
+//! Minimal dense matrix used by `quality()`.
+//!
+//! `determinant` and `inverse` go through one LU decomposition with partial
+//! pivoting — O(n³) and numerically stable. The previous implementation
+//! expanded cofactors recursively (O(n!), allocating a `Vec` per minor) and
+//! only implemented `inverse` for the 1×1 case, which limited `quality()` to
+//! exactly two rating groups.
+
 use std::ops;
 
-fn det(m: &[f64], x: usize) -> f64 {
-    if x == 1 {
-        m[0]
-    } else if x == 2 {
-        m[0] * m[3] - m[1] * m[2]
-    } else {
-        let mut d = 0.0;
-
-        for n in 0..x {
-            let ms = m
-                .iter()
-                .enumerate()
-                .skip(x)
-                .filter(|(i, _)| (i % x) != n)
-                .map(|(_, v)| *v)
-                .collect::<Vec<_>>();
-
-            d += (-1.0f64).powi(n as i32) * m[n] * det(&ms, x - 1);
-        }
-
-        d
-    }
-}
-
 #[derive(Clone, Debug)]
 pub struct Matrix {
     data: Box<[f64]>,
@@ -31,6 +15,107 @@ pub struct Matrix {
     width: usize,
 }
 
+/// LU decomposition with partial pivoting: `PA = LU`, stored compactly.
+///
+/// `lu` holds `L` below the diagonal (unit diagonal implied) and `U` on and
+/// above it. `sign` is the determinant sign contributed by row swaps, or 0.0
+/// when the matrix is singular.
+struct Lu {
+    lu: Vec<f64>,
+    perm: Vec<usize>,
+    n: usize,
+    sign: f64,
+}
+
+impl Lu {
+    fn decompose(m: &Matrix) -> Self {
+        debug_assert_eq!(m.width, m.height, "LU requires a square matrix");
+
+        let n = m.width;
+        let mut lu = m.data.to_vec();
+        let mut perm: Vec<usize> = (0..n).collect();
+        let mut sign = 1.0;
+
+        for col in 0..n {
+            // Partial pivot: take the largest-magnitude candidate to limit
+            // growth of round-off in the elimination below.
+            let mut pivot_row = col;
+            let mut pivot_max = lu[col * n + col].abs();
+
+            for row in (col + 1)..n {
+                let candidate = lu[row * n + col].abs();
+                if candidate > pivot_max {
+                    pivot_max = candidate;
+                    pivot_row = row;
+                }
+            }
+
+            if pivot_max == 0.0 {
+                sign = 0.0;
+                continue;
+            }
+
+            if pivot_row != col {
+                for k in 0..n {
+                    lu.swap(col * n + k, pivot_row * n + k);
+                }
+                perm.swap(col, pivot_row);
+                sign = -sign;
+            }
+
+            let pivot = lu[col * n + col];
+
+            for row in (col + 1)..n {
+                let factor = lu[row * n + col] / pivot;
+                lu[row * n + col] = factor;
+
+                for k in (col + 1)..n {
+                    lu[row * n + k] -= factor * lu[col * n + k];
+                }
+            }
+        }
+
+        Self { lu, perm, n, sign }
+    }
+
+    fn determinant(&self) -> f64 {
+        if self.sign == 0.0 {
+            return 0.0;
+        }
+
+        let mut det = self.sign;
+        for i in 0..self.n {
+            det *= self.lu[i * self.n + i];
+        }
+
+        det
+    }
+
+    /// Solve `Ax = b` for a single column of the identity, giving one column
+    /// of the inverse.
+    fn solve_column(&self, col: usize, out: &mut [f64]) {
+        let n = self.n;
+
+        // Forward substitution through L, applying the row permutation.
+        for i in 0..n {
+            let mut sum = if self.perm[i] == col { 1.0 } else { 0.0 };
+            for (k, &solved) in out.iter().enumerate().take(i) {
+                sum -= self.lu[i * n + k] * solved;
+            }
+            out[i] = sum;
+        }
+
+        // Back substitution through U.
+        for i in (0..n).rev() {
+            let mut sum = out[i];
+            for (k, &solved) in out.iter().enumerate().skip(i + 1) {
+                sum -= self.lu[i * n + k] * solved;
+            }
+            out[i] = sum / self.lu[i * n + i];
+        }
+    }
+}
+
 impl Matrix {
     pub fn new(height: usize, width: usize) -> Matrix {
         Matrix {
@@ -52,73 +137,59 @@ impl Matrix {
         matrix
     }
 
-    pub fn minor(&self, row_n: usize, col_n: usize) -> Matrix {
-        let mut matrix = Matrix::new(self.height - 1, self.width - 1);
-
-        let mut nr = 0;
-
-        for r in 0..self.height {
-            if r == row_n {
-                continue;
-            }
-
-            let mut nc = 0;
-
-            for c in 0..self.width {
-                if c == col_n {
-                    continue;
-                }
-
-                matrix[(nr, nc)] = self[(r, c)];
-
-                nc += 1;
-            }
-
-            nr += 1;
-        }
-
-        matrix
-    }
-
+    /// Determinant of a square matrix. The 0×0 determinant is 1 by convention
+    /// (the empty product).
+    ///
+    /// # Panics
+    ///
+    /// Panics if the matrix is not square.
     pub fn determinant(&self) -> f64 {
-        debug_assert!(self.width == self.height);
+        assert_eq!(
+            self.width, self.height,
+            "determinant requires a square matrix, got {}x{}",
+            self.height, self.width
+        );
 
-        det(&self.data, self.width)
+        if self.width == 0 {
+            return 1.0;
+        }
+
+        Lu::decompose(self).determinant()
     }
 
-    pub fn adjugate(&self) -> Matrix {
-        debug_assert!(self.width == self.height);
+    /// Matrix inverse via LU decomposition.
+    ///
+    /// # Panics
+    ///
+    /// Panics if the matrix is not square or is singular.
+    pub fn inverse(&self) -> Matrix {
+        assert_eq!(
+            self.width, self.height,
+            "inverse requires a square matrix, got {}x{}",
+            self.height, self.width
+        );
 
-        let mut matrix = Matrix::new(self.height, self.width);
+        let n = self.width;
+        let mut inverse = Matrix::new(n, n);
 
-        if matrix.height == 2 {
-            matrix[(0, 0)] = self[(1, 1)];
-            matrix[(0, 1)] = -self[(0, 1)];
-            matrix[(1, 0)] = -self[(1, 0)];
-            matrix[(1, 1)] = self[(0, 0)];
-        } else {
-            for r in 0..matrix.height {
-                for c in 0..matrix.width {
-                    let sign = if (r + c) % 2 == 0 { 1.0 } else { -1.0 };
+        if n == 0 {
+            return inverse;
+        }
 
-                    matrix[(r, c)] = self.minor(r, c).determinant() * sign;
-                }
+        let lu = Lu::decompose(self);
+        assert!(lu.sign != 0.0, "cannot invert a singular matrix");
+
+        let mut column = vec![0.0; n];
+
+        for c in 0..n {
+            lu.solve_column(c, &mut column);
+
+            for (r, &value) in column.iter().enumerate() {
+                inverse[(r, c)] = value;
             }
         }
 
-        matrix
-    }
-
-    pub fn inverse(&self) -> Matrix {
-        let mut matrix = Matrix::new(self.width, self.height);
-
-        if self.height == self.width && self.height == 1 {
-            matrix[(0, 0)] = 1.0 / self[(0, 0)];
-        } else {
-            panic!("eh, okey")
-        }
-
-        matrix
+        inverse
     }
 }
 
@@ -126,20 +197,62 @@ impl ops::Index<(usize, usize)> for Matrix {
     type Output = f64;
 
     fn index(&self, pos: (usize, usize)) -> &Self::Output {
+        debug_assert!(
+            pos.0 < self.height && pos.1 < self.width,
+            "index ({}, {}) out of bounds for {}x{} matrix",
+            pos.0,
+            pos.1,
+            self.height,
+            self.width
+        );
+
         &self.data[(self.width * pos.0) + pos.1]
     }
 }
 
 impl ops::IndexMut<(usize, usize)> for Matrix {
     fn index_mut(&mut self, pos: (usize, usize)) -> &mut Self::Output {
+        debug_assert!(
+            pos.0 < self.height && pos.1 < self.width,
+            "index ({}, {}) out of bounds for {}x{} matrix",
+            pos.0,
+            pos.1,
+            self.height,
+            self.width
+        );
+
         &mut self.data[(self.width * pos.0) + pos.1]
     }
 }
 
-impl<'a> ops::Mul<&'a Matrix> for f64 {
+fn multiply(lhs: &Matrix, rhs: &Matrix) -> Matrix {
+    assert_eq!(
+        lhs.width, rhs.height,
+        "cannot multiply {}x{} by {}x{}",
+        lhs.height, lhs.width, rhs.height, rhs.width
+    );
+
+    let mut matrix = Matrix::new(lhs.height, rhs.width);
+
+    for r in 0..matrix.height {
+        for c in 0..matrix.width {
+            let mut value = 0.0;
+
+            for x in 0..lhs.width {
+                value += lhs[(r, x)] * rhs[(x, c)];
+            }
+
+            matrix[(r, c)] = value;
+        }
+    }
+
+    matrix
+}
+
+impl ops::Mul<&Matrix> for f64 {
     type Output = Matrix;
 
-    fn mul(self, rhs: &'a Matrix) -> Matrix {
+    fn mul(self, rhs: &Matrix) -> Matrix {
         let mut matrix = Matrix::new(rhs.height, rhs.width);
 
         for r in 0..rhs.height {
@@ -152,54 +265,35 @@ impl<'a> ops::Mul<&'a Matrix> for f64 {
     }
 }
 
-impl<'a> ops::Mul<&'a Matrix> for Matrix {
+impl ops::Mul<&Matrix> for Matrix {
     type Output = Matrix;
 
-    fn mul(self, rhs: &'a Matrix) -> Matrix {
-        let mut matrix = Matrix::new(self.height, rhs.width);
-
-        for r in 0..matrix.height {
-            for c in 0..matrix.width {
-                let mut value = 0.0;
-
-                for x in 0..self.width {
-                    value += self[(r, x)] * rhs[(x, c)];
-                }
-
-                matrix[(r, c)] = value;
-            }
-        }
-
-        matrix
+    fn mul(self, rhs: &Matrix) -> Matrix {
+        multiply(&self, rhs)
     }
 }
 
-impl<'a> ops::Mul<&'a Matrix> for &'a Matrix {
+impl ops::Mul<&Matrix> for &Matrix {
     type Output = Matrix;
 
-    fn mul(self, rhs: &'a Matrix) -> Matrix {
-        let mut matrix = Matrix::new(self.height, rhs.width);
-
-        for r in 0..matrix.height {
-            for c in 0..matrix.width {
-                let mut value = 0.0;
-
-                for x in 0..self.width {
-                    value += self[(r, x)] * rhs[(x, c)];
-                }
-
-                matrix[(r, c)] = value;
-            }
-        }
-
-        matrix
+    fn mul(self, rhs: &Matrix) -> Matrix {
+        multiply(self, rhs)
     }
 }
 
-impl<'a> ops::Add<&'a Matrix> for &'a Matrix {
+impl ops::Add<&Matrix> for &Matrix {
     type Output = Matrix;
 
-    fn add(self, rhs: &'a Matrix) -> Matrix {
+    fn add(self, rhs: &Matrix) -> Matrix {
+        assert!(
+            self.height == rhs.height && self.width == rhs.width,
+            "cannot add {}x{} to {}x{}",
+            self.height,
+            self.width,
+            rhs.height,
+            rhs.width
+        );
+
         let mut matrix = Matrix::new(self.height, self.width);
 
         for r in 0..matrix.height {
@@ -211,3 +305,106 @@ impl<'a> ops::Add<&'a Matrix> for &'a Matrix {
         matrix
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn from_rows(rows: &[&[f64]]) -> Matrix {
+        let mut m = Matrix::new(rows.len(), rows[0].len());
+        for (r, row) in rows.iter().enumerate() {
+            for (c, &v) in row.iter().enumerate() {
+                m[(r, c)] = v;
+            }
+        }
+        m
+    }
+
+    #[test]
+    fn determinant_1x1() {
+        assert!((from_rows(&[&[3.0]]).determinant() - 3.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn determinant_2x2() {
+        let m = from_rows(&[&[1.0, 2.0], &[3.0, 4.0]]);
+        assert!((m.determinant() - (-2.0)).abs() < 1e-12);
+    }
+
+    #[test]
+    fn determinant_3x3() {
+        let m = from_rows(&[&[6.0, 1.0, 1.0], &[4.0, -2.0, 5.0], &[2.0, 8.0, 7.0]]);
+        assert!((m.determinant() - (-306.0)).abs() < 1e-10);
+    }
+
+    #[test]
+    fn determinant_requires_no_pivot_at_origin() {
+        // A zero in the top-left forces a row swap; the sign must follow.
+        let m = from_rows(&[&[0.0, 1.0], &[1.0, 0.0]]);
+        assert!((m.determinant() - (-1.0)).abs() < 1e-12);
+    }
+
+    #[test]
+    fn determinant_of_singular_is_zero() {
+        let m = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]);
+        assert!(m.determinant().abs() < 1e-12);
+    }
+
+    #[test]
+    fn inverse_1x1() {
+        let inv = from_rows(&[&[4.0]]).inverse();
+        assert!((inv[(0, 0)] - 0.25).abs() < 1e-12);
+    }
+
+    #[test]
+    fn inverse_times_original_is_identity() {
+        for rows in [
+            vec![vec![1.0, 2.0], vec![3.0, 4.0]],
+            vec![
+                vec![6.0, 1.0, 1.0],
+                vec![4.0, -2.0, 5.0],
+                vec![2.0, 8.0, 7.0],
+            ],
+            vec![
+                vec![2.0, 0.0, 1.0, 3.0],
+                vec![1.0, 5.0, 2.0, 0.0],
+                vec![0.0, 1.0, 4.0, 1.0],
+                vec![3.0, 2.0, 0.0, 6.0],
+            ],
+        ] {
+            let refs: Vec<&[f64]> = rows.iter().map(|r| r.as_slice()).collect();
+            let m = from_rows(&refs);
+            let product = &m * &m.inverse();
+
+            for r in 0..product.height {
+                for c in 0..product.width {
+                    let expected = if r == c { 1.0 } else { 0.0 };
+                    assert!(
+                        (product[(r, c)] - expected).abs() < 1e-9,
+                        "({r},{c}) = {} expected {expected}",
+                        product[(r, c)]
+                    );
+                }
+            }
+        }
+    }
+
+    #[test]
+    #[should_panic(expected = "singular")]
+    fn inverse_of_singular_panics() {
+        let _ = from_rows(&[&[1.0, 2.0], &[2.0, 4.0]]).inverse();
+    }
+
+    #[test]
+    fn empty_determinant_is_one() {
+        assert!((Matrix::new(0, 0).determinant() - 1.0).abs() < 1e-12);
+    }
+
+    #[test]
+    fn transpose_round_trips() {
+        let m = from_rows(&[&[1.0, 2.0, 3.0], &[4.0, 5.0, 6.0]]);
+        let t = m.transpose();
+        assert_eq!((t.height, t.width), (3, 2));
+        assert_eq!(t.transpose()[(1, 2)], m[(1, 2)]);
+    }
+}
diff --git a/tests/quality.rs b/tests/quality.rs
new file mode 100644
index 0000000..fe7c2c4
--- /dev/null
+++ b/tests/quality.rs
@@ -0,0 +1,119 @@
+//! `quality()` beyond two rating groups.
+//!
+//! The historical golden (two equal singletons) is asserted in
+//! `src/lib.rs::tests::test_quality`. These cover the N-group generalisation,
+//! which previously panicked with an out-of-bounds index at 3+ groups.
+
+use trueskill_tt::{Gaussian, quality};
+
+const BETA: f64 = 25.0 / 3.0 / 2.0;
+
+fn rating(mu: f64, sigma: f64) -> Gaussian {
+    Gaussian::from_ms(mu, sigma)
+}
+
+#[test]
+fn three_equal_groups_is_finite_and_in_range() {
+    let r = rating(25.0, 3.0);
+    let q = quality(&[&[r], &[r], &[r]], BETA);
+
+    assert!(q.is_finite(), "quality must be finite, got {q}");
+    assert!((0.0..=1.0).contains(&q), "quality out of range: {q}");
+}
+
+#[test]
+fn quality_supports_many_groups() {
+    let r = rating(25.0, 3.0);
+    for n in 2..=8 {
+        let holders: Vec<[Gaussian; 1]> = (0..n).map(|_| [r]).collect();
+        let groups: Vec<&[Gaussian]> = holders.iter().map(|g| g.as_slice()).collect();
+        let q = quality(&groups, BETA);
+        assert!(q.is_finite(), "n={n}: quality must be finite, got {q}");
+        assert!((0.0..=1.0).contains(&q), "n={n}: out of range: {q}");
+    }
+}
+
+/// Equal-strength groups are the best-matched case: introducing a skill gap
+/// must lower quality.
+#[test]
+fn imbalance_lowers_quality() {
+    let strong = rating(40.0, 3.0);
+    let average = rating(25.0, 3.0);
+
+    let balanced = quality(&[&[average], &[average], &[average]], BETA);
+    let lopsided = quality(&[&[strong], &[average], &[average]], BETA);
+
+    assert!(
+        lopsided < balanced,
+        "expected imbalanced quality {lopsided} < balanced {balanced}"
+    );
+}
+
+/// Quality is a property of the multiset of groups, not their order.
+#[test]
+fn quality_is_permutation_invariant() {
+    let a = rating(30.0, 2.0);
+    let b = rating(25.0, 3.0);
+    let c = rating(20.0, 4.0);
+
+    let forward = quality(&[&[a], &[b], &[c]], BETA);
+    let reversed = quality(&[&[c], &[b], &[a]], BETA);
+
+    assert!(
+        (forward - reversed).abs() < 1e-9,
+        "permutation changed quality: {forward} vs {reversed}"
+    );
+}
+
+#[test]
+fn multi_player_groups_work() {
+    let r = rating(25.0, 3.0);
+    let q = quality(&[&[r, r], &[r, r], &[r, r]], BETA);
+    assert!(q.is_finite());
+    assert!((0.0..=1.0).contains(&q));
+}
+
+#[test]
+fn uneven_group_sizes_work() {
+    let r = rating(25.0, 3.0);
+    let q = quality(&[&[r, r], &[r], &[r, r, r]], BETA);
+    assert!(q.is_finite(), "got {q}");
+    assert!((0.0..=1.0).contains(&q), "got {q}");
+}
+
+#[test]
+#[should_panic(expected = "at least 2 rating groups")]
+fn single_group_panics_with_clear_message() {
+    let r = rating(25.0, 3.0);
+    let _ = quality(&[&[r]], BETA);
+}
+
+#[test]
+#[should_panic(expected = "at least 2 rating groups")]
+fn zero_groups_panics_with_clear_message() {
+    let _ = quality(&[], BETA);
+}
+
+#[test]
+#[should_panic(expected = "non-empty")]
+fn empty_group_panics_with_clear_message() {
+    let r = rating(25.0, 3.0);
+    let _ = quality(&[&[r], &[]], BETA);
+}
+
+#[test]
+fn history_predict_quality_supports_three_teams() {
+    use trueskill_tt::History;
+
+    let mut h = History::default();
+    h.record_winner(&"a", &"b", 1).unwrap();
+    h.record_winner(&"b", &"c", 2).unwrap();
+    h.converge().unwrap();
+
+    let q = h.predict_quality(&[&[&"a"], &[&"b"], &[&"c"]]);
+    assert!(
+        q.is_finite(),
+        "3-team predict_quality must be finite, got {q}"
+    );
+    assert!((0.0..=1.0).contains(&q), "out of range: {q}");
+}