LibreOJ 528 「LibreOJ β Round #4」求和

N,M≤1013，看起来 O(N2/3) 不可做（视 N,M 同阶）。
当我知道怎么做的时候，wdnmd……

N∑i=1M∑j=1μ2(gcd(i,j))=min(N,M)∑d=1μ2(d)⌊Nd⌋∑i=1⌊Md⌋∑j=1[gcd(i,j)=1]=min(N,M)∑d=1μ2(d)min(⌊Nd⌋,⌊Md⌋)∑k=1μ(k)⌊Ndk⌋⌊Mdk⌋=min(N,M)∑T=1⌊NT⌋⌊MT⌋∑d|Tμ(d)μ2(Td)

考虑函数 μ∗μ2 的前缀和 n∑i=1∑d|iμ(d)μ2(id)。
设 f(n)=maxd2|nd，则 n∑i=1∑d|iμ(d)μ2(id)=n∑i=1∑d|iμ2(d)μ(id)=n∑i=1∑d|i[f(d)=1]μ(id)=n∑i=1∑d|iμ(id)∑k2|dμ(k)=n∑i=1∑k2|iμ(k)∑d|ik2μ(idk2)=n∑i=1∑k2|i[i=k2]μ(k)=⌊√n⌋∑k=1μ(k)

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

代码：

#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);
}