[CQOI2015]选数(bzoj3930 莫比乌斯反演+杜教筛+累加有上下界)

题目:

[CQOI2015]选数

 

题意:

\sum_{i1=L}^{R}\sum_{i2=L}^{R}...\sum_{in=L}^{R}[gcd(i1,i2,...,in)=k]

 

思路:

若L%k==0,那么累加的下界为 \left \lfloor \frac{L}{k} \right \rfloor ,若L%k!=0,那么累加的下界为  \left \lfloor \frac{L}{k} \right \rfloor+1 ,综合一下,下界为 \left \lfloor \frac{L-1}{i} \right \rfloor+1

我们设下界 = t 。

原式 = \sum_{i1=T}^{R/K}\sum_{i2=T}^{R/K}...\sum_{in=T}^{R/K}[gcd(i1,i2,...,in)=1] 

f(a)=\sum_{i1=T}^{R/K}\sum_{i2=T}^{R/K}...\sum_{in=T}^{R/K}[gcd(i1,i2,...,in)=a]

g(a)=\sum_{a|d}^{ }f(d)

f(a)=\sum_{a|d}^{ }\mu (\frac{d}{a})*g(d) (莫比乌斯反演)

f(1)=\sum_{1|d}^{ }\mu (d)*g(d)=\sum_{d=1}^{R/k }\mu (d)*g(d)

g(a)=\sum_{a|d}^{ }f(d)=\sum_{i1=t}^{R/d}*\sum_{i2=t}^{R/d}*...*\sum_{in=t}^{R/d}[a|gcd(i1,i2,...,in)]

把a除掉后,上届为 \left \lfloor\left \lfloor R/k \right \rfloor /d \right \rfloor ,下界为 \left \lfloor\left \lfloor (t-1) \right \rfloor /d \right \rfloor = \left \lfloor\left \lfloor (L-1)/k \right \rfloor /d \right \rfloor

因此, g(d)=(\left \lfloor \frac{R}{kd} \right \rfloor-\left \lfloor \frac{L-1}{kd} \right \rfloor)^{n}

因此,可以对 g(d) 进行数论分块,通过杜教筛求莫比乌斯函数的前缀和。

code:

#include<bits/stdc++.h>
using namespace std;

typedef long long ll;

const int MAX = 1e6;
const ll MOD = 1000000007;

ll n,k,L,H;
ll mu[MAX+10],prime[MAX+10];
bool check[MAX+10];

inline ll read()
{
    register ll x=0,t=1;register char ch=getchar();
    while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
    if(ch=='-')t=-1,ch=getchar();
    while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
    return x*t;
}

void get_sum()
{
    memset(check,false,sizeof(check));
    mu[1]=1;
    int tot=0;
    for(int i=2;i<=MAX;i++){
        if(!check[i]){
            prime[tot++]=i;
            mu[i]=-1;
        }
        for(int j=0;j<tot&&i*prime[j]<=MAX;j++){
            check[i*prime[j]]=true;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            else{
                mu[i*prime[j]]=mu[i]*(-1);
            }
        }
    }
    for(int i=2;i<=MAX;i++)
        mu[i]=mu[i-1]+mu[i];
}

map<ll,ll>dp;

ll cal(ll x)
{
    if(x<=MAX)   return mu[x];
    if(dp[x])   return dp[x];
    ll pos;
    ll tmp=0;
    for(ll i=2;i<=x;i=pos+1){
        pos=x/(x/i);
        tmp+=(pos-i+1)*cal(x/i);
    }
    dp[x]=1-tmp;
    return dp[x];
}

ll multiply(ll a, ll b)
{
	ll ans = 1;
	while (b)
	{
		if (b & 1)
		{
			ans = ((ans%MOD)*(a%MOD)) % MOD;
			b--;
		}
		b /= 2;
		a = ((a%MOD)*(a%MOD)) % MOD;
	}
	return ans;
}

int main()
{
    get_sum();
    n=read();
    k=read();
    L=read();
    H=read();
    ll ans=0;
    ll x=H/k;
    ll y=(L-1)/k;
    ll pos=0;
    ll i;
    for(i=1;i<=x;i=pos+1){
        if(i<=y)
            pos=min(x/(x/i),y/(y/i));
        else
            pos=x/(x/i);
        ans=(ans+multiply(x/i-y/i,n)*((cal(pos)-cal(i-1)+MOD)%MOD)%MOD)%MOD;
    }
    printf("%lld\n",ans);
    return 0;
}

 

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

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值