luogu P4240 毒瘤之神的考验（莫比乌斯反演+递推前缀和+数论分块）

problem

洛谷链接

solution

$$\sum _{i=1}^n\sum _{j=1}^m\phi (ij)=\sum _{i=1}^n\sum _{j=1}^m\frac{\phi (i)\phi (j)\gcd (i,j)}{\phi (\gcd (i,j))}$$$$=\sum _{d=1}^{\min (n,m)}\frac{d}{\phi (d)}\sum _{i=1}^n\sum _{j=1}^m\phi (i)\phi (j)[(i,j)=d]$$$$=\sum _{d=1}^{\min (n,m)}\frac{d}{\phi (d)}\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\phi (id)\phi (jd)[(i,j)=1]$$$$=\sum _{d=1}^{\min (n,m)}\frac{d}{\phi (d)}\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\phi (id)\phi (jd)\sum _{k|(i,j)}\mu (k)$$$$=\sum _{d=1}^{\min (n,m)}\frac{d}{\phi (d)}\sum _{k=1}^{\min (\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )}\mu (k)\sum _{i=1}^{\lfloor \frac{n}{dk}\rfloor }\phi (ikd)\sum _{j=1}^{\lfloor \frac{m}{dk}\rfloor }\phi (jkd)$$$$\Vert t=dk\iff k=\frac{t}{d}$$$$=\sum _{t=1}^{\min (n,m)}\sum _{d|t}\frac{d}{\phi (d)}\mu (\frac{t}{d})+\sum _{i=1}^{\lfloor \frac{n}{t}\rfloor }\phi (it)\sum _{j=1}^{\lfloor \frac{m}{t}\rfloor }\phi (jt)$$

令 $f(t)={\sum }_{d|t}\frac{d}{\phi (d)}\mu (\frac{t}{d})$，在 $O(n\log n)$ 的时间内预处理。

设 $g(a,b)={\sum }_{i=1}^b\phi (i\cdot a)$。
$$\sum _{i=1}^n\sum _{j=1}^m\phi (ij)\iff \sum _{i=1}^{\min (n,m)}f(i)\cdot g(i,\lfloor \frac{n}{i}\rfloor )\cdot g(i,\lfloor \frac{m}{i}\rfloor )$$注意到有用的 $g(a,b)$ 均满足 $ab\le \max (n,m)$，可以用动态二维数组预处理。

根据定义容易得递推式 $:g(a,b)=g(a,b-1)+\phi (ab)$。

设前缀和 $s(x,y,n)={\sum }_{i=1}^{n}f(i)g(i,x)g(i,y)$，则有递推式 $s(x,y,n)=s(x,y,n-1)+f(n)g(n,x)g(n,y)$。

设定一个阈值 $B$，预处理所有 $x\le B\wedge y\le B$ 的 $s(x,y,z)$，此时还满足 $z\le \frac{n}{\max (x,y)}$，同理用动态三维数组存。

询问时枚举 $i$。
对于 $\lfloor \frac{n}{i}\rfloor \le B\wedge \lfloor \frac{m}{i}\rfloor \le B$ 直接暴力求 $f(i)g(i,\lfloor \frac{n}{i}\rfloor )g(i,\lfloor \frac{m}{i}\rfloor )$（反正都是预处理过的信息）。

其余部分利用预处理的前缀和 $s(x,y,z)$ 直接整除分块计算。

假设当 $i\in [l,r]$ 时，$\lfloor \frac{n}{i}\rfloor ,\lfloor \frac{m}{i}\rfloor$ 均相同，即 $f(i)g(i,\lfloor \frac{n}{i}\rfloor )g(i,\lfloor \frac{m}{i}\rfloor )$ 相同，

则 ${\sum }_{i=l}^{r}f(i)g(i,x)g(i,y)=s(x,y,r)-s(x,y,l-1)$。

查询时间复杂度 $O(T\frac{n}{B}+T\sqrt{n})$，预处理时间 $O(n\log n+nB)$，空间复杂度 $O(nB)$。

当 $B=n^{\frac{1}{3}}$ 时复杂度较优秀。

code

#include <bits/stdc++.h>
using namespace std;
#define int long long
#define mod 998244353
#define maxn 100005
#define maxB 60
int *g[maxn], *s[maxB][maxB];
int mu[maxn], phi[maxn], prime[maxn], f[maxn];
bool vis[maxn];
int cnt;

int qkpow( int x, int y ) {
    int ans = 1;
    while( y ) {
        if( y & 1 ) ans = ans * x % mod;
        x = x * x % mod;
        y >>= 1;
    }
    return ans;
}

void sieve( int n = 1e5, int B = 50 ) {
    mu[1] = phi[1] = 1;
    for( int i = 2;i <= n;i ++ ) {
        if( ! vis[i] ) {
            prime[++ cnt] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for( int j = 1;j <= cnt and i * prime[j] <= n;j ++ ) {
            vis[i * prime[j]] = 1;
            if( i % prime[j] == 0 ) {
                mu[i * prime[j]] = 0;
                phi[i * prime[j]] = phi[i] * prime[j] % mod;
                break;
            }
            else {
                mu[i * prime[j]] = -mu[i];
                phi[i * prime[j]] = phi[i] * phi[prime[j]] % mod;
            }
        }
    }
    for( int i = 1;i <= n;i ++ )
        for( int j = i;j <= n;j += i )
            ( f[j] += i * qkpow( phi[i], mod - 2 ) % mod * mu[j / i] ) %= mod;
    for( int i = 1;i <= n;i ++ ) {
        g[i] = new int [n / i + 5];
        g[i][0] = 0;
        for( int j = 1;i * j <= n;j ++ )
            g[i][j] = ( g[i][j - 1] + phi[i * j] ) % mod;
    }
    for( int i = 1;i <= B;i ++ )
        for( int j = 1;j <= B;j ++ ) {
            s[i][j] = new int [n / max( i, j ) + 5];
            s[i][j][0] = 0;
            for( int k = 1;k * max( i, j ) <= n;k ++ )
                s[i][j][k] = ( s[i][j][k - 1] + f[k] * g[k][i] % mod * g[k][j] ) % mod;
        }
}

signed main() {
    sieve();
    int T, n, m, B = 50;
    scanf( "%lld", &T );
    while( T -- ) {
        scanf( "%lld %lld", &n, &m );
        if( n > m ) swap( n, m );
        int ans = 0;
        for( int i = 1;i <= m / B;i ++ )
            ( ans += f[i] * g[i][n / i] % mod * g[i][m / i] ) %= mod;
        for( int l = m / B + 1, r;l <= n;l = r + 1 ) {
            r = min( n / ( n / l ), m / ( m / l ) );
            ( ans += s[n / l][m / l][r] - s[n / l][m / l][l - 1] ) %= mod;
        }
        printf( "%lld\n", ( ans + mod ) % mod );
    }
    return 0;
}