Streaming Algorithm Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// A simple, practical Count-Min Sketch implementation for insertion-only streams.
5
// Guarantees: with width w = ceil(e/eps) and depth d = ceil(log(1/delta)),
6
// estimate(key) <= true_count(key) + eps * N with probability >= 1 - delta,
7
// where N is the total number of insertions.
8
9
struct CountMinSketch {
10
    int w, d;                       // width and depth
11
    vector<vector<long long>> tbl;  // counters (signed here for simplicity)
12
    vector<uint64_t> seeds;         // independent seeds per row
13
14
    // SplitMix64: fast 64-bit hash mixer
15
    static inline uint64_t splitmix64(uint64_t x) {
16
        x += 0x9e3779b97f4a7c15ULL;
17
        x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL;
18
        x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL;
19
        return x ^ (x >> 31);
20
    }
21
22
    // Construct from (eps, delta)
23
    CountMinSketch(double eps = 1e-3, double delta = 1e-6) {
24
        if (eps <= 0) eps = 1e-9;
25
        if (delta <= 0) delta = 1e-12;
26
        w = (int)ceil(exp(1.0) / eps);
27
        d = (int)ceil(log(1.0 / delta));
28
        tbl.assign(d, vector<long long>(w, 0));
29
        seeds.resize(d);
30
        std::mt19937_64 rng(123456789);
31
        for (int i = 0; i < d; ++i) seeds[i] = rng();
32
    }
33
34
    // Hash (key, row) -> [0, w)
35
    inline uint32_t index_hash(uint64_t key, int row) const {
36
        uint64_t h = splitmix64(key ^ seeds[row]);
37
        return (uint32_t)(h % (uint64_t)w);
38
    }
39
40
    // Update with increment delta >= 0 (insertion-only). For generality we allow signed.
41
    void update(uint64_t key, long long delta = 1) {
42
        for (int r = 0; r < d; ++r) {
43
            uint32_t c = index_hash(key, r);
44
            tbl[r][c] += delta; // CMS works best with nonnegative totals
45
        }
46
    }
47
48
    // Point query: min across rows
49
    long long estimate(uint64_t key) const {
50
        long long ans = LLONG_MAX;
51
        for (int r = 0; r < d; ++r) {
52
            uint32_t c = index_hash(key, r);
53
            ans = min(ans, tbl[r][c]);
54
        }
55
        return ans;
56
    }
57
};
58
59
int main() {
60
    ios::sync_with_stdio(false);
61
    cin.tie(nullptr);
62
63
    // Example: track frequencies of numbers 0..999999 with tight memory.
64
    CountMinSketch cms(1e-3, 1e-6); // eps=0.001, delta=1e-6
65
66
    // Simulate a skewed stream: key 42 is heavy.
67
    const int N = 500000;
68
    std::mt19937 rng(42);
69
    std::discrete_distribution<int> dist({1, 100}); // 2 keys: 7 is light, 42 is heavy
70
71
    long long true_7 = 0, true_42 = 0;
72
    for (int i = 0; i < N; ++i) {
73
        int pick = dist(rng);
74
        uint64_t key = (pick == 0 ? 7ULL : 42ULL);
75
        if (key == 7ULL) ++true_7; else ++true_42;
76
        cms.update(key, 1);
77
    }
78
79
    long long est7 = cms.estimate(7);
80
    long long est42 = cms.estimate(42);
81
82
    cout << "True count(7)  = " << true_7 << "\n";
83
    cout << "CMS  estimate(7)= " << est7 << "\n";
84
    cout << "True count(42) = " << true_42 << "\n";
85
    cout << "CMS  estimate(42)= " << est42 << "\n";
86
87
    // The overestimation should be bounded by eps * N (about 500 here) with high probability.
88
    return 0;
89
}
90

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Misra–Gries keeps at most k-1 counters and finds all items with frequency > n/k.
5
// The stored counts are LOWER bounds; each candidate's true count can be up to n/k higher.
6
7
template <class Key>
8
struct MisraGries {
9
    int k; // we store at most k-1 keys
10
    unordered_map<Key, long long> cnt; // candidate counters
11
    long long n_seen = 0; // total processed items
12
13
    explicit MisraGries(int k_) : k(max(2, k_)) {}
14
15
    void update(const Key& x, long long times = 1) {
16
        // Standard MG is for insertion-only (times >= 0). Apply 'times' updates one by one if small.
17
        // For simplicity, assume times >= 0 in this demo.
18
        for (long long t = 0; t < times; ++t) {
19
            ++n_seen;
20
            auto it = cnt.find(x);
21
            if (it != cnt.end()) {
22
                ++(it->second);
23
            } else if ((int)cnt.size() < k - 1) {
24
                cnt[x] = 1;
25
            } else {
26
                // Decrement all
27
                vector<Key> to_erase;
28
                for (auto &p : cnt) {
29
                    if (--p.second == 0) to_erase.push_back(p.first);
30
                }
31
                for (auto &z : to_erase) cnt.erase(z);
32
            }
33
        }
34
    }
35
36
    // Return current candidates and their lower-bound counts
37
    unordered_map<Key, long long> candidates() const { return cnt; }
38
39
    // Error bound: for any key i, true f_i - stored[i] <= n_seen / k
40
    pair<long long, long long> bounds_for(const Key& x) const {
41
        long long lb = 0;
42
        auto it = cnt.find(x);
43
        if (it != cnt.end()) lb = it->second;
44
        long long ub = lb + n_seen / k; // MG bound
45
        return {lb, ub};
46
    }
47
};
48
49
int main() {
50
    ios::sync_with_stdio(false);
51
    cin.tie(nullptr);
52
53
    MisraGries<int> mg(101); // find all items with freq > n/101 (~top-1%)
54
55
    // Simulate a stream with a few heavy hitters.
56
    std::mt19937 rng(123);
57
    vector<int> keys = {1,2,3,4,5,6,7,8,9,10};
58
    std::discrete_distribution<int> dist({1,1,1,1,1,1,1,1,40,50}); // 9 and 10 are heavy
59
60
    const int N = 200000;
61
    for (int i = 0; i < N; ++i) {
62
        int x = keys[dist(rng)];
63
        mg.update(x);
64
    }
65
66
    auto cand = mg.candidates();
67
    cout << "Number of candidates <= k-1 = " << cand.size() << "\n";
68
    // Report bounds for some keys
69
    for (int q : {9, 10, 1, 2}) {
70
        auto [lb, ub] = mg.bounds_for(q);
71
        cout << "Key " << q << ": lower_bound=" << lb << ", upper_bound<=" << ub << "\n";
72
    }
73
74
    // Note: To get exact frequencies of these candidates, perform a second pass or
75
    // maintain an auxiliary exact map only for candidates.
76
    return 0;
77
}
78

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// AMS F2 sketch: maintains R independent signed sums Z_j = sum_i g_j(i) * f_i.
5
// Query returns a robust aggregate (median) of Z_j^2 to estimate F2 = sum_i f_i^2.
6
// This demo uses a pairwise-universal hash mod 2^61-1 to generate g_j(i) in {-1, +1}.
7
8
struct PairHash61 {
9
    // Hash function: h(x) = (a*x + b) mod (2^61 - 1). We use its lowest bit for sign.
10
    static const unsigned __int128 MOD = ((unsigned __int128)1 << 61) - 1; // 2^61 - 1
11
    unsigned long long a, b; // parameters in [1, MOD-1]
12
13
    explicit PairHash61(uint64_t seed) {
14
        std::mt19937_64 rng(seed);
15
        std::uniform_int_distribution<unsigned long long> dist(1ULL, (1ULL<<61) - 2ULL);
16
        a = dist(rng);
17
        b = dist(rng);
18
    }
19
20
    static inline unsigned __int128 mul_mod(unsigned __int128 x, unsigned __int128 y) {
21
        unsigned __int128 z = x * y;
22
        // reduce z mod (2^61 - 1)
23
        unsigned __int128 lo = (z & ((((unsigned __int128)1)<<61) - 1));
24
        unsigned __int128 hi = (z >> 61);
25
        unsigned __int128 r = lo + hi;
26
        if (r >= MOD) r -= MOD;
27
        return r;
28
    }
29
30
    static inline unsigned __int128 add_mod(unsigned __int128 x, unsigned __int128 y) {
31
        unsigned __int128 r = x + y;
32
        if (r >= MOD) r -= MOD;
33
        return r;
34
    }
35
36
    inline uint64_t hash64(uint64_t x) const {
37
        unsigned __int128 X = x;
38
        unsigned __int128 A = a, B = b;
39
        unsigned __int128 t = add_mod(mul_mod(A, X), B);
40
        // return as 64-bit value (we only need LSB for sign)
41
        return (uint64_t)t;
42
    }
43
44
    inline int sign(uint64_t x) const {
45
        // map to +1/-1 using lowest bit
46
        return (hash64(x) & 1ULL) ? +1 : -1;
47
    }
48
};
49
50
struct AMSSketchF2 {
51
    int R; // number of repetitions
52
    vector<PairHash61> h; // independent hashers for signs
53
    vector<long double> Z; // accumulators (can be large in magnitude)
54
55
    explicit AMSSketchF2(int R_, uint64_t base_seed = 2024) : R(max(1, R_)) {
56
        h.reserve(R);
57
        Z.assign(R, 0.0L);
58
        std::mt19937_64 rng(base_seed);
59
        for (int j = 0; j < R; ++j) h.emplace_back(rng());
60
    }
61
62
    // Turnstile update: (key, delta), delta can be negative.
63
    void update(uint64_t key, long long delta) {
64
        for (int j = 0; j < R; ++j) {
65
            int s = h[j].sign(key);
66
            Z[j] += (long double)s * (long double)delta;
67
        }
68
    }
69
70
    // Estimate F2 via robust aggregation of Z_j^2. Median is simple and effective.
71
    long double estimate_F2() const {
72
        vector<long double> vals(R);
73
        for (int j = 0; j < R; ++j) vals[j] = Z[j] * Z[j];
74
        nth_element(vals.begin(), vals.begin() + R/2, vals.end());
75
        return vals[R/2];
76
    }
77
};
78
79
int main() {
80
    ios::sync_with_stdio(false);
81
    cin.tie(nullptr);
82
83
    // Build an AMS sketch with R repetitions. For (eps, delta) one often sets R = c * eps^{-2} * log(1/delta).
84
    int R = 81; // e.g., eps ~ 1/sqrt(R) ~ 0.11 (rough intuition); use larger R for tighter error
85
    AMSSketchF2 ams(R);
86
87
    // We'll also keep an exact map for validation in this demo (not part of the sketch itself).
88
    unordered_map<uint64_t, long long> exact;
89
90
    // Simulate a turnstile stream with insertions and deletions
91
    std::mt19937_64 rng(7);
92
    uniform_int_distribution<uint64_t> keydist(1, 1'000'000);
93
    for (int t = 0; t < 200000; ++t) {
94
        uint64_t k = keydist(rng);
95
        int delta = (rng() & 1) ? +1 : -1; // random +/- 1
96
        ams.update(k, delta);
97
        exact[k] += delta;
98
        if (exact[k] == 0) exact.erase(k);
99
    }
100
101
    // Exact F2 for comparison (requires memory and is only for demonstration)
102
    long double F2_exact = 0.0L;
103
    for (auto &p : exact) F2_exact += (long double)p.second * (long double)p.second;
104
105
    long double F2_est = ams.estimate_F2();
106
107
    cout.setf(std::ios::fixed); cout << setprecision(2);
108
    cout << "Exact F2      = " << (double)F2_exact << "\n";
109
    cout << "AMS estimate  = " << (double)F2_est << "\n";
110
    cout << "Rel. error    = " << (double)(fabsl(F2_est - F2_exact) / max(1.0L, F2_exact)) << "\n";
111
    return 0;
112
}
113