bzoj3930 [CQOI2015]選數

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) f ( n ) 表示 gcd=n g c d = n 的方案數， g(n) g ( n ) 表示 gcd|n g c d | n 的方案數，顯然有

g(n)=∑n|df(d) g ( n ) = ∑ n | d f ( d )

反演可得 f(n)=∑n|dg(d)μ(dn) f ( n ) = ∑ n | d g ( d ) μ ( d n )

容易發現g(d)的求法非常簡單，即 g(d)=(⌊Hd⌋−⌊L−1d⌋)n g ( d ) = ( ⌊ H d ⌋ − ⌊ L − 1 d ⌋ ) n

那麼有 f(k)=∑k|dg(d)μ(dk)=∑d=1⌊Hk⌋μ(d)(⌊Hkd⌋−⌊L−1kd⌋)n f ( k ) = ∑ k | d g ( d ) μ ( d k ) = ∑ d = 1 ⌊ H k ⌋ μ ( d ) ( ⌊ H k d ⌋ − ⌊ L − 1 k d ⌋ ) n

看到下取整就可以分塊，μ可以杜教篩，注意到 ⌊⌊nk⌋d⌋ ⌊ ⌊ n k ⌋ d ⌋ 實際上是和 ⌊nkd⌋ ⌊ n k d ⌋ 等價的，是以分塊的時候怎麼求右端點十分随意

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=;
const int N=;

std:: map <LL,LL> map;

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

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

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

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

void solve(int n,int k,int l,int h) {
    LL ans=;
    for (int i=,j;i<=h;i=j+) {
        j=std:: min((l/i)?(l/(l/i)):MOD,h/(h/i));
        ans=(ans+(get_mu(j)-get_mu(i-))*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-=)/=k;
    solve(n,k,l,h);
    return ;
}