[JZOJ5054]统计

题目大意

给定$n,k$，求满足以下条件的整数数组$a$的数量：
$\bullet\ $数组$a$的大小为$k$
$\bullet\ \forall i\in[1,k],a_i\in[1,n]$
$\bullet\ \forall i\in[1,k),a_i\le a_{i+1}$
$\bullet\ \gcd_{i=1}^k\{a_i\}=1$
一个测试点$T$组数据。答案对$10^9+7$取模。

$T\le 5,1\le n\le10^9,1\le k\le10^3$

---

题目分析

一眼莫比乌斯反演，于是设$F(x)$表示数值在$[1,x]$范围内，长度为$k$的不降序列个数。
根据莫比乌斯反演，显然答案是

$$\sum_{i=1}^n\mu(i)F\left(\lfloor\frac n i\rfloor\right)$$
这个对$\lfloor\frac n i\rfloor$分块，使用杜教筛来求$\mu(x)$的前缀和就可以了。
至于$F(x)$，由挡板原理得$F(x)={k+x-1\choose x-1}$，这个怎么算呢？我预处理$10^7$以内的阶乘以及阶乘逆元，对于该范围内的组合数我直接算，否则我暴力计算（反正一个长度为$k$的连乘）。
至于时间复杂度，低于$\mathrm O(n^{2\over 3}k)$就是了。

---

代码实现

#include <iostream>
#include <cstdio>

using namespace std;

const int P=1000000007;
const int K=1000;
const int N=10000000;
const int V=60000;
const int A=32000;

int pri[N+5],mu[N+5],sum[N+5],f[N+5],fact[N+5],invf[N+5];
int vis[A],S[A];
int T,n,k,cnt;

int quick_power(int x,int y)
{
    int ret=1;
    for (;y;y>>=1,x=1ll*x*x%P) if (y&1) ret=1ll*ret*x%P;
    return ret;
}

inline int getid(int x){return n/x;}

void pre()
{
    fact[0]=1;
    for (int i=1;i<=N;++i) fact[i]=1ll*fact[i-1]*i%P;
    invf[N]=quick_power(fact[N],P-2);
    for (int i=N;i>=1;--i) invf[i-1]=1ll*invf[i]*i%P;
    sum[1]=mu[1]=1,f[1]=1;
    for (int i=2;i<=N;++i)
    {
        if (!f[i]) pri[++pri[0]]=f[i]=i,mu[i]=-1;
        for (int j=1,x;j<=pri[0];++j)
        {
            if (1ll*pri[j]*i>N) break;
            f[x=pri[j]*i]=pri[j];
            mu[x]=f[x]==f[i]?0:-mu[i];
            if (!(i%pri[j])) break;
        }
        sum[i]=sum[i-1]+mu[i];
    }
}

int sieve(int n)
{
    if (n<=N) return sum[n];
    int id=getid(n);
    if (vis[id]==cnt) return S[id];
    vis[id]=cnt;
    int ret=0;
    for (int st=2,en,x;st<=n;st=en+1) x=n/st,en=n/x,(ret+=1ll*sieve(x)*(en-st+1)%P)%=P;
    return S[id]=ret=(P-ret+1)%P;
}

inline int F(int x)
{
    if (k+x-1<=N) return 1ll*fact[k+x-1]*invf[x-1]%P*invf[k]%P;
    int ret=invf[k];
    for (int i=x;i<=k+x-1;++i) ret=1ll*ret*i%P;
    return ret;
}

int calc()
{
    int ret=0;
    for (int st=1,en,x;st<=n;st=en+1)
    {
        x=n/st,en=n/x;
        (ret+=1ll*(((++cnt,sieve(en))-(++cnt,sieve(st-1))+P)%P)*F(x)%P)%=P;
    }
    return ret;
}

int main()
{
    freopen("count.in","r",stdin),freopen("count.out","w",stdout);
    for (pre(),scanf("%d",&T);T--;printf("%d\n",calc())) scanf("%d%d",&n,&k);
    fclose(stdin),fclose(stdout);
    return 0;
}