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 
狄利克雷卷积+函数推演+欧拉筛+分块
posted @ 2017-07-13 21:15 Doggu 阅读(...) 评论(...) 编辑 收藏