bzoj3930 [CQOI2015]选数

最新推荐文章于 2018-12-01 08:48:00 发布

原创最新推荐文章于 2018-12-01 08:48:00 发布 · 202 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#bzoj

c++ 同时被 3 个专栏收录

1073 篇文章

订阅专栏

莫比乌斯反演

20 篇文章

订阅专栏

杜教筛

9 篇文章

订阅专栏

本文探讨了一个数学问题：如何计算从指定区间内选取特定数量整数时，这些整数的最大公约数恰好为K的不同方案的数量。通过定义两个函数f(n)和g(n)，分别表示最大公约数等于n和最大公约数能整除n的方案数，并利用反演公式和杜教筛法进行高效计算。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

Description

我们知道，从区间[L,H]（L和H为整数）中选取N个整数，总共有(H-L+1)^N种方案。小z很好奇这样选出的数的最大公约数的规律，他决定对每种方案选出的N个整数都求一次最大公约数，以便进一步研究。然而他很快发现工作量太大了，于是向你寻求帮助。你的任务很简单，小z会告诉你一个整数K，你需要回答他最大公约数刚好为K的选取方案有多少个。由于方案数较大，你只需要输出其除以1000000007的余数即可。

对于100%的数据，1≤N,K≤10^9，1≤L≤H≤10^9，H-L≤10^5

Solution

可能是我之前做过太多暴力推的柿子了，反而忘了反演的本质是啥

设 $f(n)$ 表示 $gcd=n$ 的方案数， $g(n)$ 表示 $gcd|n$ 的方案数，显然有

g (n) = \sum n | d f (d)

$g(n)=\sum_{n|d}f(d)$ 反演可得

f (n) = \sum n | d g (d) μ (d n)

$f(n)=\sum_{n|d}g(d)\mu(\frac{d}{n})$
容易发现g(d)的求法非常简单，即

g (d) = (⌊ H d ⌋ - ⌊ L - 1 d ⌋) n

$g(d)=(\lfloor\frac{H}{d}\rfloor-\lfloor\frac{L-1}{d}\rfloor)^n$
那么有

f (k) = \sum k | d g (d) μ (d k) = \sum d = 1 ⌊ H k ⌋ μ (d) (⌊ H k d ⌋ - ⌊ L - 1 k d ⌋) n

$f(k)=\sum_{k|d}g(d)\mu(\frac{d}{k})=\sum_{d=1}^{\lfloor\frac{H}{k}\rfloor}\mu(d)(\lfloor\frac{H}{kd}\rfloor-\lfloor\frac{L-1}{kd}\rfloor)^n$
看到下取整就可以分块，μ可以杜教筛，注意到

⌊⌊nk⌋d⌋⌊⌊nk⌋d⌋ $\lfloor\frac{\lfloor\frac{n}{k}\rfloor}{d}\rfloor$ 实际上是和

⌊nkd⌋⌊nkd⌋ $\lfloor\frac{n}{kd}\rfloor$ 等价的，因此分块的时候怎么求右端点十分随意

Code

#include <stdio.h>
#include <string.h>
#include <map>
#define rep(i,st,ed) for (int i=st;i<=ed;++i)

typedef long long LL;
const int MOD=1000000007;
const int N=5000000;

std:: map <LL,LL> map;

int prime[N/5];
short mu[N+5];
bool not_prime[N+5];

LL ksm(LL x,LL dep) {
    LL ret=1;
    while (dep) {
        if (dep&1) ret=ret*x%MOD;
        dep/=2; x=x*x%MOD;
    }
    return ret;
}

void pre_work() {
    mu[1]=1;
    rep(i,2,N) {
        if (!not_prime[i]) {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=prime[0]&&i*prime[j]<=N;j++) {
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    rep(i,1,N) mu[i]+=mu[i-1];
}

LL get_mu(LL x) {
    if (x<=N) return mu[x];
    if (map[x]) return map[x];
    LL ret=1;
    for (int i=2,j;i<=x;i=j+1) {
        j=x/(x/i);
        ret=(ret+MOD-(LL)(j-i+1)*get_mu(x/i)%MOD)%MOD;
    }
    return map[x]=ret;
}

void solve(int n,int k,int l,int h) {
    LL ans=0;
    for (int i=1,j;i<=h;i=j+1) {
        j=std:: min((l/i)?(l/(l/i)):MOD,h/(h/i));
        ans=(ans+(get_mu(j)-get_mu(i-1))*ksm(h/i-l/i,n)%MOD)%MOD;
    }
    printf("%lld\n", (ans+MOD)%MOD);
}

int main(void) {
    pre_work();
    int n,k,l,h; scanf("%d%d%d%d",&n,&k,&l,&h);
    h/=k; (l-=1)/=k;
    solve(n,k,l,h);
    return 0;
}