3309: DZY Loves Math
Time Limit: 20 Sec Memory Limit: 512 MB
Submit: 761 Solved: 401
[Submit][Status][Discuss]
Description
对于正整数n，定义f(n)为n所含质因子的最大幂指数。例如f(1960)=f(2^3 * 5^1 * 7^2)=3, f(10007)=1, f(1)=0。
给定正整数a,b，求sigma(sigma(f(gcd(i,j)))) (i=1..a, j=1..b)。
Input
第一行一个数T，表示询问数。
接下来T行，每行两个数a,b，表示一个询问。
Output
对于每一个询问，输出一行一个非负整数作为回答。
Sample Input
10000
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957
Sample Output
14225956593420
4332838845846
15400094813
HINT
【数据规模】
T<=10000
1<=a,b<=10^7
Source
分析：
还是莫比乌斯反演...然而公式还是很好推...都是套路...但是我不会求f(i)TAT...
Σ(1<=i<=n) Σ(1<=j<=m) f(gcd(i,j)) 依旧令n<=m
=Σ(1<=i<=n) Σ(i|j&&j<=n) μ(j/i)*(n/j)*(m/j)
=Σ(1<=j<=n) (n/j)*(m/j) Σ(i|j) f(i)*μ(j/i)
看ij不顺眼...还是换成x和d吧...
=Σ(1<=x<=n) (n/x)*(m/x) Σ(d|x) f(d)*μ(x/d)
然后怎么办，用树状数组维护前缀和？？感觉还是会T的很惨...直觉告诉我这题复杂度绝对是T*(sqrt(n)+sqrt(m))的......
再仔细分析一下：
我们发现当x/d存在平方因子也就是f(x/d)>1的时候μ(x/d)=0...
但是x和d是分不开的...质因数分解一下？？
x=p1^a1*p2^a2*p3^a3*......*pt^at
d=p1^b1*p2^b2*p3^b3*......*pt^bt
然后呢？额，还是借鉴了一下别人的代码TAT...
因为μ函数改变了f的正负，所以有些东西可以相消，所以我们可以发现对于x的a数组可以分为以下两种情况：
No.1
不存在ai≠aj，那么当前x对ans做出的贡献就是(-1)^(t+1)
证明如下：
当f有贡献的时候当且仅当μ不为0，也就是说如果设ai=aj=a，那么当且仅当f(d)=a||a-1的时候，首先考虑f(d)=a-1的时候：
只有一种情况就是所有的质因子次数都去a-1的时候，此时ans+=(a-1)*(-1)^t...
当f(d)=a的时候，那么对ans的贡献就是a*((-1)^0*C(t,0)+(-1)^1*C(t,1)+......+(-1)^(t-1)*C(t,t-1))
根据杨辉三角我们可以得到这个式子的结果是a*(-1)^(t+1)...
那么加起来ans=(-1)^(t+1)
No.2
存在ai≠aj的时候，根据组合数组合一下就可以发现加起来成了0(懒得证了...自己推一推吧...)
也就是说我们只需要计算ai=aj=a的时候的(-1)^(t+1)，这个问题就可以在线性筛的时候解决...
代码：
#include<algorithm>
#include<iostream>
#include<cstring>
#include<cstdio>
//by NeighThorn
using namespace std;
//大鹏一日同风起，扶摇直上九万里 const int maxn=+; int a,b,cas,cnt,f[maxn],num[maxn],vis[maxn],lala[maxn],prime[maxn]; long long ans; signed main(void){
cnt=;
for(int i=;i<=;i++){
if(!vis[i])
prime[++cnt]=i,lala[i]=num[i]=,f[i]=;
for(int j=;j<=cnt&&i*prime[j]<=;j++){
int x=i*prime[j];vis[x]=;
if(i%prime[j]==){
lala[x]=lala[i];num[x]=num[i]+;
if(lala[i]==)
f[x]=;
else
f[x]=(num[lala[x]]==num[x]?-f[lala[x]]:);
break;
}
lala[x]=i,num[x]=,f[x]=(num[i]==?-f[i]:);
}
}
for(int i=;i<=;i++)
f[i]+=f[i-];
scanf("%d",&cas);
while(cas--){
scanf("%d%d",&a,&b);
if(a>b)
swap(a,b);
ans=;
for(int i=,r;i<=a;i=r+){
r=min(a/(a/i),b/(b/i));
ans+=(long long)(f[r]-f[i-])*(a/i)*(b/i);
}
printf("%lld\n",ans);
}
return ;
}//Cap ou pas cap. Cap.
By NeighThorn