[ bzoj2820] YY的GCD

[ bzoj2820] YY的GCD

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分解質因數爲a₁^b₁a₂^b₂...a_m^b_m  (a₁,a₂...a_m爲彼此不一樣的質數)

　　　　　　　　　　　　若b₁,b₂...b_n中有大於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(n²)級別，難以接受。

轉換思路，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];

你要頹代碼麼？珍重OI節操

實際複雜度是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 }

Update 6/13