bzoj 3930: [CQOI2015]选数 （杜教筛+反演）

题目描述

传送门

题目大意：从区间[L,H]（L和H为整数）中选取N个整数，总共有(H-L+1)^N种方案。求从区间[L,H]中选取N个整数，使其最大公约数刚好为K的选取方案有多少个。

题解

$f(d)$表示gcd为d的方案有多少个。
$F(d)$表示d|gcd的方案有多少个。
$F(d)=({\lfloor {r\over d} \rfloor}-{\lfloor {{l-1}\over d} \rfloor})^N$
那么根据莫比乌斯反演的公式$f(d)=\sum\limits_{d|n} \mu({n \over d}) *F(n)$
转换一下枚举的方式

$$f(d)=\sum\limits_{i=1}^{\lfloor {r\over d} \rfloor} \mu(i)*({\lfloor {r\over d*i} \rfloor}-{\lfloor {{l-1}\over d*i} \rfloor})^N$$
这个式子只要预处理前缀和就可以在 $O(\sqrt {\lfloor {r\over d} \rfloor}+\sqrt{\lfloor {{l-1}\over d} \rfloor})$ 的时间内求解。
由于r的范围很大，所有不能直接进行前缀和的预处理，我们需要用到杜教筛，因为杜教筛会对计算过的值进行记录，所有避免了重复计算，而且我们可以用线性筛处理 $n^{2\over 3}$的 $\mu$,所以可以得到非常优越的时间复杂度。

代码

#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
#include<map>
#define LL long long 
#define N 1000003
#define p 1000000007
using namespace std;
int t,n,m,k,l,r,n1,prime[N],pd[N],mu[N];
LL ans;
map<int,LL> sum;
void init()
{
    mu[1]=1;
    for (int i=2;i<=n1;i++) {
        if (!pd[i]) {
            prime[++prime[0]]=i;
            mu[i]=-1;
        }
        for (int j=1;j<=prime[0];j++) {
            if (prime[j]*i>n1) break;
            pd[prime[j]*i]=1;
            if (i%prime[j]==0) {
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
    }
    for (int i=1;i<=n1;i++) mu[i]=mu[i-1]+mu[i];
}
LL get_sum(int n)
{
    if (n<=n1) return mu[n];
    if (sum[n]) return sum[n];
    LL tmp=1;
    for (int i=2,j;i<=n;i=j+1) {
        if (n/i) j=n/(n/i);
        else j=n;
        tmp-=(LL)(j-i+1)*get_sum(n/i)%p;
        tmp%=p;
    }
    return sum[n]=tmp;
}
LL quickpow(LL num,int x)
{
    LL base=(num%p+p)%p; LL ans=1;
    while (x) {
        if (x&1) ans=ans*base%p;
        x>>=1;
        base=base*base%p;
    }
    return ans;
}
int main()
{
    scanf("%d%d%d%d",&t,&k,&l,&r);
    n=(l-1)/k; m=r/k;
    n1=ceil(pow(m,2.0/3.0)); 
    init(); 
    for (int i=1,j;i<=m;i=j+1) {
     if (n/i) j=min(n/(n/i),m/(m/i));
     else j=m/(m/i);
     //cout<<j<<" "<<get_sum(j)<<" "<<get_sum(i-1)<<" "<<quickpow(m/i-n/i,k)<<endl;
     ans=(ans+(get_sum(j)-get_sum(i-1))*quickpow(m/i-n/i,t)%p)%p;
    }
    printf("%lld\n",(ans%p+p)%p);
}