51Nod 1584 加权约数和

51Nod 1584 加权约数和


要求:\(\sum_{i=1}^n\sum_{j=1}^n\max(i,j)\sigma_1(ij)\)

枚举\(\max(i,j)\),式子变为:\(2\sum_{i=1}^ni\sum_{j=1}^i\sigma_1(ij)-\sum_{i=1}^ni\sigma_1(i^2)\)

这里不加证明地给出结论:\(\sigma_1^k(ij)=\sum_{p|i}\sum_{q|j}(p\frac{j}{q})^k[\gcd(p,q)=1]\)

先看右边部分,可以线性筛\(\sigma_1(n^2)\)

(筛约数个数和的方法:\(\sigma_1(n=\prod_{i=1}^mp_i^{e_i})=\prod_{i=1}^m(1+p_i+\cdots+p_i^{e_i})\),筛平方的约数个数和类似,事实上任意\(k\)次方的都能筛)

或者也可以推式子:

\(\sum_{i=1}^ni\sigma_1(i^2)\)

\(=\sum_{i=1}^ni\sum_{p|i}\sum_{q|i}p\frac iq\sum_{d|p,d|q}\mu(d)\)

\(=\sum_{i=1}^ni\sum_{p|i}\sum_{q|i}p\frac iq\sum_{d|p,d|q}\mu(d)\)

\(=\sum_{i=1}^ni\sum_{d|i}\mu(d)\sum_{d|p}^ip\sum_{d|q}^i\frac iq\)

\(=\sum_{i=1}^ni\sum_{d|i}\mu(d)\sum_{p|\frac id}^idp\sum_{q|\frac id}^i\frac i{dq}\)

\(=\sum_{i=1}^ni\sum_{d|i}\mu(d)d\sigma_1(\frac id)^2\)

左边部分类似,直接写结论了:

\(\sum_{i=1}^ni\sum_{j=1}^i\sigma_1(ij)=\sum_{i=1}i\sum_{d|i}\mu(\frac id)(\frac id)\sigma_1(d)S_{\sigma_1}(d)\)

其中\(S_f(n)=\sum_{i=1}^nf(i)\)


右边用线性筛:

#include<bits/stdc++.h>
#define il inline
#define vd void
#define mod 1000000007
typedef long long ll;
il ll gi(){
    ll x=0,f=1;
    char ch=getchar();
    while(!isdigit(ch))f^=ch=='-',ch=getchar();
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return f?x:-x;
}
int pr[1000010],P,sd[1000010],_sd[1000010],sum[1000010],ssd[1000010],mu[1000010];
int sd2[1000010],_sd2[1000010],ssd2[1000010],sum2[1000010];
bool yes[1000010];
int ans[1000010];
int main(){
#ifdef XZZSB
    freopen("in.in","r",stdin);
    freopen("out.out","w",stdout);
#endif
    int N=1000000;
    sd[1]=sd2[1]=1;mu[1]=1;
    for(int i=2;i<=N;++i){
        if(!yes[i])pr[++P]=i,sd[i]=sum[i]=i+1,_sd[i]=1,mu[i]=mod-1,sd2[i]=sum2[i]=(1+i+1ll*i*i)%mod,_sd2[i]=1;
        for(int j=1;j<=P&&i*pr[j]<=N;++j){
            yes[i*pr[j]]=1;
            if(i%pr[j]==0){
                sum[i*pr[j]]=(1ll*sum[i]*pr[j]+1)%mod;
                sd[i*pr[j]]=1ll*_sd[i]*sum[i*pr[j]]%mod;
                _sd[i*pr[j]]=_sd[i];
                sum2[i*pr[j]]=(1ll*sum2[i]*pr[j]%mod*pr[j]+pr[j]+1)%mod;
                sd2[i*pr[j]]=1ll*_sd2[i]*sum2[i*pr[j]]%mod;
                _sd2[i*pr[j]]=_sd2[i];
                mu[i*pr[j]]=0;
                break;
            }
            sum[i*pr[j]]=1+pr[j];
            sd[i*pr[j]]=1ll*sd[i]*sum[i*pr[j]]%mod;
            _sd[i*pr[j]]=sd[i];
            sum2[i*pr[j]]=(1ll*pr[j]*pr[j]+pr[j]+1)%mod;
            sd2[i*pr[j]]=1ll*sd2[i]*sum2[i*pr[j]]%mod;
            _sd2[i*pr[j]]=sd2[i];
            mu[i*pr[j]]=mod-mu[i];
        }
    }
    for(int i=1;i<=N;++i)ssd[i]=(ssd[i-1]+sd[i])%mod;
    for(int i=1;i<=N;++i)ssd2[i]=(ssd2[i-1]+1ll*i*sd2[i])%mod;
    for(int i=1;i<=N;++i)
        for(int j=i;j<=N;j+=i)
            if(mu[i])ans[j]=(ans[j]+2ll*j*mu[i]%mod*i%mod*sd[j/i]%mod*ssd[j/i])%mod;
    for(int i=1;i<=N;++i)ans[i]=(ans[i]+ans[i-1])%mod;
    int T=gi(),n;
    for(int i=1;i<=T;++i)n=gi(),printf("Case #%d: %d\n",i,(ans[n]-ssd2[n]+mod)%mod);
    return 0;
}

右边用推式子:

#include<bits/stdc++.h>
#define il inline
#define vd void
#define mod 1000000007
typedef long long ll;
il ll gi(){
    ll x=0,f=1;
    char ch=getchar();
    while(!isdigit(ch))f^=ch=='-',ch=getchar();
    while(isdigit(ch))x=x*10+ch-'0',ch=getchar();
    return f?x:-x;
}
int pr[1000010],P,sd[1000010],_sd[1000010],sum[1000010],ssd[1000010],mu[1000010];
bool yes[1000010];
int ans[1000010];
int main(){
#ifdef XZZSB
    freopen("in.in","r",stdin);
    freopen("out.out","w",stdout);
#endif
    int N=1000000;
    sd[1]=1;mu[1]=1;
    for(int i=2;i<=N;++i){
        if(!yes[i])pr[++P]=i,sd[i]=i+1,_sd[i]=1,sum[i]=i+1,mu[i]=mod-1;
        for(int j=1;j<=P&&i*pr[j]<=N;++j){
            yes[i*pr[j]]=1;
            if(i%pr[j]==0){
                sum[i*pr[j]]=(1ll*sum[i]*pr[j]+1)%mod;
                sd[i*pr[j]]=1ll*_sd[i]*sum[i*pr[j]]%mod;
                _sd[i*pr[j]]=_sd[i];
                mu[i*pr[j]]=0;
                break;
            }
            sum[i*pr[j]]=1+pr[j];
            sd[i*pr[j]]=1ll*sd[i]*sum[i*pr[j]]%mod;
            _sd[i*pr[j]]=sd[i];
            mu[i*pr[j]]=mod-mu[i];
        }
    }
    for(int i=1;i<=N;++i)ssd[i]=(ssd[i-1]+sd[i])%mod;
    for(int i=1;i<=N;++i)
        for(int j=i;j<=N;j+=i){
            ans[j]=(ans[j]+2ll*j*mu[j/i]%mod*(j/i)%mod*sd[i]%mod*ssd[i])%mod;
            ans[j]=(ans[j]-1ll*j*mu[i]%mod*i%mod*sd[j/i]%mod*sd[j/i]%mod+mod)%mod;
        }
    for(int i=1;i<=N;++i)ans[i]=(ans[i]+ans[i-1])%mod;
    int T=gi(),n;
    for(int i=1;i<=T;++i)n=gi(),printf("Case #%d: %d\n",i,ans[n]);
    return 0;
}

转载于:https://www.cnblogs.com/xzz_233/p/11142149.html

  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值