lyndon_benchmark.cpp

#include <algorithm>
#include <cassert>
#include <chrono>
#include <cstdint>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

struct AlgoStats {
    uint64_t extend_calls = 0;
    uint64_t fast_returns = 0;
    uint64_t reuses = 0;
    uint64_t while_iterations = 0;
    uint64_t suffix_links = 0;
    uint64_t char_cmps = 0;
    uint64_t assignments = 0;
    uint64_t cleanup = 0;

    uint64_t window_hits = 0;
    uint64_t reuse_from_next = 0;
    uint64_t reuse_from_prev_lce = 0;
    uint64_t fast_return_success = 0;
    uint64_t fast_return_to_rhs = 0;
    uint64_t rhs_fallback_total = 0;
    uint64_t fast_return_lce_total = 0;
    uint64_t eq_plce_branch = 0;
    uint64_t gt_plce_branch = 0;
    uint64_t explicit_extend_calls = 0;
    uint64_t explicit_match_cmps = 0;
    uint64_t explicit_mismatch_cmps = 0;
};

struct BenchStatsAvg {
    double extend_calls = 0.0;
    double fast_returns = 0.0;
    double reuses = 0.0;
    double while_iterations = 0.0;
    double suffix_links = 0.0;
    double char_cmps = 0.0;
    double assignments = 0.0;
    double cleanup = 0.0;
    double window_hits = 0.0;
    double reuse_from_next = 0.0;
    double reuse_from_prev_lce = 0.0;
    double fast_return_success = 0.0;
    double fast_return_to_rhs = 0.0;
    double rhs_fallback_total = 0.0;
    double fast_return_lce_total = 0.0;
    double eq_plce_branch = 0.0;
    double gt_plce_branch = 0.0;
    double explicit_extend_calls = 0.0;
    double explicit_match_cmps = 0.0;
    double explicit_mismatch_cmps = 0.0;
};

struct BenchResult {
    int n = 0;
    double nss_us = 0.0;
    double ngs_us = 0.0;
    double recover_nss_us = 0.0;
    double recover_ngs_us = 0.0;
};

struct BenchProfileResult : BenchResult {
    BenchStatsAvg nss_stats;
    BenchStatsAvg ngs_stats;
};

using Clock = std::chrono::high_resolution_clock;

static inline int t_nss_left(const unsigned char *text, int x, int n) {
    return (x >= 0) ? ((x < n) ? (int)text[x] : -1) : -1;
}
static inline int t_nss_right(const unsigned char *text, int x, int n) {
    return (x < n) ? ((x >= 0) ? (int)text[x] : -1) : -1;
}
static inline int t_ngs_left(const unsigned char *text, int x, int n) {
    return (x >= 0) ? ((x < n) ? (int)text[x] : 256) : 256;
}
static inline int t_ngs_right(const unsigned char *text, int x, int n) {
    return (x < n) ? ((x >= 0) ? (int)text[x] : 256) : 256;
}

static inline int first_mismatch_u64(const unsigned char *a, const unsigned char *b) {
    for (int k = 0; k < 8; ++k) {
        if (a[k] != b[k]) return k;
    }
    return 8;
}

static inline int extend_lce_nss_plain(const unsigned char *text,
                                       const int *nss,
                                       const int *nlce,
                                       const int *plce,
                                       int n,
                                       int l,
                                       int r,
                                       int known,
                                       int &d,
                                       int &rhs) {
    if (r + known < rhs) {
        known = (nss[l - d] == r - d) ? nlce[l - d] : plce[r - d];
        if (r + known < rhs) return known;
        known = rhs - r;
    }

    while (true) {
        int li = l + known;
        int ri = r + known;
        while (li >= 0 && li + 8 <= n && ri + 8 <= n) {
            uint64_t wa = 0, wb = 0;
            std::memcpy(&wa, text + li, sizeof(uint64_t));
            std::memcpy(&wb, text + ri, sizeof(uint64_t));
            if (wa == wb) {
                known += 8;
                li += 8;
                ri += 8;
            } else {
                known += first_mismatch_u64(text + li, text + ri);
                rhs = r + known;
                d = r - l;
                return known;
            }
        }
        int cl = t_nss_left(text, li, n);
        int cr = t_nss_right(text, ri, n);
        if (cl != cr) break;
        ++known;
    }
    rhs = r + known;
    d = r - l;
    return known;
}

static inline int extend_lce_ngs_plain(const unsigned char *text,
                                       const int *ngs,
                                       const int *nlce,
                                       const int *plce,
                                       int n,
                                       int l,
                                       int r,
                                       int known,
                                       int &d,
                                       int &rhs) {
    if (r + known < rhs) {
        known = (ngs[l - d] == r - d) ? nlce[l - d] : plce[r - d];
        if (r + known < rhs) return known;
        known = rhs - r;
    }

    while (true) {
        int li = l + known;
        int ri = r + known;
        while (li >= 0 && li + 8 <= n && ri + 8 <= n) {
            uint64_t wa = 0, wb = 0;
            std::memcpy(&wa, text + li, sizeof(uint64_t));
            std::memcpy(&wb, text + ri, sizeof(uint64_t));
            if (wa == wb) {
                known += 8;
                li += 8;
                ri += 8;
            } else {
                known += first_mismatch_u64(text + li, text + ri);
                rhs = r + known;
                d = r - l;
                return known;
            }
        }
        int cl = t_ngs_left(text, li, n);
        int cr = t_ngs_right(text, ri, n);
        if (cl != cr) break;
        ++known;
    }
    rhs = r + known;
    d = r - l;
    return known;
}

void lce_nss_plain(const char *text_in, int *nss, int *nlce, int *pss, int *plce, int n) {
    const auto *text = reinterpret_cast<const unsigned char *>(text_in);
    int d = 0;
    int rhs = -1;
    for (int r = 0; r < n; ++r) {
        int l = r - 1;
        int lce = extend_lce_nss_plain(text, nss, nlce, plce, n, l, r, 0, d, rhs);
        while (t_nss_left(text, l + lce, n) > t_nss_right(text, r + lce, n)) {
            nss[l] = r;
            nlce[l] = lce;
            if (lce == plce[l]) lce = extend_lce_nss_plain(text, nss, nlce, plce, n, pss[l], r, lce, d, rhs);
            else if (lce > plce[l]) lce = plce[l];
            l = pss[l];
        }
        pss[r] = l;
        plce[r] = lce;
    }
    int l = n - 1;
    while (l >= 0) {
        nss[l] = n;
        nlce[l] = 0;
        l = pss[l];
    }
}

void lce_ngs_plain(const char *text_in, int *ngs, int *nlce, int *pgs, int *plce, int n) {
    const auto *text = reinterpret_cast<const unsigned char *>(text_in);
    int d = 0;
    int rhs = -1;
    for (int r = 0; r < n; ++r) {
        int l = r - 1;
        int lce = extend_lce_ngs_plain(text, ngs, nlce, plce, n, l, r, 0, d, rhs);
        while (t_ngs_left(text, l + lce, n) < t_ngs_right(text, r + lce, n)) {
            ngs[l] = r;
            nlce[l] = lce;
            if (lce == plce[l]) lce = extend_lce_ngs_plain(text, ngs, nlce, plce, n, pgs[l], r, lce, d, rhs);
            else if (lce > plce[l]) lce = plce[l];
            l = pgs[l];
        }
        pgs[r] = l;
        plce[r] = lce;
    }
    int l = n - 1;
    while (l >= 0) {
        ngs[l] = n;
        nlce[l] = 0;
        l = pgs[l];
    }
}

static inline int extend_lce_nss_profile(const unsigned char *text,
                                         const int *nss,
                                         const int *nlce,
                                         const int *plce,
                                         int n,
                                         int l,
                                         int r,
                                         int known,
                                         int &d,
                                         int &rhs,
                                         AlgoStats &st) {
    ++st.extend_calls;
    if (r + known < rhs) {
        ++st.window_hits;
        ++st.reuses;
        if (nss[l - d] == r - d) {
            ++st.reuse_from_next;
            known = nlce[l - d];
        } else {
            ++st.reuse_from_prev_lce;
            known = plce[r - d];
        }
        if (r + known < rhs) {
            ++st.fast_returns;
            ++st.fast_return_success;
            st.fast_return_lce_total += (uint64_t)known;
            return known;
        }
        ++st.fast_return_to_rhs;
        known = rhs - r;
        st.rhs_fallback_total += (uint64_t)known;
    }

    ++st.explicit_extend_calls;
    while (true) {
        const int li = l + known;
        const int ri = r + known;
        const int cl = t_nss_left(text, li, n);
        const int cr = t_nss_right(text, ri, n);
        ++st.char_cmps;
        if (cl != cr) {
            ++st.explicit_mismatch_cmps;
            break;
        }
        ++st.explicit_match_cmps;
        ++known;
    }
    rhs = r + known;
    d = r - l;
    return known;
}

static inline int extend_lce_ngs_profile(const unsigned char *text,
                                         const int *ngs,
                                         const int *nlce,
                                         const int *plce,
                                         int n,
                                         int l,
                                         int r,
                                         int known,
                                         int &d,
                                         int &rhs,
                                         AlgoStats &st) {
    ++st.extend_calls;
    if (r + known < rhs) {
        ++st.window_hits;
        ++st.reuses;
        if (ngs[l - d] == r - d) {
            ++st.reuse_from_next;
            known = nlce[l - d];
        } else {
            ++st.reuse_from_prev_lce;
            known = plce[r - d];
        }
        if (r + known < rhs) {
            ++st.fast_returns;
            ++st.fast_return_success;
            st.fast_return_lce_total += (uint64_t)known;
            return known;
        }
        ++st.fast_return_to_rhs;
        known = rhs - r;
        st.rhs_fallback_total += (uint64_t)known;
    }

    ++st.explicit_extend_calls;
    while (true) {
        const int li = l + known;
        const int ri = r + known;
        const int cl = t_ngs_left(text, li, n);
        const int cr = t_ngs_right(text, ri, n);
        ++st.char_cmps;
        if (cl != cr) {
            ++st.explicit_mismatch_cmps;
            break;
        }
        ++st.explicit_match_cmps;
        ++known;
    }
    rhs = r + known;
    d = r - l;
    return known;
}

void lce_nss_profile(const char *text_in, int *nss, int *nlce, int *pss, int *plce, int n, AlgoStats &st) {
    const auto *text = reinterpret_cast<const unsigned char *>(text_in);
    int d = 0;
    int rhs = -1;
    for (int r = 0; r < n; ++r) {
        int l = r - 1;
        int lce = extend_lce_nss_profile(text, nss, nlce, plce, n, l, r, 0, d, rhs, st);
        while (t_nss_left(text, l + lce, n) > t_nss_right(text, r + lce, n)) {
            ++st.while_iterations;
            ++st.assignments;
            nss[l] = r;
            nlce[l] = lce;
            if (lce == plce[l]) {
                ++st.eq_plce_branch;
                lce = extend_lce_nss_profile(text, nss, nlce, plce, n, pss[l], r, lce, d, rhs, st);
            } else if (lce > plce[l]) {
                ++st.gt_plce_branch;
                lce = plce[l];
            }
            l = pss[l];
            ++st.suffix_links;
        }
        pss[r] = l;
        plce[r] = lce;
    }
    int l = n - 1;
    while (l >= 0) {
        ++st.cleanup;
        nss[l] = n;
        nlce[l] = 0;
        l = pss[l];
    }
}

void lce_ngs_profile(const char *text_in, int *ngs, int *nlce, int *pgs, int *plce, int n, AlgoStats &st) {
    const auto *text = reinterpret_cast<const unsigned char *>(text_in);
    int d = 0;
    int rhs = -1;
    for (int r = 0; r < n; ++r) {
        int l = r - 1;
        int lce = extend_lce_ngs_profile(text, ngs, nlce, plce, n, l, r, 0, d, rhs, st);
        while (t_ngs_left(text, l + lce, n) < t_ngs_right(text, r + lce, n)) {
            ++st.while_iterations;
            ++st.assignments;
            ngs[l] = r;
            nlce[l] = lce;
            if (lce == plce[l]) {
                ++st.eq_plce_branch;
                lce = extend_lce_ngs_profile(text, ngs, nlce, plce, n, pgs[l], r, lce, d, rhs, st);
            } else if (lce > plce[l]) {
                ++st.gt_plce_branch;
                lce = plce[l];
            }
            l = pgs[l];
            ++st.suffix_links;
        }
        pgs[r] = l;
        plce[r] = lce;
    }
    int l = n - 1;
    while (l >= 0) {
        ++st.cleanup;
        ngs[l] = n;
        nlce[l] = 0;
        l = pgs[l];
    }
}

void recover_lyndon_array(const int *nss, int *lambda, int n) {
    for (int i = 0; i < n; ++i) lambda[i] = nss[i] - i;
}
void recover_inverse_lyndon_array(const int *ngs, const int *nlce, int *lambda_inv, int n) {
    for (int i = 0; i < n; ++i) lambda_inv[i] = ngs[i] - i + nlce[i];
}

bool is_lyndon(const char *text, int i, int m) {
    for (int k = 1; k < m; ++k) {
        const int sufflen = m - k;
        bool decided = false;
        for (int j = 0; j < sufflen; ++j) {
            const unsigned char cw = (unsigned char)text[i + j];
            const unsigned char cs = (unsigned char)text[i + k + j];
            if (cw < cs) { decided = true; break; }
            if (cw > cs) return false;
        }
        if (!decided) return false;
    }
    return true;
}

bool is_inverse_lyndon(const char *text, int i, int m) {
    for (int k = 1; k < m; ++k) {
        const int sufflen = m - k;
        for (int j = 0; j < sufflen; ++j) {
            const unsigned char cs = (unsigned char)text[i + k + j];
            const unsigned char cw = (unsigned char)text[i + j];
            if (cs < cw) break;
            if (cs > cw) return false;
        }
    }
    return true;
}

void bruteforce_lyndon_array(const char *text, int *lambda, int n) {
    for (int i = 0; i < n; ++i) {
        lambda[i] = 1;
        for (int m = n - i; m > 1; --m) {
            if (is_lyndon(text, i, m)) {
                lambda[i] = m;
                break;
            }
        }
    }
}

void bruteforce_inverse_lyndon_array(const char *text, int *lambda_inv, int n) {
    for (int i = 0; i < n; ++i) {
        lambda_inv[i] = 1;
        for (int m = n - i; m > 1; --m) {
            if (is_inverse_lyndon(text, i, m)) {
                lambda_inv[i] = m;
                break;
            }
        }
    }
}

bool full_verify(const char *text, int n, bool verbose) {
    std::vector<int> nss(n), nlce_nss(n), pss(n), plce_nss(n), lambda(n), lambda_bf(n);
    std::vector<int> ngs(n), nlce_ngs(n), pgs(n), plce_ngs(n), lambda_inv(n), lambda_inv_bf(n);
    lce_nss_plain(text, nss.data(), nlce_nss.data(), pss.data(), plce_nss.data(), n);
    lce_ngs_plain(text, ngs.data(), nlce_ngs.data(), pgs.data(), plce_ngs.data(), n);
    recover_lyndon_array(nss.data(), lambda.data(), n);
    recover_inverse_lyndon_array(ngs.data(), nlce_ngs.data(), lambda_inv.data(), n);
    bruteforce_lyndon_array(text, lambda_bf.data(), n);
    bruteforce_inverse_lyndon_array(text, lambda_inv_bf.data(), n);
    for (int i = 0; i < n; ++i) {
        if (lambda[i] != lambda_bf[i] || lambda_inv[i] != lambda_inv_bf[i]) {
            if (verbose) std::cerr << "Verification failed at i=" << i << "\n";
            return false;
        }
    }
    return true;
}

static inline void add_stats(BenchStatsAvg &dst, const AlgoStats &src) {
    dst.extend_calls += src.extend_calls;
    dst.fast_returns += src.fast_returns;
    dst.reuses += src.reuses;
    dst.while_iterations += src.while_iterations;
    dst.suffix_links += src.suffix_links;
    dst.char_cmps += src.char_cmps;
    dst.assignments += src.assignments;
    dst.cleanup += src.cleanup;
    dst.window_hits += src.window_hits;
    dst.reuse_from_next += src.reuse_from_next;
    dst.reuse_from_prev_lce += src.reuse_from_prev_lce;
    dst.fast_return_success += src.fast_return_success;
    dst.fast_return_to_rhs += src.fast_return_to_rhs;
    dst.rhs_fallback_total += src.rhs_fallback_total;
    dst.fast_return_lce_total += src.fast_return_lce_total;
    dst.eq_plce_branch += src.eq_plce_branch;
    dst.gt_plce_branch += src.gt_plce_branch;
    dst.explicit_extend_calls += src.explicit_extend_calls;
    dst.explicit_match_cmps += src.explicit_match_cmps;
    dst.explicit_mismatch_cmps += src.explicit_mismatch_cmps;
}

static inline void scale_stats(BenchStatsAvg &st, double inv_runs) {
    st.extend_calls *= inv_runs;
    st.fast_returns *= inv_runs;
    st.reuses *= inv_runs;
    st.while_iterations *= inv_runs;
    st.suffix_links *= inv_runs;
    st.char_cmps *= inv_runs;
    st.assignments *= inv_runs;
    st.cleanup *= inv_runs;
    st.window_hits *= inv_runs;
    st.reuse_from_next *= inv_runs;
    st.reuse_from_prev_lce *= inv_runs;
    st.fast_return_success *= inv_runs;
    st.fast_return_to_rhs *= inv_runs;
    st.rhs_fallback_total *= inv_runs;
    st.fast_return_lce_total *= inv_runs;
    st.eq_plce_branch *= inv_runs;
    st.gt_plce_branch *= inv_runs;
    st.explicit_extend_calls *= inv_runs;
    st.explicit_match_cmps *= inv_runs;
    st.explicit_mismatch_cmps *= inv_runs;
}

BenchResult benchmark_plain(const char *text, int n, int runs) {
    std::vector<int> nss(n), nlce_nss(n), pss(n), plce_nss(n), lambda(n);
    std::vector<int> ngs(n), nlce_ngs(n), pgs(n), plce_ngs(n), lambda_inv(n);

    double tnss = 0.0, tngs = 0.0, trnss = 0.0, trngs = 0.0;
    for (int r = 0; r < runs; ++r) {
        if ((r & 1) == 0) {
            auto s1 = Clock::now();
            lce_nss_plain(text, nss.data(), nlce_nss.data(), pss.data(), plce_nss.data(), n);
            auto e1 = Clock::now();
            tnss += std::chrono::duration_cast<std::chrono::microseconds>(e1 - s1).count();

            auto s1r = Clock::now();
            recover_lyndon_array(nss.data(), lambda.data(), n);
            auto e1r = Clock::now();
            trnss += std::chrono::duration_cast<std::chrono::microseconds>(e1r - s1r).count();

            auto s2 = Clock::now();
            lce_ngs_plain(text, ngs.data(), nlce_ngs.data(), pgs.data(), plce_ngs.data(), n);
            auto e2 = Clock::now();
            tngs += std::chrono::duration_cast<std::chrono::microseconds>(e2 - s2).count();

            auto s2r = Clock::now();
            recover_inverse_lyndon_array(ngs.data(), nlce_ngs.data(), lambda_inv.data(), n);
            auto e2r = Clock::now();
            trngs += std::chrono::duration_cast<std::chrono::microseconds>(e2r - s2r).count();
        } else {
            auto s2 = Clock::now();
            lce_ngs_plain(text, ngs.data(), nlce_ngs.data(), pgs.data(), plce_ngs.data(), n);
            auto e2 = Clock::now();
            tngs += std::chrono::duration_cast<std::chrono::microseconds>(e2 - s2).count();

            auto s2r = Clock::now();
            recover_inverse_lyndon_array(ngs.data(), nlce_ngs.data(), lambda_inv.data(), n);
            auto e2r = Clock::now();
            trngs += std::chrono::duration_cast<std::chrono::microseconds>(e2r - s2r).count();

            auto s1 = Clock::now();
            lce_nss_plain(text, nss.data(), nlce_nss.data(), pss.data(), plce_nss.data(), n);
            auto e1 = Clock::now();
            tnss += std::chrono::duration_cast<std::chrono::microseconds>(e1 - s1).count();

            auto s1r = Clock::now();
            recover_lyndon_array(nss.data(), lambda.data(), n);
            auto e1r = Clock::now();
            trnss += std::chrono::duration_cast<std::chrono::microseconds>(e1r - s1r).count();
        }
    }
    return {n, tnss / runs, tngs / runs, trnss / runs, trngs / runs};
}

BenchProfileResult benchmark_profile(const char *text, int n, int runs) {
    std::vector<int> nss(n), nlce_nss(n), pss(n), plce_nss(n), lambda(n);
    std::vector<int> ngs(n), nlce_ngs(n), pgs(n), plce_ngs(n), lambda_inv(n);

    BenchProfileResult res;
    res.n = n;
    for (int r = 0; r < runs; ++r) {
        AlgoStats st_nss, st_ngs;
        if ((r & 1) == 0) {
            auto s1 = Clock::now();
            lce_nss_profile(text, nss.data(), nlce_nss.data(), pss.data(), plce_nss.data(), n, st_nss);
            auto e1 = Clock::now();
            res.nss_us += std::chrono::duration_cast<std::chrono::microseconds>(e1 - s1).count();

            auto s1r = Clock::now();
            recover_lyndon_array(nss.data(), lambda.data(), n);
            auto e1r = Clock::now();
            res.recover_nss_us += std::chrono::duration_cast<std::chrono::microseconds>(e1r - s1r).count();

            auto s2 = Clock::now();
            lce_ngs_profile(text, ngs.data(), nlce_ngs.data(), pgs.data(), plce_ngs.data(), n, st_ngs);
            auto e2 = Clock::now();
            res.ngs_us += std::chrono::duration_cast<std::chrono::microseconds>(e2 - s2).count();

            auto s2r = Clock::now();
            recover_inverse_lyndon_array(ngs.data(), nlce_ngs.data(), lambda_inv.data(), n);
            auto e2r = Clock::now();
            res.recover_ngs_us += std::chrono::duration_cast<std::chrono::microseconds>(e2r - s2r).count();
        } else {
            auto s2 = Clock::now();
            lce_ngs_profile(text, ngs.data(), nlce_ngs.data(), pgs.data(), plce_ngs.data(), n, st_ngs);
            auto e2 = Clock::now();
            res.ngs_us += std::chrono::duration_cast<std::chrono::microseconds>(e2 - s2).count();

            auto s2r = Clock::now();
            recover_inverse_lyndon_array(ngs.data(), nlce_ngs.data(), lambda_inv.data(), n);
            auto e2r = Clock::now();
            res.recover_ngs_us += std::chrono::duration_cast<std::chrono::microseconds>(e2r - s2r).count();

            auto s1 = Clock::now();
            lce_nss_profile(text, nss.data(), nlce_nss.data(), pss.data(), plce_nss.data(), n, st_nss);
            auto e1 = Clock::now();
            res.nss_us += std::chrono::duration_cast<std::chrono::microseconds>(e1 - s1).count();

            auto s1r = Clock::now();
            recover_lyndon_array(nss.data(), lambda.data(), n);
            auto e1r = Clock::now();
            res.recover_nss_us += std::chrono::duration_cast<std::chrono::microseconds>(e1r - s1r).count();
        }
        add_stats(res.nss_stats, st_nss);
        add_stats(res.ngs_stats, st_ngs);
    }
    const double inv_runs = 1.0 / runs;
    res.nss_us *= inv_runs;
    res.ngs_us *= inv_runs;
    res.recover_nss_us *= inv_runs;
    res.recover_ngs_us *= inv_runs;
    scale_stats(res.nss_stats, inv_runs);
    scale_stats(res.ngs_stats, inv_runs);
    return res;
}

std::string generate_random(int n, int sigma, uint64_t seed) {
    std::mt19937_64 gen(seed);
    std::uniform_int_distribution<int> dist(0, sigma - 1);
    std::string s(n, 'a');
    for (int i = 0; i < n; ++i) s[i] = char('a' + dist(gen));
    return s;
}

std::string load_file(const std::string &path, int maxlen) {
    std::ifstream f(path, std::ios::binary);
    if (!f) throw std::runtime_error("Cannot open file: " + path);
    f.seekg(0, std::ios::end);
    std::streamoff sz = f.tellg();
    f.seekg(0, std::ios::beg);
    const int n = (int)std::min<std::streamoff>(sz, maxlen);
    std::string s(n, '\0');
    f.read(&s[0], n);
    return s;
}

std::string make_repetitive_constant(int n) {
    return std::string(n, 'a');
}
std::string make_repetitive_periodic(int n) {
    std::string s(n, 'a');
    for (int i = 0; i < n; ++i) s[i] = (i & 1) ? 'b' : 'a';
    return s;
}
std::string make_long_border(int n, int pct) {
    int b = std::max(1, (n * pct) / 100);
    if (2 * b > n) b = n / 2;
    std::string pref = generate_random(b, 4, 1000 + pct + n);
    std::string mid = generate_random(std::max(0, n - 2 * b), 4, 2000 + pct + n);
    return pref + mid + pref;
}
std::string make_quasi_monotone(int n, int noise_pct) {
    std::string s(n, 'a');
    for (int i = 0; i < n; ++i) s[i] = char('a' + std::min(25, i / std::max(1, n / 26)));
    std::mt19937_64 gen(3000 + noise_pct + n);
    std::uniform_int_distribution<int> pos(0, n - 1);
    std::uniform_int_distribution<int> chr(0, 25);
    int flips = (n * noise_pct) / 100;
    for (int k = 0; k < flips; ++k) s[pos(gen)] = char('a' + chr(gen));
    return s;
}

void print_profile_csv_row(const std::string &dataset, const std::string &param, const BenchProfileResult &r) {
    auto rd = [](double x) -> long long { return (long long)std::llround(x); };
    std::cout << dataset << ',' << param << ',' << r.n << ','
              << std::fixed << std::setprecision(1)
              << r.nss_us << ',' << r.ngs_us << ',' << r.recover_nss_us << ',' << r.recover_ngs_us << ','
              << std::setprecision(4) << (r.nss_us > 0.0 ? r.ngs_us / r.nss_us : 0.0)
              << ',' << rd(r.nss_stats.extend_calls) << ',' << rd(r.ngs_stats.extend_calls)
              << ',' << rd(r.nss_stats.fast_returns) << ',' << rd(r.ngs_stats.fast_returns)
              << ',' << rd(r.nss_stats.reuses) << ',' << rd(r.ngs_stats.reuses)
              << ',' << rd(r.nss_stats.while_iterations) << ',' << rd(r.ngs_stats.while_iterations)
              << ',' << rd(r.nss_stats.suffix_links) << ',' << rd(r.ngs_stats.suffix_links)
              << ',' << rd(r.nss_stats.char_cmps) << ',' << rd(r.ngs_stats.char_cmps)
              << ',' << rd(r.nss_stats.assignments) << ',' << rd(r.ngs_stats.assignments)
              << ',' << rd(r.nss_stats.cleanup) << ',' << rd(r.ngs_stats.cleanup)
              << ',' << rd(r.nss_stats.window_hits) << ',' << rd(r.ngs_stats.window_hits)
              << ',' << rd(r.nss_stats.reuse_from_next) << ',' << rd(r.ngs_stats.reuse_from_next)
              << ',' << rd(r.nss_stats.reuse_from_prev_lce) << ',' << rd(r.ngs_stats.reuse_from_prev_lce)
              << ',' << rd(r.nss_stats.fast_return_success) << ',' << rd(r.ngs_stats.fast_return_success)
              << ',' << rd(r.nss_stats.fast_return_to_rhs) << ',' << rd(r.ngs_stats.fast_return_to_rhs)
              << ',' << rd(r.nss_stats.rhs_fallback_total) << ',' << rd(r.ngs_stats.rhs_fallback_total)
              << ',' << rd(r.nss_stats.fast_return_lce_total) << ',' << rd(r.ngs_stats.fast_return_lce_total)
              << ',' << rd(r.nss_stats.eq_plce_branch) << ',' << rd(r.ngs_stats.eq_plce_branch)
              << ',' << rd(r.nss_stats.gt_plce_branch) << ',' << rd(r.ngs_stats.gt_plce_branch)
              << ',' << rd(r.nss_stats.explicit_extend_calls) << ',' << rd(r.ngs_stats.explicit_extend_calls)
              << ',' << rd(r.nss_stats.explicit_match_cmps) << ',' << rd(r.ngs_stats.explicit_match_cmps)
              << ',' << rd(r.nss_stats.explicit_mismatch_cmps) << ',' << rd(r.ngs_stats.explicit_mismatch_cmps)
              << '\n';
}

void print_plain_csv_row(const std::string &dataset, const std::string &param, const BenchResult &r) {
    std::cout << dataset << ',' << param << ',' << r.n << ','
              << std::fixed << std::setprecision(1)
              << r.nss_us << ',' << r.ngs_us << ',' << r.recover_nss_us << ',' << r.recover_ngs_us << ','
              << std::setprecision(4) << (r.nss_us > 0.0 ? r.ngs_us / r.nss_us : 0.0) << '\n';
}

void print_border_profile_csv_row(const std::string &family, const BenchProfileResult &r) {
    auto rd = [](double x) -> long long { return (long long)std::llround(x); };
    std::cout << family << ',' << r.n << ','
              << std::fixed << std::setprecision(1)
              << r.nss_us << ',' << r.ngs_us << ','
              << std::setprecision(4) << (r.nss_us > 0.0 ? r.ngs_us / r.nss_us : 0.0)
              << ',' << rd(r.nss_stats.char_cmps) << ',' << rd(r.ngs_stats.char_cmps)
              << ',' << rd(r.nss_stats.reuses) << ',' << rd(r.ngs_stats.reuses)
              << ',' << rd(r.nss_stats.explicit_extend_calls) << ',' << rd(r.ngs_stats.explicit_extend_calls)
              << '\n';
}

void usage() {
    std::cerr << "Usage:\n"
              << "  lyndon_benchmark verify <n> [iters]\n"
              << "  lyndon_benchmark random <n> [runs]\n"
              << "  lyndon_benchmark file <path> [maxlen] [runs]\n"
              << "  lyndon_benchmark bench [runs]\n"
              << "  lyndon_benchmark bench_profile [runs]\n"
              << "  lyndon_benchmark bench_struct [runs]\n"
              << "  lyndon_benchmark bench_struct_profile [runs]\n"
              << "  lyndon_benchmark bench_border [runs]\n"
              << "  lyndon_benchmark bench_border_profile [runs]\n"
              << "  lyndon_benchmark bench_files <runs> <maxlen> <file1> [file2 ...]\n"
              << "  lyndon_benchmark bench_files_profile <runs> <maxlen> <file1> [file2 ...]\n";
}

int main(int argc, char **argv) {
    if (argc < 2) {
        usage();
        return 1;
    }
    const std::string mode = argv[1];

    if (mode == "verify") {
        int n = (argc >= 3) ? std::stoi(argv[2]) : 2000;
        int iters = (argc >= 4) ? std::stoi(argv[3]) : 50;
        for (int sigma : {2, 4, 26}) {
            for (int t = 0; t < iters; ++t) {
                auto s = generate_random(n, sigma, 42 + sigma * 1000ULL + t);
                if (!full_verify(s.data(), (int)s.size(), false)) {
                    std::cerr << "Verification failed for sigma=" << sigma << " iter=" << t << "\n";
                    return 1;
                }
            }
        }
        std::cout << "All verifications passed.\n";
        return 0;
    }

    if (mode == "random") {
        if (argc < 3) {
            usage();
            return 1;
        }
        int n = std::stoi(argv[2]);
        int runs = (argc >= 4) ? std::stoi(argv[3]) : 5;
        auto s = generate_random(n, 26, 12345);
        auto r = benchmark_plain(s.data(), (int)s.size(), runs);
        std::cout << "n=" << n << " runs=" << runs << "\n";
        std::cout << "LCE-NSS: " << std::fixed << std::setprecision(1) << r.nss_us << " us\n";
        std::cout << "LCE-NGS: " << r.ngs_us << " us\n";
        std::cout << "Ratio:   " << std::setprecision(4) << (r.nss_us > 0.0 ? r.ngs_us / r.nss_us : 0.0) << "\n";
        return 0;
    }

    if (mode == "file") {
        if (argc < 3) {
            usage();
            return 1;
        }
        int maxlen = (argc >= 4) ? std::stoi(argv[3]) : 100000000;
        int runs = (argc >= 5) ? std::stoi(argv[4]) : 5;
        auto s = load_file(argv[2], maxlen);
        auto r = benchmark_plain(s.data(), (int)s.size(), runs);
        std::cout << "file=" << argv[2] << " n=" << s.size() << " runs=" << runs << "\n";
        std::cout << "LCE-NSS: " << std::fixed << std::setprecision(1) << r.nss_us << " us\n";
        std::cout << "LCE-NGS: " << r.ngs_us << " us\n";
        std::cout << "Ratio:   " << std::setprecision(4) << (r.nss_us > 0.0 ? r.ngs_us / r.nss_us : 0.0) << "\n";
        return 0;
    }

    if (mode == "bench" || mode == "bench_profile") {
        int runs = (argc >= 3) ? std::stoi(argv[2]) : 5;
        const int sizes[] = {1000, 5000, 10000, 50000, 100000, 500000, 1000000, 2000000, 5000000};
        const int sigmas[] = {2, 4, 26};
        if (mode == "bench") {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio\n";
            for (int sigma : sigmas) {
                for (int n : sizes) {
                    auto s = generate_random(n, sigma, 12345 + sigma);
                    print_plain_csv_row("random", "sigma=" + std::to_string(sigma), benchmark_plain(s.data(), n, runs));
                }
            }
        } else {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio,"
                         "nss_extend_calls,ngs_extend_calls,nss_fast_returns,ngs_fast_returns,nss_reuses,ngs_reuses,"
                         "nss_while_iters,ngs_while_iters,nss_suffix_links,ngs_suffix_links,nss_char_cmps,ngs_char_cmps,"
                         "nss_assignments,ngs_assignments,nss_cleanup,ngs_cleanup,"
                         "nss_window_hits,ngs_window_hits,nss_reuse_from_next,ngs_reuse_from_next,"
                         "nss_reuse_from_prev_lce,ngs_reuse_from_prev_lce,nss_fast_return_success,ngs_fast_return_success,"
                         "nss_fast_return_to_rhs,ngs_fast_return_to_rhs,nss_rhs_fallback_total,ngs_rhs_fallback_total,"
                         "nss_fast_return_lce_total,ngs_fast_return_lce_total,nss_eq_plce_branch,ngs_eq_plce_branch,"
                         "nss_gt_plce_branch,ngs_gt_plce_branch,nss_explicit_extend_calls,ngs_explicit_extend_calls,"
                         "nss_explicit_match_cmps,ngs_explicit_match_cmps,nss_explicit_mismatch_cmps,ngs_explicit_mismatch_cmps\n";
            for (int sigma : sigmas) {
                for (int n : sizes) {
                    auto s = generate_random(n, sigma, 12345 + sigma);
                    print_profile_csv_row("random", "sigma=" + std::to_string(sigma), benchmark_profile(s.data(), n, runs));
                }
            }
        }
        return 0;
    }

    if (mode == "bench_struct" || mode == "bench_struct_profile") {
        int runs = (argc >= 3) ? std::stoi(argv[2]) : 5;
        const int sizes[] = {1000, 5000, 10000, 50000, 100000, 500000, 1000000, 2000000, 5000000};
        if (mode == "bench_struct") {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio\n";
            for (int n : sizes) {
                auto s0 = make_repetitive_constant(n);
                print_plain_csv_row("repetitive_constant", "a", benchmark_plain(s0.data(), n, runs));
                auto s1 = make_repetitive_periodic(n);
                print_plain_csv_row("repetitive_periodic", "ab", benchmark_plain(s1.data(), n, runs));
                auto s2 = make_long_border(n, 25);
                print_plain_csv_row("long_border", "25pct", benchmark_plain(s2.data(), n, runs));
                auto s3 = make_long_border(n, 40);
                print_plain_csv_row("long_border", "40pct", benchmark_plain(s3.data(), n, runs));
                auto s4 = make_quasi_monotone(n, 1);
                print_plain_csv_row("quasi_monotone", "noise1pct", benchmark_plain(s4.data(), n, runs));
                auto s5 = make_quasi_monotone(n, 5);
                print_plain_csv_row("quasi_monotone", "noise5pct", benchmark_plain(s5.data(), n, runs));
            }
        } else {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio,"
                         "nss_extend_calls,ngs_extend_calls,nss_fast_returns,ngs_fast_returns,nss_reuses,ngs_reuses,"
                         "nss_while_iters,ngs_while_iters,nss_suffix_links,ngs_suffix_links,nss_char_cmps,ngs_char_cmps,"
                         "nss_assignments,ngs_assignments,nss_cleanup,ngs_cleanup,"
                         "nss_window_hits,ngs_window_hits,nss_reuse_from_next,ngs_reuse_from_next,"
                         "nss_reuse_from_prev_lce,ngs_reuse_from_prev_lce,nss_fast_return_success,ngs_fast_return_success,"
                         "nss_fast_return_to_rhs,ngs_fast_return_to_rhs,nss_rhs_fallback_total,ngs_rhs_fallback_total,"
                         "nss_fast_return_lce_total,ngs_fast_return_lce_total,nss_eq_plce_branch,ngs_eq_plce_branch,"
                         "nss_gt_plce_branch,ngs_gt_plce_branch,nss_explicit_extend_calls,ngs_explicit_extend_calls,"
                         "nss_explicit_match_cmps,ngs_explicit_match_cmps,nss_explicit_mismatch_cmps,ngs_explicit_mismatch_cmps\n";
            for (int n : sizes) {
                auto s0 = make_repetitive_constant(n);
                print_profile_csv_row("repetitive_constant", "a", benchmark_profile(s0.data(), n, runs));
                auto s1 = make_repetitive_periodic(n);
                print_profile_csv_row("repetitive_periodic", "ab", benchmark_profile(s1.data(), n, runs));
                auto s2 = make_long_border(n, 25);
                print_profile_csv_row("long_border", "25pct", benchmark_profile(s2.data(), n, runs));
                auto s3 = make_long_border(n, 40);
                print_profile_csv_row("long_border", "40pct", benchmark_profile(s3.data(), n, runs));
                auto s4 = make_quasi_monotone(n, 1);
                print_profile_csv_row("quasi_monotone", "noise1pct", benchmark_profile(s4.data(), n, runs));
                auto s5 = make_quasi_monotone(n, 5);
                print_profile_csv_row("quasi_monotone", "noise5pct", benchmark_profile(s5.data(), n, runs));
            }
        }
        return 0;
    }

    if (mode == "bench_border" || mode == "bench_border_profile") {
        int runs = (argc >= 3) ? std::stoi(argv[2]) : 3;
        const int sizes[] = {100000, 500000, 1000000};
        if (mode == "bench_border") {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio\n";
            for (int n : sizes) {
                auto s25 = make_long_border(n, 25);
                print_plain_csv_row("border_heavy", "25pct", benchmark_plain(s25.data(), n, runs));
                auto s40 = make_long_border(n, 40);
                print_plain_csv_row("border_heavy", "40pct", benchmark_plain(s40.data(), n, runs));
            }
        } else {
            std::cout << "family,n,nss_us,ngs_us,ratio,nss_char_cmps,ngs_char_cmps,nss_reuses,ngs_reuses,nss_explicit_extend_calls,ngs_explicit_extend_calls\n";
            for (int n : sizes) {
                auto s25 = make_long_border(n, 25);
                print_border_profile_csv_row("25pct_border", benchmark_profile(s25.data(), n, runs));
                auto s40 = make_long_border(n, 40);
                print_border_profile_csv_row("40pct_border", benchmark_profile(s40.data(), n, runs));
            }
        }
        return 0;
    }

    if (mode == "bench_files" || mode == "bench_files_profile") {
        if (argc < 5) {
            usage();
            return 1;
        }
        int runs = std::stoi(argv[2]);
        int maxlen = std::stoi(argv[3]);
        if (mode == "bench_files") {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio\n";
            for (int i = 4; i < argc; ++i) {
                auto s = load_file(argv[i], maxlen);
                print_plain_csv_row("file", argv[i], benchmark_plain(s.data(), (int)s.size(), runs));
            }
        } else {
            std::cout << "dataset,param,n,lce_nss_us,lce_ngs_us,recover_nss_us,recover_ngs_us,ratio,"
                         "nss_extend_calls,ngs_extend_calls,nss_fast_returns,ngs_fast_returns,nss_reuses,ngs_reuses,"
                         "nss_while_iters,ngs_while_iters,nss_suffix_links,ngs_suffix_links,nss_char_cmps,ngs_char_cmps,"
                         "nss_assignments,ngs_assignments,nss_cleanup,ngs_cleanup,"
                         "nss_window_hits,ngs_window_hits,nss_reuse_from_next,ngs_reuse_from_next,"
                         "nss_reuse_from_prev_lce,ngs_reuse_from_prev_lce,nss_fast_return_success,ngs_fast_return_success,"
                         "nss_fast_return_to_rhs,ngs_fast_return_to_rhs,nss_rhs_fallback_total,ngs_rhs_fallback_total,"
                         "nss_fast_return_lce_total,ngs_fast_return_lce_total,nss_eq_plce_branch,ngs_eq_plce_branch,"
                         "nss_gt_plce_branch,ngs_gt_plce_branch,nss_explicit_extend_calls,ngs_explicit_extend_calls,"
                         "nss_explicit_match_cmps,ngs_explicit_match_cmps,nss_explicit_mismatch_cmps,ngs_explicit_mismatch_cmps\n";
            for (int i = 4; i < argc; ++i) {
                auto s = load_file(argv[i], maxlen);
                print_profile_csv_row("file", argv[i], benchmark_profile(s.data(), (int)s.size(), runs));
            }
        }
        return 0;
    }

    usage();
    return 1;
}