BZOJ 3309: DZY Loves Math

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

4
10000
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957

Sample Output

35793453939901
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),这个问题就可以在线性筛的时候解决...

代码:

 1 #include<algorithm>
 2 #include<iostream>
 3 #include<cstring>
 4 #include<cstdio>
 5 //by NeighThorn
 6 using namespace std;
 7 //大鹏一日同风起,扶摇直上九万里 
 8 
 9 const int maxn=10000000+5;
10 
11 int a,b,cas,cnt,f[maxn],num[maxn],vis[maxn],lala[maxn],prime[maxn];
12 
13 long long ans;
14 
15 signed main(void){
16     cnt=0;
17     for(int i=2;i<=10000000;i++){
18         if(!vis[i])
19             prime[++cnt]=i,lala[i]=num[i]=1,f[i]=1;
20         for(int j=1;j<=cnt&&i*prime[j]<=10000000;j++){
21             int x=i*prime[j];vis[x]=1;
22             if(i%prime[j]==0){
23                 lala[x]=lala[i];num[x]=num[i]+1;
24                 if(lala[i]==1)
25                     f[x]=1;
26                 else
27                     f[x]=(num[lala[x]]==num[x]?-f[lala[x]]:0);
28                 break;
29             }
30             lala[x]=i,num[x]=1,f[x]=(num[i]==1?-f[i]:0);
31         }
32     }
33     for(int i=1;i<=10000000;i++)
34         f[i]+=f[i-1];
35     scanf("%d",&cas);
36     while(cas--){
37         scanf("%d%d",&a,&b);
38         if(a>b)
39             swap(a,b);
40         ans=0;
41         for(int i=1,r;i<=a;i=r+1){
42             r=min(a/(a/i),b/(b/i));
43             ans+=(long long)(f[r]-f[i-1])*(a/i)*(b/i);
44         }
45         printf("%lld\n",ans);
46     }
47     return 0;
48 }//Cap ou pas cap. Cap. 
View Code

By NeighThorn

 

posted @ 2016-12-26 12:03  NeighThorn  阅读(185)  评论(0编辑  收藏  举报