library

This documentation is automatically generated by online-judge-tools/verification-helper

---

Project maintained by tko919 Hosted on GitHub Pages — Theme by mattgraham

Prime Sum
(Math/lucydp.hpp)

• View this file on GitHub
• Last update: 2025-06-29 02:34:27+00:00
• Include: #include "Math/lucydp.hpp"
• Link: View error logs on GitHub Actions

使い方

PrimeSum(ll n): テンプレートには

• T: 返り値の型
• T (*F)(ll): 完全乗法的関数 $f$ のprefix sumを返す関数 を指定。

T operator[](ll x): $\sum_{p \leq x:\mbox{prime}} f(p)$ を出力。 $x=\lfloor n/d \rfloor$ で表される必要がある。

Depends on

• Prime Sieve (Math/sieve.hpp)

Verified with

• Verify/LC_counting_primes.test.cpp
• Verify/YUKI_1781.test.cpp

Code

#pragma once
#include "Math/sieve.hpp"


template <typename T, T (*F)(ll)> struct LucyDP {
    ll N, SQ, sz;
    vector<ll> quo;
    vector<T> dat;
    static inline pair<vector<ll>, vector<int>> gen(ll n) {
        ll sq = sqrtl(n);
        vector<ll> quo(sq + n / (sq + 1));
        iota(ALL(quo), 1);
        for (ll i = sq, x = n / (sq + 1); x; x--, i++)
            quo[i] = n / x;
        auto ps = sieve(sq);
        return {quo, ps};
    }
    LucyDP() {}
    LucyDP(ll n, vector<ll> &_quo, vector<int> &ps)
        : N(n), SQ(sqrtl(N)), sz(SQ + n / (SQ + 1)), quo(_quo), dat(sz) {
        rep(i, 0, sz) dat[i] = F(quo[i]) - 1;
        for (auto &p : ps) {
            T coe = dat[p - 1] - dat[p - 2];
            for (int i = sz - 1;; i--) {
                if (quo[i] < ll(p) * p)
                    break;
                dat[i] -= (dat[idx(quo[i] / p)] - dat[p - 2]) * coe;
            }
        }
    }
    T operator[](ll x) {
        return dat[idx(x)];
    }

  private:
    int idx(ll x) const {
        if (x <= SQ)
            return x - 1;
        else
            return sz - N / x;
    }
};

/**
 * @brief Prime Sum
 * @docs docs/primesum.md
 */

#line 2 "Math/sieve.hpp"

template <int L = 101010101> vector<int> sieve(int N) {
    bitset<L> isp;
    int n, sq = ceil(sqrt(N));
    for (int z = 1; z <= 5; z += 4) {
        for (int y = z; y <= sq; y += 6) {
            for (int x = 1; x <= sq and (n = 4 * x * x + y * y) <= N; ++x) {
                isp[n].flip();
            }
            for (int x = y + 1; x <= sq and (n = 3 * x * x - y * y) <= N;
                 x += 2) {
                isp[n].flip();
            }
        }
    }
    for (int z = 2; z <= 4; z += 2) {
        for (int y = z; y <= sq; y += 6) {
            for (int x = 1; x <= sq and (n = 3 * x * x + y * y) <= N; x += 2) {
                isp[n].flip();
            }
            for (int x = y + 1; x <= sq and (n = 3 * x * x - y * y) <= N;
                 x += 2) {
                isp[n].flip();
            }
        }
    }
    for (int y = 3; y <= sq; y += 6) {
        for (int z = 1; z <= 2; ++z) {
            for (int x = z; x <= sq and (n = 4 * x * x + y * y) <= N; x += 3) {
                isp[n].flip();
            }
        }
    }
    for (int n = 5; n <= sq; ++n)
        if (isp[n]) {
            for (int k = n * n; k <= N; k += n * n) {
                isp[k] = false;
            }
        }
    isp[2] = isp[3] = true;

    vector<int> ret;
    for (int i = 2; i <= N; i++)
        if (isp[i]) {
            ret.push_back(i);
        }
    return ret;
}

/**
 * @brief Prime Sieve
 */
#line 3 "Math/lucydp.hpp"

template <typename T, T (*F)(ll)> struct LucyDP {
    ll N, SQ, sz;
    vector<ll> quo;
    vector<T> dat;
    static inline pair<vector<ll>, vector<int>> gen(ll n) {
        ll sq = sqrtl(n);
        vector<ll> quo(sq + n / (sq + 1));
        iota(ALL(quo), 1);
        for (ll i = sq, x = n / (sq + 1); x; x--, i++)
            quo[i] = n / x;
        auto ps = sieve(sq);
        return {quo, ps};
    }
    LucyDP() {}
    LucyDP(ll n, vector<ll> &_quo, vector<int> &ps)
        : N(n), SQ(sqrtl(N)), sz(SQ + n / (SQ + 1)), quo(_quo), dat(sz) {
        rep(i, 0, sz) dat[i] = F(quo[i]) - 1;
        for (auto &p : ps) {
            T coe = dat[p - 1] - dat[p - 2];
            for (int i = sz - 1;; i--) {
                if (quo[i] < ll(p) * p)
                    break;
                dat[i] -= (dat[idx(quo[i] / p)] - dat[p - 2]) * coe;
            }
        }
    }
    T operator[](ll x) {
        return dat[idx(x)];
    }

  private:
    int idx(ll x) const {
        if (x <= SQ)
            return x - 1;
        else
            return sz - N / x;
    }
};

/**
 * @brief Prime Sum
 * @docs docs/primesum.md
 */

Back to top page