problem
洛谷链接
solution
∑i=1n∑j=1mφ(ij)=∑i=1n∑j=1mφ(i)φ(j)gcd(i,j)φ(gcd(i,j))\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)=\sum_{i=1}^n\sum_{j=1}^m\frac{\varphi(i)\varphi(j)\gcd(i,j)}{\varphi{(\gcd(i,j))}}i=1∑nj=1∑mφ(ij)=i=1∑nj=1∑mφ(gcd(i,j))φ(i)φ(j)gcd(i,j)=∑d=1min(n,m)dφ(d)∑i=1n∑j=1mφ(i)φ(j)[(i,j)=d]=\sum_{d=1}^{\min(n,m)}\frac{d}{\varphi(d)}\sum_{i=1}^n\sum_{j=1}^m\varphi(i)\varphi(j)\big[(i,j)=d\big]=d=1∑min(n,m)φ(d)di=1∑nj=1∑mφ(i)φ(j)[(i,j)=d]=∑d=1min(n,m)dφ(d)∑i=1⌊nd⌋∑j=1⌊md⌋φ(id)φ(jd)[(i,j)=1]=\sum_{d=1}^{\min(n,m)}\frac{d}{\varphi(d)}\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}\varphi(id)\varphi(jd)\big[(i,j)=1\big]=d=1∑min(n,m)φ(d)di=1∑⌊dn⌋j=1∑⌊dm⌋φ(id)φ(jd)[(i,j)=1]=∑d=1min(n,m)dφ(d)∑i=1⌊nd⌋∑j=1⌊md⌋φ(id)φ(jd)∑k∣(i,j)μ(k)=\sum_{d=1}^{\min(n,m)}\frac{d}{\varphi(d)}\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac m d\rfloor}\varphi(id)\varphi(jd)\sum_{k|(i,j)}\mu(k)=d=1∑min(n,m)φ(d)di=1∑⌊dn⌋j=1∑⌊dm⌋φ(id)φ(jd)k∣(i,j)∑μ(k)=∑d=1min(n,m)dφ(d)∑k=1min(⌊nd⌋,⌊md⌋)μ(k)∑i=1⌊ndk⌋φ(ikd)∑j=1⌊mdk⌋φ(jkd)=\sum_{d=1}^{\min(n,m)}\frac{d}{\varphi(d)}\sum_{k=1}^{\min(\lfloor\frac nd\rfloor,\lfloor\frac m d\rfloor)}\mu(k)\sum_{i=1}^{\lfloor\frac{n}{dk}\rfloor}\varphi(ikd)\sum_{j=1}^{\lfloor\frac{m}{dk}\rfloor}\varphi(jkd)=d=1∑min(n,m)φ(d)dk=1∑min(⌊dn⌋,⌊dm⌋)μ(k)i=1∑⌊dkn⌋φ(ikd)j=1∑⌊dkm⌋φ(jkd)∥t=dk⇔k=td\Bigg\lVert{t=dk\Leftrightarrow k=\frac{t}{d}}∥∥∥∥∥t=dk⇔k=dt=∑t=1min(n,m)∑d∣tdφ(d)μ(td)+∑i=1⌊nt⌋φ(it)∑j=1⌊mt⌋φ(jt)=\sum_{t=1}^{\min(n,m)}\sum_{d|t}\frac{d}{\varphi(d)}\mu(\frac td)+\sum_{i=1}^{\lfloor\frac nt\rfloor}\varphi(it)\sum_{j=1}^{\lfloor\frac mt\rfloor}\varphi(jt) =t=1∑min(n,m)d∣t∑φ(d)dμ(dt)+i=1∑⌊tn⌋φ(it)j=1∑⌊tm⌋φ(jt)
令 f(t)=∑d∣tdφ(d)μ(td)f(t)=\sum_{d|t}\frac{d}{\varphi(d)}\mu(\frac td)f(t)=∑d∣tφ(d)dμ(dt),在 O(nlogn)O(n\log n)O(nlogn) 的时间内预处理。
设 g(a,b)=∑i=1bφ(i⋅a)g(a,b)=\sum_{i=1}^b\varphi(i·a)g(a,b)=∑i=1bφ(i⋅a)。
∑i=1n∑j=1mφ(ij)⇔∑i=1min(n,m)f(i)⋅g(i,⌊ni⌋)⋅g(i,⌊mi⌋)\sum_{i=1}^n\sum_{j=1}^m\varphi(ij)\Leftrightarrow \sum_{i=1}^{\min(n,m)}f(i)·g(i,\lfloor\frac ni\rfloor)·g(i,\lfloor\frac mi\rfloor) i=1∑nj=1∑mφ(ij)⇔i=1∑min(n,m)f(i)⋅g(i,⌊in⌋)⋅g(i,⌊im⌋)注意到有用的 g(a,b)g(a,b)g(a,b) 均满足 ab≤max(n,m)ab\le \max(n,m)ab≤max(n,m),可以用动态二维数组预处理。
根据定义容易得递推式 :g(a,b)=g(a,b−1)+φ(ab):g(a,b)=g(a,b-1)+\varphi(ab):g(a,b)=g(a,b−1)+φ(ab)。
设前缀和 s(x,y,n)=∑i=1nf(i)g(i,x)g(i,y)s(x,y,n)=\sum_{i=1}^nf(i)g(i,x)g(i,y)s(x,y,n)=∑i=1nf(i)g(i,x)g(i,y),则有递推式 s(x,y,n)=s(x,y,n−1)+f(n)g(n,x)g(n,y)s(x,y,n)=s(x,y,n-1)+f(n)g(n,x)g(n,y)s(x,y,n)=s(x,y,n−1)+f(n)g(n,x)g(n,y)。
设定一个阈值 BBB,预处理所有 x≤B∧y≤Bx\le B\wedge y\le Bx≤B∧y≤B 的 s(x,y,z)s(x,y,z)s(x,y,z),此时还满足 z≤nmax(x,y)z\le \frac{n}{\max(x,y)}z≤max(x,y)n,同理用动态三维数组存。
询问时枚举 iii。
对于 ⌊ni⌋≤B∧⌊mi⌋≤B\lfloor\frac n i\rfloor\le B\wedge\lfloor\frac m i\rfloor\le B⌊in⌋≤B∧⌊im⌋≤B 直接暴力求 f(i)g(i,⌊ni⌋)g(i,⌊mi⌋)f(i)g(i,\lfloor\frac ni\rfloor)g(i,\lfloor\frac m i\rfloor)f(i)g(i,⌊in⌋)g(i,⌊im⌋)(反正都是预处理过的信息)。
其余部分利用预处理的前缀和 s(x,y,z)s(x,y,z)s(x,y,z) 直接整除分块计算。
假设当 i∈[l,r]i\in[l,r]i∈[l,r] 时,⌊ni⌋,⌊mi⌋\lfloor\frac n i\rfloor,\lfloor\frac m i\rfloor⌊in⌋,⌊im⌋ 均相同,即 f(i)g(i,⌊ni⌋)g(i,⌊mi⌋)f(i)g(i,\lfloor\frac n i\rfloor)g(i,\lfloor\frac m i\rfloor)f(i)g(i,⌊in⌋)g(i,⌊im⌋) 相同,
则 ∑i=lrf(i)g(i,x)g(i,y)=s(x,y,r)−s(x,y,l−1)\sum_{i=l}^rf(i)g(i,x)g(i,y)=s(x,y,r)-s(x,y,l-1)∑i=lrf(i)g(i,x)g(i,y)=s(x,y,r)−s(x,y,l−1)。
查询时间复杂度 O(TnB+Tn)O(T\frac nB+T\sqrt{n})O(TBn+Tn),预处理时间 O(nlogn+nB)O(n\log n+nB)O(nlogn+nB),空间复杂度 O(nB)O(nB)O(nB)。
当 B=n13B=n^{\frac 1 3}B=n31 时复杂度较优秀。
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;
}