np_algorithms.cpp

// np_algorithms.cpp
// =====================================================================
// NP-Complete / NP-Hard 問題的實戰解法示範（C++17）
//
// 涵蓋三大類技巧：
//   A. 精確解（指數時間，但小規模實用）
//      1. N-Queens          —— 回溯 + 剪枝
//      2. DPLL SAT Solver   —— 單元傳播 + 回溯（現代 CDCL 的雛形）
//      3. Subset Sum        —— 偽多項式 DP + Meet-in-the-Middle
//      4. 0/1 Knapsack      —— 偽多項式 DP
//      5. Graph Coloring    —— 回溯 + 剪枝
//      6. Hamiltonian Path  —— 位元壓縮 DP  O(n^2 * 2^n)
//      7. TSP (Held-Karp)   —— 位元壓縮 DP  O(n^2 * 2^n)
//      8. Vertex Cover      —— 分支限界（參數化 FPT 思想）
//   B. 近似演算法（多項式時間 + 品質保證）
//      9. Vertex Cover 2-近似（取極大匹配）
//     10. Set Cover 貪婪 ln(n)-近似
//     11. Bin Packing First-Fit Decreasing（<= 11/9 OPT + 1）
//     12. Metric TSP：最近鄰啟發 + 2-opt 局部搜尋
//   C. 特例回到 P（結構讓問題變簡單）
//     13. 2-SAT —— 用 SCC（Tarjan）在線性時間求解
//
// Compile: g++ -std=c++17 -O2 -Wall -Wextra -o np_algorithms np_algorithms.cpp
// Run:     ./np_algorithms
// =====================================================================

#include <iostream>
#include <vector>
#include <algorithm>
#include <numeric>
#include <cmath>
#include <climits>
#include <random>
#include <functional>
#include <iomanip>

using namespace std;

// =====================================================================
// 1. N-Queens —— 回溯 + 剪枝（用三個布林陣列 O(1) 檢查衝突）
// =====================================================================
class NQueens {
    int n;
    long long solutions = 0;
    vector<bool> col, diag1, diag2;
    void solve(int row) {
        if (row == n) { ++solutions; return; }
        for (int c = 0; c < n; ++c) {
            if (col[c] || diag1[row + c] || diag2[row - c + n - 1]) continue;
            col[c] = diag1[row + c] = diag2[row - c + n - 1] = true;
            solve(row + 1);                       // 只往合法分支走（剪枝）
            col[c] = diag1[row + c] = diag2[row - c + n - 1] = false;
        }
    }
public:
    long long count(int n_) {
        n = n_; solutions = 0;
        col.assign(n, false);
        diag1.assign(2 * n - 1, false);
        diag2.assign(2 * n - 1, false);
        solve(0);
        return solutions;
    }
};

// =====================================================================
// 2. DPLL SAT Solver —— CNF 可滿足性判定
//    文字編碼：變數 v 的正文字 = +v，負文字 = -v（v 由 1 起算）
// =====================================================================
class DPLL {
    int numVars;
    vector<vector<int>> clauses;

    // 對部分賦值 assign（0=未定, 1=真, -1=假）做單元傳播
    // 回傳：1 = 全部子句已滿足, -1 = 出現空子句（矛盾）, 0 = 未定
    int propagate(vector<int>& assign) {
        bool changed = true;
        while (changed) {
            changed = false;
            bool allSat = true;
            for (const auto& cl : clauses) {
                int unassigned = 0, lastLit = 0;
                bool sat = false;
                for (int lit : cl) {
                    int v = abs(lit), want = lit > 0 ? 1 : -1;
                    if (assign[v] == want) { sat = true; break; }
                    if (assign[v] == 0) { ++unassigned; lastLit = lit; }
                }
                if (sat) continue;
                allSat = false;
                if (unassigned == 0) return -1;          // 空子句 -> 矛盾
                if (unassigned == 1) {                   // 單元子句 -> 強制賦值
                    assign[abs(lastLit)] = lastLit > 0 ? 1 : -1;
                    changed = true;
                }
            }
            if (allSat) return 1;
        }
        return 0;
    }

    bool dfs(vector<int> assign) {
        int st = propagate(assign);
        if (st == 1) { model = assign; return true; }
        if (st == -1) return false;
        int v = 1;
        while (assign[v] != 0) ++v;                      // 挑第一個未賦值變數
        for (int val : {1, -1}) {                        // 先試真、再試假
            assign[v] = val;
            if (dfs(assign)) return true;
        }
        return false;
    }
public:
    vector<int> model;                                   // 求出的賦值
    bool solve(int nv, const vector<vector<int>>& cls) {
        numVars = nv; clauses = cls;
        return dfs(vector<int>(nv + 1, 0));
    }
};

// =====================================================================
// 3. Subset Sum
// =====================================================================
// 3a. 偽多項式 DP：O(n * target)
bool subsetSumDP(const vector<int>& a, int target) {
    vector<bool> dp(target + 1, false);
    dp[0] = true;
    for (int x : a)
        for (int s = target; s >= x; --s)
            if (dp[s - x]) dp[s] = true;
    return dp[target];
}

// 3b. Meet-in-the-Middle：O(2^(n/2) * n)，值域大時的利器
bool subsetSumMITM(const vector<int>& a, long long target) {
    int n = (int)a.size(), h = n / 2;
    auto enumerate = [](const vector<int>& part) {
        vector<long long> sums;
        int m = (int)part.size();
        for (int mask = 0; mask < (1 << m); ++mask) {
            long long s = 0;
            for (int i = 0; i < m; ++i)
                if (mask & (1 << i)) s += part[i];
            sums.push_back(s);
        }
        return sums;
    };
    vector<int> L(a.begin(), a.begin() + h), R(a.begin() + h, a.end());
    vector<long long> sl = enumerate(L), sr = enumerate(R);
    sort(sr.begin(), sr.end());
    for (long long s : sl)
        if (binary_search(sr.begin(), sr.end(), target - s)) return true;
    return false;
}

// =====================================================================
// 4. 0/1 Knapsack —— 偽多項式 DP：O(n * W)
// =====================================================================
int knapsack01(const vector<int>& w, const vector<int>& v, int W) {
    vector<int> dp(W + 1, 0);
    for (size_t i = 0; i < w.size(); ++i)
        for (int cap = W; cap >= w[i]; --cap)
            dp[cap] = max(dp[cap], dp[cap - w[i]] + v[i]);
    return dp[W];
}

// =====================================================================
// 5. Graph Coloring —— 回溯：能否用 k 色著色
// =====================================================================
class GraphColoring {
    const vector<vector<int>>& adj;
    int n, k;
    vector<int> color;
    bool ok(int u, int c) {
        for (int v : adj[u]) if (color[v] == c) return false;
        return true;
    }
    bool dfs(int u) {
        if (u == n) return true;
        for (int c = 0; c < k; ++c) {
            if (!ok(u, c)) continue;                    // 剪枝：衝突不往下走
            color[u] = c;
            if (dfs(u + 1)) return true;
            color[u] = -1;
        }
        return false;
    }
public:
    GraphColoring(const vector<vector<int>>& g) : adj(g), n((int)g.size()) {}
    bool colorable(int k_) {
        k = k_; color.assign(n, -1);
        return dfs(0);
    }
};

// =====================================================================
// 6. Hamiltonian Path —— 位元壓縮 DP：O(n^2 * 2^n)
//    dp[mask][v] = 走過集合 mask、目前停在 v，是否可行
// =====================================================================
bool hamiltonianPath(const vector<vector<int>>& adj) {
    int n = (int)adj.size();
    if (n == 0) return false;
    if (n == 1) return true;
    vector<vector<char>> dp(1 << n, vector<char>(n, 0));
    for (int v = 0; v < n; ++v) dp[1 << v][v] = 1;      // 任意起點
    for (int mask = 1; mask < (1 << n); ++mask)
        for (int v = 0; v < n; ++v) {
            if (!dp[mask][v]) continue;
            for (int u : adj[v])
                if (!(mask & (1 << u)))
                    dp[mask | (1 << u)][u] = 1;
        }
    for (int v = 0; v < n; ++v)
        if (dp[(1 << n) - 1][v]) return true;
    return false;
}

// =====================================================================
// 7. TSP —— Held-Karp 位元壓縮 DP：O(n^2 * 2^n)，n<=20 實用
// =====================================================================
long long tspHeldKarp(const vector<vector<int>>& dist) {
    int n = (int)dist.size();
    const long long INF = LLONG_MAX / 4;
    vector<vector<long long>> dp(1 << n, vector<long long>(n, INF));
    dp[1][0] = 0;                                        // 從城市 0 出發
    for (int mask = 1; mask < (1 << n); ++mask) {
        if (!(mask & 1)) continue;                       // 必含起點
        for (int v = 0; v < n; ++v) {
            if (dp[mask][v] == INF) continue;
            for (int u = 0; u < n; ++u)
                if (!(mask & (1 << u)))
                    dp[mask | (1 << u)][u] =
                        min(dp[mask | (1 << u)][u], dp[mask][v] + dist[v][u]);
        }
    }
    long long best = INF;
    for (int v = 1; v < n; ++v)
        if (dp[(1 << n) - 1][v] != INF)
            best = min(best, dp[(1 << n) - 1][v] + dist[v][0]);  // 回到起點
    return best;
}

// =====================================================================
// 8. Vertex Cover —— 分支限界（FPT 思想）：cover 大小 <= k 是否可行？
//    關鍵觀察：任一條邊 (u,v)，u 或 v 至少要選一個 -> 二分支
// =====================================================================
bool vertexCoverK(vector<vector<int>> adj, int k) {
    // 找任一條還存在的邊
    int u = -1, v = -1;
    for (int i = 0; i < (int)adj.size() && u == -1; ++i)
        if (!adj[i].empty()) { u = i; v = adj[i][0]; }
    if (u == -1) return true;                            // 沒邊了 -> 已覆蓋
    if (k == 0) return false;                            // 還有邊但配額用完

    auto removeVertex = [](vector<vector<int>> g, int x) {
        for (int y : g[x]) {
            auto& lst = g[y];
            lst.erase(remove(lst.begin(), lst.end(), x), lst.end());
        }
        g[x].clear();
        return g;
    };
    return vertexCoverK(removeVertex(adj, u), k - 1)     // 選 u
        || vertexCoverK(removeVertex(adj, v), k - 1);    // 選 v
}

// =====================================================================
// 9. Vertex Cover 2-近似 —— 取「極大匹配」的兩端點
//    保證：|解| <= 2 * OPT（每條匹配邊至少要 1 個點，我們取了 2 個）
// =====================================================================
vector<int> vertexCover2Approx(int n, vector<pair<int,int>> edges) {
    vector<bool> inCover(n, false);
    vector<int> cover;
    for (auto [u, v] : edges) {
        if (inCover[u] || inCover[v]) continue;          // 邊已被覆蓋
        inCover[u] = inCover[v] = true;                  // 兩端點都收
        cover.push_back(u); cover.push_back(v);
    }
    return cover;
}

// =====================================================================
// 10. Set Cover 貪婪近似 —— 每輪挑「新覆蓋元素最多」的集合
//     保證：|解| <= ln(n) * OPT + 1（幾乎是能做到的最佳近似比）
// =====================================================================
vector<int> setCoverGreedy(int universe, const vector<vector<int>>& sets) {
    vector<bool> covered(universe, false);
    int remaining = universe;
    vector<int> chosen;
    while (remaining > 0) {
        int best = -1, bestGain = 0;
        for (int i = 0; i < (int)sets.size(); ++i) {
            int gain = 0;
            for (int x : sets[i]) if (!covered[x]) ++gain;
            if (gain > bestGain) { bestGain = gain; best = i; }
        }
        if (best == -1) break;                           // 無法覆蓋
        chosen.push_back(best);
        for (int x : sets[best])
            if (!covered[x]) { covered[x] = true; --remaining; }
    }
    return chosen;
}

// =====================================================================
// 11. Bin Packing —— First-Fit Decreasing
//     保證：使用箱數 <= 11/9 * OPT + 1
// =====================================================================
int binPackingFFD(vector<double> items, double capacity) {
    sort(items.rbegin(), items.rend());                  // 由大到小
    vector<double> bins;                                 // 每箱剩餘空間
    for (double it : items) {
        bool placed = false;
        for (double& room : bins)
            if (room >= it - 1e-12) { room -= it; placed = true; break; }
        if (!placed) bins.push_back(capacity - it);      // 開新箱
    }
    return (int)bins.size();
}

// =====================================================================
// 12. Metric TSP 啟發式 —— 最近鄰建構 + 2-opt 改善
// =====================================================================
static double tourLength(const vector<int>& tour,
                         const vector<vector<double>>& d) {
    double len = 0;
    int n = (int)tour.size();
    for (int i = 0; i < n; ++i) len += d[tour[i]][tour[(i + 1) % n]];
    return len;
}

vector<int> tspNearestNeighbor(const vector<vector<double>>& d) {
    int n = (int)d.size();
    vector<bool> used(n, false);
    vector<int> tour = {0};
    used[0] = true;
    for (int step = 1; step < n; ++step) {
        int cur = tour.back(), best = -1;
        for (int v = 0; v < n; ++v)
            if (!used[v] && (best == -1 || d[cur][v] < d[cur][best])) best = v;
        used[best] = true;
        tour.push_back(best);
    }
    return tour;
}

// 2-opt：反覆嘗試「反轉一段路徑」，只要變短就採用
void tsp2opt(vector<int>& tour, const vector<vector<double>>& d) {
    int n = (int)tour.size();
    bool improved = true;
    while (improved) {
        improved = false;
        for (int i = 0; i < n - 1; ++i)
            for (int j = i + 2; j < n; ++j) {
                if (i == 0 && j == n - 1) continue;      // 同一條邊
                int a = tour[i], b = tour[i + 1];
                int c = tour[j], e = tour[(j + 1) % n];
                double delta = d[a][c] + d[b][e] - d[a][b] - d[c][e];
                if (delta < -1e-10) {
                    reverse(tour.begin() + i + 1, tour.begin() + j + 1);
                    improved = true;
                }
            }
    }
}

// =====================================================================
// 13. 2-SAT —— 特例回到 P：蘊含圖 + Tarjan SCC，O(n + m)
//     子句 (a ∨ b) 等價於 ¬a→b 與 ¬b→a
// =====================================================================
class TwoSAT {
    int n;                                // 變數個數
    vector<vector<int>> g;                // 蘊含圖，節點 2v=正、2v+1=負
    vector<int> comp, order;
    vector<bool> visited;

    void dfs1(int u) {
        visited[u] = true;
        for (int v : g[u]) if (!visited[v]) dfs1(v);
        order.push_back(u);
    }
    void dfs2(int u, int c, const vector<vector<int>>& gr) {
        comp[u] = c;
        for (int v : gr[u]) if (comp[v] == -1) dfs2(v, c, gr);
    }
public:
    vector<bool> value;                   // 求出的賦值
    TwoSAT(int vars) : n(vars), g(2 * vars) {}

    // 加入子句 (x_i = vi) ∨ (x_j = vj)，vi/vj 為 true 表示正文字
    void addClause(int i, bool vi, int j, bool vj) {
        int a = 2 * i + (vi ? 0 : 1);     // 文字 a
        int b = 2 * j + (vj ? 0 : 1);     // 文字 b
        g[a ^ 1].push_back(b);            // ¬a -> b
        g[b ^ 1].push_back(a);            // ¬b -> a
    }

    bool solve() {                        // Kosaraju 兩遍 DFS 求 SCC
        int N = 2 * n;
        visited.assign(N, false);
        order.clear();
        for (int i = 0; i < N; ++i) if (!visited[i]) dfs1(i);
        vector<vector<int>> gr(N);
        for (int u = 0; u < N; ++u) for (int v : g[u]) gr[v].push_back(u);
        comp.assign(N, -1);
        int c = 0;
        for (int i = N - 1; i >= 0; --i)
            if (comp[order[i]] == -1) dfs2(order[i], c++, gr);
        value.assign(n, false);
        for (int v = 0; v < n; ++v) {
            if (comp[2 * v] == comp[2 * v + 1]) return false; // x 與 ¬x 同分量
            value[v] = comp[2 * v] > comp[2 * v + 1];
        }
        return true;
    }
};

// =====================================================================
// Demonstrations
// =====================================================================
static void hr(const string& title) {
    cout << "\n" << string(62, '-') << "\n" << title << "\n";
}

int main() {
    cout << fixed << setprecision(3);

    // ---- 1. N-Queens ----
    hr("1. N-Queens (backtracking + pruning)");
    NQueens nq;
    for (int n : {4, 8, 10})
        cout << "  n=" << n << " -> " << nq.count(n) << " solutions\n";

    // ---- 2. DPLL SAT ----
    hr("2. DPLL SAT solver");
    {
        // (x1 v x2) & (~x1 v x3) & (~x2 v ~x3) & (x1 v x3)
        DPLL s;
        bool sat = s.solve(3, {{1,2},{-1,3},{-2,-3},{1,3}});
        cout << "  formula A: " << (sat ? "SAT" : "UNSAT");
        if (sat) {
            cout << "  model:";
            for (int v = 1; v <= 3; ++v)
                cout << " x" << v << "=" << (s.model[v] == 1);
        }
        cout << "\n";
        // 不可滿足：(x1) & (~x1)
        DPLL s2;
        cout << "  formula B (x & ~x): "
             << (s2.solve(1, {{1},{-1}}) ? "SAT" : "UNSAT") << "\n";
    }

    // ---- 3. Subset Sum ----
    hr("3. Subset Sum (DP + meet-in-the-middle)");
    {
        vector<int> a = {3, 34, 4, 12, 5, 2};
        cout << "  {3,34,4,12,5,2} sum=9  : DP=" << subsetSumDP(a, 9)
             << " MITM=" << subsetSumMITM(a, 9) << "\n";
        cout << "  {3,34,4,12,5,2} sum=30 : DP=" << subsetSumDP(a, 30)
             << " MITM=" << subsetSumMITM(a, 30) << "\n";
    }

    // ---- 4. Knapsack ----
    hr("4. 0/1 Knapsack (pseudo-polynomial DP)");
    {
        vector<int> w = {10, 20, 30}, v = {60, 100, 120};
        cout << "  W=50, weights{10,20,30}, values{60,100,120} -> best value = "
             << knapsack01(w, v, 50) << " (expect 220)\n";
    }

    // ---- 5. Graph Coloring ----
    hr("5. Graph Coloring (backtracking)");
    {
        // K4（完全圖）需要 4 色；C5（奇環）需要 3 色
        vector<vector<int>> k4 = {{1,2,3},{0,2,3},{0,1,3},{0,1,2}};
        vector<vector<int>> c5 = {{1,4},{0,2},{1,3},{2,4},{3,0}};
        GraphColoring g1(k4), g2(c5);
        cout << "  K4 with 3 colors: " << (g1.colorable(3) ? "yes" : "no")
             << " | with 4: " << (g1.colorable(4) ? "yes" : "no") << "\n";
        cout << "  C5 with 2 colors: " << (g2.colorable(2) ? "yes" : "no")
             << " | with 3: " << (g2.colorable(3) ? "yes" : "no") << "\n";
    }

    // ---- 6. Hamiltonian Path ----
    hr("6. Hamiltonian Path (bitmask DP)");
    {
        // path 圖 0-1-2-3 有；星狀圖（中心 0）沒有
        vector<vector<int>> path = {{1},{0,2},{1,3},{2}};
        vector<vector<int>> star = {{1,2,3},{0},{0},{0}};
        cout << "  path graph: " << (hamiltonianPath(path) ? "yes" : "no")
             << " | star graph: " << (hamiltonianPath(star) ? "yes" : "no") << "\n";
    }

    // ---- 7. TSP Held-Karp ----
    hr("7. TSP exact (Held-Karp bitmask DP)");
    {
        vector<vector<int>> d = {
            {0, 10, 15, 20},
            {10, 0, 35, 25},
            {15, 35, 0, 30},
            {20, 25, 30, 0}};
        cout << "  4-city classic -> optimal tour = " << tspHeldKarp(d)
             << " (expect 80)\n";
    }

    // ---- 8. Vertex Cover exact ----
    hr("8. Vertex Cover k-branching (FPT style)");
    {
        // 三角形 + 一條尾巴: 0-1,1-2,2-0,2-3
        vector<vector<int>> g = {{1,2},{0,2},{0,1,3},{2}};
        cout << "  triangle+tail: k=1 " << (vertexCoverK(g,1) ? "yes" : "no")
             << " | k=2 " << (vertexCoverK(g,2) ? "yes" : "no") << "\n";
    }

    // ---- 9. Vertex Cover 2-approx ----
    hr("9. Vertex Cover 2-approximation (maximal matching)");
    {
        vector<pair<int,int>> edges = {{0,1},{0,2},{1,2},{2,3},{3,4}};
        auto cover = vertexCover2Approx(5, edges);
        cout << "  cover = {";
        for (size_t i = 0; i < cover.size(); ++i)
            cout << cover[i] << (i + 1 < cover.size() ? "," : "");
        cout << "} size=" << cover.size() << " (OPT=2: {2,?} -> bound 2*OPT=4)\n";
    }

    // ---- 10. Set Cover greedy ----
    hr("10. Set Cover greedy approximation");
    {
        // 宇宙 {0..4}；S0={0,1,2}, S1={1,3}, S2={2,3}, S3={3,4}
        vector<vector<int>> sets = {{0,1,2},{1,3},{2,3},{3,4}};
        auto chosen = setCoverGreedy(5, sets);
        cout << "  chosen sets:";
        for (int i : chosen) cout << " S" << i;
        cout << "  (count=" << chosen.size() << ")\n";
    }

    // ---- 11. Bin Packing FFD ----
    hr("11. Bin Packing First-Fit Decreasing");
    {
        vector<double> items = {0.42, 0.25, 0.27, 0.07, 0.72, 0.86, 0.09, 0.44, 0.50, 0.68};
        cout << "  10 items, capacity 1.0 -> bins used = "
             << binPackingFFD(items, 1.0)
             << " (lower bound = ceil(sum) = "
             << (int)ceil(accumulate(items.begin(), items.end(), 0.0)) << ")\n";
    }

    // ---- 12. TSP heuristic: NN + 2-opt vs exact ----
    hr("12. Metric TSP: nearest-neighbor + 2-opt vs Held-Karp");
    {
        mt19937 rng(2026);
        int n = 12;
        vector<double> xs(n), ys(n);
        for (int i = 0; i < n; ++i) {
            xs[i] = (double)(rng() % 1000);
            ys[i] = (double)(rng() % 1000);
        }
        vector<vector<double>> d(n, vector<double>(n, 0));
        vector<vector<int>> di(n, vector<int>(n, 0));
        for (int i = 0; i < n; ++i)
            for (int j = 0; j < n; ++j) {
                d[i][j] = hypot(xs[i] - xs[j], ys[i] - ys[j]);
                di[i][j] = (int)llround(d[i][j] * 1000);   // 整數版給 Held-Karp
            }
        auto tour = tspNearestNeighbor(d);
        double nnLen = tourLength(tour, d);
        tsp2opt(tour, d);
        double optLen2 = tourLength(tour, d);
        double exact = (double)tspHeldKarp(di) / 1000.0;
        cout << "  nearest neighbor : " << nnLen  << "\n";
        cout << "  after 2-opt      : " << optLen2 << "\n";
        cout << "  exact (Held-Karp): " << exact << "\n";
        cout << "  2-opt gap        : "
             << 100.0 * (optLen2 - exact) / exact << "%\n";
    }

    // ---- 13. 2-SAT ----
    hr("13. 2-SAT via SCC (a P-time special case)");
    {
        // (x0 v x1) & (~x0 v x1) & (~x1 v x2) —— 可滿足
        TwoSAT ts(3);
        ts.addClause(0, true, 1, true);
        ts.addClause(0, false, 1, true);
        ts.addClause(1, false, 2, true);
        bool ok = ts.solve();
        cout << "  formula C: " << (ok ? "SAT" : "UNSAT");
        if (ok) {
            cout << "  model:";
            for (int v = 0; v < 3; ++v) cout << " x" << v << "=" << ts.value[v];
        }
        cout << "\n";
        // (x0 v x0) & (~x0 v ~x0) —— 矛盾
        TwoSAT ts2(1);
        ts2.addClause(0, true, 0, true);
        ts2.addClause(0, false, 0, false);
        cout << "  formula D (x & ~x): " << (ts2.solve() ? "SAT" : "UNSAT") << "\n";
    }

    cout << "\nAll demos finished.\n";
    return 0;
}