BZOJ 3930 [CQOI2015]选数 分块+前缀和+玄

BZOJ 3930 [CQOI2015]选数

Solution

题目要求:

a1=LRa2=LRaN=LR[gcd(a1,a2,,aN)=K] ∑ a 1 = L R ∑ a 2 = L R ⋯ ∑ a N = L R [ g c d ( a 1 , a 2 , ⋯ , a N ) = K ]

设要求的答案为 f(K) f ( K ) ,表示

a1=LRa2=LRaN=LR[gcd(a1,a2,,aN)=K] ∑ a 1 = L R ∑ a 2 = L R ⋯ ∑ a N = L R [ g c d ( a 1 , a 2 , ⋯ , a N ) = K ]

考虑构造一个函数 F(K) F ( K ) ,表示
a1=LRa2=LRaN=LR[K|gcd(a1,a2,,aN)] ∑ a 1 = L R ∑ a 2 = L R ⋯ ∑ a N = L R [ K | g c d ( a 1 , a 2 , ⋯ , a N ) ]

那么 f,F f , F 之间有如下关系
F(n)=Rn|df(d) F ( n ) = ∑ n | d R f ( d )

进行莫比乌斯反演
得到 f(n)=Rn|dμ(dn)F(d) f ( n ) = ∑ n | d R μ ( d n ) F ( d )

F(n)=(RnL1n)N F ( n ) = ( ⌊ R n ⌋ − ⌊ L − 1 n ⌋ ) N

所以原式化为:

f(K)=K|dRμ(dK)(RdL1d)N f ( K ) = ∑ K | d R μ ( d K ) ( ⌊ R d ⌋ − ⌊ L − 1 d ⌋ ) N

T=dK T = d K ,枚举 T T ,得:

f(K)=T=1RKμ(T)(RKTL1KT)N

用一个分块+前缀和搞定
看上去好像做完了的样子
但是若 K=1,R=1e9 K = 1 , R = 1 e 9 就会神奇的发现
前缀和爆了
显然只能开到 1e7 1 e 7 (某些奇怪的OJ只有128MB内存)
但是要求你求到 1e9 1 e 9

所以说需要一些神奇的方法(玄学到自己都不清不白)
因为 d|iμ(d)=[i==1] ∑ d | i μ ( d ) = [ i == 1 ]
可得

i=1nd|iμ(d)=1 ∑ i = 1 n ∑ d | i μ ( d ) = 1

枚举 d d (约数)可得
d=1ni=1ndμ(d)=1

等价于
i=1,d=1id<=ndμ(i) ∑ i = 1 , d = 1 i ∗ d <= n d ∗ μ ( i )

因为 i,j i , j 等价,所以
i=1,d=1id<=niμ(d) ∑ i = 1 , d = 1 i ∗ d <= n i ∗ μ ( d )

等价于
d=1ni=1ndμ(i)=1 ∑ d = 1 n ∑ i = 1 ⌊ n d ⌋ μ ( i ) = 1

等价于枚举 d d (次数)
d=1ni=1ndμ(i)=1


sum(n)=i=1nμ(i) s u m ( n ) = ∑ i = 1 n μ ( i )

那么

d=1nsum(nd)=1 ∑ d = 1 n s u m ( ⌊ n d ⌋ ) = 1

sum(n)=1d=2nsum(nd) s u m ( n ) = 1 − ∑ d = 2 n s u m ( ⌊ n d ⌋ )

这样我们就可以递归的将 sum(n) s u m ( n ) 求出
而且将前 1e6 1 e 6 sum s u m 打表出来,直接代入,可极大加速

详细代码如下:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int INF = 1000000010;
const int N = 1000005;
const int mod = 1000000007;
unordered_map <int,long long> mu_sum;
int n,k,l,r;
ll pri[N],tot,mu[N];
ll s[N];
bool mark[N];
void get() {
    mu[1]=1;
    for(int i=2;i<=1000000;++i) {
        if(!mark[i]) {
            pri[++tot]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=tot && i*pri[j]<=1000000;++j) {
            mark[i*pri[j]]=1;
            if(i%pri[j]==0) {break;}
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=1000000;++i) {
        mu[i]+=mu[i-1];
        s[i]=s[i-1]+mu[i];
    }
}
ll qpow(ll a,ll b) {
    ll ans=1;
    while(b) {
        if(b&1) {
            ans*=a;
            ans%=mod;
        }
        b>>=1;
        a*=a;
        a%=mod;
    }
    return ans;
}
ll get_mu(int x) {
    if(x<=1000000) return mu[x];
    if(mu_sum.find(x)!=mu_sum.end())
        return mu_sum[x];
    int pos;
    ll ans=1;
    for(int i=1;i<=x;i=pos+1) {
        pos=x/(x/i);
        if(x/i-1) ans-=(get_mu(pos)-get_mu(i-1))*(x/i-1);
    }
    return mu_sum[x]=ans;
}
int main() {
    scanf("%d%d%d%d",&n,&k,&l,&r);
    r/=k;l=(l-1)/k;
    get();
    ll ans=0;
    int pos=0;
    for(int i=1;i<=r;i=pos+1) {
        pos=min(r/(r/i),l/i?l/(l/i):INF);
        ans+=(get_mu(pos)-get_mu(i-1))*qpow((r/i)-(l/i),n);
        ans%=mod;
    }
    printf("%lld\n",(ans+mod)%mod);
    return 0;
}
  • 1
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

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

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

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

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

抵扣说明:

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

余额充值