[ 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，累加對其他數的貢獻，實際復雜度會下降不少。

講不明白，上代碼，還是直接頹代碼比較直接痛快。

for(int i=1,j=prime[1];i<=nop;j=prime[++i])
        for(int k=j;k<=maxn;k+=j)
            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：本人碼風自帶防抄，對拍等等可以，提交請謹慎。

#include<cstdio>
#include<cmath>
using namespace std;
const int maxn=1e7;
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;
int main(){
    for(int i=1;i<=maxn;++i)miu[i]=1;
    for(int i=2;i<=maxn;++i){
        if(!notprime[i]){nop++;prime[nop]=i;miu[i]=-1;}
        for(int j=1;j<=nop&&i*prime[j]<=maxn;++j){
            notprime[i*prime[j]]=1;
            if(i%prime[j]){miu[i*prime[j]]=-miu[i];continue;}
            miu[i*prime[j]]=0;
            break;
        }
    }
    for(int i=1,j=prime[1];i<=nop;j=prime[++i])
        for(int k=j;k<=maxn;k+=j)
            f[k]+=miu[k/j];
    for(int i=1;i<=maxn;++i)qz[i]=qz[i-1]+f[i];
    scanf("%d",&t);
    while(t--){
        scanf("%d%d",&m,&n);sda=0;sdb=0;if(m>n)sqr=m,m=n,n=sqr;
        sqr=sqrt(n);ans=0;
        for(int i=1;i<=sqr;++i)ca[++sda]=n/i,enda[sda]=i;
        for(int i=sqr;i;--i)ca[++sda]=i,enda[sda]=n/i;
        sqr=sqrt(m);
        for(int i=1;i<=sqr;++i)cb[++sdb]=m/i,endb[sdb]=i;
        for(int i=sqr;i;--i)cb[++sdb]=i,endb[sdb]=m/i;
        while(enda[sda-1]>=m)sda--;enda[sda]=sqr=m;
        while(sda&&sdb){
            if(enda[sda-1]==endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--,sda--;
            else if(enda[sda-1]<endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--;
            else ans+=(qz[sqr]-qz[enda[sda-1]])*ca[sda]*cb[sdb],sqr=enda[sda-1],sda--;
        }
        printf("%lld\n",ans);
    }
}

 update 6/13:依據Dybala_zdy大神的優化.

先讀入全部詢問，根據其最大值確定線性篩等預處理的范圍，優化至5700ms。

#include<cstdio>
#include<cmath>
using namespace std;
#define maxnnn 10000000
inline int max(int a,int b){return a>b?a:b;}
int maxn=0,tm[10005],tn[10005];
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;
int main(){scanf("%d",&t);
    for(int i=1;i<=t;++i){
        scanf("%d%d",&tm[i],&tn[i]);
        maxn=max(maxn,max(tm[i],tn[i]));
    }
    for(int i=1;i<=maxn;++i)miu[i]=1;
    for(int i=2;i<=maxn;++i){
        if(!notprime[i]){nop++;prime[nop]=i;miu[i]=-1;}
        for(int j=1;j<=nop&&i*prime[j]<=maxn;++j){
            notprime[i*prime[j]]=1;
            if(i%prime[j]){miu[i*prime[j]]=-miu[i];continue;}
            miu[i*prime[j]]=0;
            break;
        }
    }
    for(int i=1,j=prime[1];i<=nop;j=prime[++i])
        for(int k=j;k<=maxn;k+=j)
            f[k]+=miu[k/j];
    for(int i=1;i<=maxn;++i)qz[i]=qz[i-1]+f[i];
    for(int ii=1;ii<=t;++ii){
        m=tm[ii];n=tn[ii];sda=0;sdb=0;if(m>n)sqr=m,m=n,n=sqr;
        sqr=sqrt(n);ans=0;
        for(int i=1;i<=sqr;++i)ca[++sda]=n/i,enda[sda]=i;
        for(int i=sqr;i;--i)ca[++sda]=i,enda[sda]=n/i;
        sqr=sqrt(m);
        for(int i=1;i<=sqr;++i)cb[++sdb]=m/i,endb[sdb]=i;
        for(int i=sqr;i;--i)cb[++sdb]=i,endb[sdb]=m/i;
        while(enda[sda-1]>=m)sda--;enda[sda]=sqr=m;
        while(sda&&sdb){
            if(enda[sda-1]==endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--,sda--;
            else if(enda[sda-1]<endb[sdb-1])ans+=(qz[sqr]-qz[endb[sdb-1]])*ca[sda]*cb[sdb],sqr=endb[sdb-1],sdb--;
            else ans+=(qz[sqr]-qz[enda[sda-1]])*ca[sda]*cb[sdb],sqr=enda[sda-1],sda--;
        }
        printf("%lld\n",ans);
    }
}