diff --git a/content/graph/WeightedMatching.h b/content/graph/WeightedMatching.h
index 17f0b5eba..e02d0de6f 100644
--- a/content/graph/WeightedMatching.h
+++ b/content/graph/WeightedMatching.h
@@ -1,47 +1,40 @@
 /**
- * Author: Benjamin Qi, chilli
- * Date: 2020-04-04
+ * Author: Jeremy Lim, Joshua Andersson
+ * Date: 2026-07-09
  * License: CC0
- * Source: https://github.com/bqi343/USACO/blob/master/Implementations/content/graphs%20(12)/Matching/Hungarian.h
+ * Source: https://github.com/scipy/scipy/blob/main/scipy/optimize/rectangular_lsap/rectangular_lsap.cpp
  * Description: Given a weighted bipartite graph, matches every node on
  * the left with a node on the right such that no
  * nodes are in two matchings and the sum of the edge weights is minimal. Takes
- * cost[N][M], where cost[i][j] = cost for L[i] to be matched with R[j] and
+ * cost[N][M], where cost[i][j] = cost for L[i] to be matched with R[j], and
  * returns (min cost, match), where L[i] is matched with
  * R[match[i]]. Negate costs for max cost. Requires $N \le M$.
  * Time: O(N^2M)
- * Status: Tested on kattis:cordonbleu, stress-tested
+ * Status: Tested on kattis:cordonbleu, kattis:engaging, stress-tested
  */
 #pragma once
 
-pair<int, vi> hungarian(const vector<vi> &a) {
-	if (a.empty()) return {0, {}};
-	int n = sz(a) + 1, m = sz(a[0]) + 1;
-	vi u(n), v(m), p(m), ans(n - 1);
-	rep(i,1,n) {
-		p[0] = i;
-		int j0 = 0; // add "dummy" worker 0
-		vi dist(m, INT_MAX), pre(m, -1);
-		vector<bool> done(m + 1);
-		do { // dijkstra
-			done[j0] = true;
-			int i0 = p[j0], j1, delta = INT_MAX;
-			rep(j,1,m) if (!done[j]) {
-				auto cur = a[i0 - 1][j - 1] - u[i0] - v[j];
-				if (cur < dist[j]) dist[j] = cur, pre[j] = j0;
-				if (dist[j] < delta) delta = dist[j], j1 = j;
+template <typename T> pair<T, vector<int>>
+weighted_matching(vector<vector<T>> &C) {
+	int i = sz(C), m = i ? sz(C[0]) : 0, c, s, r;
+	vector<T> dist(m), potential(m);
+	vi row_match(i), col_match(m, -1), cols(m), prev(m);
+	T d, nd, cost = 0;
+	while (i--) {
+		rep(c,0,m) dist[c] = C[i][c], cols[c] = c, prev[c] = i;
+		for (s = 0;;) {
+			rep(j,s,m) {
+				c = cols[j], nd = dist[c] - potential[c];
+				if (j == s || d > nd) d = nd, swap(cols[s], cols[j]);
 			}
-			rep(j,0,m) {
-				if (done[j]) u[p[j]] += delta, v[j] -= delta;
-				else dist[j] -= delta;
-			}
-			j0 = j1;
-		} while (p[j0]);
-		while (j0) { // update alternating path
-			int j1 = pre[j0];
-			p[j0] = p[j1], j0 = j1;
+			if (!~(r = col_match[c = cols[s++]])) break;
+			rep(j,0,m) if (dist[j] > (nd = C[r][j]-C[r][c]+dist[c]))
+				dist[j] = nd, prev[j] = r;
 		}
+		for (cost += dist[c]; s--;)
+			potential[cols[s]] = dist[cols[s]] - d;
+		for (; r != i; swap(c, row_match[r]))
+			r = col_match[c] = prev[c];
 	}
-	rep(j,1,m) if (p[j]) ans[p[j] - 1] = j - 1;
-	return {-v[0], ans}; // min cost
+	return {cost, row_match};
 }
diff --git a/stress-tests/graph/WeightedMatching.cpp b/stress-tests/graph/WeightedMatching.cpp
index 01f295550..352fe5686 100644
--- a/stress-tests/graph/WeightedMatching.cpp
+++ b/stress-tests/graph/WeightedMatching.cpp
@@ -8,7 +8,7 @@
 
 void test(int N, int mxCost, int iters) {
 	for (int it = 0; it < iters; it++) {
-		int n = randRange(0, N), m = randRange(0, N);
+		int n = randIncl(0, N), m = randIncl(0, N);
 		if (n > m)
 			swap(n, m);
 
@@ -29,7 +29,7 @@ void test(int N, int mxCost, int iters) {
 		}
 		mcmf.setpi(s);
 		auto maxflow = mcmf.maxflow(s, t);
-		auto matching = hungarian(cost);
+		auto matching = weighted_matching(cost);
 		assert(maxflow.first == n);
 		assert(maxflow.second == matching.first);
 		int matchSum = 0;
@@ -40,7 +40,6 @@ void test(int N, int mxCost, int iters) {
 			used.insert(matching.second[i]);
 		}
 		assert(matchSum == matching.first);
-		return;
 	}
 }
 signed main() {