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

4
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957

Sample Output

35793453939901
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

狄利克雷卷积+函数推演+欧拉筛+分块

转载于:https://www.cnblogs.com/Doggu/p/bzoj3309.html