OI

min25 筛

简述

时隔一年之后再学一个新的筛法是不是没救了啊。

主要原理

min25 筛的主要思路来源于埃式筛法,通过神奇的 DP 数组进行预处理之后,再用很短的时间就可以求出整个前缀和了。最终的主要思路是快速计算 \(f(p^k)\) 来表示要算的积性函数 \(h(x)\),所以 min25 筛求解 \(h(x)\) 可能需要计算多个不同的 \(f(p^k)\) 前缀和来进行对 \(h(x)\) 的组合表示。

首先需要用线性筛算出足够多的 \(pre[n] = \sum_{i = 1}^n prime_i^k\)。假设我们现在要求的一个 \(f(p^k)\) 为 \(f(x) = \sum_{i = 1} [i \in prime] i^k\),那么我们可以来设一个 DP (别问我为什么要这么干,真的不知道咋想出来的)。

定义 \(g(n, i)\) 为前 \(n\) 个「质数或者是最小值因子大于 \(prime_i\) 的数」的 \(k\) 次方之和,可以理解为埃式筛法筛完 \(i\) 轮后剩下的局部解。最后 \(\sum_{i = 1}^n [i \in primes] i^k = g(n, \sqrt{n})\)。

乍一看如果直接 DP 的话我们的时间复杂度是 \(\Theta(n \sqrt n)\) 的,但是我们发现其实这个 DP 的转移不需要依赖于 \(n – 1\),而是从 \(\lfloor \frac{n}{primes_i} \rfloor\) 转移。所以最后的层数最多 \(\sqrt n\)层,然而时间复杂度是 \( \Theta(\frac{n^{\frac{3}{4}}}{\log_2 n}) \),我不会算(zzt 的论文写了证明,神仙可以去看看)。

如何进行转移呢?如果假设这个 \(primes_i^2 > n\) 了,那么显然这一轮埃筛只筛了空气出来;否则,我们需要认真思考下怎样去筛掉那些最小值因子恰好为 \(primes_i\) 的数的 \(k\) 次方。

首先,我们已经知道要筛掉的数是最小值因子恰好为 \(primes_i\) 的数的 \(k\) 次方。假设为 \(x^k\),那么可以被表示为 \(primes_i^k x’^k\)。前面那个 \(primes_i^k\) 可以直接提出来,关键是怎么去算 \(x’^k\)。可以考虑直接拿 \(n\) 去暴力除掉一个 \(primes_i\),拿 \(g(\lfloor \frac{n}{primes_i} \rfloor, i – 1)\) 即可。

不过要注意一点,最小值因子恰好为 \(primes_i\) 的数的 \(k\) 次方中,虽然除去了 \(primes_i)\) 得到 \(g(\lfloor \frac{n}{primes_i} \rfloor, i – 1)\),但是也要注意保持条件「最小值因子恰好为 \(primes_i\) 的数」,我们需要从 \(g(\lfloor \frac{n}{primes_i} \rfloor, i – 1)\) 去掉前 \(i – 1\) 个质数的 \(k\) 次方来保持答案的正确性。

整理得到:

\[ g(n, i) = \begin{cases} g(n, i – 1), primes_i^2 > n \\  g(n, i – 1) – primes_i^k (g(\lfloor \frac{n}{primes_i} \rfloor, i – 1) – pre[i – 1]) \end{cases} \]

所以你会发现,那个 \(n\) 其实就被优化成了 \(\sqrt n\)。我们只需要存 \(\sqrt n\) 个即可。

算完这个之后,考虑算最后的答案。设计任务 \(solve(n, nid)\),考虑当前要算的数的最小值因子编号 \(nid\),然后先考虑把质数的答案算干净,为 \( g(n, maxPrimeId) – pre[nid – 1] \)。再考虑算合数的贡献,枚举下一个质因子 \(primes_j\),因为是积性函数,所以从 \(solve(\frac{n}{primes_j}, j + 1)\) 拿到之后直接乘上一个 \(f(primes_j^k)\)。

总结为:

\[ solve(n, nid) = g(n, maxPrimeId) – pre[nid – 1] + \sum_{j = i}^{primes_j^2 \leq n} \sum_{k = 1}^{primes_j^{k + 1} \leq n} (solve(\lfloor \frac{n}{primes_j^k} \rfloor, nid + 1) * f(primes_j^k) + f(primes_j^{k + 1}) ) \]

代码

// min25.cpp
// This code file is a template of min25 sieve with rich comments for easier understanding;
// This code is based on the solution of problem https://www.luogu.com.cn/problem/P5325;

// Mathematical background:
// The goal is to get {\sum_{i = 1}^n f(i), i \in [1, 10^10]}, which f(p^k) = p^k(p^k - 1).
// Of course, the function f is multiplicative function.

#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int MAX_N = 2e5 + 200, mod = 1e9 + 7;

// primes[] : used by linear sieve;
// primePre[] : \sum_{i = 0} [i \in prime] i;
// primePre[] : \sum_{i = 0} [i \in prime] i^2;
// idx[2][] : locations;
// facts[] & ftot : factor stored;
// g[2][] : the original g[][] with space saved;

int primes[MAX_N], tot, primePre[MAX_N], primePre2[MAX_N], idx[2][MAX_N], ftot, g[2][MAX_N];
ll facts[MAX_N];
bool vis[MAX_N];
ll n;

// helper function to locate the factor address;

struct locator
{
    int &operator[](const ll &x) { return x < MAX_N ? idx[0][x] : idx[1][n / x]; }
} loc;

int fpow(int bas, int tim)
{
    int ret = 1;
    while (tim)
    {
        if (tim & 1)
            ret = 1LL * ret * bas % mod;
        bas = 1LL * bas * bas % mod;
        tim >>= 1;
    }
    return ret;
}

const int inv2 = fpow(2, mod - 2), inv3 = fpow(3, mod - 2);

// preprocess;

void sieve()
{
    for (int i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            primes[++tot] = i, primePre[tot] = (0LL + primePre[tot - 1] + i) % mod,
            primePre2[tot] = (0LL + primePre2[tot - 1] + 1LL * i * i % mod) % mod;
        for (int j = 1; j <= tot && 1LL * i * primes[j] < MAX_N; j++)
        {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0)
                break;
        }
    }
}

int sieve(ll nbound, int nid)
{
    if (primes[nid] >= nbound)
        return 0;
    // answer would be splited into 2 parts;
    // the \sum_{i = 1} f(i) [i \in primes] part, and the rest part;
    int id = loc[nbound], nppart = 0, ppart = (0LL + g[1][id] + mod - g[0][id] + mod - (0LL + primePre2[nid] + mod - primePre[nid]) % mod) % mod;
    for (int i = nid + 1; i <= tot && 1LL * primes[i] * primes[i] <= nbound; i++)
        for (ll j = 1, acc = primes[i]; acc <= nbound; j++, acc *= primes[i])
            nppart = (0LL + nppart + 1LL * acc % mod * ((acc - 1) % mod) % mod * ((0LL + sieve(nbound / acc, i) + (j != 1)) % mod) % mod) % mod;
    return (0LL + nppart + ppart) % mod;
}

int main()
{
    scanf("%lld", &n), sieve();
    for (ll i = 1, gx; i <= n; i = gx + 1)
    {
        gx = n / (n / i);
        int id = ++ftot;
        facts[id] = n / i;
        // x^1;
        g[0][ftot] = (1LL * facts[id] % mod * ((facts[id] + 1) % mod) % mod * inv2 % mod + mod - 1) % mod;
        // x^2;
        g[1][ftot] = (1LL * facts[id] % mod * ((facts[id] + 1) % mod) % mod * inv2 % mod * (2LL * facts[id] % mod + 1) % mod * inv3 % mod + mod - 1) % mod;
        loc[facts[id]] = id;
    }
    // do dp;
    for (int i = 1; i <= tot; i++)
        for (int j = 1; j <= ftot && 1LL * primes[i] * primes[i] <= facts[j]; j++)
        {
            int pos = loc[facts[j] / primes[i]];
            g[0][j] = (0LL + g[0][j] + mod - 1LL * primes[i] * ((0LL + g[0][pos] + mod - primePre[i - 1]) % mod) % mod) % mod;
            g[1][j] = (0LL + g[1][j] + mod - 1LL * primes[i] * primes[i] % mod * ((0LL + g[1][pos] + mod - primePre2[i - 1]) % mod) % mod) % mod;
        }
    printf("%lld\n", (1LL + sieve(n, 0)) % mod);
    return 0;
}

例题 A:经典积性函数筛 – BZOJ 3944 Sum

先考虑欧拉函数 \(\varphi(x)\) 的前缀和。发现 \(\varphi(p^k) = p^k – p^{k – 1} = p^{k – 1}(p – 1)\),所以我们只需要维护质数积和即可。

再考虑莫比乌斯函数 \(\mu(x)\),其中发现合数肯定都是 \(0\),唯一贡献的只有 \(x = \prod p_i\),所以只需要算区间内的质数个数即可,因为只有 \(\mu(p^k), k = 1\) 有贡献。

// BZ3944.cpp
#include <bits/stdc++.h>

using namespace std;

typedef long long ll;

const int MAX_N = 2e5 + 200;

int T, primes[MAX_N], tot, idx[2][MAX_N];
ll g[2][MAX_N], pre[2][MAX_N], n, facts[MAX_N], ftot;
bool vis[MAX_N];

struct locator
{
    int &operator[](const ll &x) { return x < MAX_N ? idx[0][x] : idx[1][n / x]; }
} loc;

void sieve()
{
    for (int i = 2; i < MAX_N; i++)
    {
        if (!vis[i])
            primes[++tot] = i, pre[0][tot] = pre[0][tot - 1] + i, pre[1][tot] = -tot;
        for (int j = 1; j <= tot && 1LL * i * primes[j] < MAX_N; j++)
        {
            vis[i * primes[j]] = true;
            if (i % primes[j] == 0)
                break;
        }
    }
}

ll sieve_phi(ll nbound, int nid)
{
    if (primes[nid] >= nbound || nbound <= 1)
        return 0;
    int id = loc[nbound];
    ll ret = g[0][id] - g[1][id] - pre[0][nid] + nid;
    for (int i = nid + 1; i <= tot && 1LL * primes[i] * primes[i] <= nbound; i++)
        for (ll acc = primes[i] * primes[i], phi_acc = primes[i] - 1; acc <= nbound;)
        {
            ret += 1LL * phi_acc * sieve_phi(nbound / (acc / primes[i]), i) + phi_acc * primes[i];
            acc *= primes[i], phi_acc *= primes[i];
        }
    return ret;
}

ll sieve_mu(ll nbound, int nid)
{
    if (primes[nid] >= nbound || nbound <= 1)
        return 0;
    int id = loc[nbound];
    ll ret = -(g[1][id] - nid);
    for (int i = nid + 1; i <= tot && 1LL * primes[i] * primes[i] <= nbound; i++)
        ret -= sieve_mu(nbound / primes[i], i);
    return ret;
}

void solve()
{
    ftot = 0;
    for (ll i = 1, gx; i <= n; i = gx + 1)
    {
        gx = n / (n / i);
        int id = ++ftot;
        facts[id] = n / i;
        // prime sum & number;
        g[0][ftot] = (facts[id] + 2) * (facts[id] - 1) / 2, g[1][ftot] = facts[id] - 1;
        // prime num;
        loc[facts[id]] = id;
    }
    for (int i = 1; i <= tot; i++)
        for (int j = 1; j <= ftot && 1LL * primes[i] * primes[i] <= facts[j]; j++)
        {
            int pos = loc[facts[j] / primes[i]];
            g[0][j] -= 1LL * primes[i] * (0LL + g[0][pos] - pre[0][i - 1]);
            g[1][j] -= g[1][pos] - i + 1;
        }
    if (n == 0)
        puts("0 0");
    else
        printf("%lld %lld\n", sieve_phi(n, 0) + 1, sieve_mu(n, 0) + 1);
}

int main()
{
    scanf("%d", &T), sieve();
    while (T--)
        scanf("%lld", &n), solve();
    return 0;
}

 


kal0rona

http://kaloronahuang.com

江西师大附中全机房最弱