JZOJ 4583 求和

首先为了简洁 n=N,m=M。

然后 μ(i)μ(j)σ(ij) 这个式子，当 i 和 j 都没有平方质因子时， μ(i)μ(j)σ(ij)=μ(i)μ(j)σ(igcd(i,j))σ(jgcd(i,j))σ(gcd2(i,j))=μ(i)μ(j)σ(i)σ(j)σ(gcd2(i,j))σ2(gcd(i,j))

当 i 或 j 有平方质因子时，显然 μ(i)μ(j)=0，所以没有贡献。
于是 n∑i=1m∑j=1μ(i)μ(j)σ(i)σ(j)σ(gcd2(i,j))σ2(gcd(i,j))=min(n,m)∑d=1σ(gcd2(i,j))σ2(gcd(i,j))⋅∑id≤n∑jd≤m[gcd(i,j)=1]μ(id)μ(jd)σ(id)σ(jd)=min(n,m)∑d=1σ(gcd2(i,j))σ2(gcd(i,j))⋅∑id≤n∑jd≤m∑k|gcd(i,j)μ(k)μ(id)μ(jd)σ(id)σ(jd)=min(n,m)∑d=1σ(gcd2(i,j))σ2(gcd(i,j))⋅∑dk≤min(n,m)μ(k)∑idk≤n∑jdk≤mμ(idk)μ(jdk)σ(idk)σ(jdk)=min(n,m)∑T=1f(n,T)f(m,T)g(T)

其中 f(n,m)=∑im≤nμ(im)σ(im),g(n)=∑d|nσ(d2)σ2(d)μ(nd)。
求出所有的 f(n,T),f(m,T),g(T) 都是 O(nlogn) 的（认为 n,m 同阶），于是线性筛出 μ(i),σ(i),σ(i2) 这些东西就可以做了。
然鹅我处处求逆元常数过大（又没法递推）……于是要吸氧才能过……

代码：

#pragma GCC optimize (2)
#include <cstdio>
#include <algorithm>
using namespace std;
const int N = 5e5;
const long long mod = 998244353;
inline long long fpow(long long a,long long b)
{
    long long ret =  1;
    for(;b;b >>= 1)
        (b & 1) && (ret = ret * a % mod),a = a * a % mod;
    return ret;
}
int n,m;
int vis[N + 5],cnt,prime[N + 5];
long long sigma[N + 5],mn[N + 5],sqr[N + 5],sqrmn[N + 5],mu[N + 5];
long long f[2][N + 5],g[N + 5],t[N + 5];
long long ans;
int main()
{
    mu[1] = sigma[1] = mn[1] = sqr[1] = sqrmn[1] = 1;
    for(register int i = 2;i <= N;++i)
    {
        if(!vis[i])
            mu[prime[++cnt] = i] = mod - 1,sigma[i] = 1 + i,mn[i] = 1 + i,sqr[i] = (1 + i + (long long)i * i % mod) % mod,sqrmn[i] = (1 + i + (long long)i * i % mod) % mod;
        for(register int j = 1;j <= cnt && i * prime[j] <= N;++j)
        {
            vis[i * prime[j]] = 1;
            if(!(i % prime[j]))
            {
                sigma[i * prime[j]] = sigma[i] * fpow(mn[i],mod - 2) % mod * (mn[i] * prime[j] % mod + 1) % mod,mn[i * prime[j]] = (mn[i] * prime[j] % mod + 1) % mod;
                sqr[i * prime[j]] = sqr[i] * fpow(sqrmn[i],mod - 2) % mod * ((sqrmn[i] * prime[j] % mod + 1) % mod * prime[j] % mod + 1) % mod,sqrmn[i * prime[j]] = ((sqrmn[i] * prime[j] % mod + 1) % mod * prime[j] % mod + 1) % mod;
                break;
            }
            else
            {
                sigma[i * prime[j]] = sigma[i] * (prime[j] + 1) % mod,mn[i * prime[j]] = prime[j] + 1;
                sqr[i * prime[j]] = sqr[i] * (((long long)prime[j] * prime[j] + prime[j] + 1) % mod) % mod,sqrmn[i * prime[j]] = ((long long)prime[j] * prime[j] % mod + prime[j] + 1) % mod;
                mu[i * prime[j]] = mod - mu[i];
            }
        }
    }
    scanf("%d%d",&n,&m);
    if(n > m)
        swap(n,m);
    for(register int i = 1;i <= n;++i)
        for(register int j = i,k = 1;j <= n;j += i,++k)
            g[j] = (g[j] + sqr[i] * fpow(sigma[i] * sigma[i] % mod,mod - 2) % mod * mu[k] % mod) % mod;
    for(register int i = 1;i <= n;++i)
    {
        for(register int j = i;j <= n;j += i)
            f[0][i] = (f[0][i] + mu[j] * sigma[j] % mod) % mod;
        for(register int j = i;j <= m;j += i)
            f[1][i] = (f[1][i] + mu[j] * sigma[j] % mod) % mod;
        ans = (ans + f[0][i] * f[1][i] % mod * g[i] % mod) % mod;
    }
    printf("%lld\n",ans);
}