2018 ACM-ICPC 沈阳网络赛 G. Spare Tire （容斥）

最新推荐文章于 2019-10-09 20:58:22 发布

兜率工

最新推荐文章于 2019-10-09 20:58:22 发布

阅读量155

点赞数

分类专栏：组合数学

本文链接：https://blog.csdn.net/qq_40679299/article/details/82891447

版权

组合数学专栏收录该内容

41 篇文章 2 订阅

订阅专栏

根据题意打表可以得到通项公式 n*(n+1)

转化题意，即求在1-n中与m互质的数套用上述公式的和。

那么这就是最简单的容斥模型。看这个：https://blog.csdn.net/qq_40679299/article/details/82873829

将求个数改为求和即可，原理跟上面一样，

计算平方和的时候直接套用公式就好。

附代码：

#include<bits/stdc++.h>

using namespace std;

typedef long long ll;

const ll mod=1e9+7;

int p[50];
int cnt;

ll exgcd(ll a,ll b,ll &x,ll &y){
    if(a%b==0){
        x=0ll;y=1ll;
        return b;
    }
    ll v,tx,ty;
    v=exgcd(b,a%b,tx,ty);
    x=ty;
    y=tx-a/b*ty;
    return v;
}

ll inv(ll a,ll p){
    if(!a) return 0ll;
    ll x,y;
    exgcd(a,p,x,y);
    x=(x%p+p)%p;
    return x;
}

ll cal(ll n){
    ll res=0;
    for(int i=1;i<(1<<cnt);i++){
        int t=i,k=0;
        ll tmp=1;
        int len=0;
        while(t){
            if(t&1){
                tmp=tmp*p[k]%mod;
                len++;
            }
            t>>=1;
            k++;
        }
        ll f=tmp*tmp%mod*(n/tmp)%mod*(n/tmp+1)%mod*(n/tmp*2+1)%mod*inv(6,mod)%mod+(n/tmp)*(1+n/tmp)%mod*inv(2,mod)%mod*tmp%mod;
        if(len&1) res+=(f%mod);
        else res-=(f%mod);
        res=(res+mod)%mod;
    }
    return ((n*(n+1)%mod*(2*n+1)%mod*inv(6,mod)%mod+(n*(1+n)%mod)*inv(2,mod)%mod)%mod-res+mod)%mod;
}

int main(){
    ll n,m;
    while(scanf("%lld%lld",&n,&m)!=EOF){
        cnt=0;
        for(int i=2;i*i<=m;i++){
            if(m%i) continue;
            while(m%i==0) m/=i;
            p[cnt++]=i;
        }
        if(m!=1) p[cnt++]=m;
        printf("%lld\n",cal(n));
    }
    return 0;
}