min25筛

这也是个用于求积性函数前缀和的筛法。

对于积性函数 f(n),满足这样几个要求时可以用min25来做。

  • n \in prime 时, f(n) 可以使用关于 n 的多项式表示。

  • n \in prime 时,f(n^k) 要能够快速计算答案。

预处理 Σf(p)

在min25筛的过程中,我们需要预处理多项式的 x^k 相关的这个东西

\sum_{i=2}^n i^k [i \in prime]

那么记所有质数的集合为 P,第 j 个素数为 P_jx 的最小质因子为 minp(x),我们需要预处理这样一个东西

g(n,j) = \sum_{i=2}^n i^k [i \in P ~ or ~ minp(i) > P_j]

举个例子吧

\begin{aligned} g(n,0) &= 2^k + 3^k + 4^k + 5^k + 6^k + 7^k + \cdots + n^k \\ g(2n,1) &= 2^k + 3^k + 5^k + 7^k + 9^k + 11^k + \cdots + (2n-1)^k \\ g(6n,2) &= 2^k + 3^k + 5^k + 7^k + 11^k + \cdots + (6n-5)^k + (6n-1)^k \end{aligned}

相当于在这个求和中,g(n,j) 将前 j 个素数的倍数全部去掉了,只留下了那些素数的 k 次方和。

\begin{aligned} g(2n,0) – g(2n,1) &= 4^k + 6^k + 8^k + 10^k + 12^k + \cdots + (2n)^k \\ &= 2^k \left( 2^k + 3^k + 4^k + 5^k + 6^k + \cdots + n^k \right) \\ \\ g(6n,1) – g(6n,2) &= 9^k + 15^k + 21^k + 27^k + 33^k + \cdots + (6n-3)^k \\ &= 3^k \left( 3^k + 5^k + 7^k + 9^k + 11^k + \cdots + (2n-1)^k \right) \end{aligned}

那么从 g(n,j-1) 转移到 g(n,j) 的过程中,我们考虑去掉 p_j^k 的倍数 \sum (p_j \prod p_{j+i}^{e_i})^k

差不多就是有这样一个转移

g(n,j) = g(n,j-1) – p_j^k \left( g\left(\left\lfloor\frac{n}{p_j}\right\rfloor, j-1\right) – \sum_{x=1}^{j-1} p_x^k \right)

为什么要减去前 j-1 个素数 k 次方和的转移呢?因为 (p_x p_j)^k 在这之前都被筛掉了啊。

如果 p_j^2 > n,那肯定没有东西能被删掉了。所以有这样一个转移

g(n,j) = \begin{cases} g(n,j-1), &p_j^2 \gt n \\ g(n,j-1) – p_j^k \left( g(\lfloor n / p_j\rfloor, j-1) – \sum_{x=1}^{j-1} p_x^k \right), &p_j^2 \le n \end{cases}

我们需要一些转换,才能把上面那个东西从递归变成迭代。

我们会发现,我们可能需要经常性的计算 g(\lfloor n/x\rfloor, j) 的值,那我们就一起维护了吧。

这里是预处理的代码。

// euler sieve, spk for sum p^k from p_1 ~ p_pcn
int pri[MAXN], pcn, spk[MAXN];
int m, Sqrt, id1[MAXN], id2[MAXN];
long long gk[MAXN], w[MAXN];

inline long long sum_ik(lld n); // for i=2..n, sum i^k

void init()
{
    m = 0, Sqrt = sqrt(n);

    for (long long L = 1, R; L <= n; L = R+1)
    {
        R = n / (n / L), w[++m] = n / L;
        (w[m] <= Sqrt ? id1[w[m]] : id2[R]) = m;
        gk[m] = sum_ik(w[m]);
    }

    for (int j = 1; j <= pcn; j++) // p_pcn^2 should be greater than n
    {
        for (int i = 1; i <= m && pri[j] <= w[i]/pri[j]; i++)
        {
            int d = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[n/(w[i]/pri[j])];
            gk[i] -= (spk[j] - spk[j-1]) * (gk[d] - spk[j-1]);
        }
    }
}

其中,id1[i] 表示数论分块中整除值 \lfloor n/l \rfloor = i 的块编号,id2[n/i] 表示整除值 \lfloor n / l \rfloor = i 的块编号,w[s] 表示数论分块第 s 个块的整除值。这样比较节省空间。

然后考虑下面那个迭代。首先,w_i 就是 g 函数的第一个参数。

p_j^2 > w_i,也就是第 i 个块整除值的时候,没有必要继续处理了。内层的 i=1 \cdots m 相当于从大往小枚举块,也就不用利用滚动数组什么的记录上一轮 j-1 时的具体值了。

最后,init执行完时,gk[i] 就是第 i 个整除值 w_i 时,g(w_i, pcn) 的值了,也就是我们想要的东西了。

由于那个多项式可能会有好几个 p^x,所以都维护出来吧。

求 Σf(i)

我们刚才去掉了合数处的点值,那么我们还要根据积性函数的性质把那些点值加回来。

那么我们现在构造这样一个函数

S(n,j) = \sum_{x=2}^n f(x)~[ minp(x) \ge p_j ]

例如

\begin{aligned} S(20,3) &= f(5) + f(7) + f(11) + f(13) + f(17) + f(19) \\ S(20,2) – S(20,3) &= f(3) + f(9) + f(15) \\ S(20,1) – S(20,2) &= f(2) + f(4) + f(6) + \cdots + f(18) + f(20) \end{aligned}

我们需要从最大的质数开始,一个个往回补上点值。

S(n,j) = \sum_{k \ge j} \sum_{e} f(p_k^e)\left(1 + S\left(\left\lfloor n/p_k^e\right\rfloor, k+1\right)\right)

考虑 p_k^e 的倍数,显然当 e=0 的时候不能算上贡献;当 e>0 的时候,\sum p_k^e \prod p_{k+i}^{e_i} 的贡献由 f(p_k^e) 乘上 \sum f(\prod p_{k+i}^{e_i}) 得到。为什么还有个 1+ 呢?S(x,k) 里面可没有 p_j 的东西,不加就算不上 \sum f(p_k^e) 的贡献了呀。

那我们整理一下这个鬼一样的式子。

\begin{aligned} S(n, j) &= \sum_{k \ge j} f(p_k) + \sum_{k \ge j} \sum_{e} \left( f(p_k^e) S(\lfloor n/p_k^e\rfloor, k+1) + f(p_k^{e+1}) \right) \\ &= g(n, |P|) – \sum_{i=1}^{j-1} f(p_i) + \sum_{k \ge j} \sum_{e} \left( f(p_k^e) S(\lfloor n/p_k^e\rfloor) + f(p_k^{e+1}) \right) \end{aligned}

所以这一部分的代码可以写成

long long S(long long x, int y)
{
    if (x <= 1 || pri[y] > x) return 0;
    int d = x <= Sqrt ? id1[x] : id2[n/x];
    long long ans = g(x) - g(pri[y]-1);

    for (int i = y; i <= pcn && pri[i] <= x / pri[i]; i++)
    {
        long long p1 = pri[i], p2 = p1 * p1;
        for (int e = 1; p1 <= x / pri[i]; ++e, p1 = p2, p2 *= pri[i])
            ans += S(x/p1, i+1) * f(i,e) + f(i,e+1);
    }

    return ans;
}

其中 g(n)\sum_{i=1}^n f(i) [i \in prime],也就是上面式子中的 g(n,|P|) 可以从上面预处理的式子中得到。而 f(i,e) 就是 f(p_i^e)

最后你要求的 \sum_{i=1}^n f(i) = S(n,1) + 1

举例来说,我已知积性函数 f(p^e) = C_{e+A-1}^{A-1},那么 f(p)=A,我可以预处理上面预处理部分 sp0, g0 等内容,那么 g(x) - g(pri[y]-1) 应该写为 A * (g0[d]-(y-1)),而 f(i,e) 则写为 C(e+A-1, A-1)。对应例题为 HDU6537。

举例来说,我已知积性函数 f(p^e) = Ae+1,那么 f(p)=A+1,我可以预处理 sp0, g0 等内容,那么 g(x) - g(pri[y]-1) 应该写为 (A+1) * (g0[d]-(y-1)),而 f(i,e) 则写为 A*e+1。对应具体例题为 SPOJ DIVCNT3。

再举例来说,我已知积性函数 f(p^e) = p \oplus e,那么可以在 p=2f(2)=3,然后其他时候有 f(p)=p-1,预处理 sp0, sp1, g0, g1 等内容,那么 g(x) - g(pri[y]-1) 应该写为 g1[d]-g0[d]-sp1[y-1]+y-1,并特判 if (y == 1) ans += 2;,而 f(i,e) 则写为 pri[i]^e。对应具体例题为 LOJ6053。

后记

所以说这是不是很像埃式筛啊,所以它又被称为拓展埃式筛。

据说时间复杂度是

O\left( \frac{n^{\frac{3}{4}}}{\log n} \right)

咱也不会证,咱也不敢证。

2019杭电多校第三场 E. Easy Math Problem

题目想要对于给定 nk

\sum_{i=1}^n \sum_{j=1}^n \gcd(i,j)^k ~ lcm(i,j) ~[\gcd(i,j) \in prime]

那我们考虑一下枚举某个素数 d = \gcd(i,j) 的情况吧。

\sum_{d \in P} d^{k-1} \sum_{i=1}^n \sum_{j=1}^n ij ~ [\gcd(i,j) = d]

我们先来考虑一下

f(N) = \sum_{i=1}^N \sum_{j=1}^N ij [\gcd(i,j)=1]

先根据对称性,可以得到

\sum_{i=1}^N \sum_{j=1}^i ij [\gcd(i,j)=1] = \frac{f(N)+1}{2}

再根据非常显然的 \gcd(i,j)=1 \Leftrightarrow \gcd(i,j-i)=1,可以配对

\sum_{j=1}^i j [\gcd(i,j)=1] = \frac{i \varphi(i) + [i=1]}{2}

那么就有

\begin{aligned} \frac{f(N)+1}{2} &= \sum_{i=1}^N \sum_{j=1}^i ij [\gcd(i,j)=1] \\ &= \frac{\sum_{i=1}^N i^2 \varphi(i) +1}{2} \end{aligned}

所以说刚才求的 f(N) 就相当于对 i^2 \varphi(i) 求了个前缀和。根据杜教筛一般套路,给配上个 g(n) = n^2 就可以筛筛筛了。

f(n) = \frac{n(n+1)(2n+1)}{6} – \sum_{i=2}^n i^2 f\left(\left\lfloor\frac{n}{i}\right\rfloor\right)

再回到原式,其实相当于求

\sum_{d \in P} d^{k+1} \sum_{i=1}^{\lfloor n/d\rfloor} i^2 \varphi(i)

那么我们对 n 进行数论分块,就可以获得答案了。考虑块 [L,R] 内前面一个东西贡献的值,其实相当于

\sum_{i=L}^R i^{k+1} [i \in prime]

那么把它变成一个前缀和的差,构造一个积性函数

h(p) = p^{k+1}; h(p^e) = 0, e>1

然后考虑用min25筛,就可以做了。

由于min25筛需要快速求出 p^k 的前缀和,那么就得上拉格朗日插值一下。

推荐看 这篇文章 学习min25。我也要好好花时间学习一下啦 QAQ

#include <bits/stdc++.h>
using namespace std;
typedef long long lld;
const int MAXN = 7e6+5;
const int MOD = 1e9+7;
lld N, K, ans;

lld fpow(lld a, lld k) {
    lld b = 1;
    while (k) {
        if (k & 1) b = b * a % MOD;
        a = a * a % MOD;
        k >>= 1;
    }
    return b;
}

namespace interpolation
{
    typedef vector<lld> Poly;

    Poly& operator*=(Poly &a, const Poly &b) {
        Poly c(a.size()+b.size()-1, 0);
        for (int i = 0; i < a.size(); i++)
            for (int j = 0; j < b.size(); j++)
                (c[i+j] += a[i] * b[j] % MOD) %= MOD;
        return a = c;
    }

    Poly& operator*=(Poly &a, lld qwq) {
        for (lld &ai : a)
            ai = ai * qwq % MOD;
        return a;
    }

    Poly& operator+=(Poly &a, const Poly &b) {
        a.resize(max(a.size(), b.size()));
        for (int i = 0; i < b.size(); i++)
            (a[i] += b[i]) %= MOD;
        return a;
    }

    Poly Lagrange(const Poly &X, const Poly &FX) {
        Poly c;
        int n = X.size();
        for (int i = 0; i < n; i++) {
            Poly x({FX[i]});
            for (int j = 0; j < n; j++)
                if (j != i)
                    x *= Poly({-X[j],1}),
                    x *= fpow(X[i]-X[j], MOD-2);
            c += x;
        }
        return c;
    }

    lld evaluate(const Poly &a, lld x) {
        lld xx = 1, ans = 0;
        for (int i = 0; i < a.size(); i++)
            (ans += a[i] * xx % MOD) %= MOD,
            xx = xx * x % MOD;
        return ans;
    }

    Poly init_lag() {
        Poly a, b;
        lld last = 0;
        for (int i = 0; i <= K+2; i++) {
            a.push_back(i);
            (last += fpow(i, K+1)) %= MOD;
            b.push_back(last);
        }
        return Lagrange(a, b);
    }
}

namespace euler_seive
{
    int pri[MAXN], _B1[MAXN], pcn;
    bool isnp[MAXN];
    unordered_map<lld, int> _B2;
    const lld inv6 = fpow(6, MOD-2);

    void init_prime()
    {
        _B1[1] = 1;

        for (int i = 2; i < MAXN; i++)
        {
            if (!isnp[i])
            {
                pri[++pcn] = i;
                _B1[i] = i-1;
            }

            for (int j = 1; j <= pcn; j++)
            {
                if (1ll * i * pri[j] >= MAXN) break;
                isnp[i * pri[j]] = true;
                if (i % pri[j] == 0) {
                    _B1[i * pri[j]] = _B1[i] * pri[j];
                    break;
                } else {
                    _B1[i * pri[j]] = _B1[i] * (pri[j]-1);
                }
            }
        }

        for (int i = 2; i < MAXN; i++)
            _B1[i] = (_B1[i-1] + 1ll * _B1[i] * i % MOD * i % MOD) % MOD;
    }

    inline lld sum2(lld a) { a %= MOD; return a * (a+1) % MOD * (2*a+1) % MOD * inv6 % MOD; }

    int sumId2Phi(lld n)
    {
        if (n < MAXN) return _B1[n];
        else if (_B2.count(n)) return _B2[n];
        lld ans = n % MOD * (n%MOD+1) % MOD * ((MOD+1)/2) % MOD;
        ans = ans * ans % MOD; // \sum n^3
        for (lld l = 2, r; l <= n; l = r+1)
            r = n / (n / l),
            ans -= (sum2(r) - sum2(l-1)) * sumId2Phi(n/l) % MOD;
        ans = (ans % MOD + MOD) % MOD;
        return _B2[n] = ans;
    }
}

namespace min25_seive
{
    using euler_seive::pri;
    using euler_seive::pcn;
    interpolation::Poly _sumik_f;
    int _sumpk[MAXN], m, Sqrt;
    int id1[MAXN], id2[MAXN], tot;
    lld g[MAXN], w[MAXN];

    inline lld sumik(lld n) {
        return interpolation::evaluate(_sumik_f, n % MOD) - 1;
    }

    inline void init()
    {
        for (int i = 1; i <= pcn; i++)
            _sumpk[i] = (_sumpk[i-1] + fpow(pri[i], K+1)) % MOD;
        _sumik_f = interpolation::init_lag();

        m = 0;
        Sqrt = sqrt(N);

        for (lld l = 1, r; l <= N; l = r+1)
        {
            r = N / (N / l), w[++m] = N / l;
            (w[m] <= Sqrt ? id1[w[m]] : id2[N/w[m]]) = m;
            g[m] = sumik(w[m]);
        }

        tot = lower_bound(pri+1, pri+1+pcn, Sqrt) - pri;
        for (int j = 1; j <= tot; j++)
        {
            for (int i = 1; i <= m && pri[j] <= w[i]/pri[j]; i++)
            {
                int k = (w[i]/pri[j]<=Sqrt) ? id1[w[i]/pri[j]] : id2[N/(w[i]/pri[j])];
                (g[i] -= (_sumpk[j] - _sumpk[j-1]) * (g[k] - _sumpk[j-1]) % MOD) %= MOD;
            }
        }
    }

    inline lld sumpk(lld n) {
        return g[n <= Sqrt ? id1[n] : id2[N/n]];
    }
}

void solve()
{
    scanf("%lld %lld", &N, &K);
    min25_seive::init();

    ans = 0;
    for (lld l = 1, r, s = 0; l <= N; l = r+1)
    {
        r = N / (N / l);
        lld ss = min25_seive::sumpk(r);
        ans += euler_seive::sumId2Phi(N/l) * (ss-s) % MOD;
        s = ss;
    }

    printf("%lld\n", (ans % MOD + MOD) % MOD);
}

int main()
{
    euler_seive::init_prime();
    int T;
    scanf("%d", &T);
    while (T--) solve();
    return 0;
}

我好菜啊,反演不会推,min25也不会写,一道题写了四个小时到处是bug,到处TLE和MLE。

std好快啊。qwq

2019杭电多校第一场 K. Function

题目想要求

\sum_{i=1}^n \gcd\left(\left\lfloor\sqrt[3]i\right\rfloor,i\right) mod~998244353

首先可以想到,利用 \lfloor\sqrt[3]{i}\rfloor 的不同值进行分块。我们先记 m=\lfloor\sqrt[3]{n}\rfloor, 那么原式变为

\sum_{c=1}^{m} \sum_{i=c^3}^{\min((c+1)^3-1, n)} \gcd(c,i)

看着那个min比较碍眼,不如考虑当 c 的情况,那么有

H(c)=\sum_{i=c^3}^{(c+1)^3-1} \gcd(c, i) = \sum_{d|c} d \sum_{i=c^3}^{(c+1)^3-1} \left[\gcd(c,i)=d\right]

看起来有个什么类似于前缀和的东西?再设

f(n,d) = \sum_{i=1}^{n} \left[ \gcd(d,i) = 1 \right]

根据莫比乌斯函数的容斥意义容易得到

f(n,d) = \sum_{g|d} \mu(g) \left\lfloor\frac{n}{g}\right\rfloor

回到上面那个式子

\begin{aligned} H(c) &= \sum_{d|c} d \times \left( f\left(\left\lfloor\frac{(c+1)^3-1}{d}\right\rfloor,\frac{c}{d}\right) – f\left(\left\lfloor\frac{c^3-1}{d}\right\rfloor,\frac{c}{d}\right) \right) \\ &= \sum_{d|c} d \sum_{g|\frac{c}{d}} \mu(g) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{gkd=c} d ~\mu(g) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{kt=c} \left( \sum_{gd=t} d ~\mu(g) \right) \left( \left\lfloor\frac{(c+1)^3-1}{gd}\right\rfloor – \left\lfloor\frac{c^3-1}{gd}\right\rfloor \right) \\ &= \sum_{kt=c} \varphi(t) \left( \left\lfloor\frac{(c+1)^3-1}{t}\right\rfloor – \left\lfloor\frac{c^3-1}{t}\right\rfloor \right) \end{aligned}

好像可以发现 c | (c+1)^3-1,那后面那个玩意又有点碍眼,看看好像只有 t=1 的时候不太一样?不如我们……

\begin{aligned} H(c) &= \sum_{kt=c} \varphi(t) \left( \frac{(c+1)^3-1}{t} – \left\lfloor\frac{c^3}{t} \right\rfloor + 1 \right) \\ &= \sum_{kt=c} \varphi(t)~k~(3c+3) + \sum_{kt=c} \varphi(t)k \\ &= c + 3(c+1) (\varphi * id)(c) \end{aligned}

那我们继续看后面那部分吧

\begin{aligned} H'(m) &= \sum_{d|m} d \sum_{i=m^3}^n [\gcd(m,i) = d] \\ &= \sum_{d|m} d \sum_{g|\frac{m}{d}} \mu(g) \left( \left\lfloor\frac{n}{gd}\right\rfloor – \left\lfloor\frac{m^3-1}{gd}\right\rfloor \right) \\ &= \sum_{t|m} \varphi(t) \left( \left\lfloor\frac{n}{t}\right\rfloor – \left\lfloor\frac{m^3-1}{t}\right\rfloor \right) \end{aligned}

再回去看看最开始要的东西

H'(m) + \sum_{i=1}^{m-1} H(i)

不过 H(i) 看起来可以线性筛再前缀和?

好像 O(n^{\frac{1}{6}}) 单组询问的复杂度可以接受?Let’s go~

#include <bits/stdc++.h>
typedef long long lld;
using namespace std;
const size_t MAXN = 1e7+5;
const int MOD = 998244353;
bool isnp[MAXN];
int pcn, pri[MAXN/10];
lld g[MAXN], s[MAXN];
int phi[MAXN], q[MAXN];
int mod;

void init_prime()
{
    g[1] = q[1] = phi[1] = 1;

    for (int i = 2; i < MAXN; i++)
    {
        if (!isnp[i])
        {
            g[i] = s[i] = 2*i-1;
            q[i] = i-1;
            phi[i] = i-1;
            pri[pcn++] = i;
        }

        for (int j = 0; j < pcn; j++)
        {
            int k = i * pri[j];
            if (k >= MAXN) break;
            isnp[k] = true;

            if (i % pri[j] == 0)
            {
                s[k] = (s[i] + q[i] + 0ll) * pri[j];
                q[k] = q[i] * pri[j];
                g[k] = g[i] / s[i] * s[k];
                phi[k] = phi[i] * pri[j];
                break;
            }
            else
            {
                phi[k] = phi[i] * (pri[j] - 1);
                q[k] = q[pri[j]];
                s[k] = s[pri[j]];
                g[k] = g[i] * s[k];
            }
        }
    }

    for (int i = 1; i < MAXN; i++)
    {
        g[i] = (g[i-1] + i + (3 * i + 3) * (g[i] % MOD)) % MOD;
    }
}

int two(__int128 n)
{
    if (n < 8) return 1;
    int l = 1, r = 1e9+7, m;
    while (r - l > 1)
    {
        m = (l + r) >> 1;
        if (m > n / m / m)
            r = m;
        else
            l = m;
    }
    return r-1;
}

void solve()
{
    __int128 n = 0;
    char buf[30];
    scanf("%s", buf);
    for (int i = 0; buf[i]; i++)
        n = n * 10 + (buf[i] - '0');
    int m = two(n);
    lld ans = g[m-1];
    __int128 qwq = m; qwq *= m; qwq *= m; qwq -= 1;
    for (int i = 1; i * i <= m; i++)
    {
        if (m % i == 0)
        {
            ans += phi[i] * ((n / i % MOD) - (qwq / i % MOD)) % MOD;
            if (m != i * i)
                ans += phi[m / i] * (n / (m / i) % MOD - (qwq / (m / i) % MOD)) % MOD;
        }
    }
    ans %= MOD;
    if (ans < 0) ans += MOD;
    printf("%lld\n", ans);
}

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