Polynomial/PrefixSumH.hpp

• View this file on GitHub
• Last update: 2023-04-12 00:49:31+08:00
• Include: #include "Polynomial/PrefixSumH.hpp"

Depends on

• Polynomial/Ntt.hpp
• Polynomial/Poly.hpp
• Template/Power.hpp
• Template/Template.hpp

Code

#pragma once

#include "../Template/Template.hpp"
#include "Poly.hpp"
template <class T>
T PrefixSumH(i64 n) {
    /*
        \sum_{i=0}^n frac{1}{i}
        Time Complexity: O(\sqrt(n)\log n)
    */
    assert(n != 0);
    if (n == 1) return 1;
    using poly = Poly<T>;
    i64 v = 1;
    while (v * v < n) v *= 2;
    T iv = T(v).inverse();
    std::array<poly, 3> G = {poly({1, v + 1}), poly({1, 1}), poly({1, v + 1})};
    for (i64 d = 1; d != v; d <<= 1ll) {
        std::array<poly, 3> G1, G2, G3;
        for (int i = 0; i < 3; i++) {
            G1[i] = G[i].shift(T(d) * iv);
            G2[i] = G[i].shift(T(d * v + v) * iv);
            G3[i] = G[i].shift(T(d * v + d + v) * iv);
        }

        for (int t = 0; t <= d; t++) {
            G[1][t] = G1[1][t] * G[0][t] + G[1][t] * G1[2][t];
            G2[1][t] = G2[1][t] * G3[2][t] + G3[1][t] * G2[0][t];
            G[0][t] *= G1[0][t];
            G[2][t] *= G1[2][t];
            G2[0][t] *= G3[0][t];
            G2[2][t] *= G3[2][t];
        }
        for (int i = 0; i < 3; i++) {
            copy(begin(G2[i].a), end(G2[i].a) - 1, back_inserter(G[i].a));
        }
    }
    using M = std::array<std::array<T, 2>, 2>;

    M res = {std::array<T, 2>({1, 0}), std::array<T, 2>({0, 1})};
    auto M_multiply = [&](const M &a, const M &b) -> M {
        M res;
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                for (int k = 0; k < 2; k++) res[i][j] += a[i][k] * b[k][j];
        return res;
    };
    i64 i = 0;
    while (i + v <= n) {
        M now = {std::array<T, 2>({G[0][i / v], 0}),
                 std::array<T, 2>({G[1][i / v], G[2][i / v]})};
        res = M_multiply(now, res);
        i += v;
    }
    while (i < n) {
        M now = {std::array<T, 2>({i + 1, 0}), std::array<T, 2>({1, i + 1})};
        i++;
        res = M_multiply(now, res);
    }
    return res[1][0] / res[0][0];
}

#line 2 "Polynomial/PrefixSumH.hpp"

#line 2 "Template/Template.hpp"

#include <bits/stdc++.h>

using i64 = std::int64_t;
#line 2 "Polynomial/Poly.hpp"

#line 2 "Polynomial/Ntt.hpp"

#line 1 "Template/Power.hpp"
template <class T>
T power(T a, int b) {
    T res = 1;
    for (; b; b /= 2, a *= a) {
        if (b % 2) {
            res *= a;
        }
    }
    return res;
}
#line 5 "Polynomial/Ntt.hpp"
template <class Z>
struct NTT {
    std::vector<int> rev;
    std::vector<Z> roots{0, 1};

    static constexpr int getRoot() {
        auto _mod = Z::get_mod();
        using u64 = uint64_t;
        u64 ds[32] = {};
        int idx = 0;
        u64 m = _mod - 1;
        for (u64 i = 2; i * i <= m; ++i) {
            if (m % i == 0) {
                ds[idx++] = i;
                while (m % i == 0) m /= i;
            }
        }
        if (m != 1) ds[idx++] = m;

        int _pr = 2;
        for (;;) {
            int flg = 1;
            for (int i = 0; i < idx; ++i) {
                u64 a = _pr, b = (_mod - 1) / ds[i], r = 1;
                for (; b; a = a * a % _mod, b /= 2) {
                    if (b % 2 == 1) r = r * a % _mod;
                }
                if (r == 1) {
                    flg = 0;
                    break;
                }
            }
            if (flg == 1) break;
            ++_pr;
        }
        return _pr;
    };

    static constexpr int rt = getRoot();

    void dft(std::vector<Z> &a) {
        int n = a.size();

        if (int(rev.size()) != n) {
            int k = __builtin_ctz(n) - 1;
            rev.resize(n);
            for (int i = 0; i < n; i++) {
                rev[i] = rev[i >> 1] >> 1 | (i & 1) << k;
            }
        }

        for (int i = 0; i < n; i++) {
            if (rev[i] < i) {
                std::swap(a[i], a[rev[i]]);
            }
        }
        if (int(roots.size()) < n) {
            int k = __builtin_ctz(roots.size());
            roots.resize(n);
            while ((1 << k) < n) {
                Z e = power(Z(rt), (Z::get_mod() - 1) >> (k + 1));
                for (int i = 1 << (k - 1); i < (1 << k); i++) {
                    roots[2 * i] = roots[i];
                    roots[2 * i + 1] = roots[i] * e;
                }
                k++;
            }
        }
        for (int k = 1; k < n; k *= 2) {
            for (int i = 0; i < n; i += 2 * k) {
                for (int j = 0; j < k; j++) {
                    Z u = a[i + j];
                    Z v = a[i + j + k] * roots[k + j];
                    a[i + j] = u + v;
                    a[i + j + k] = u - v;
                }
            }
        }
    }
    void idft(std::vector<Z> &a) {
        int n = a.size();
        reverse(a.begin() + 1, a.end());
        dft(a);
        Z inv = (1 - Z::get_mod()) / n;
        for (int i = 0; i < n; i++) {
            a[i] *= inv;
        }
    }
    std::vector<Z> multiply(std::vector<Z> a, std::vector<Z> b) {
        int sz = 1, tot = a.size() + b.size() - 1;

        if (tot <= 20) {
            std::vector<Z> ret(tot);
            for (size_t i = 0; i < a.size(); i++)
                for (size_t j = 0; j < b.size(); j++) ret[i + j] += a[i] * b[j];
            return ret;
        }

        while (sz < tot) {
            sz *= 2;
        }

        a.resize(sz), b.resize(sz);
        dft(a), dft(b);

        for (int i = 0; i < sz; ++i) {
            a[i] = a[i] * b[i];
        }

        idft(a);
        a.resize(tot);
        return a;
    }
};
#line 4 "Polynomial/Poly.hpp"

template <class Z>
struct Poly {
    std::vector<Z> a;
    Poly() {}
    Poly(int sz, Z val) { a.assign(sz, val); }
    Poly(const std::vector<Z> &a) : a(a) {}
    Poly(const std::initializer_list<Z> &a) : a(a) {}
    int size() const { return a.size(); }
    void resize(int n) { a.resize(n); }
    Z operator[](int idx) const {
        if (idx < size()) {
            return a[idx];
        } else {
            return 0;
        }
    }
    Z &operator[](int idx) { return a[idx]; }
    Poly mulxk(int k) const {
        auto b = a;
        b.insert(b.begin(), k, 0);
        return Poly(b);
    }
    Poly modxk(int k) const {
        k = std::min(k, size());
        return Poly(std::vector<Z>(a.begin(), a.begin() + k));
    }
    Poly divxk(int k) const {
        if (size() <= k) {
            return Poly();
        }
        return Poly(std::vector<Z>(a.begin() + k, a.end()));
    }
    friend Poly operator+(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] + b[i];
        }
        return Poly(res);
    }
    friend Poly operator-(const Poly &a, const Poly &b) {
        std::vector<Z> res(std::max(a.size(), b.size()));
        for (int i = 0; i < int(res.size()); i++) {
            res[i] = a[i] - b[i];
        }
        return Poly(res);
    }

    friend Poly operator*(Poly a, Poly b) {
        if (a.size() == 0 || b.size() == 0) {
            return Poly();
        }
        static NTT<Z> ntt;
        return ntt.multiply(a.a, b.a);
    }
    friend Poly operator*(Z a, Poly b) {
        for (int i = 0; i < int(b.size()); i++) {
            b[i] *= a;
        }
        return b;
    }
    friend Poly operator*(Poly a, Z b) {
        for (int i = 0; i < int(a.size()); i++) {
            a[i] *= b;
        }
        return a;
    }
    Poly &operator+=(Poly b) { return (*this) = (*this) + b; }
    Poly &operator-=(Poly b) { return (*this) = (*this) - b; }
    Poly &operator*=(Poly b) { return (*this) = (*this) * b; }
    Poly operator/(const Poly &r) const { return Poly(this->a) /= r; }
    Poly rev(int deg = -1) const {
        Poly ret(this->a);
        if (deg != -1) ret.a.resize(deg, Z(0));
        reverse(begin(ret.a), end(ret.a));
        return ret;
    }
    Poly &operator/=(const Poly &r) {
        if (this->size() < r.size()) {
            this->a.clear();
            return *this;
        }
        int n = this->size() - r.size() + 1;
        return *this = (rev().modxk(n) * r.rev().inv(n)).modxk(n).rev(n);
    }
    Poly deriv() const {
        if (a.empty()) {
            return Poly();
        }
        std::vector<Z> res(size() - 1);
        for (int i = 0; i < size() - 1; ++i) {
            res[i] = Z(i + 1) * a[i + 1];
        }
        return Poly(res);
    }
    Poly integr() const {
        std::vector<Z> res(size() + 1);
        for (int i = 0; i < size(); ++i) {
            res[i + 1] = a[i] / (i + 1);
        }
        return Poly(res);
    }
    Poly inv(int m) const {
        Poly x{a[0].inverse()};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{2} - modxk(k) * x)).modxk(k);
        }
        return x.modxk(m);
    }
    Poly log(int m) const { return (deriv() * inv(m)).integr().modxk(m); }
    Poly exp(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x * (Poly{1} - x.log(k) + modxk(k))).modxk(k);
        }
        return x.modxk(m);
    }
    Poly pow(int k, int m) const {
        int i = 0;
        while (i < size() && a[i].get() == 0) {
            i++;
        }
        if (i == size() || 1LL * i * k >= m) {
            return Poly(std::vector<Z>(m));
        }
        Z v = a[i];
        auto f = divxk(i) * v.inverse();
        return (f.log(m - i * k) * k).exp(m - i * k).mulxk(i * k) * power(v, k);
    }
    Poly sqrt(int m) const {
        Poly x{1};
        int k = 1;
        while (k < m) {
            k *= 2;
            x = (x + (modxk(k) * x.inv(k)).modxk(k)) * ((Z::get_mod() + 1) / 2);
        }
        return x.modxk(m);
    }
    Poly mulT(Poly b) const {
        if (b.size() == 0) {
            return Poly();
        }
        int n = b.size();
        reverse(b.a.begin(), b.a.end());
        return ((*this) * b).divxk(n - 1);
    }
    std::vector<Z> eval(std::vector<Z> x) const {
        if (size() == 0) {
            return std::vector<Z>(x.size(), 0);
        }
        const int n = std::max(int(x.size()), size());
        std::vector<Poly> q(4 * n);
        std::vector<Z> ans(x.size());
        x.resize(n);
        std::function<void(int, int, int)> build = [&](int p, int l, int r) {
            if (r - l == 1) {
                q[p] = Poly{1, -x[l]};
            } else {
                int m = (l + r) / 2;
                build(2 * p, l, m);
                build(2 * p + 1, m, r);
                q[p] = q[2 * p] * q[2 * p + 1];
            }
        };
        build(1, 0, n);
        std::function<void(int, int, int, const Poly &)> work =
            [&](int p, int l, int r, const Poly &num) {
                if (r - l == 1) {
                    if (l < int(ans.size())) {
                        ans[l] = num[0];
                    }
                } else {
                    int m = (l + r) / 2;
                    work(2 * p, l, m, num.mulT(q[2 * p + 1]).modxk(m - l));
                    work(2 * p + 1, m, r, num.mulT(q[2 * p]).modxk(r - m));
                }
            };
        work(1, 0, n, mulT(q[1].inv(n)));
        return ans;
    }
    Poly inter(const Poly &y) const {
        std::vector<Poly> Q(a.size() * 4), P(a.size() * 4);
        std::function<void(int, int, int)> build = [&](int p, int l, int r) {
            int m = (l + r) >> 1;
            if (l == r) {
                Q[p] = Poly{-a[m], Z(1)};
            } else {
                build(p * 2, l, m);
                build(p * 2 + 1, m + 1, r);
                Q[p] = Q[p * 2] * Q[p * 2 + 1];
            }
        };
        build(1, 0, a.size() - 1);
        Poly f;
        f.a.resize((int)(Q[1].size()) - 1);
        for (int i = 0; i + 1 < Q[1].size(); i += 1)
            f[i] = (Q[1][i + 1] * (i + 1));
        Poly g = f.eval(a);
        std::function<void(int, int, int)> work = [&](int p, int l, int r) {
            int m = (l + r) >> 1;
            if (l == r) {
                P[p].a.push_back(y[m] * power(g[m], Z::get_mod() - 2));
                return;
            }
            work(p * 2, l, m);
            work(p * 2 + 1, m + 1, r);
            P[p].a.resize(r - l + 1);
            Poly A = P[p * 2] * Q[p * 2 + 1];
            Poly B = P[p * 2 + 1] * Q[p * 2];
            for (int i = 0; i <= r - l; i++) P[p][i] = (A[i] + B[i]);
        };
        work(1, 0, a.size() - 1);
        return P[1];
    }
    Poly shift(Z t, int m = -1) {
        /*
            input: y(0) , y(1) , ... , y(n-1)
            output: y(t) , t(t+1) , ... ,y (t+m-1)
            ## m = -1 => m = n
        */
        if (m == -1) m = this->size();
        i64 T = t.get();
        int k = (int)(this->size()) - 1;
        T %= Z::get_mod();

        if (T <= k) {
            Poly ret(m, 0);
            int ptr = 0;
            for (i64 i = T; i <= k && ptr < m; i++) ret[ptr++] = a[i];
            if (k + 1 < T + m) {
                auto suf = shift(k + 1, m - ptr);
                for (int i = k + 1; i < T + m; i++)
                    ret[ptr++] = suf[i - (k + 1)];
            }
            return ret;
        }
        if (T + m > Z::get_mod()) {
            auto pref = shift(T, Z::get_mod() - T);
            auto suf = shift(0, m - pref.size());
            copy(begin(suf.a), end(suf.a), back_inserter(pref.a));
            return pref;
        }

        Poly finv(k + 1, 1), d(k + 1, 0);
        for (int i = 2; i <= k; i++) finv[k] *= i;
        finv[k] = Z(1) / finv[k];
        for (int i = k; i >= 1; i--) finv[i - 1] = finv[i] * i;
        for (int i = 0; i <= k; i++) {
            d[i] = finv[i] * finv[k - i] * a[i];
            if ((k - i) & 1) d[i] = -d[i];
        }

        Poly h(m + k, 0);
        for (int i = 0; i < m + k; i++) {
            h[i] = Z(1) / (T - k + i);
        }

        auto dh = d * h;
        Poly ret(m, 0);
        Z cur = T;
        for (int i = 1; i <= k; i++) cur *= T - i;
        for (int i = 0; i < m; i++) {
            ret[i] = cur * dh[k + i];
            cur *= T + i + 1;
            cur *= h[i];
        }
        return ret;
    }
};
#line 5 "Polynomial/PrefixSumH.hpp"
template <class T>
T PrefixSumH(i64 n) {
    /*
        \sum_{i=0}^n frac{1}{i}
        Time Complexity: O(\sqrt(n)\log n)
    */
    assert(n != 0);
    if (n == 1) return 1;
    using poly = Poly<T>;
    i64 v = 1;
    while (v * v < n) v *= 2;
    T iv = T(v).inverse();
    std::array<poly, 3> G = {poly({1, v + 1}), poly({1, 1}), poly({1, v + 1})};
    for (i64 d = 1; d != v; d <<= 1ll) {
        std::array<poly, 3> G1, G2, G3;
        for (int i = 0; i < 3; i++) {
            G1[i] = G[i].shift(T(d) * iv);
            G2[i] = G[i].shift(T(d * v + v) * iv);
            G3[i] = G[i].shift(T(d * v + d + v) * iv);
        }

        for (int t = 0; t <= d; t++) {
            G[1][t] = G1[1][t] * G[0][t] + G[1][t] * G1[2][t];
            G2[1][t] = G2[1][t] * G3[2][t] + G3[1][t] * G2[0][t];
            G[0][t] *= G1[0][t];
            G[2][t] *= G1[2][t];
            G2[0][t] *= G3[0][t];
            G2[2][t] *= G3[2][t];
        }
        for (int i = 0; i < 3; i++) {
            copy(begin(G2[i].a), end(G2[i].a) - 1, back_inserter(G[i].a));
        }
    }
    using M = std::array<std::array<T, 2>, 2>;

    M res = {std::array<T, 2>({1, 0}), std::array<T, 2>({0, 1})};
    auto M_multiply = [&](const M &a, const M &b) -> M {
        M res;
        for (int i = 0; i < 2; i++)
            for (int j = 0; j < 2; j++)
                for (int k = 0; k < 2; k++) res[i][j] += a[i][k] * b[k][j];
        return res;
    };
    i64 i = 0;
    while (i + v <= n) {
        M now = {std::array<T, 2>({G[0][i / v], 0}),
                 std::array<T, 2>({G[1][i / v], G[2][i / v]})};
        res = M_multiply(now, res);
        i += v;
    }
    while (i < n) {
        M now = {std::array<T, 2>({i + 1, 0}), std::array<T, 2>({1, i + 1})};
        i++;
        res = M_multiply(now, res);
    }
    return res[1][0] / res[0][0];
}

Back to top page