【JZOJ5054】【GDOI2017模拟二试4.12】统计

Description

给定n,k,求满足一下条件的整数数组a[]的数量：
1.a[]中共有k个元素;
2. a[i] ∈ [1,n];
3. ∀i∈[1,k),a[i]≤a[i+1];
4、gcd(a1,a2…ak)=1
答案可能很大，请mod(109+7)后输出

Data Constraint

对所有数据，我们有：
T≤5
n≤10^9,k≤1000

Solution

这是道反演题，而我在考场上只打出了前35分和i=2的杜教筛。
我们设出f[i]表示gcd为i的数量，g[i]表示gcd为i的倍数的数量。

g[d]=Ck⌊n/d⌋+k−1=∑d|kf[k]

f[d]=∑i=1n/dg[i∗d]∗μ(i),f[1]=∑i=1ng[i]∗μ(i)

我们发现g只有N−−√种取值，所以果断上个杜教筛。最后提醒一下：由于组合数计算，直接计算的话是O(N−−√∗K∗T)，明显爆炸,所以我们可以预处理出前10^7的阶乘以及对应的逆元，这样在小于107时可以O(1)计算，而大于10^7的数又十分少。

Code

#include<iostream>#include<cmath>#include<cstring>#include<cstdio>#include<algorithm>#define ll long longusing namespace std;const int maxn=1e7,mo=1e9+7;ll c[maxn+5],c1[maxn+5];int p[maxn+5],bz[maxn+5],d[maxn],h[maxn],g[maxn];ll n,m,i,t,j,k,l,x,y,z,ans;int hash(int x){    int t=x%maxn;    while (h[t]!=x && h[t]) t=(t+1)%maxn;    return t;}ll dg(ll n){    if (n<=maxn) return p[n];    ll x=1,t,k=hash(n),i=2;    if (h[k]) return g[k];h[k]=n;    while (i<n){        t=n/(n/i);        x-=(t-i+1)*dg(n/i)%mo;i=t+1;    }    x%=mo;    g[k]=x;return x;}ll mi(ll x,ll y){    if (y==1) return x;    ll t=mi(x,y/2);    if (y%2)return t*t%mo*x%mo;return t*t%mo;}int main(){    freopen("count.in","r",stdin);freopen("count.out","w",stdout);    p[1]=1;    for (i=2;i<=maxn;i++){        if (!bz[i]) d[++d[0]]=i,p[i]=-1;        for (j=1;j<=d[0];j++){            if (i*d[j]>maxn) break;            bz[i*d[j]]=1;            if (i%d[j]) p[i*d[j]]=-p[i];            else break;        }    }c[1]=1;    for (i=2;i<=maxn;i++)        p[i]=(p[i-1]+p[i])%mo,c[i]=c[i-1]*i%mo;    c1[maxn]=mi(c[maxn],mo-2);    for (i=maxn-1;i>=0;i--)        c1[i]=c1[i+1]*(i+1)%mo;    scanf("%d",&l);    while (l){        scanf("%d%d",&n,&m);        i=1;ans=0;        while (i<=n){            t=n/(n/i);x=c1[m];            if (n/i+m-1<=maxn) x=x*c[n/i+m-1]%mo*c1[n/i-1]%mo;            else{                for (j=n/i;j<=n/i+m-1;j++) x=x*j%mo;            }            ans=ans+x*(dg(t)-dg(i-1))%mo;i=t+1;        }        ans%=mo;ans=(ans+mo)%mo;        printf("%lld\n",ans);        l--;    }}