complete implemenation of bsb: need more debugging

Makefile

@@ -1,6 +1,7 @@
 SHELL = /bin/bash
 
-CXX = g++
+# CXX = g++
+CXX = clang++
 
 SRC_DIR = src
 GEN_DIR = generator
@@ -223,4 +224,4 @@ clean-tests: clean-inputs clean-outputs
 .PHONY: clean
 clean:
 	@echo "Cleaning up $(BUILD_DIR)/..."
-	@rm -rf $(BUILD_DIR)
+	@rm -rf $(BUILD_DIR)

src/bsb/bsb.hpp

@@ -0,0 +1,248 @@
+#ifndef BSB_HPP
+#define BSB_HPP
+
+#include <algorithm>
+#include <boost/heap/fibonacci_heap.hpp>
+#include <cassert>
+#include <cmath>
+#include <list>
+#include <map>
+#include <optional>
+#include <tuple>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+namespace GMS {
+
+using namespace std;
+using ll = unsigned long long;
+
+template <typename KEY, typename VALUE>
+class L33 {
+   private:
+    using KV = pair<KEY, VALUE>;
+    using KV_LIST = vector<KV>;
+    using K_LIST = vector<KEY>;
+
+    using BLOCK = list<KV>;
+    using BLOCK_LIST = list<BLOCK>;
+
+    using ITEM = pair<VALUE, typename BLOCK_LIST::iterator>;
+    using ITEM_LIST = list<ITEM>;
+
+    const ll M;
+    const VALUE B;
+
+    BLOCK_LIST blocks;
+
+    ITEM_LIST D0, D1;
+
+    unordered_map<
+        KEY, tuple<bool, typename ITEM_LIST::iterator,
+                   typename BLOCK_LIST::iterator, typename BLOCK::iterator> >
+        ptr;
+
+    map<VALUE, typename ITEM_LIST::iterator> UB;
+
+    void Delete(const KEY& a, optional<VALUE> b = std::nullopt) {
+        assert(ptr.count(a));
+
+        auto [d01, it_ILIST, it_BLIST, it_B] = ptr[a];
+        ptr.erase(a);
+
+        assert(!b || it_B->second == b.value());
+        it_BLIST->erase(it_B);
+
+        if (it_BLIST->size() == 0) {
+            blocks.erase(it_BLIST);
+
+            if (d01 == 0)
+                D0.erase(it_ILIST);
+            else {
+                UB.erase(it_ILIST->first);
+                D1.erase(it_ILIST);
+
+                if (D1.empty()) {
+                    auto it_BLIST = blocks.emplace(blocks.end());
+                    auto it_ILIST = D1.emplace(D1.end(), B, it_BLIST);
+                    UB[B] = it_ILIST;
+                }
+            }
+        }
+    }
+
+    void split(bool d01, typename ITEM_LIST::iterator it_ILIST) {
+        auto it_BLIST = it_ILIST->second;
+
+        BLOCK& B = *it_BLIST;
+        if (B.size() <= M) return;
+
+        const ll L = B.size() / 2;
+        KV_LIST S;
+        S.reserve(B.size());
+
+        for (auto& item : B) {
+            ptr.erase(item.first);
+            S.emplace_back(std::move(item));
+        }
+
+        auto nth = S.begin();
+        advance(nth, L - 1);
+
+        nth_element(S.begin(), nth, S.end());
+        advance(nth, 1);
+
+        BLOCK X, Y;
+        auto it = S.begin();
+        while (it != nth) {
+            X.emplace_back(std::move(*it)), ++it;
+        }
+        while (it != S.end()) Y.emplace_back(std::move(*it)), ++it;
+
+        VALUE ubx = Y.front().second;
+        for (auto& [a, b] : Y) ubx = min(ubx, b);
+
+        auto it_BLIST_X = blocks.emplace(it_BLIST, std::move(X));
+        auto it_BLIST_Y = blocks.emplace(it_BLIST, std::move(Y));
+        blocks.erase(it_BLIST);
+
+        auto it_ILIST_Y = it_ILIST;
+        it_ILIST->second = it_BLIST_Y;
+        auto it_ILIST_X = (d01 ? D0 : D1).insert(it_ILIST_Y, {ubx, it_BLIST_X});
+
+        if (d01) UB[ubx] = it_ILIST_X;
+
+        for (auto it = it_BLIST_X->begin(); it != it_BLIST_X->end(); it++) {
+            const auto [a, b] = *it;
+            ptr[a] = {d01, it_ILIST_X, it_BLIST_X, it};
+        }
+
+        for (auto it = it_BLIST_Y->begin(); it != it_BLIST_Y->end(); it++) {
+            const auto [a, b] = *it;
+            ptr[a] = {d01, it_ILIST_Y, it_BLIST_Y, it};
+        }
+
+        split(d01, it_ILIST_X);
+        split(d01, it_ILIST_Y);
+    }
+
+   public:
+    L33(ll _M, VALUE _B) : M(_M), B(_B) {
+        auto it_BLIST = blocks.emplace(blocks.end());
+        auto it_ILIST = D1.emplace(D1.end(), B, it_BLIST);
+        UB[B] = it_ILIST;
+    }
+
+    void Insert(const KEY& a, const VALUE& b) {
+        if (ptr.count(a)) {
+            auto [d01, it_ILIST, it_BLIST, it_B] = ptr[a];
+            if (it_B->second <= b) return;
+
+            Delete(a);
+        }
+
+        auto it_ILIST = UB.lower_bound(b)->second;
+        auto it_BLIST = it_ILIST->second;
+        auto it_B = it_BLIST->emplace(it_BLIST->end(), a, b);
+
+        ptr[a] = {1, it_ILIST, it_BLIST, it_B};
+
+        split(1, it_ILIST);
+    }
+
+    void BatchPrepend(const KV_LIST& L) {
+        if (L.size() == 0) return;
+
+        BLOCK tmp;
+        optional<VALUE> u = nullopt;
+        for (auto [a, b] : L) {
+            if (ptr.count(a)) {
+                auto [d01, it_ILIST, it_BLIST, it_B] = ptr[a];
+                if (it_B->second <= b) continue;
+
+                Delete(a);
+            }
+
+            if (!u)
+                u = b;
+            else
+                u = max(u.value(), b);
+
+            tmp.push_back({a, b});
+        }
+
+        if (u == nullopt) return;
+
+        auto it_BLIST = blocks.emplace(blocks.end(), std::move(tmp));
+        auto it_ILIST = D0.emplace(D0.begin(), u.value(), it_BLIST);
+
+        split(0, it_ILIST);
+    }
+
+    pair<VALUE, K_LIST> Pull() {
+        BLOCK S0, S1;
+
+        auto it = D0.begin();
+        while (S0.size() < M + 1 && it != D0.end()) {
+            S0.splice(S0.end(), *it->second);
+            blocks.erase(it->second);
+
+            it = D0.erase(it);
+        }
+
+        it = D1.begin();
+        while (S1.size() < M + 1 && it != D1.end()) {
+            S1.splice(S0.end(), *it->second);
+            blocks.erase(it->second);
+
+            UB.erase(it->first);
+            it = D1.erase(it);
+        }
+
+        KV_LIST S;
+        S.reserve(S0.size() + S1.size());
+
+        for (auto& item : S0) S.emplace_back(std::move(item));
+        for (auto& item : S1) S.emplace_back(std::move(item));
+
+        for (auto [a, b] : S) ptr.erase(a);
+
+        if (D1.empty()) {
+            auto it_BLIST = blocks.emplace(blocks.end());
+            auto it_ILIST = D1.emplace(D1.end(), B, it_BLIST);
+            UB[B] = it_ILIST;
+        }
+
+        if (S.size() <= M) {
+            K_LIST s;
+            s.reserve(S.size());
+            for (auto [a, b] : S) s.push_back(a);
+
+            return {B, s};
+        }
+
+        auto nth = S.begin();
+        std::advance(nth, M - 1);
+        std::nth_element(S.begin(), nth, S.end());
+
+        KV_LIST R(make_move_iterator(nth + 1), make_move_iterator(S.end()));
+        S.resize(M);
+
+        VALUE x = R.front().second;
+        for (auto [a, b] : R) x = min(x, b);
+
+        BatchPrepend(R);
+
+        K_LIST s;
+        s.reserve(S.size());
+        for (auto [a, b] : S) s.push_back(a);
+
+        return {x, s};
+    }
+    bool empty() const { return ptr.empty(); }
+};
+
+}  // namespace GMS
+
+#endif

src/bsb/main.cpp

@@ -1,95 +1,171 @@
 #include <cassert>
+#include <cmath>
 #include <iomanip>
 #include <iostream>
+#include <iterator>
 #include <limits>
+#include <optional>
+#include <unordered_map>
+#include <unordered_set>
 #include <utility>
 #include <vector>
 
-#include <cmath>
+#include "bsb.hpp"
-#include <optional>
-#include <unordered_set>
-#include <unordered_map>
-
-#include <boost/heap/fibonacci_heap.hpp>
 
 using namespace std;
+using namespace GMS;
+
+using DIST = tuple<double, ll, ll>;
+using DIST_LIST = vector<DIST>;
+
+using Item = pair<DIST, ll>;
+using Fiboheap =
+    boost::heap::fibonacci_heap<Item, boost::heap::compare<greater<Item> > >;
 
-using ll = long long;
 using Edge = tuple<ll, ll, double>;
 using ADJ_LIST = vector<vector<pair<ll, double> > >;
 
-using Item = pair<double, ll>;
-using Fiboheap = boost::heap::fibonacci_heap<Item, boost::heap::compare<greater<Item> > >;
-
 const double INF = numeric_limits<double>::infinity();
-pair<vector<ll>, vector<ll> > FindPivots(double B, vector<ll> S, const ADJ_LIST& adj, vector<double>& dist, ll k, ll t)
-{
-    unordered_set<ll> U;
-    for(auto i: S) U.insert(i);
 
-    unordered_set<ll> bef = U;
+bool relex(ll s, DIST& ds, double len, ll i, DIST& di) {
-    for(ll tr = 1; tr <= k; tr++)
+    DIST tmp = {get<0>(ds) + len, get<1>(ds) + 1, s};
-    {
+
+    if (di >= tmp) {
+        di = tmp;
+        return true;
+    }
+
+    return false;
+}
+pair<vector<ll>, vector<ll> > FindPivots(DIST B, vector<ll> S,
+                                         const ADJ_LIST& adj, DIST_LIST& dist,
+                                         const ll k, const ll t) {
+    unordered_set<ll> W;
+    for (auto i : S) W.insert(i);
+
+    unordered_set<ll> bef = W;
+    for (ll tr = 1; tr <= k; tr++) {
         unordered_set<ll> st;
-        for(auto s:bef)
+        for (auto s : bef) {
-        {
+            for (auto [i, len] : adj[s]) {
-            for(auto [i, len] : adj[s])
+                if (relex(s, dist[s], len, i, dist[i])) {
-            {
+                    if (dist[i] < B) st.insert(i);
-                if(dist[i] >= dist[s] + len)
-                {
-
                 }
             }
         }
 
-        for(auto i:st) U.insert(i);
+        for (auto i : st) W.insert(i);
         swap(st, bef);
+
+        if (W.size() > k * S.size()) return {S, vector<ll>(W.begin(), W.end())};
     }
+
+    unordered_map<ll, ll> par, cnt;
+
+    function<int(int)> root = [&W, &par, &cnt, &dist, &root](int s) -> int {
+        int p = get<2>(dist[s]);
+        if (!W.count(p)) {
+            cnt[s]++;
+            return s;
+        }
+
+        if (par.count(s)) {
+            cnt[par[s]]++;
+            return par[s];
+        }
+
+        par[s] = root(p);
+        cnt[par[s]]++;
+        return par[s];
+    };
+
+    for (auto i : W) root(i);
+
+    vector<ll> P;
+    for (auto [i, c] : cnt)
+        if (c >= k) P.push_back(i);
+
+    return {P, vector<ll>(W.begin(), W.end())};
 }
 
-pair<double, vector<ll> > BMSSP(ll l, vector<ll> S, double B, const ADJ_LIST& adj, vector<double>& dist, ll k, ll t)
+pair<DIST, unordered_set<ll> > BMSSP(ll l, vector<ll> S, DIST B, const ADJ_LIST& adj,
-{
+                              DIST_LIST& dist, const ll k, const ll t) {
-    assert(S.size() <= 1LL<<(l*t));
+    assert(S.size() <= (1uLL << (l * t)));
 
-    if(l == 0)
+    if (l == 0) {
-    {
         ll s = S.front();
 
         vector<ll> U;
-        Fiboheap pq; unordered_map<ll, optional<Fiboheap::handle_type> > handle; 
+        Fiboheap pq;
-        
+        unordered_map<ll, optional<Fiboheap::handle_type> > handle;
+
         handle[s] = pq.push({dist[s], s});
-        
+
-        while(!pq.empty() && U.size() < k+1)
+        while (!pq.empty() && U.size() < k + 1) {
-        {
+            auto [du, u] = pq.top();
-            auto [du, u] = pq.top(); pq.pop();
+            pq.pop();
             handle[u] = nullopt;
-            
+
             U.push_back(u);
 
-            for(auto [v, len] : adj[u])
+            for (auto [v, len] : adj[u]) {
-            {
+                if (relex(u, dist[u], len, v, dist[v]) && dist[v] < B &&
-                if(dist[u]+len >= dist[v] && dist[u]+len < B && (!handle.count(v) || handle[v] != nullopt))
+                    (!handle.count(v) || handle[v] != nullopt)) {
-                {
+                    if (handle.count(v))
-                    dist[v] = dist[u] + len;
+                        pq.increase(handle[v].value(), {dist[v], v});
-                    if(handle.count(v)) pq.increase(handle[v].value(), {dist[v], v});
+                    else
-                    else handle[v] = pq.push({dist[v], v});
+                        handle[v] = pq.push({dist[v], v});
                 }
             }
         }
-        
-        if(U.size() <= k) return {B, U};
-        double B_ = dist[U.front()];
-        for(auto i:U) B_ = max(B_, dist[i]);
 
-        return {B_, U};
+        DIST B_ = dist[U.back()];
+
+        unordered_set<ll> u(make_move_iterator(U.begin()), make_move_iterator(U.end()));
+
+        if (U.size() <= k) return {B, u};
+        return {B_, u};
     }
 
-    double B_ = B;
-    vector<ll> U = S;
-
     auto [P, W] = FindPivots(B, S, adj, dist, k, t);
 
+    const ll M = (1uLL << ((l - 1) * t));
+    L33<ll, DIST> D(M, B);
+    for (auto i : P) D.Insert(i, dist[i]);
+
+    const ll LIMIT = k * (1uLL << (l * t));
+
+    unordered_set<ll> U; DIST B_;
+    while (U.size() < LIMIT && !D.empty()) {
+        auto [Bi, Si] = D.Pull();
+        auto [Bi_, Ui] = BMSSP(l - 1, Si, Bi, adj, dist, k, t);
+        B_ = Bi_;
+
+        
+        unordered_set<ll> K;
+
+        for(auto u : Ui) for(auto [v, len] : adj[u])
+        {
+            if(relex(u, dist[u], len, v, dist[v]))
+            {
+                if(Bi <= dist[v] && dist[v] < B) D.Insert(v, dist[v]);
+                else if(Bi_ <= dist[v] && dist[v] < Bi) K.insert(v);
+            }
+        }
+        for(auto s:Si) if(Bi_ <= dist[s] && dist[s] < Bi)
+            K.insert(s);
+
+        U.merge(Ui);
+        
+        vector<pair<ll, DIST> > K_; K_.reserve(K.size());
+        for(auto k:K) K_.push_back({k, dist[k]});
+
+        D.BatchPrepend(K_);
+    }
+
+    B_ = min(B_, B);
+    for(auto s:W) if(dist[s] < B_) U.insert(s);
+
     return {B_, U};
 }
 
@@ -112,41 +188,46 @@ int main() {
     vector<vector<ll> > renum(n + 1);
     ll c = 0;
 
-    auto get = [&c, &renum](ll x) -> ll {
+    auto get_renum = [&c, &renum](ll x) -> ll {
         renum[x].push_back(++c);
         return c;
     };
 
-    ADJ_LIST adj(2*m + 1);
+    ADJ_LIST adj(2 * m + 1);
 
     for (auto [s, e, len] : edges) {
-        adj[get(s)].emplace_back(get(e), len);
+        adj[get_renum(s)].emplace_back(get_renum(e), len);
     }
 
-    assert(c == 2*m);
+    assert(c == 2 * m);
 
     for (ll i = 1; i <= n; i++) {
         if (renum[i].size() <= 1) continue;
 
         ll k = renum[i].size();
         for (ll s = 0; s < k; s++) {
-            adj[renum[i][s]].emplace_back(renum[i][(s+1)%k], 0.0);
+            adj[renum[i][s]].emplace_back(renum[i][(s + 1) % k], 0.0);
         }
     }
 
-    const ll k = floor(pow(log2(2*m), 1.0/3));
+    const ll k = floor(pow(log2(2 * m), 1.0 / 3));
-    const ll t = floor(pow(log2(2*m), 2.0/3));
+    const ll t = floor(pow(log2(2 * m), 2.0 / 3));
 
-    const ll l = ceil((log2(2*m)/t));
+    const ll l = ceil((log2(2 * m) / t));
 
     const int s = 1;
-    vector<double> dist(2*m+1, INF);
-    dist[s] = 1;
 
-    auto [B_, U] = BMSSP(l, {s}, INF, adj, dist, k, t);
+    // distance,
+    DIST_LIST dist(2 * m + 1, {INF, INF, -1});
 
-    vector<double> ans(n+1, INF);
+    dist[s] = {0, 0, -1};
-    for(ll i=1;i<=n;i++) if(!renum[i].empty()) ans[i] = dist[renum[i].front()];
+
+    const DIST Binf = {INF, -1, -1};
+    auto [B_, U] = BMSSP(l, {s}, Binf, adj, dist, k, t);
+
+    vector<double> ans(n + 1, INF);
+    for (ll i = 1; i <= n; i++)
+        if (!renum[i].empty()) ans[i] = get<0>(dist[renum[i].front()]);
 
     cout << fixed << setprecision(15);
     for (ll i = 1; i <= n; i++) cout << ans[i] << '\n';

src/dijk/main.cpp

@@ -9,7 +9,8 @@ using namespace std;
 using ll = long long;
 
 using Item = pair<double, ll>;
-using Fiboheap = boost::heap::fibonacci_heap<Item, boost::heap::compare<greater<Item> > >;
+using Fiboheap =
+    boost::heap::fibonacci_heap<Item, boost::heap::compare<greater<Item> > >;
 
 const double INF = numeric_limits<double>::infinity();