BZOJ3309: DZY Loves Math

传送门

神题！（似乎对我这样的低智选手来说那道题都是神题。

你的莫比乌斯怎么这么熟练啊：

$ans=\sum\limits_{i=1}^{N} f(i) \times g(i)$

$f(d)=\sum\limits_{i=1}^{N}\sum\limits_{j=1}^{M} gcd(i,j)==d$

$f(d)=\sum\limits_{d|k} \frac{N}{k} \frac{M}{k} \mu (\frac{k}{d})$

$ans=\sum\limits_{i=1}^{N} \sum\limits_{i|k} \frac{N}{k} \frac{M}{k} \mu (\frac{k}{i}) \times g(i)$

观察到$i$与$N,M$关系不大，提到前面

$ans=\sum\limits_{k=1}^{N}\frac{N}{k}\frac{M}{k}\sum\limits_{i|k}g(i)\mu(\frac{k}{i})$

现在考虑怎么在快速的求出后面那一段。

剩下的这一段比较玄学，我胡说一下..

我们只考虑$\mu(\frac{k}{i}) \not= 0$的那一段。

先放结论，如果$k$的所有质因子的质数全部相等的话，设$i$表示不同的质因子数，那么结果就是$(-1)^{i+1}$，否则就是$0$。

我们从整体上来分析这个式子，对于这个式子的$g(i)$讲，我们不妨设$k$的质因数的指数的最大值为$p$，那么$g(i)$的取值只有可能为$p$或$p-1$，否则$\mu(\frac{k}{i})$的值为$0$。

那么如果$k$的各个质因子的指数不相等的话，一定能拆分成两堆相等的方案，使其和为$0$。

如果相等的话，那么会存在一个单独的$a-1$与一个$-a$互相抵消为$-1$,因为质因数的数量影响着那个数的正负，所以结果就和$i$相关。

那么这个东西可以考虑用线筛搞出来。

我们在线筛的时候记录两个值，一个是$low(i)$代表$i$的最小的质因子的指数，$pow(i)$代表$i$的最小的质因子的质数幂。

具体参考代码吧，自己推推能明白。

//BZOJ 3309
//by Cydiater
//2017.2.21
#include <iostream>
#include <queue>
#include <map>
#include <cstring>
#include <string>
#include <algorithm>
#include <cstdlib>
#include <cstdio>
#include <iomanip>
#include <ctime>
#include <cmath>
#include <bitset>
#include <set>
#include <vector>
#include <complex>
using namespace std;
#define ll long long
#define up(i,j,n)	for(int i=j;i<=n;i++)
#define down(i,j,n)	for(int i=j;i>=n;i--)
#define cmax(a,b)	a=max(a,b)
#define cmin(a,b)	a=min(a,b)
const ll MAXN=1e7+5;
const ll LIM=1e7;
inline ll read(){
char ch=getchar();ll x=0,f=1;
while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
while(ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
return x*f;
}
int g[MAXN],prime[MAXN],cnt=0,lw[MAXN],pw[MAXN],T,N,M;
bool vis[MAXN];
namespace solution{
void Prepare(){
g[1]=0;
up(i,2,LIM){
if(!vis[i]){
g[i]=1;lw[i]=1;pw[i]=i;
prime[++cnt]=i;
}
up(j,1,cnt){
ll ip=i*prime[j];
if(ip>LIM)break;
vis[ip]=1;
if(i%prime[j]){
lw[ip]=1;pw[ip]=prime[j];
g[ip]=(lw[i]==1?-g[i]:0);
}else{
ll tmp=i/pw[i];
lw[ip]=lw[i]+1;pw[ip]=pw[i]*prime[j];
if(tmp==1)	g[ip]=1;
else 		g[ip]=(lw[tmp]==lw[ip]?-g[tmp]:0);
break;
}
}
}
up(i,2,LIM)g[i]+=g[i-1];
//cout<<"Time has passed:"<<1.0*clock()/1000<<"s!"<<endl;
}
void Solve(){
T=read();
while(T--){
N=read();M=read();
if(N>M)swap(N,M);
ll j,ans=0;
for(ll i=1;i<=N;i=j+1){
j=min(N/(N/i),M/(M/i));
ans+=(ll)(g[j]-g[i-1])*(ll)(N/i)*(ll)(M/i);
}
printf("%lld\n",ans);
}
//cout<<"Time has passed:"<<1.0*clock()/1000<<"s!"<<endl;
}
}
int main(){
//freopen("input.in","r",stdin);
using namespace solution;
Prepare();
Solve();
return 0;
}