Portal

Description

給出\(n,k,L,R(\leq10^9)\)，求從\([L,R]\)中選出\(n\)個可相同有順序的數使得其gcd為\(k\)的方案數。

Solution

記\(f(x)\)表示gcd為\(x\)時的方案數，那麽我們要求的就是\(f(k)\)。設\(F(x)=\sum_{x|d}f(d)\)表示gcd為\(x\)的倍數時的方案數，即\(F(x)=(?\dfrac{R}{x}?-?\dfrac{L-1}{x}?)^n\)。於是我們得到
\[\begin{align*} f(k) &= \sum_{k|d}\mu(\frac{k}{d})F(d) \&= \sum_{i=1}^{+∞}\mu(i)(?\frac{R}{ik}?-?\frac{L-1}{ik}?)^n \&= \sum_{i=1}^{+∞}\mu(i)(?\frac{?\frac{R}{k}?}{i}?-?\frac{?\frac{L-1}{k}?}{i}?)^n \end{align*}\]

Code

//「CQOI2015」選數
#include <cstdio>
#include <map>
typedef long long lint;
using std::map;
inline int min(int x,int y) {return x<y?x:y;}
const int P=1e9+7;
const int N=2e6+10;
int n0;
int prCnt,pr[N]; bool prNot[N];
int mu[N],sum0[N];
void init(int n)
{
    mu[1
 
]=1;
    for(int i=2;i<=n;i++)
    {
        if(!prNot[i]) pr[++prCnt]=i,mu[i]=-1;
        for(int j=1;j<=prCnt;j++)
        {
            int x=i*pr[j];
            if(x>n) break;
            prNot[x]=true;
            if(i%pr[j]) mu[x]=-mu[i];
            else break; 
        }
    }
    for
 
(int i=1;i<=n;i++) sum0[i]=sum0[i-1]+mu[i];
}
map<int,int> S;
map<int,bool> get;
int sum(int x)
{
    if(x<=n0) return sum0[x];
    if(get[x]) return S[x];
    int r=1;
    for(int L=2,R;L<=x;L=R+1)
    {
        int v=x/L; R=x/v;
        r=(r-1LL*(R-L+1)*sum(v)%P+P)%P;
    }
    get[x]=true; return S[x]=r;
}
int pow(int x,int y)
{
    int r=1,t=x;
    for(int i=y;i;i>>=1,t=(1LL*t*t)%P) if(i&1) r=(1LL*r*t)%P;
    return r;
}
int main()
{
    int n,k,a,b; init(n0=2e6);
    scanf("%d%d%d%d",&n,&k,&a,&b);
    a=(a-1)/k,b=b/k;
    lint ans=0;
    for(int L=1,R;L<=b;L=R+1)
    {
        int v1=b/L,v2=a/L; R=v2?min(b/v1,a/v2):b/v1;
        ans=(ans+1LL*(sum(R)-sum(L-1)+P)*pow(v1-v2,n)%P)%P;
    }
    printf("%lld\n",ans);
    return 0;   
}

P.S.

因為\(\mu(x)\)的部分和有可能是負數，所以運算之前要加一個\(10^9+7\)變成正數...

LibreOJ2095 - 「CQOI2015」選數