BZOJ 3309 DZY Loves Math
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
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957
Sample Output
14225956593420
4332838845846
15400094813
HINT
【数据规模】
T<=10000
1<=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 #include <cstdio> 11 #include <algorithm> 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
