【杜教筛】【GDOI2018Day1模拟4.17】呼吸决定

题意

求$\sum_{i=1}^n μ(i)*i^m$
m<=2e5
n<=1e9

题解

设$pi= μ(i)*i^m |si=\sum_{i=1}^n pi$，现在尝试让p卷一个g（狄利克雷卷积）

$\sum_{i=1}^n (s*g)(i)=\sum_{i=1}^n \sum_{d|i} p(\frac{i}{d})g(d)$

$=\sum_{d=1}^n g(d) \sum_{d|i}p(i/d)$
$=\sum_{d=1}^n g(d) \sum_{i=1}^{\frac{n}{d}}p(i)$
$=\sum_{d=1}^n g(d) *s(\frac{n}{d})$

我们现在要求的是s(n)，那么考虑一下拆一下式子

$\sum_{d=1}^n g(d) *s(\frac{n}{d})=g(1)*s(n)-\sum_{d=2}^n g(d) *s(\frac{n}{d})$

那么如果我们的$\sum_{d=1}^n g(d) *s(\frac{n}{d})$是一个可以快速计算的东西（比如说常数之类的），那么我们现在就可以先预处理出前面的一部分s，然后对后面的递归处理

如果一段g的计算是O(1)的那么复杂度在阈值取$n^{\frac{2}{3}}$的时候复杂度是$n^{\frac{2}{3}}$

比如我们如果要求的是μ的前缀和，那么我们直接卷一个1，上式左部分答案就是1，一段g的和也可以O（1）求

然后这题我们卷的是 $i^m$，所以g是要用拉格朗日插值求的

但是如果我们直接插值的话复杂度就变成了 $n^{\frac{2}{3}}*m$

但是事实上我们可以把所有要用的连续一段的d的和预处理出来
有一个很有用的东西
$\frac{\frac{n}{a}}{b}=\frac{n}{ab}$（下取整），所以说我们完全可以把每一个有用的比阈值大的d求出来，这样的d的个数显然是n/阈值级别的

贴代码

#include<iostream>
#include<algorithm>
#include<cstdio>
#include<cstring>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
#define fo1(i,b,a) for(i=b;i>=a;i--)
#define ll long long

using namespace std;

const int md=998244353,mt=4e6;

int f[mt+5],st[mt+5],s[mt+5],miu[mt+5];
int pre[200005],suf[200005],yu[200005],Inv[mt+5];
int sb[200005];
bool bz[mt+5],bb;
int i,j,k,l,m,n,x,y,tp,mi,ans,be,ed;

int quickmi(int x,int y){
    int t1=1;
    while (y){
        if ((y & 1)==1) t1=1ll*t1*x%md;
        x=1ll*x*x%md;
        y=y/2;
    }
    return t1;
}
void ge_yu(){
    Inv[1]=1;
    for (int i=2;i<mt;i++) Inv[i]=(ll)Inv[md%i]*(md-md/i)%md;
    Inv[0]=1;
    for (int i=1;i<mt;i++) Inv[i]=(ll)Inv[i-1]*Inv[i]%md;
    for (int i=1,sig=(m&1)?1:-1;i<=m+2;i++,sig=-sig) yu[i]=(ll)sig*Inv[m+2-i]*Inv[i-1]%md;
}
ll power(int x){
    if (x<=mt) return f[x];
    if (bb) return sb[n/x];
    pre[0]=1; suf[m+3]=1;
    fo(i,1,m+2) pre[i]=1ll*pre[i-1]*(x-i)%md;
    fo1(i,m+2,1) suf[i]=1ll*suf[i+1]*(x-i)%md;
    ll ans=1ll*f[1]*suf[2]%md*yu[1]%md;
    fo(i,2,m+2){
        ans=(ans+1ll*pre[i-1]*suf[i+1]%md*f[i]%md*yu[i])%md;
    }
    return (md+ans)%md;
}
int solve(int x){
    int pic=0,be=2,ed;
    if (x<=mt) return s[x]; else
    if (bz[n/x]==true) return st[n/x]; else{
        while (be<=x){
            ed=x/(x/be);
            pic=(pic+(power(ed)-power(be-1)+md)%md*solve(x/be)%md)%md;
            be=ed+1;
        }
    }
    st[n/x]=(md+1-pic)%md; bz[n/x]=true;
    return st[n/x];
}
int main(){
    freopen("t3.in","r",stdin);
//  freopen("calc.out","w",stdout);
    scanf("%d%d",&n,&m);
    mi=min(n,mt);
    ge_yu();
    miu[1]=1; s[1]=1; f[1]=1;
    fo(i,2,mi){
        if (bz[i]==false){
            miu[i]=-1; st[++tp]=i;
        }
        be=quickmi(i,m);
        s[i]=(s[i-1]+1ll*miu[i]*be%md+md)%md;
        f[i]=(f[i-1]+be)%md;
        fo(j,1,tp){
            if (i*st[j]>mi) break;
            bz[i*st[j]]=true;
            if (i%st[j]==0){
                miu[i*st[j]]=0; break;
            }
            miu[i*st[j]]=-miu[i];
        }
    }
    memset(st,0,sizeof(st));
    memset(bz,false,sizeof(bz));
    be=1;
    while (be<=n){
        ed=n/(n/be);
        if (ed<=200000) sb[ed]=power(n/be); else break;
        be=ed+1;
    }
    bb=true;
    ans=solve(n)%md;
    printf("%d\n",ans);
    return 0;
}