DZY Loves Math
万年没写莫反都快不会了。。认真推一次式子吧 /kel
首先把题面写下来：
[sum_{i=1}^nsum_{j=1}^m f(gcd(i,j))
]
其中 $ f $ 为素因数次幂最高的数的指数。
首先，按照莫反的套路枚举 $ gcd $ 的值
[sum_{i=1}^n sum_{j=1}^m sum_{d=1}^n [gcd(i,j) = d] f(d)
]
然后交换循环次序
[sum_{d=1}^n f(d) sum_{i=1}^n sum_{j=1}^m [gcd(i,j) = d]
]
接着把 $ gcd = d $ 变成 $ = 1 $ 来莫反
[sum_{d=1}^n f(d) sum_{i=1}^{lfloorfrac{n}{d}
floor} sum_{j=1}^{lfloor frac m d
floor } [gcd(i,j) = 1]
]
于是把 $ epsilon = I imes mu $ 带进去
[sum_{d=1}^n f(d)sum_{i=1}^{lfloorfrac{n}{d}
floor} sum_{j=1}^{lfloor frac m d
floor } sum_{k | gcd(i,j)} mu(k)
]
下面把 $ sum k $ 提前，于是有
[sum_{d=1} ^ n f(d) sum_{k=1}^{lfloor frac n d
floor} lfloor frac n {dk}
floor lfloor frac n {dk}
floor mu(k)
]
然后就是莫反的经典套路，换 $ T $ 法。设 $ T = dk $ ，于是问题变成了枚举 $ d,k $
[sum_{T=1}^n lfloor frac n T
floor lfloor frac m T
floor sum_{d|T} f(d) mu(frac T k)
]
其实后面那一堆显然是一个狄利克雷卷积，写成卷积形式就是
[sum_{T=1}^{n} lfloor frac n T
floor lfloor frac m T
floor (f imes mu)(T)
]
好，如果我们求出了 $ f imes mu $ 和它的前缀和就可以数论分块做了。对于 $ f imes mu $ 的做法，可以线筛（但是略有麻烦），于是我们也可以用前几天学的 这个 的做法来 $ O(nloglog n) $ 做。
它并不是一个标准的狄利克雷卷积？但是它卷了 $ mu $ ，不妨设 $ g = f imes mu $ 那么有 $ f = g imes I $。
把 $ f = g imes I $ 写开就是 $ f(x) = sum_{d | x} g(d) $ 这就是标准的逆卷积了，可以套用那里面的第三种结论。
于是我们就有了一种 $ O(nloglog n + Tsqrt n) $ 的做法，可以通过这个题。
#include "iostream"
#include "algorithm"
#include "cstring"
#include "cstdio"
#include "vector"
#include "map"
using namespace std;
#define MAXN 10000006
typedef long long ll;
int n , m;
int pri[MAXN] , en , mm[MAXN]; // mx : 最大次数 mm : 最小质数的次数
long long mx[MAXN];
void sieve( ) {
for( int i = 2 ; i < MAXN ; ++ i ) {
if( !pri[i] ) pri[++ en] = i , mx[i] = mm[i] = 1;
for( int j = 1 ; j <= en && pri[j] * i < MAXN ; ++ j ) {
pri[i * pri[j]] = 1;
if( i % pri[j] == 0 ) {
mm[i * pri[j]] = mm[i] + 1 , mx[i * pri[j]] = max( mx[i] , mm[i] + 1ll );
break;
}
mm[i * pri[j]] = 1 , mx[i * pri[j]] = max( mx[i] , 1ll );
}
}
for( int i = en ; i ; -- i )
for( int j = ( MAXN - 1 ) / pri[i] ; j ; -- j ) {
mx[j * pri[i]] -= mx[j];
}
for( int i = 1 ; i < MAXN ; ++ i ) mx[i] += mx[i - 1];
}
long long ans;
int main() {
sieve();
int T;cin >> T;
while( T-- ) {
scanf("%d%d",&n,&m); ans = 0;
int x = min( n , m );
for( int l = 1 , r ; l <= x ; l = r + 1) {
r = min( n / ( n / l ) , m / ( m / l ) );
ans += 1ll * ( mx[r] - mx[l - 1] ) * ( n / l ) * ( m / l );
}
printf("%lld
",ans);
}
}