梁祝之套马的汉子,sbs官网直播,快乐蹲坑
time limit : 3000 ms
description
神犇yy虐完数论后给傻×kac出了一题给定n, m,求1<=x<=n, 1<=y<=m且gcd(x, y)为质数的(x, y)有多少对kac这种傻×必然不会了，于是向你来请教……
多组输入
input
第一行一个整数t 表述数据组数接下来t行，每行两个正整数，表示n, m
output
t行，每行一个整数表示第i组数据的结果
sample input
2
10 10
100 100
sample output
30
2791
hint
t = 10000 n, m <= 10000000
莫比乌斯函数，莫比乌斯反演，欧拉线性筛，数论分块（分段），约数个数
大概说一说吧我也是第一次做这种题，会的跳过直接进正题
莫比乌斯函数：k分解质因数为a1b1a2b2...ambm (a1,a2...am为彼此不同的质数)
若b1,b2...bn中有大于1的项，则μ(k)=0
否则μ(k)=(-1)m
莫比乌斯反演：对于已知函数f和g，若满足
f(n)=σ g(d) (其中d是n的约数)
则有
g(n)=σ μ(d)*f(n/d) (其中d是n的约数)
常见形式为σ μ(d)= [n] (其中d是n的约数，[n]等同于代码中的（n==1?1:0）)
欧拉线性筛：可以在线性复杂度内求出素数，欧拉函数，莫比乌斯函数，详见代码
约数个数：n的约数最多只有2√n个。
对于一个约数a，p/a也是它的约数，除恰好为√n外成对存在。
如果a<√n，那么p/a>√n。小于√n的显然最多有√n个，那么最多就有√n对即2√n个。
数论分块:求和σ时,随着参数递增而函数值变化频率很低时,直接求每个函数值的出现次数×函数值。
如：给出k，求σ(k/i) (/为整除，i<=k)。
k=135时，只要68<=i<=135，(k/i)就是1，每次都除显然浪费时间，故求出68与135这两个端点，乘以(k/i)的值1即可。
结合约数个数那一条，这个问题可以从o(n)降低为o(2√n)求解
话说这个题目名竟然没有被和谐
首先，并不难列出最基础的式子：
其中p为质数，那么我们可以通过枚举质数来求解，假设n<=m，否则swap。
将p提出，同时x和y的含义发生变化，变为（是p的几倍），这样x和y的枚举上界就缩小了p倍。
根据莫比乌斯反演：σ μ(d)= [n] (其中d是n的约数，[n]等同于代码中的（n==1?1:0）)
x和y能整除gcd(x,y)，而gcd(x,y)能整除d，则x和y都能整除d。
改变枚举顺序，先枚举d,那么x，y都是d的倍数，则再次改变x,y的含义(pd的倍数)。
d的范围自然是min(n/p,m/p),这样就可以把μ(d)提到外层。
里面两层的σ的累加值既然是1，那么就是枚举几次就是几咯(|_ a_|表示a向下取整)
可以发现σ后面的式子与p*d关联挺大，所以转换思路枚举k=p*d
(不要问我是怎么想到的，不然我也不会是一个在这里写博客的蒟蒻了)
(从这个式子往下省略下取整符号，我太弱了找不到它！所有除号都表示整除，为了方便区分加上了括号)
其中(n/k)*(m/k)项与p无关，可以提出来。
看起来好像化简到头了，为了放松放松脑子，我们先不推公式了，我们来初始化它。
σ μ(k/p)只与k有关而与m，n无关，可以预处理出来，设它是f(k)=σ μ(k/p),其中p是小于等于k的质数。
看起来通过枚举k及k以下的p，复杂度是o(n2)级别，难以接受。
转换思路，k不能枚举那就枚举p呗。
通过枚举每个质数p，累加对其他数的贡献，实际复杂度会下降不少。
讲不明白，上代码，还是直接颓代码比较直接痛快。
1 for(int i=1,j=prime[1];i<=nop;j=prime[++i]) 2 for(int k=j;k<=maxn;k+=j) 3 f[k]+=miu[k/j];
实际复杂度是o(n log n),可以接受。
那么f函数就搞定了，原算式更简单了。
这样我们就可以o(n)处理每个询问了，可是面对数据范围还是只有部分分。
内心：就这样算了，考场上肯定就码暴力了。
不行，不屈的衡中学子我们要追求卓越。
我们可以发现当k趋近于n时，(n/k)与(m/k)两项变化很慢,k增大1时这两项很可能都不变,考虑分段处理。
而每次都在变化的f(k),我们可以用前缀和来求区间和,存在totf里，即totf(k)=σf(i) (i<=k)
(n/k)最多2√n段，(m/k)也是这个级别，所以一共最多4√n段。
每次询问的复杂度为√n级别，可以接受。
呼，完成啦！
最终10个测试点8900ms，比较充裕。
附码量很小（脑量很大）的ac代码。
ps：分段那部分打麻烦了，其实码量可以更小，详见其他神犇
（我这个蒟蒻总不能比他们打的好吧，其实是我故意的）
ps：本人码风自带防抄，对拍等等可以，提交请谨慎。
我要评论 1 #include<cstdio>
2 #include<cmath>
3 using namespace std;
4 const int maxn=1e7;
5 int prime[maxn],notprime[maxn+5],nop,n,m,miu[maxn+5],t,enda[6666],endb[6666],ca[6666],cb[6666],sda,sdb,sqr;long long f[maxn+5],qz[maxn+5],ans;
6 int main(){
7 for(int i=1;i<=maxn;++i)miu[i]=1;
8 for(int i=2;i<=maxn;++i){
9 if(!notprime[i]){nop++;prime[nop]=i;miu[i]=-1;}
10 for(int j=1;j<=nop&&i*prime[j]<=maxn;++j){
11 notprime[i*prime[j]]=1;
12 if(i%prime[j]){miu[i*prime[j]]=-miu[i];continue;}
13 miu[i*prime[j]]=0;
14 break;
15 }
16 }
17 for(int i=1,j=prime[1];i<=nop;j=prime[++i])
18 for(int k=j;k<=maxn;k+=j)
19 f[k]+=miu[k/j];
20 for(int i=1;i<=maxn;++i)qz[i]=qz[i-1]+f[i];
21 scanf("%d",&t);
22 while(t--){
23 scanf("%d%d",&m,&n);sda=0;sdb=0;if(m>n)sqr=m,m=n,n=sqr;
24 sqr=sqrt(n);ans=0;
25 for(int i=1;i<=sqr;++i)ca[++sda]=n/i,enda[sda]=i;
26 for(int i=sqr;i;--i)ca[++sda]=i,enda[sda]=n/i;
27 sqr=sqrt(m);
28 for(int i=1;i<=sqr;++i)cb[++sdb]=m/i,endb[sdb]=i;
29 for(int i=sqr;i;--i)cb[++sdb]=i,endb[sdb]=m/i;
30 while(enda[sda-1]>=m)sda--;enda[sda]=sqr=m;
31 while(sda&&sdb){
32 if(enda[sda-1]==endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--,sda--;
33 else if(enda[sda-1]<endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--;
34 else ans+=(qz[sqr]-qz[enda[sda-1]])*ca[sda]*cb[sdb],sqr=enda[sda-1],sda--;
35 }
36 printf("%lld\n",ans);
37 }
38 }
update 6/13:依据dybala_zdy大神的优化.
先读入全部询问，根据其最大值确定线性筛等预处理的范围，优化至5700ms。
1 #include<cstdio>
2 #include<cmath>
3 using namespace std;
4 #define maxnnn 10000000
5 inline int max(int a,int b){return a>b?a:b;}
6 int maxn=0,tm[10005],tn[10005];
7 int prime[maxnnn],notprime[maxnnn+5],nop,n,m,miu[maxnnn+5],t,enda[6666],endb[6666],ca[6666],cb[6666],sda,sdb,sqr;long long f[maxnnn+5],qz[maxnnn+5],ans;
8 int main(){scanf("%d",&t);
9 for(int i=1;i<=t;++i){
10 scanf("%d%d",&tm[i],&tn[i]);
11 maxn=max(maxn,max(tm[i],tn[i]));
12 }
13 for(int i=1;i<=maxn;++i)miu[i]=1;
14 for(int i=2;i<=maxn;++i){
15 if(!notprime[i]){nop++;prime[nop]=i;miu[i]=-1;}
16 for(int j=1;j<=nop&&i*prime[j]<=maxn;++j){
17 notprime[i*prime[j]]=1;
18 if(i%prime[j]){miu[i*prime[j]]=-miu[i];continue;}
19 miu[i*prime[j]]=0;
20 break;
21 }
22 }
23 for(int i=1,j=prime[1];i<=nop;j=prime[++i])
24 for(int k=j;k<=maxn;k+=j)
25 f[k]+=miu[k/j];
26 for(int i=1;i<=maxn;++i)qz[i]=qz[i-1]+f[i];
27 for(int ii=1;ii<=t;++ii){
28 m=tm[ii];n=tn[ii];sda=0;sdb=0;if(m>n)sqr=m,m=n,n=sqr;
29 sqr=sqrt(n);ans=0;
30 for(int i=1;i<=sqr;++i)ca[++sda]=n/i,enda[sda]=i;
31 for(int i=sqr;i;--i)ca[++sda]=i,enda[sda]=n/i;
32 sqr=sqrt(m);
33 for(int i=1;i<=sqr;++i)cb[++sdb]=m/i,endb[sdb]=i;
34 for(int i=sqr;i;--i)cb[++sdb]=i,endb[sdb]=m/i;
35 while(enda[sda-1]>=m)sda--;enda[sda]=sqr=m;
36 while(sda&&sdb){
37 if(enda[sda-1]==endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--,sda--;
38 else if(enda[sda-1]<endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--;
39 else ans+=(qz[sqr]-qz[enda[sda-1]])*ca[sda]*cb[sdb],sqr=enda[sda-1],sda--;
40 }
41 printf("%lld\n",ans);
42 }
43 }
如对本文有疑问，请在下面进行留言讨论，广大热心网友会与你互动！！ 点击进行留言回复
如何在没有core文件的情况下用dmesg+addr2line定位段错误
用QT制作3D点云显示器——QtDataVisualization
网友评论