【牛客小白月赛8-J】好的序列

mrsrz 4月 14, 2020

整除分块，拉格朗日插值，杜教筛。

题目大意：

给定 $m,k$，定义 $f(m)$ 表示长度为 $k$ 的、每个数大小不超过 $m$ 的、所有数的 $\gcd$ 为 $1$ 的不同序列个数。求下面的式子的值。

$\sum_{i=1}^n f(i)$

链接

题解：

长度为 $k$ 的，每个数大小不超过 $m$ 的序列个数为 $m^k$。我们考虑枚举 $\gcd$，并容斥得到 $f(m)$ 的值。

$f(m)=\sum_{i=1}^m \mu(i)\left(\lfloor\frac m i\rfloor\right)^k$

于是我们要求的就是：

$\begin{align} \sum_{m=1}^n \sum_{i=1}^m \mu(i)\left(\lfloor\frac m i\rfloor\right)^k&=\sum_{i=1}^n \mu(i) \sum_{m=i}^n\left(\lfloor\frac m i\rfloor\right)^k\\ &=\sum_{i=1}^n \mu(i) \sum_{m=1}^n\left(\lfloor\frac m i\rfloor\right)^k\\ &=\sum_{i=1}^n\mu(i)\left(g(\lfloor\frac {n+1} i\rfloor-1)\cdot i+ ((n+1)\bmod i)\cdot(\lfloor\frac{n+1}i \rfloor)^k\right)\\ &=\sum_{i=1}^n\mu(i)\left(g(\lfloor\frac {n+1} i\rfloor-1)\cdot i+ (n+1-\lfloor\frac{n+1}i\rfloor\cdot i)\cdot(\lfloor\frac{n+1}i \rfloor)^k\right) \end{align}$

其中 $g(n)=\sum\limits_{i=1}^n i^k$

我们把这些项拆开来，分别使用整除分块计算即可。

这里涉及到 $\mu(i)$ 的前缀和以及 $i\cdot mu(i)$ 的前缀和的计算，需要使用杜教筛。

考虑如何求 $g(n)$。~~一个比较无脑的办法是多项式多点求值。~~

预处理连续点值后可以用拉格朗日插值在 $O(k)$ 内插出一个点处的值。

直接在整除分块里做拉格朗日插值的话，将会有 $O(k\sqrt n)$ 的复杂度，不能接受。

我们可以预处理 $S$ 以内的前缀 $k$ 次方和，对于小于等于 $S$ 的直接 $O(1)$ 查询答案，对于大的则使用拉格朗日插值获取答案。

理论上取 $S=\sqrt {nk}$ 可以较好地将这部分的复杂度平衡至 $O(\sqrt{nk})$，实际上当然是能预处理越多越好。

那么时间复杂度为 $O(\sqrt[3]{n^2}+\sqrt{nk})$。

Code：

#include<cstdio>
#include<cstring>
#include<unordered_map>
typedef long long LL;
const int md=998244353,N=1e5+7;
int n,k,X[10000008],fac[N],iv[N],ans,pre[N],suf[N];
int P[10000008],pri[10000001/10],tot,mu[10000008],imu[10000008];
bool vis[10000008];
std::unordered_map<int,int>Mu,IMu;
inline void upd(int&a){a+=a>>31&md;}
inline int pow(int a,int b){
    int ret=1;
    for(;b;b>>=1,a=(LL)a*a%md)if(b&1)ret=(LL)ret*a%md;
    return ret;
}
inline void init(const int n){
    P[1]=mu[1]=1;
    for(int i=2;i<=n;++i){
        if(!vis[i])pri[++tot]=i,P[i]=pow(i,k),mu[i]=-1;
        for(int j=1;j<=tot&&i*pri[j]<=n;++j){
            vis[i*pri[j]]=1,P[i*pri[j]]=(LL)P[i]*P[pri[j]]%md;
            if(i%pri[j]==0)break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;++i)upd(X[i]=P[i]+X[i-1]-md);
}
inline int query(int x){
    if(x<=10000000)return X[x];
    int ret=0;
    for(int i=1;i<=k+3;++i)
    pre[i]=pre[i-1]*(LL)(x-i+md)%md;
    for(int i=k+3;i;--i)
    suf[i]=suf[i+1]*(LL)(x-i+md)%md;
    for(int i=1;i<=k+3;++i){
        int y=X[i]*(LL)pre[i-1]%md*suf[i+1]%md*iv[i-1]%md*iv[k+3-i]%md;
        if((k+3-i)&1)upd(y=-y);
        upd(ret+=y-md);
    }
    return ret;
}
inline int sumd(int x){return(x+1LL)*x/2%md;}
int getmu(int n){
    if(n<=10000000)return mu[n];
    if(Mu.count(n))return Mu[n];
    int&S=Mu[n];S=0;
    for(unsigned l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        S=(S+(r-l+1LL)*getmu(n/l))%md;
    }
    upd(S=1-S);
    return S;
}
int getimu(int n){
    if(n<=10000000)return imu[n];
    if(IMu.count(n))return IMu[n];
    int&S=IMu[n];S=0;
    for(int l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        S=(S+((LL)sumd(r)-sumd(l-1)+md)*getimu(n/l))%md;
    }
    upd(S=1-S);
    return S;
}
inline int solve(int t){
    return(t*(LL)query((n+1)/t-1)+(LL)((n+1)%t)*pow((n+1)/t,k))%md;
}
int main(){
    scanf("%d%d",&n,&k);
    init(10000000);
    pre[0]=suf[k+4]=1;
    for(int i=*fac=1;i<=k+3;++i)fac[i]=(LL)i*fac[i-1]%md;
    iv[k+3]=pow(fac[k+3],md-2);
    for(int i=k+2;~i;--i)iv[i]=(i+1LL)*iv[i+1]%md;
    if(n<=5000000){
        for(int i=1;i<=n;++i)
        ans=(ans+(LL)(mu[i]+md)*solve(i))%md;
        printf("%d\n",ans);
    }else{
        for(int i=1;i<=5000000;++i)
        ans=(ans+(LL)(mu[i]+md)*solve(i))%md;
        ++n;
        for(int i=1;i<=10000000;++i)
        imu[i]=(imu[i-1]+(LL)i*mu[i]+md)%md,mu[i]=(mu[i]+mu[i-1]+md)%md;
        for(int l=5000001,r;l<n;l=r+1){
            r=n/(n/l);
            if(r==n)--r;
            int sum_mu=getmu(r)-getmu(l-1);upd(sum_mu);
            int sum_imu=getimu(r)-getimu(l-1);upd(sum_imu);
            ans=(ans+(LL)sum_imu*query(n/l-1)%md)%md;
            int t=pow(n/l,k),x=0;
            x=(x+(LL)n*t%md*sum_mu)%md;
            x=(x-(LL)(n/l)*t%md*sum_imu%md+md)%md;
            upd(ans+=x-md);
        }
        printf("%d\n",ans);
    }
    return 0;
}