3309: DZY Loves Math
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
Sample Output
HINT
【数据规模】
T<=100001<=a,b<=10^7这道题我足足写了一白板。因为这实在是一道好题。
我们要求∑∑f(gcd(i,j))，可以枚举d=gcd(i,j)
则有∑f(d)∑∑[gcd(i,j)==d]=∑f(d)∑∑e(gcd(i,j))=∑f(d)∑∑∑[k|i][k|j]mu(k)
即为∑f(d)∑mu(k)(n/dk)(m/dk)
令T=dk，则d|T,k=T/d，故有∑(n/T)(m/T)∑f(d)mu(T/d)。设g(T)=∑f(d)mu(T/d)。
如果我们能够O（n）算出g(T)及其前缀和，那么按照套路便能直接用分块。O（sqrt（n））单次询问。
现在看如何预处理出g(T)=∑f(d)mu(T/d)：先唯一分解出T和d，再分析mu(T/d)与f(d)的性质。然后找最大项数目，分类讨论，再进一步分类讨论即得解。这个分类讨论当时证的很复杂。有篇博客所言甚善：
如果存在ai≠aj(i≠j)，那答案只取决于最大的幂，对于其余部分通过组合数的性质，取得的数字个数为奇数和偶数的数量相等，因此和都是0，故如果存在ai≠aj(i≠j)，则g(T)=0
如果所有的a值都相等，我们假设对于任意选取方案，f值都不变，那么由于选取奇数个元素和偶数个元素的方案数相等，和仍然为0
但是有一种选取方案的f值=a-1 即d=p1*p2*p3*...pk时，因此要把那个减去，根据莫比乌斯函数性质，减去的是(-1)^k，因为是减去所以最终结果为(-1)^(k+1)，故如果不存在ai≠aj，则g(T)=(-1)^(k+1)
然后，可以使用欧拉筛O（n）预处理。但是怎样处理呢？考虑到欧拉筛中i*prime[j]的最小质因子是prime[j]。于是我们可以直接使用i的最小质因子与j比较。所以对于每一个数，需要存g[i]和i的最小质因子的幂指数a[i]。
然后，当i%prime[j]==0，我们怎么做呢？常见的做法是把i中所有的prime[j]都用除法除去。但事实上，这样效率太慢。考虑到欧拉筛中i*prime[j]的最小质因子是prime[j]（可用反证法证明），那么i的最小质因子也是prime[j]。如果我们存下minpa[i]表示i最小质因数的a[i]次方,那么直接除去即可。
眼见为实。代码如下：
1 /************************************************************** 2 Problem: 3309 3 User: Doggu 4 Language: C++ 5 Result: Accepted 6 Time:11000 ms 7 Memory:166836 kb 8 ****************************************************************/ 9 10 #include11 #include 12 const int N = 1e7+5;13 int minpa[N], a[N], g[N], prime[N], ptot;14 bool vis[N];15 void EULER(int n) {16 for( int i = 2; i <= n; i++ ) {17 if(!vis[i]) prime[++ptot]=i, minpa[i]=i, a[i]=1, g[i]=1;18 for( int j = 1; j <= ptot; j++ ) {19 if((long long)i*prime[j]>n) break;20 vis[i*prime[j]]=1;21 minpa[i*prime[j]]=prime[j];22 a[i*prime[j]]=1;23 g[i*prime[j]]=(a[i]==1?-g[i]:0);24 if(i%prime[j]==0) {25 minpa[i*prime[j]]=minpa[i]*prime[j];26 a[i*prime[j]]=a[i]+1;27 int temp=i/minpa[i];28 if(temp==1) g[i*prime[j]]=1;29 else g[i*prime[j]]=(a[temp]==a[i*prime[j]]?-g[temp]:0);30 break;31 }32 }33 }34 for( int i = 1; i <= n; i++ ) g[i]+=g[i-1];35 }36 int main() {37 int T, n, m;38 EULER(1e7);39 scanf("%d",&T);40 while(T--) {41 scanf("%d%d",&n,&m);42 long long ans=0;43 if(n>m) std::swap(n,m);44 for( int a = 1, ed; a <= n; a = ed + 1 ) {45 ed=std::min(n/(n/a),m/(m/a));46 ans+=(long long)(g[ed]-g[a-1])*(n/a)*(m/a);47 }48 printf("%lld\n",ans);49 }50 }51