LibreOJ 528 「LibreOJ β Round #4」求和

发表于 | 分类于 |

\(N,M \le 10^{13}\)，看起来 \(O(N^{2/3})\) 不可做（视 \(N,M\) 同阶）。
当我知道怎么做的时候，wdnmd……

\[ \begin{align*} \sum\limits_{i=1}^N\sum\limits_{j=1}^M \mu^2(\gcd(i,j)) &=\sum\limits_{d=1}^{\min(N,M)} \mu^2(d) \sum\limits_{i=1}^{\lfloor\frac Nd\rfloor}\sum\limits_{j=1}^{\lfloor\frac Md\rfloor}[\gcd(i,j)=1] \\ &=\sum\limits_{d=1}^{\min(N,M)} \mu^2(d) \sum\limits_{k=1}^{\min(\lfloor\frac Nd\rfloor,\lfloor\frac Md\rfloor)} \mu(k)\lfloor\frac N{dk}\rfloor\lfloor\frac M{dk}\rfloor \\ &=\sum\limits_{T=1}^{\min(N,M)} \lfloor\frac NT\rfloor\lfloor\frac MT\rfloor \sum\limits_{d|T} \mu(d)\mu^2(\frac Td) \end{align*} \]

考虑函数 \(\mu * \mu^2\) 的前缀和 \(\sum\limits_{i=1}^n \sum\limits_{d|i} \mu(d)\mu^2(\frac id)\)。
设 \(f(n) = \max\limits_{d^2|n} d\)，则 \[ \begin{align*} \sum\limits_{i=1}^n\sum\limits_{d|i}\mu(d)\mu^2(\frac id) &=\sum\limits_{i=1}^n\sum\limits_{d|i}\mu^2(d)\mu(\frac id) \\ &=\sum\limits_{i=1}^n\sum\limits_{d|i}[f(d)=1]\mu(\frac id) \\ &=\sum\limits_{i=1}^n\sum\limits_{d|i}\mu(\frac id)\sum\limits_{k^2|d}\mu(k) \\ &=\sum\limits_{i=1}^n\sum\limits_{k^2|i}\mu(k)\sum\limits_{d|\frac i{k^2}} \mu(\frac i{dk^2}) \\ &=\sum\limits_{i=1}^n\sum\limits_{k^2|i}[i=k^2]\mu(k) \\ &=\sum\limits_{k=1}^{\lfloor\sqrt n\rfloor} \mu(k) \end{align*} \]

……
然后就变成一道大水题了……

代码： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
const long long N = 1e13;
const int MX = 1e7;
const int mod = 998244353;
long long n,m;
int vis[MX + 5],cnt,prime[MX + 5],mu[MX + 5];
int ans;
int main()
{
    mu[1] = 1;
    for(register int i = 2;i <= MX;++i)
    {
        if(!vis[i])
            mu[prime[++cnt] = i] = mod - 1;
        for(register int j = 1;j <= cnt && i * prime[j] <= MX;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
                break;
            mu[i * prime[j]] = mod - mu[i];
        }
        mu[i] = (mu[i] + mu[i - 1]) % mod;
    }
    scanf("%lld%lld",&n,&m);
    for(register long long l = 1,r;l <= min(n,m);l = r + 1)
    {
        r = min(n / (n / l),m / (m / l));
        ans = (ans + (n / l % mod) * (m / l % mod) % mod * (mu[(int)sqrt(r)] - mu[(int)sqrt(l - 1)] + mod) % mod) % mod;
    }
    printf("%d\n",ans);
}