loj#6229. 这是一道简单的数学题 (??反演+杜教筛)

题意：给定\(n\le 10^9\)，求：\(F(n)=\sum_{i=1}^n\sum_{j=1}^i\frac{\mathrm{lcm}(i,j)}{\mathrm{gcd}(i,j)}\)，对1e9+7取模

推式子：

\(F(n)=\sum_{i=1}^n\sum_{j=1}^i\frac{\mathrm{lcm}(i,j)}{\mathrm{gcd}(i,j)}\)

\(=\sum_{i=1}^n\sum_{j=1}^i\frac{ij}{\gcd^2(i,j)}\)

\(=\sum_{p=1}^n\frac1{p^2}\sum_{i=1}^n\sum_{j=1}^iij[\gcd(i,j)=p]\)

\(=\sum_{p=1}^n\sum_{i=1}^{n/p}\sum_{j=1}^{i/p}ij[\gcd(i,j)=1]\)

\(=\sum_{p=1}^n\sum_{i=1}^{n/p}i\sum_{j=1}^{i/p}j[\gcd(i,j)=1]\)

根据一个经典式子：\(\sum_{i=1}^ni[\gcd(i,n)=1]=\frac {[n=1]+n\varphi(n)}{2}\)

\(=\sum_{p=1}^n\sum_{i=1}^{n/p}i\frac{[i=1]+i\varphi(i)}{2}\)

\(=\frac {n+\sum_{i=1}^n\lfloor\frac n i\rfloor i^2\varphi(i)}2\)

现在我们考虑求\(\sum_{i=1}^n\lfloor\frac n i\rfloor i^2\varphi(i)\)，暴力线性筛是\(O(n)\)的，显然会T，考虑用杜教筛优化

设\(f(i)=i^2\varphi(i),S(n)=\sum_{i=1}^nf(i)\)

原式=\(\sum_{i=1}^n\lfloor\frac n i\rfloor f(i)\)

由于\(\lfloor\frac ni\rfloor\)的取值只有\(O(\sqrt n)\)种，所以可以计算每一块比下一块多多少，然后计算下前缀和

这块可能不太好理解，例如n=11时

ans=11f(1)+5f(2)+3f(3)+2(f(4)+f(5))+1(f(6)+...+f(11))，替换成S则有

ans=(11-5)S(1)+(5-3)S(2)+(3-2)S(3)+(2-1)S(5)+S(11)

直接上杜教筛就行了，因为这些取值在计算S(11)时候都要用到，并且每个n/x只有一种，开个数组记忆化即可

代码里枚举的顺序不太一样，是从S大到S小一样，两块的差就相当于某一块的大小

#include <cstdio>
using namespace std;

const int p = 1000000007, fuck = 1000000;

bool vis[fuck + 233];
int prime[fuck], tot;
int phi[fuck + 233];

int qpow(int x, int y)
{
    int res = 1;
    while (y > 0)
    {
        if (y & 1) res = res * (long long)x % p;
        x = x * (long long)x % p;
        y >>= 1;
    }
    return res;
}

int inv2 = qpow(2, p - 2), inv6 = qpow(6, p - 2);

int s1(int x) { x %= p; return x* (long long)(x + 1) % p * inv2 % p; }
int s2(int x) { x %= p; return x * (long long)(x + 1) % p * (x * 2 + 1) % p * inv6 % p; }
int s3(int x) { return qpow(s1(x), 2); }

bool count[4000010];
int ans[400010], n;

int chuans(int x)
{
    if (x <= fuck) { return phi[x]; }
    if (count[n / x]) return ans[n / x];
    count[n / x] = 1;
    int res = s3(x);
    for (int i = 2, j; i <= x; i = j + 1)
    {
        j = x / (x / i);
        res = ((res - (s2(j) - s2(i - 1)) * (long long)chuans(x / i) % p) % p + p) % p;
    }
    return ans[n / x] = res;
}

int main()
{
    phi[1] = 1;
    for (int i = 2; i <= fuck; i++)
    {
        if (vis[i] == false) prime[++tot] = i, phi[i] = i - 1;
        for (int j = 1; j <= tot && i * prime[j] <= fuck; j++)
        {
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0) { phi[i * prime[j]] = phi[i] * prime[j]; break; }
            else phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
        phi[i] = (phi[i - 1] + phi[i] * (long long)i % p * i % p) % p;
    }
    scanf("%d", &n);
    int sum = 0;
    for (int i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        sum = (sum + (j - i + 1) * (long long)chuans(n / i) % p) % p;
    }
    sum = (sum + n) % p * (long long)inv2 % p;
    printf("%d\n", sum);
    return 0;
}

Min_25筛表示他需要重复计算根N次，就GG了

原文地址：https://www.cnblogs.com/oier/p/10294169.html