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

Jul 30, 2019 由 小羊

题目想要对于给定 $n$ 和 $k$ 求

$\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