zcmimi's blog
查看原题

点击跳转

前置知识:

  1. 杜教筛(包括狄利克雷卷积)

  2. 数论分块

  3. 欧拉函数 或 莫比乌斯函数

(会杜教筛的大佬上述都会吧)

欧拉函数卷积推导

\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j)

根据 \sum\limits_{i|n}\varphi(n)=n (1 * \varphi = Id)

= \sum_{i=1}^n\sum_{j=1}^n ij \sum_{k|i,k|j} \varphi(k)

调换枚举顺序

= \sum_{k=1}^n \varphi(k) \sum_{k|i,k|j} ij \\ = \sum_{k=1}^n \varphi(k) \sum_{k|i}i \sum_{k|j}j \\ = \sum_{k=1}^n \varphi(k) (\sum_{i=1}^{\left \lfloor\frac nk\right\rfloor}ki)^2 \\ = \sum_{k=1}^n \varphi(k) k^2(\sum_{i=1}^{\left \lfloor\frac nk\right\rfloor}i)^2

莫比乌斯函数推导

\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j) \\ =\sum_{d=1}^n d \sum_{i=1}^n\sum_{j=1}^n ij [gcd(i,j)=d] \\ =\sum_{d=1}^n d^3 \sum_{i=1}^{\left \lfloor\frac nd\right\rfloor}\sum_{j=1}^{\left \lfloor\frac nd\right\rfloor} ij [gcd(i,j)=1] \\ =\sum_{d=1}^n d^3 \sum_{i=1}^{\left \lfloor\frac nd\right\rfloor}\sum_{j=1}^{\left \lfloor\frac nd\right\rfloor} ij \sum_{k|i,k|j}\mu(k) \\ =\sum_{d=1}^n d^3 \sum_{k=1}^{\left \lfloor\frac nd\right\rfloor}\mu(k) \sum_{k|i,k|j} ij \\ =\sum_{d=1}^n d^3 \sum_{k=1}^{\left \lfloor\frac nd\right\rfloor}\mu(k) k^2(\sum_{k|i}i)^2

为了方便,设sum(T)=\sum_{i=1}^T i

=\sum_{d=1}^n d^3 \sum_{k=1}^{\left \lfloor\frac nd\right\rfloor}\mu(k)k^2sum(\left \lfloor \frac n{dk}\right \rfloor)^2

T=dk

=\sum_{T=1}^n sum(\left \lfloor \frac nT\right \rfloor)^2 \sum_{d|T}d^3 \frac{T^2}{d^2}\mu(\frac Td) \\ =\sum_{T=1}^n sum(\left \lfloor \frac nT\right \rfloor)^2 T^2\sum_{d|T}d\mu(\frac Td)

其中: \sum_{d|T}d\mu(\frac Td)=\varphi(T)(\mu * Id = \varphi)

=\sum_{T=1}^n sum(\left \lfloor \frac nT\right \rfloor)^2 T^2 \varphi(T)

上述两种方法推导出的最终结果是一样的

接下来:

\sum_{k=1}^n \varphi(k)用杜教筛求出前缀和

k^2(\sum_{i=1}^{\left \lfloor\frac nk\right\rfloor}i)^2可以用数论分块求解

#include<bits/stdc++.h>
typedef long long ll;
const int N=5000000,P=1000007;
struct hash{
    int cnt=0,head[P];
    struct edge{int to,nxt,w;}e[P];
    void add(int x,int y,int w){e[++cnt].to=y;e[cnt].nxt=head[x];head[x]=cnt;e[cnt].w=w;}
    void ins(int x,int v){add(x%P,x,v);}
    int ask(int x){for(int i=head[x%P];i;i=e[i].nxt)if(e[i].to==x)return e[i].w;return -1;}
}S;
int pri[N+11],phi[N+11],p,inv,lim;
bool b[N+11];
int sum(ll n){n%=p;return (1ll*n*(n+1)>>1)%p;}
int SUM(ll n){n%=p;return 1ll*n*(n+1)%p*(n+n+1)%p*inv%p;}
il void init(){
    phi[1]=1;
    int cnt=0;
    for(int i=2;i<=lim;++i){
        if(!b[i])pri[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt&&i*pri[j]<=lim;++j){
            b[i*pri[j]]=1;
            if(i%pri[j])phi[i*pri[j]]=phi[i]*(pri[j]-1)%p;
            else{phi[i*pri[j]]=phi[i]*pri[j]%p;break;}
        }
    }
    for(int i=1;i<=lim;++i)phi[i]=(phi[i-1]+1ll*i*i%p*phi[i]%p)%p;
}
int PHI(ll n){
    if(n<=lim)return phi[n];
    int t=S.ask(n);if(~t)return t;
    int res=sum(n);res=1ll*res*res%p;
    for(ll l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        int s=(SUM(r)-SUM(l-1)+p)%p;
        res=(res-1ll*PHI(n/l)*s%p+p)%p;
    }
    S.ins(n,res);
    return res;
}
int pw(int x,int b){
    int ans=1;
    while(b){
        if(b&1)ans=1ll*ans*x%p;
        b>>=1;x=1ll*x*x%p;
    }
    return ans;
}
int main(){
    ll n;
    scanf("%lld%lld",&p,&n);
    inv=pw(6,p-2),lim=pow(n,2.0/3);
    init();
    int res=0,s,phi;
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        s=sum(n/l),phi=(PHI(r)-PHI(l-1)+p)%p;
        res=(res+1ll*s*s%p*phi%p)%p;
    }
    printf("%d\n",res);
}

下方是更加玄学优化的代码

#include<bits/stdc++.h>
typedef long long ll;
const int N=5000000;
int b[N+11],pri[N+11],phi[N+11],p,inv,lim,s[N],S[N],sqr,ans[N];ll n;
int id(ll x){return x<=sqr?x:sqr+n/x;}
int PHI(ll n){
    int x=id(n);
    if(n<=lim)return phi[n];
    if(ans[x])return ans[x];
    int res=s[x];
    for(register ll l=2,r,i;l<=n;l=r+1)
        i=n/l,r=n/i,
        res=(res-1ll*PHI(i)*(S[id(r)]-S[id(l-1)]))%p;
    return ans[x]=res;
}
int pw(int x,int b){
    int ans=1;
    while(b){
        if(b&1)ans=1ll*ans*x%p;
        b>>=1;x=1ll*x*x%p;
    }
    return ans;
}
int main(){
    scanf("%d%lld",&p,&n);
    inv=pw(6,p-2),lim=pow(n,2.0/3);sqr=sqrt(n);
    phi[1]=1;
    int cnt=0,x,res=0;
    for(register int i=2;i<=lim;++i){
        if(!b[i])b[i]=pri[++cnt]=i,phi[i]=i-1;
        for(register int j=1,y;j<=cnt&&pri[j]<=b[i]&&(y=pri[j]*i)<=lim;++j)
            b[y]=pri[j],phi[y]=i%pri[j]?phi[i]*phi[pri[j]]:phi[i]*pri[j];
    }
    for(register int i=1;i<=lim;++i)phi[i]=(phi[i-1]+1ll*i*i%p*phi[i]%p)%p;
    for(register ll l=1,r,i;l<=n;l=r+1)
        i=n/l,r=n/i,
        x=id(i),
        i%=p,
        S[x]=1ll*inv*i%p*(i+1)%p*(i+i+1)%p,
        s[x]=((1ll*i*(i+1))>>1)%p,s[x]=1ll*s[x]*s[x]%p;    
    for(register ll l=1,r,i;l<=n;l=r+1)
        i=n/l,r=n/i,
        res=(res+1ll*s[id(i)]*(PHI(r)-PHI(l-1)))%p;
    if(res<0)res+=p;
    printf("%d\n",res);
}
LG 3768 简单的数学题
comment评论
Search
search