1. 程式人生 > >BZOJ 4176: Lucas的數論(莫比烏斯反演+杜教篩)

BZOJ 4176: Lucas的數論(莫比烏斯反演+杜教篩)

alt 思路 難了 lol cst 每一個 r+ 去年 ons

Description

去年的Lucas非常喜歡數論題,但是一年以後的Lucas卻不那麽喜歡了。
在整理以前的試題時,發現了這樣一道題目“求Sigma(f(i)),其中1<=i<=N”,其中 表示i的約數個數。他現在長大了,題目也變難了。
求如下表達式的值:
技術分享圖片

其中 表示ij的約數個數。
他發現答案有點大,只需要輸出模1000000007的值。
Input

第一行一個整數n。
Output

一行一個整數ans,表示答案模1000000007的值。
Sample Input
2
Sample Output
8
HINT

對於100%的數據n <= 10^9。

解題思路

  好神仙的一道題。首先有一個結論就是\(\sigma_0(ij)=\sum\limits_{a\mid i}\sum\limits_{b\mid j}[gcd(a,b)=1]\)

。證明的話就是考慮每一個質數的貢獻,發現左右兩邊相等。有了這個結論就可以推式子了。
\[ ans=\sum\limits_{i=1}^n \sum\limits_{j=1}^n \sigma_0(ij)\]
\[ ans=\sum\limits_{i=1}^n \sum\limits_{j=1}^n\sum\limits_{a\mid i}\sum\limits_{b\mid j}[gcd(a,b)=1]\]

  莫比烏斯反演得:
\[ans=\sum\limits_{d=1}^n \mu(d)\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{n}{d}}\sum\limits_{a\mid i}\sum\limits_{b\mid j}1\]


  把後面的式子合並一下:
  \[ans=\sum\limits_{d=1}^n\mu(d)(\sum\limits_{i=1}^{\frac{n}{d}}\frac{n}{id})^2\]

  發現這樣就可以分塊套分塊去做了,\(\sum\mu\)用杜教篩即可。分塊第一層可以看做枚舉\(n/d\),第二層就是裏面的\(n/(id)\)

代碼

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<map>

using namespace std;
const int N=4000005;
const int MOD=1e9+7;
typedef long long LL;

int n,sum[N],prime[N],cnt,miu[N],ans;
bool vis[N];
map<int,int> mp;

inline int Add(int x,int y){
    x+=y; return (x>=MOD)?(x-MOD):x;
}
inline int Sub(int x,int y){
    x-=y; return x<0?x+MOD:x;
}

inline void prework(int n){
    miu[1]=1;
    for(int i=2;i<=n;i++){
        if(!vis[i]) prime[++cnt]=i,miu[i]=-1;
        for(int j=1;j<=cnt && 1ll*prime[j]*i<=n;j++){
            vis[i*prime[j]]=1;
            if(!(i%prime[j])) break;
            miu[i*prime[j]]=-miu[i];
        }
    }
    for(int i=1;i<=n;i++) sum[i]=sum[i-1]+miu[i];
}

int calc(int x){
    if(x<=N-5) return sum[x];
    if(mp.count(x)) return mp[x];
    int ret=1;
    for(int l=2,r;l<=x;l=r+1){
        r=x/(x/l);
        ret-=(r-l+1)*calc(x/l);
    } 
    return mp[x]=ret;
}

int calc2(int x){
    int ret=0;
    for(int l=1,r;l<=x;l=r+1){
        r=x/(x/l);
        ret=Add(ret,1ll*(r-l+1)*(x/l)%MOD);
    }
    return 1ll*ret*ret%MOD;
}

int main(){
    scanf("%d",&n); prework(4000000);
    for(int l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ans=Add(ans,1ll*Sub(calc(r),calc(l-1))*calc2(n/l)%MOD);
    }
    printf("%d\n",ans);
    return 0;
}

BZOJ 4176: Lucas的數論(莫比烏斯反演+杜教篩)