fork download
copy 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3.  
  4. using ll = long long;
  5. const ll INF = (1LL << 60); // A very large number (safe for sums up to ~1e18)
  6.  
  7. // ===================================================================
  8. // This file contains a collection of algorithms to solve the 
  9. // "Assignment Problem" (and related matching problems). 
  10. // 
  11. // The main workhorse here is the Hungarian Algorithm (also called 
  12. // Kuhn-Munkres algorithm). 
  13. // 
  14. // Each function is ready to be used as a "black box". 
  15. // Read the comments above each one to understand:
  16. //   - What it solves
  17. //   - What input it expects
  18. //   - What it returns
  19. //   - Time complexity
  20. //   - Important constraints / assumptions
  21. // 
  22. // Technical terms explained simply:
  23. //   - "Assignment Problem": You have N workers and M jobs. 
  24. //     Each worker must be assigned to exactly one job (and each job 
  25. //     to at most one worker). You want the total cost to be minimum.
  26. //   - "Bipartite Graph": A graph with two sets of nodes (Left side = workers, 
  27. //     Right side = jobs). Edges only go from Left to Right.
  28. //   - "Matching": A set of edges where no two edges share a node.
  29. //     (e.g., each worker gets a unique job).
  30. //   - "Potentials" (u, v): Internal numbers used by the Hungarian 
  31. //     algorithm to guide the search. You don't need to understand them 
  32. //     to use the functions.
  33. //   - "Bitmask DP": Dynamic Programming where a "bitmask" (an integer 
  34. //     like 1011) represents a set of jobs that have been taken.
  35. // ===================================================================
  36.  
  37. // ===================================================================
  38. // 1) Hungarian Algorithm (Square Matrix) - Minimum Cost
  39. //    This is the standard O(n^3) algorithm.
  40. //    It finds the optimal assignment for a square cost matrix.
  41. // ===================================================================
  42.  
  43. // 1.1) Solve the assignment problem for a square matrix.
  44. //      Parameters:
  45. //        - cost: a 2D vector of size n x n. 
  46. //                cost[i][j] = cost of assigning worker i to job j.
  47. //      Returns:
  48. //        - a pair <total_cost, assignment>.
  49. //          - total_cost (long long): the minimum total cost.
  50. //          - assignment (vector<int>): a vector of size n.
  51. //            assignment[i] = j means worker i is assigned to job j.
  52. //            (Indices are 0-based).
  53. //      Time complexity: O(n^3), where n = number of rows = number of cols.
  54. //      Constraints:
  55. //        - Matrix must be square (n rows, n columns).
  56. //        - Costs can be negative, zero, or positive.
  57. //        - n should be >= 1. If n == 0, returns {0, {}}.
  58. //      Note: 
  59. //        - This function uses internal arrays "u", "v", "p", "way".
  60. //          You do NOT need to understand these to use the function.
  61. //        - The algorithm guarantees the optimal assignment.
  62. pair<ll, vector<int>> hungarianMinCost(const vector<vector<ll>>& cost) {
  63.     int n = (int)cost.size();
  64.     if (n == 0) return {0, {}};
  65.  
  66.     // Internal variables (do not worry about their meaning)
  67.     vector<ll> u(n + 1), v(n + 1);
  68.     vector<int> p(n + 1), way(n + 1);
  69.  
  70.     for (int i = 1; i <= n; i++) {
  71.         p[0] = i;
  72.         int j0 = 0;
  73.         vector<ll> minv(n + 1, INF);
  74.         vector<char> used(n + 1, false);
  75.         do {
  76.             used[j0] = true;
  77.             int i0 = p[j0];
  78.             ll delta = INF;
  79.             int j1 = 0;
  80.             for (int j = 1; j <= n; j++) {
  81.                 if (!used[j]) {
  82.                     ll cur = cost[i0 - 1][j - 1] - u[i0] - v[j];
  83.                     if (cur < minv[j]) {
  84.                         minv[j] = cur;
  85.                         way[j] = j0;
  86.                     }
  87.                     if (minv[j] < delta) {
  88.                         delta = minv[j];
  89.                         j1 = j;
  90.                     }
  91.                 }
  92.             }
  93.             for (int j = 0; j <= n; j++) {
  94.                 if (used[j]) {
  95.                     u[p[j]] += delta;
  96.                     v[j] -= delta;
  97.                 } else {
  98.                     minv[j] -= delta;
  99.                 }
  100.             }
  101.             j0 = j1;
  102.         } while (p[j0] != 0);
  103.  
  104.         do {
  105.             int j1 = way[j0];
  106.             p[j0] = p[j1];
  107.             j0 = j1;
  108.         } while (j0);
  109.     }
  110.  
  111.     vector<int> assignment(n);
  112.     for (int j = 1; j <= n; j++) {
  113.         if (p[j] > 0) {
  114.             assignment[p[j] - 1] = j - 1;
  115.         }
  116.     }
  117.     ll totalCost = -v[0]; // The total minimum cost
  118.     return {totalCost, assignment};
  119. }
  120.  
  121. // ===================================================================
  122. // 2) Hungarian Algorithm (Square Matrix) - Maximum Cost
  123. //    To maximize the total cost, we just negate all costs and run 
  124. //    the minimum cost version.
  125. // ===================================================================
  126.  
  127. // 2.1) Solve the assignment problem for a square matrix to MAXIMIZE cost.
  128. //      Parameters:
  129. //        - cost: a 2D vector of size n x n (the profit matrix).
  130. //      Returns:
  131. //        - a pair <max_profit, assignment>.
  132. //          - max_profit (long long): the maximum total profit.
  133. //          - assignment (vector<int>): optimal assignment.
  134. //      Time complexity: O(n^3).
  135. //      Constraints: Same as min cost version (square matrix).
  136. //      Note: 
  137. //        - Internally, it multiplies costs by -1 and calls the min-cost 
  138. //          Hungarian. So if costs are up to 1e9, the negated values are 
  139. //          safe within 64-bit.
  140. pair<ll, vector<int>> hungarianMaxCost(const vector<vector<ll>>& cost) {
  141.     int n = (int)cost.size();
  142.     vector<vector<ll>> negCost(n, vector<ll>(n));
  143.     for (int i = 0; i < n; i++) {
  144.         for (int j = 0; j < n; j++) {
  145.             negCost[i][j] = -cost[i][j];
  146.         }
  147.     }
  148.     auto res = hungarianMinCost(negCost);
  149.     return {-res.first, res.second}; // Negate the total cost back to positive
  150. }
  151.  
  152. // ===================================================================
  153. // 3) Hungarian Algorithm (Rectangular Matrix) - Minimum Cost
  154. //    In many problems, we have N workers and M jobs, where N <= M.
  155. //    We only need to assign each worker to a unique job.
  156. //    This function handles that case in O(N^2 * M).
  157. // ===================================================================
  158.  
  159. // 3.1) Solve assignment for rectangular matrix (N rows, M cols) with N <= M.
  160. //      Parameters:
  161. //        - cost: a 2D vector of size n x m. (n = rows, m = cols).
  162. //                Assumes n <= m (more jobs than workers).
  163. //      Returns:
  164. //        - a pair <total_cost, assignment>.
  165. //          - total_cost (ll): minimum cost to assign every row to a unique column.
  166. //          - assignment (vector<int>): size n. assignment[i] = j.
  167. //      Time complexity: O(n^2 * m). 
  168. //      Constraints:
  169. //        - n <= m (we cannot assign more workers than jobs).
  170. //        - If n > m, swap rows/cols (or transpose the matrix) before calling.
  171. //        - Costs can be negative.
  172. //      Note:
  173. //        - This is the standard CP-algorithms implementation adapted for 
  174. //          rectangular matrices. 
  175. //        - If you have more workers than jobs (n > m), you can call 
  176. //          hungarianMinCostRectangularTransposed (provided below).
  177. pair<ll, vector<int>> hungarianMinCostRectangular(const vector<vector<ll>>& cost) {
  178.     int n = (int)cost.size();    // rows (workers)
  179.     int m = (int)cost[0].size(); // cols (jobs)
  180.     if (n > m) {
  181.         // If workers > jobs, we cannot assign everyone.
  182.         // You must handle this case separately (e.g., add dummy jobs).
  183.         // This implementation assumes n <= m.
  184.         throw invalid_argument("Number of rows must be <= number of columns.");
  185.     }
  186.  
  187.     vector<ll> u(n + 1), v(m + 1);
  188.     vector<int> p(m + 1), way(m + 1);
  189.  
  190.     for (int i = 1; i <= n; i++) {
  191.         p[0] = i;
  192.         int j0 = 0;
  193.         vector<ll> minv(m + 1, INF);
  194.         vector<char> used(m + 1, false);
  195.         do {
  196.             used[j0] = true;
  197.             int i0 = p[j0];
  198.             ll delta = INF;
  199.             int j1 = 0;
  200.             for (int j = 1; j <= m; j++) {
  201.                 if (!used[j]) {
  202.                     ll cur = cost[i0 - 1][j - 1] - u[i0] - v[j];
  203.                     if (cur < minv[j]) {
  204.                         minv[j] = cur;
  205.                         way[j] = j0;
  206.                     }
  207.                     if (minv[j] < delta) {
  208.                         delta = minv[j];
  209.                         j1 = j;
  210.                     }
  211.                 }
  212.             }
  213.             for (int j = 0; j <= m; j++) {
  214.                 if (used[j]) {
  215.                     u[p[j]] += delta;
  216.                     v[j] -= delta;
  217.                 } else {
  218.                     minv[j] -= delta;
  219.                 }
  220.             }
  221.             j0 = j1;
  222.         } while (p[j0] != 0);
  223.  
  224.         do {
  225.             int j1 = way[j0];
  226.             p[j0] = p[j1];
  227.             j0 = j1;
  228.         } while (j0);
  229.     }
  230.  
  231.     vector<int> assignment(n);
  232.     for (int j = 1; j <= m; j++) {
  233.         if (p[j] > 0 && p[j] <= n) {
  234.             assignment[p[j] - 1] = j - 1;
  235.         }
  236.     }
  237.     ll totalCost = -v[0];
  238.     return {totalCost, assignment};
  239. }
  240.  
  241. // 3.2) Helper to handle the case where rows > cols (more workers than jobs).
  242. //      Parameters:
  243. //        - cost: n x m matrix where n > m.
  244. //      Returns: same as 3.1.
  245. //      How it works: 
  246. //        - It transposes the matrix (swaps rows and columns) so that the 
  247. //          number of rows becomes the smaller dimension, then calls 3.1.
  248. //        - Then it translates the assignment back to the original order.
  249. //      Time complexity: O(m^2 * n).
  250. pair<ll, vector<int>> hungarianMinCostRectangularTransposed(const vector<vector<ll>>& cost) {
  251.     int n = (int)cost.size();    // rows
  252.     int m = (int)cost[0].size(); // cols
  253.     if (n <= m) {
  254.         // If already n <= m, just call the normal one.
  255.         return hungarianMinCostRectangular(cost);
  256.     }
  257.     // Transpose: new matrix has size m x n
  258.     vector<vector<ll>> transCost(m, vector<ll>(n));
  259.     for (int i = 0; i < n; i++) {
  260.         for (int j = 0; j < m; j++) {
  261.             transCost[j][i] = cost[i][j];
  262.         }
  263.     }
  264.     auto res = hungarianMinCostRectangular(transCost); // res.assignment is of size m
  265.     // Now we need to map back.
  266.     // transCost assignment: transAssignment[row_j] = col_i.
  267.     // This means original job 'row_j' is assigned to original worker 'col_i'.
  268.     vector<int> originalAssignment(n, -1);
  269.     for (int j = 0; j < m; j++) {
  270.         int worker = res.second[j];
  271.         originalAssignment[worker] = j;
  272.     }
  273.     return {res.first, originalAssignment};
  274. }
  275.  
  276. // 3.3) Maximum Cost version for Rectangular matrices.
  277. //      Same logic: negate costs and call the min-cost version.
  278. pair<ll, vector<int>> hungarianMaxCostRectangular(const vector<vector<ll>>& cost) {
  279.     int n = (int)cost.size();
  280.     int m = (int)cost[0].size();
  281.     vector<vector<ll>> negCost(n, vector<ll>(m));
  282.     for (int i = 0; i < n; i++) {
  283.         for (int j = 0; j < m; j++) {
  284.             negCost[i][j] = -cost[i][j];
  285.         }
  286.     }
  287.     auto res = hungarianMinCostRectangular(negCost); // handles n <= m
  288.     return {-res.first, res.second};
  289. }
  290.  
  291. // ===================================================================
  292. // 4) Bitmask DP (Dynamic Programming) for Assignment 
  293. //    This is the "easy" way to solve assignment when N (or M) is small.
  294. //    It is VERY common in ECPC/ACPC when N <= 20.
  295. //    It is not as fast as Hungarian for N=1000, but it is much easier 
  296. //    to modify for extra constraints (like "must pick exactly K items").
  297. // ===================================================================
  298.  
  299. // 4.1) Assignment using Bitmask DP (Minimum Cost).
  300. //      Parameters:
  301. //        - cost: n x m matrix. 
  302. //        - We assign workers (rows) to jobs (cols) one by one.
  303. //        - n is the number of workers.
  304. //        - m is the number of jobs. (Usually n <= m, but if n > m, we can 
  305. //          flip the loops or handle it).
  306. //      Returns:
  307. //        - minimum total cost to assign all workers to unique jobs.
  308. //        - If impossible (n > m), returns INF.
  309. //      Time complexity: O(n * 2^m) or O(m * 2^n). 
  310. //        - Choose the smaller dimension for the bitmask to be efficient.
  311. //      Constraints:
  312. //        - The dimension used for the bitmask (which is the number of jobs 
  313. //          or workers) must be <= 20 (or 22 with optimization).
  314. //        - Works with negative costs.
  315. //      Note: 
  316. //        - This returns ONLY the total cost, not the assignment.
  317. //        - If you need the assignment, you can store the "parent" choice 
  318. //          in a separate array.
  319. //        - This is a "black box" for small N only.
  320. ll assignmentBitmaskDP(const vector<vector<ll>>& cost) {
  321.     int n = (int)cost.size(); // workers
  322.     int m = (int)cost[0].size(); // jobs
  323.  
  324.     if (n > m) return INF; // Not enough jobs for all workers.
  325.  
  326.     // DP over subsets of jobs.
  327.     // dp[mask] = minimum cost to assign the first k workers 
  328.     // (where k = number of set bits in mask) to the jobs in mask.
  329.     vector<ll> dp(1 << m, INF);
  330.     dp[0] = 0;
  331.  
  332.     for (int mask = 0; mask < (1 << m); mask++) {
  333.         int worker = __builtin_popcount(mask); // how many workers assigned so far
  334.         if (worker == n) continue; // all workers assigned
  335.  
  336.         for (int j = 0; j < m; j++) {
  337.             if (!(mask & (1 << j))) {
  338.                 int newMask = mask | (1 << j);
  339.                 dp[newMask] = min(dp[newMask], dp[mask] + cost[worker][j]);
  340.             }
  341.         }
  342.     }
  343.  
  344.     ll ans = INF;
  345.     for (int mask = 0; mask < (1 << m); mask++) {
  346.         if (__builtin_popcount(mask) == n) {
  347.             ans = min(ans, dp[mask]);
  348.         }
  349.     }
  350.     return ans;
  351. }
  352.  
  353. // 4.2) Assignment Bitmask DP (Maximum Cost).
  354. //      Just negate the costs inside.
  355. ll assignmentBitmaskMaxDP(const vector<vector<ll>>& cost) {
  356.     int n = (int)cost.size();
  357.     int m = (int)cost[0].size();
  358.     vector<vector<ll>> negCost(n, vector<ll>(m));
  359.     for (int i = 0; i < n; i++) {
  360.         for (int j = 0; j < m; j++) {
  361.             negCost[i][j] = -cost[i][j];
  362.         }
  363.     }
  364.     return -assignmentBitmaskDP(negCost);
  365. }
  366.  
  367. // ===================================================================
  368. // 5) Kuhn's Algorithm (Maximum Bipartite Matching)
  369. //    This finds the maximum number of pairs we can match in a 
  370. //    bipartite graph. 
  371. //    It is used as a helper for the "Minimize Maximum Cost" trick.
  372. // ===================================================================
  373.  
  374. // 5.1) Kuhn's Algorithm (DFS-based) to find maximum matching.
  375. //      Parameters:
  376. //        - adj: adjacency list of the left side (size n).
  377. //               adj[i] contains the list of right-side nodes (0-based) 
  378. //               that left node i can connect to.
  379. //        - n: number of nodes on the left.
  380. //        - m: number of nodes on the right.
  381. //      Returns:
  382. //        - The maximum number of edges in the matching.
  383. //      Time complexity: O(n * E) where E is the total number of edges.
  384. //        - In practice, very fast for sparse graphs.
  385. //      Constraints:
  386. //        - Works for any bipartite graph.
  387. //        - If you have both sides <= 500, it runs fine.
  388. //      Note:
  389. //        - This function does NOT return the actual matching edges, 
  390. //          only the count. (You can modify it to return the matching array 
  391. //          if needed, but for the binary-search trick below, we only 
  392. //          need the count).
  393. int kuhnMatchingCount(const vector<vector<int>>& adj, int n, int m) {
  394.     vector<int> matchR(m, -1); // matchR[j] = which left node is matched to right j
  395.     vector<int> visited;
  396.  
  397.     function<bool(int)> dfs = [&](int u) {
  398.         for (int v : adj[u]) {
  399.             if (visited[v]) continue;
  400.             visited[v] = 1;
  401.             if (matchR[v] == -1 || dfs(matchR[v])) {
  402.                 matchR[v] = u;
  403.                 return true;
  404.             }
  405.         }
  406.         return false;
  407.     };
  408.  
  409.     int matching = 0;
  410.     for (int u = 0; u < n; u++) {
  411.         visited.assign(m, 0);
  412.         if (dfs(u)) matching++;
  413.     }
  414.     return matching;
  415. }
  416.  
  417. // ===================================================================
  418. // 6) Advanced Trick: Minimize the Maximum Cost (Minimax Assignment)
  419. //    Problem: We want to assign every worker to a unique job, but we 
  420. //    want the maximum cost among all selected edges to be as small 
  421. //    as possible (instead of minimizing the sum).
  422. //    This is solved by Binary Search on the answer + Kuhn Matching.
  423. // ===================================================================
  424.  
  425. // 6.1) Check if we can assign all workers with every edge cost <= limit.
  426. //      Parameters:
  427. //        - cost: n x m matrix (n <= m).
  428. //        - limit: the maximum allowed cost for any chosen edge.
  429. //      Returns:
  430. //        - true if a perfect matching exists using only edges with cost <= limit.
  431. //      Time complexity: O(n * E) for each call.
  432. //      Note:
  433. //        - This is a feasibility check used inside binary search.
  434. bool canAssignWithMaxCost(const vector<vector<ll>>& cost, ll limit) {
  435.     int n = (int)cost.size();
  436.     int m = (int)cost[0].size();
  437.     if (n > m) return false; // cannot assign more workers than jobs.
  438.  
  439.     vector<vector<int>> adj(n);
  440.     for (int i = 0; i < n; i++) {
  441.         for (int j = 0; j < m; j++) {
  442.             if (cost[i][j] <= limit) {
  443.                 adj[i].push_back(j);
  444.             }
  445.         }
  446.     }
  447.     int maxMatch = kuhnMatchingCount(adj, n, m);
  448.     return maxMatch == n; // we matched every worker
  449. }
  450.  
  451. // 6.2) Find the minimum possible maximum cost.
  452. //      Parameters:
  453. //        - cost: n x m matrix.
  454. //      Returns:
  455. //        - The minimum value X such that we can assign all workers using 
  456. //          only edges with cost <= X.
  457. //      Time complexity: O(log(range) * n * E).
  458. //      Constraints:
  459. //        - Costs can be negative. We handle that by taking the min/max 
  460. //          of the matrix as the binary search bounds.
  461. //        - Assumes at least one valid assignment exists (if not, returns INF).
  462. ll minMaxAssignmentCost(const vector<vector<ll>>& cost) {
  463.     int n = (int)cost.size();
  464.     int m = (int)cost[0].size();
  465.     if (n > m) return INF;
  466.  
  467.     ll low = INF, high = -INF;
  468.     for (int i = 0; i < n; i++) {
  469.         for (int j = 0; j < m; j++) {
  470.             low = min(low, cost[i][j]);
  471.             high = max(high, cost[i][j]);
  472.         }
  473.     }
  474.  
  475.     // If n == 0, return 0.
  476.     if (n == 0) return 0;
  477.  
  478.     ll ans = high;
  479.     while (low <= high) {
  480.         ll mid = low + (high - low) / 2;
  481.         if (canAssignWithMaxCost(cost, mid)) {
  482.             ans = mid;
  483.             high = mid - 1;
  484.         } else {
  485.             low = mid + 1;
  486.         }
  487.     }
  488.     return ans;
  489. }
  490.  
  491. // ===================================================================
  492. // 7) TRICKS & EXTRA PATTERNS (FULLY IMPLEMENTED FUNCTIONS)
  493. //    These are the actual standalone functions for the common 
  494. //    assignment‑problem variations that appear in ECPC / ACPC.
  495. // ===================================================================
  496.  
  497. // ===================================================================
  498. // 7.1) Convert any rectangular matrix into a SQUARE matrix by adding 
  499. //      dummy rows or columns with cost ZERO.
  500. //      This lets you use the standard square Hungarian algorithm 
  501. //      (hungarianMinCost) even when N != M.
  502. // ===================================================================
  503.  
  504. // PURPOSE:
  505. //   Takes an N x M cost matrix and returns a square matrix (size K x K)
  506. //   where K = max(N, M). 
  507. //   - If N > M, it adds (N - M) dummy columns (jobs) with cost 0.
  508. //   - If M > N, it adds (M - N) dummy rows (workers) with cost 0.
  509. //   - If N == M, it returns a copy of the original.
  510. //
  511. // INPUT:
  512. //   cost : a 2D vector (N rows, M columns) of long long.
  513. //
  514. // OUTPUT:
  515. //   Returns a square vector<vector<ll>> of size K x K.
  516. //   Dummy rows/columns are placed at the end.
  517. //
  518. // TIME COMPLEXITY:
  519. //   O(N * M) to copy the original, plus O(K^2) overall.
  520. //
  521. // CONSTRAINTS / PRECONDITIONS:
  522. //   - None. Works for any N, M >= 0.
  523. //   - Costs can be negative, zero, or positive.
  524. //
  525. // NOTES:
  526. //   - After getting the square matrix, you can call 
  527. //     `hungarianMinCost(squareMatrix)` to get the optimal assignment.
  528. //   - If the original matrix had more rows (workers) than columns (jobs),
  529. //     adding dummy jobs (cost 0) means that extra workers will be 
  530. //     assigned to dummy jobs, i.e., they do no real work.
  531. //   - If the original had more columns (jobs) than rows (workers),
  532. //     adding dummy workers (cost 0) means that extra jobs will be 
  533. //     assigned to dummy workers, i.e., they remain unassigned.
  534. vector<vector<ll>> makeSquareMatrixForAssignment(const vector<vector<ll>>& cost) {
  535.     int n = (int)cost.size();
  536.     if (n == 0) return {}; // empty matrix
  537.     int m = (int)cost[0].size();
  538.  
  539.     int k = max(n, m);
  540.     vector<vector<ll>> sq(k, vector<ll>(k, 0));
  541.  
  542.     for (int i = 0; i < n; i++) {
  543.         for (int j = 0; j < m; j++) {
  544.             sq[i][j] = cost[i][j];
  545.         }
  546.     }
  547.     // The remaining entries are already 0 (dummy rows/columns).
  548.     return sq;
  549. }
  550.  
  551. // ===================================================================
  552. // 7.2) Assignment where you are ALLOWED to leave some workers 
  553. //      unassigned, but you pay a fixed PENALTY for each unassigned worker.
  554. //      This is extremely common in contest problems.
  555. // ===================================================================
  556.  
  557. // PURPOSE:
  558. //   Solves the assignment problem for N workers and M jobs.
  559. //   Each worker must be assigned to at most ONE job.
  560. //   If a worker is NOT assigned to any real job, you pay a fixed 
  561. //   penalty `penalty` for that worker.
  562. //   The goal is to minimize (total assignment cost + total penalties).
  563. //
  564. // INPUT:
  565. //   cost    : N x M matrix (long long). cost[i][j] is the cost of 
  566. //             assigning worker i to job j.
  567. //   penalty : a long long value (the cost per unassigned worker).
  568. //             Can be negative (if you really don't want to assign someone),
  569. //             but usually it is a positive number.
  570. //
  571. // OUTPUT:
  572. //   Returns a pair <total_cost, assignment>.
  573. //   - total_cost (long long): the minimum total cost (assignment + penalties).
  574. //   - assignment (vector<int>): size N. 
  575. //       assignment[i] = j (where 0 <= j < M) means worker i is assigned 
  576. //                        to real job j.
  577. //       assignment[i] = M + k (where k >= 0) means worker i is assigned 
  578. //                        to a dummy job, i.e., worker i is unassigned.
  579. //
  580. // TIME COMPLEXITY:
  581. //   O(N^2 * (M + N)) if you use the rectangular Hungarian, 
  582. //   but typically you will use hungarianMinCostRectangular which is 
  583. //   O(N^2 * M') where M' = M + N.
  584. //
  585. // CONSTRAINTS / PRECONDITIONS:
  586. //   - The number of real jobs M can be less than, equal to, or greater 
  587. //     than N. The function creates dummy jobs so that there are always 
  588. //     enough jobs for every worker.
  589. //   - All costs must fit in long long.
  590. //   - The returned assignment indices >= M indicate unassigned workers.
  591. //
  592. // NOTES:
  593. //   - Internally, it creates a new matrix of size N x (M + N). 
  594. //     The last N columns are dummy jobs, all having cost = penalty.
  595. //     It then calls hungarianMinCostRectangular on this matrix.
  596. //   - If you just want to allow unassigned workers without penalty, 
  597. //     set penalty = 0.
  598. pair<ll, vector<int>> assignmentWithUnassignedPenalty(
  599.     const vector<vector<ll>>& cost,
  600.     ll penalty
  601. ) {
  602.     int n = (int)cost.size();    // workers
  603.     int m = (int)cost[0].size(); // real jobs
  604.  
  605.     // New matrix: n workers, (m + n) jobs (n dummy jobs)
  606.     vector<vector<ll>> newCost(n, vector<ll>(m + n));
  607.  
  608.     for (int i = 0; i < n; i++) {
  609.         for (int j = 0; j < m; j++) {
  610.             newCost[i][j] = cost[i][j];
  611.         }
  612.         for (int j = m; j < m + n; j++) {
  613.             newCost[i][j] = penalty; // dummy job cost = penalty
  614.         }
  615.     }
  616.  
  617.     // Use the rectangular Hungarian (n rows, m+n columns, and n <= m+n)
  618.     auto res = hungarianMinCostRectangular(newCost);
  619.     return res; // assignment indices go up to m+n-1
  620. }
  621.  
  622. // ===================================================================
  623. // 7.3) Count the number of perfect matchings (total number of ways 
  624. //      to assign all workers to unique jobs) using Bitmask DP.
  625. //      This is useful when the number of jobs (or workers) is small 
  626. //      (typically <= 20).
  627. // ===================================================================
  628.  
  629. // PURPOSE:
  630. //   Given an N x M matrix indicating which edges are allowed, 
  631. //   count how many ways we can assign every worker (0..N-1) to a 
  632. //   distinct job (0..M-1) using only allowed edges.
  633. //   Every worker must get exactly one job, and no two workers share a job.
  634. //
  635. // INPUT:
  636. //   allowed : N x M matrix of integers (or booleans).
  637. //             allowed[i][j] = 1 (or true) means worker i CAN take job j.
  638. //             allowed[i][j] = 0 (or false) means it is forbidden.
  639. //
  640. // OUTPUT:
  641. //   Returns a long long integer: the total number of valid assignments.
  642. //   If N > M, returns 0 immediately (not enough jobs).
  643. //   If the count exceeds 2^63 - 1, it will overflow; use __int128 
  644. //   if you need bigger numbers, but that is rare for N <= 20.
  645. //
  646. // TIME COMPLEXITY:
  647. //   O(N * 2^M) where M is the number of columns (jobs).
  648. //   Therefore, M must be small (<= 20 for typical time limits).
  649. //
  650. // CONSTRAINTS / PRECONDITIONS:
  651. //   - M (number of jobs) should be <= 20 (or 22 with optimizations).
  652. //   - N <= M, otherwise 0 is returned.
  653. //   - Works for any allowed matrix (no need for costs).
  654. //
  655. // NOTES:
  656. //   - This uses the standard subset DP: dp[mask] stores the number of 
  657. //     ways to assign the first `popcount(mask)` workers to the jobs 
  658. //     represented by `mask`.
  659. //   - If you need to count matchings where not all workers must be assigned,
  660. //     you can sum dp[mask] over all masks (but this function counts 
  661. //     perfect matchings only, i.e., all N workers assigned).
  662. long long countPerfectMatchings(const vector<vector<int>>& allowed) {
  663.     int n = (int)allowed.size();    // workers
  664.     if (n == 0) return 1;           // empty case
  665.     int m = (int)allowed[0].size(); // jobs
  666.  
  667.     if (n > m) return 0; // cannot assign all workers
  668.  
  669.     vector<long long> dp(1 << m, 0);
  670.     dp[0] = 1;
  671.  
  672.     for (int mask = 0; mask < (1 << m); mask++) {
  673.         int worker = __builtin_popcount(mask); // how many workers assigned so far
  674.         if (worker == n) continue; // all workers are already assigned
  675.  
  676.         for (int j = 0; j < m; j++) {
  677.             if (mask & (1 << j)) continue; // job j already taken
  678.             if (!allowed[worker][j]) continue; // edge is forbidden
  679.  
  680.             int newMask = mask | (1 << j);
  681.             dp[newMask] += dp[mask];
  682.         }
  683.     }
  684.  
  685.     long long ans = 0;
  686.     for (int mask = 0; mask < (1 << m); mask++) {
  687.         if (__builtin_popcount(mask) == n) {
  688.             ans += dp[mask];
  689.         }
  690.     }
  691.     return ans;
  692. }
  693.  
  694. // ===================================================================
  695. // 7.4) Minimum cost assignment with FORBIDDEN edges.
  696. //      Some edges (i, j) are not allowed to be chosen.
  697. //      This function builds the cost matrix by setting forbidden 
  698. //      edges to INF, runs the Hungarian algorithm, and returns the result.
  699. //      If no complete assignment exists, it returns INF and an empty 
  700. //      assignment vector.
  701. // ===================================================================
  702.  
  703. // PURPOSE:
  704. //   Solves the assignment problem (N workers, M jobs, N <= M) where 
  705. //   certain edges are FORBIDDEN (cannot be used). 
  706. //   It returns the minimum total cost using only allowed edges, 
  707. //   or signals impossibility if no perfect assignment exists.
  708. //
  709. // INPUT:
  710. //   cost    : N x M matrix (long long) containing the regular costs.
  711. //   allowed : N x M matrix of integers (or booleans).
  712. //             allowed[i][j] = 1 (or true)  -> edge is allowed.
  713. //             allowed[i][j] = 0 (or false) -> edge is forbidden.
  714. //
  715. // OUTPUT:
  716. //   Returns a pair <total_cost, assignment>.
  717. //   - If a valid assignment exists:
  718. //       total_cost is the minimum cost.
  719. //       assignment is a vector of size N (assignment[i] = j).
  720. //   - If NO valid assignment exists (because forbidden edges block it):
  721. //       total_cost = INF (the global constant).
  722. //       assignment is an empty vector ({}).
  723. //
  724. // TIME COMPLEXITY:
  725. //   O(N^2 * M) (calls the rectangular Hungarian internally).
  726. //
  727. // CONSTRAINTS / PRECONDITIONS:
  728. //   - The input matrix must have N <= M (otherwise impossible to assign 
  729. //     all workers – returns {INF, {}}).
  730. //   - Costs can be negative, zero, or positive.
  731. //   - The INF constant is defined globally (usually 4e18 or (1LL<<60)).
  732. //   - The total cost of a valid assignment must be < INF/2 to be 
  733. //     considered valid.
  734. //
  735. // NOTES:
  736. //   - Internally, it builds a new matrix where forbidden edges are 
  737. //     replaced with INF (a very large number). 
  738. //   - It then calls `hungarianMinCostRectangular`. 
  739. //   - If the returned cost is >= INF/2, we conclude that no feasible 
  740. //     assignment exists.
  741. pair<ll, vector<int>> assignmentWithForbiddenEdges(
  742.     const vector<vector<ll>>& cost,
  743.     const vector<vector<int>>& allowed
  744. ) {
  745.     int n = (int)cost.size();
  746.     if (n == 0) return {0, {}};
  747.     int m = (int)cost[0].size();
  748.  
  749.     if (n > m) {
  750.         return {INF, {}}; // not enough jobs for all workers
  751.     }
  752.  
  753.     vector<vector<ll>> modifiedCost = cost;
  754.     for (int i = 0; i < n; i++) {
  755.         for (int j = 0; j < m; j++) {
  756.             if (!allowed[i][j]) {
  757.                 modifiedCost[i][j] = INF; // forbid this edge
  758.             }
  759.         }
  760.     }
  761.  
  762.     auto res = hungarianMinCostRectangular(modifiedCost);
  763.  
  764.     // If the total cost is too large, it means the algorithm was forced 
  765.     // to pick at least one forbidden edge (or the problem is infeasible).
  766.     if (res.first >= INF / 2) {
  767.         return {INF, {}};
  768.     }
  769.     return res;
  770. }
  771.  
  772. // ===================================================================
  773. // End of Section 7 functions.
  774. // ===================================================================
  775.  
  776. // ===================================================================
  777. // main() with example usage (you can ignore this part)
  778. // ===================================================================
  779.  
  780. int main() {
  781.     ios::sync_with_stdio(false);
  782.     cin.tie(nullptr);
  783.  
  784.     // Example 1: Square matrix min cost
  785.     vector<vector<ll>> cost1 = {
  786.         {4, 1, 3},
  787.         {2, 0, 5},
  788.         {3, 2, 2}
  789.     };
  790.     auto res1_min = hungarianMinCost(cost1);
  791.     cout << "Min Cost: " << res1_min.first << "\n";
  792.     cout << "Assignment: ";
  793.     for (int x : res1_min.second) cout << x << " ";
  794.     cout << "\n";
  795.  
  796.     // Example 2: Rectangular matrix (3 workers, 5 jobs)
  797.     vector<vector<ll>> cost2 = {
  798.         {1, 2, 3, 4, 5},
  799.         {5, 4, 3, 2, 1},
  800.         {2, 3, 4, 5, 6}
  801.     };
  802.     auto res2_rect = hungarianMinCostRectangular(cost2);
  803.     cout << "Rect Min Cost: " << res2_rect.first << "\n";
  804.  
  805.     // Example 3: Minimize Maximum Cost
  806.     vector<vector<ll>> cost3 = {
  807.         {1, 100, 100},
  808.         {100, 1, 100},
  809.         {100, 100, 1}
  810.     };
  811.     cout << "Min Max Cost: " << minMaxAssignmentCost(cost3) << "\n"; // Output: 1
  812.  
  813.     // Example 4: Bitmask DP (small N)
  814.     vector<vector<ll>> cost4 = {
  815.         {10, 20, 30},
  816.         {40, 50, 60},
  817.         {70, 80, 90}
  818.     };
  819.     cout << "Bitmask DP Min: " << assignmentBitmaskDP(cost4) << "\n";
  820.  
  821.     // Example 5: Penalty for unassigned workers
  822.     vector<vector<ll>> cost5 = {
  823.         {1, 2},
  824.         {3, 4},
  825.         {5, 6} // 3 workers, 2 jobs
  826.     };
  827.     ll penalty = 10;
  828.     auto res5 = assignmentWithUnassignedPenalty(cost5, penalty);
  829.     cout << "With penalty cost: " << res5.first << "\n";
  830.     cout << "Assignment: ";
  831.     for (int x : res5.second) cout << x << " ";
  832.     cout << "\n";
  833.  
  834.     // Example 6: Make square matrix and use standard square Hungarian
  835.     auto sq = makeSquareMatrixForAssignment(cost5); // 3x3
  836.     auto res6 = hungarianMinCost(sq);
  837.     cout << "Square Hungarian on padded matrix: " << res6.first << "\n";
  838.  
  839.     // Example 7: Count perfect matchings (only allowed edges)
  840.     vector<vector<int>> allowed = {
  841.         {1, 1},
  842.         {1, 0}, // worker 1 cannot take job 1
  843.         {1, 1}
  844.     };
  845.     long long ways = countPerfectMatchings(allowed);
  846.     cout << "Number of perfect matchings: " << ways << "\n";
  847.  
  848.     // Example 8: Forbidden edges
  849.     auto res7 = assignmentWithForbiddenEdges(cost5, allowed);
  850.     if (res7.first == INF) {
  851.         cout << "No feasible assignment!\n";
  852.     } else {
  853.         cout << "With forbidden edges, cost: " << res7.first << "\n";
  854.     }
  855.  
  856.     return 0;
  857. }
Success #stdin #stdout 0s 5324KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Min Cost: 5
Assignment: 1 0 2 
Rect Min Cost: 5
Min Max Cost: 1
Bitmask DP Min: 150
With penalty cost: 15
Assignment: 0 1 2 
Square Hungarian on padded matrix: 5
Number of perfect matchings: 0
No feasible assignment!
https://ideone.com/VYLaTS
language:
C++ (gcc 8.3)
created:
1 day ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!