【校内模拟】calc（多项式快速插值）（多项式多点求值）（杜教筛）

其实感慨还是挺多的，从某种意义上来说这道题算是考场AC了吧，至于那个卡空间的10pts，感觉数组实现会过的轻松一点，vector的话空间常数略大。。。但是我是一个vector多项式党

从来没有想过自己能够在考场上写出超过5.5KB的正确代码并AC，而且是精简不下去的那种。对自己来说算是一个里程碑吧，虽然说这次的基本上都是板子，但是有一堆板子是考场现推的，比如邻项快速插值的加速写法。

题解：

求$S(n)={\sum }_{i=1}^n\mu (i)i^m$，其中$n\le 1e9,m\le 2e5$。

设$f(n)=\mu (n)n^m$，则令$g(i)=i^m$，有$$\begin{gathered} (f\ast g)=ϵ \\ g(1)S(n)=[n==1]-\sum _{i=2}^{n}g(i)S(\lfloor \frac{n}{i}\rfloor ) \end{gathered}$$

现在要求$g$的前缀和，显然是关于$i$的$m+1$次多项式，快速插值之后把要求值的点搞出来多点求值即可。

也不难写，但是如果用任意点值的通用插值的话，会慢一倍，亲测，但是其他人写出来好像是慢若干倍。。。

不难想，不难写，老板来一份多项式全家桶

---

代码：

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const

using std::cerr;
using std::cout;
#define fi first
#define se second

cs int mod=998244353;
inline int add(int a,int b){a+=b-mod;return a+(a>>31&mod);}
inline int dec(int a,int b){a-=b;return a+(a>>31&mod);}
inline int mul(int a,int b){ll r=(ll)a*b;return r>=mod?r%mod:r;}
inline int power(int a,int b,int res=1){
	for(;b;b>>=1,a=mul(a,a))(b&1)&&(res=mul(res,a));
	return res;
}
inline void Inc(int &a,int b){a+=b-mod;a+=a>>31&mod;}
inline void Dec(int &a,int b){a-=b;a+=a>>31&mod;}
inline void Mul(int &a,int b){a=mul(a,b);}

typedef std::vector<int> Poly;

cs int bit=19,SIZE=1<<bit|1;

int r[SIZE],*w[bit+1];
inline void init_NTT(){
	for(int re i=1;i<=bit;++i)w[i]=new int[1<<i-1];
	int wn=power(3,mod-1>>bit);w[bit][0]=1;
	for(int re i=1;i<(1<<bit-1);++i)w[bit][i]=mul(w[bit][i-1],wn);
	for(int re i=bit-1;i;--i)
	for(int re j=0;j<(1<<i-1);++j)w[i][j]=w[i+1][j<<1];
}
inline void NTT(int *A,int len,int typ){
	for(int re i=1;i<len;++i)if(i<r[i])std::swap(A[i],A[r[i]]);
	for(int re i=1,d=1;i<len;i<<=1,++d)
	for(int re j=0;j<len;j+=i<<1){
		if(i<8){
			for(int re k=0;k<i;++k){
				int &t1=A[j+k],&t2=A[i+j+k],t=mul(t2,w[d][k]);
				t2=dec(t1,t),Inc(t1,t);
			}
		}
		else {
			for(int re k=0;k<i;k+=8){
#define work(p)\
{	\
	int &t1=A[j+k+p],&t2=A[i+j+k+p],t=mul(t2,w[d][k+p]); \
	t2=dec(t1,t),Inc(t1,t);	\
}
				work(0);work(1);work(2);work(3);
				work(4);work(5);work(6);work(7);
#undef work 
			}
		}
	}
	if(typ==-1){
		std::reverse(A+1,A+len);
		for(int re i=0,inv=power(len,mod-2);i<len;++i)Mul(A[i],inv);
	}
}
inline void NTT(Poly &A,int len,int typ){NTT(&A[0],len,typ);}
inline void init_rev(int l){for(int re i=1;i<l;++i)r[i]=r[i>>1]>>1|((i&1)?l>>1:0);}

inline Poly operator+(cs Poly &a,cs Poly &b){
	Poly c=a;if(a.size()<b.size())c.resize(b.size());
	for(int re i=0,li=b.size();i<li;++i)Inc(c[i],b[i]);
	return c;
}

inline Poly operator-(cs Poly &a,cs Poly &b){
	Poly c=a;if(a.size()<b.size())c.resize(b.size());
	for(int re i=0,li=b.size();i<li;++i)Dec(c[i],b[i]);
	return c;
}

inline Poly operator*(Poly a,Poly b){
	int n=a.size(),m=b.size(),deg=n+m-1,l=1;
	while(l<deg)l<<=1;init_rev(l);
	a.resize(l),NTT(a,l,1);
	b.resize(l),NTT(b,l,1);
	for(int re i=0;i<l;++i)Mul(a[i],b[i]);
	NTT(a,l,-1);return Poly(a.begin(),a.begin()+deg);
}

inline Poly Inv(cs Poly &a,int lim){
	int n=a.size();Poly c,b(1,power(a[0],mod-2));
	for(int re l=4;(l>>2)<lim;l<<=1){
		init_rev(l);c.resize(l>>1);
		for(int re i=0;i<(l>>1);++i)c[i]=i<n?a[i]:0;
		c.resize(l),NTT(c,l,1);
		b.resize(l),NTT(b,l,1);
		for(int re i=0;i<l;++i)Mul(b[i],dec(2,mul(b[i],c[i])));
		NTT(b,l,-1),b.resize(l>>1);
	}b.resize(lim);
	return b;
}

inline Poly operator/(Poly a,Poly b){
	if(a.size()<b.size())return Poly(1,0);
	int deg=a.size()-b.size()+1;
	std::reverse(a.begin(),a.end());
	std::reverse(b.begin(),b.end());
	b=Inv(b,deg);a=a*b;a.resize(deg);
	std::reverse(a.begin(),a.end());
	return a;
}

inline Poly operator%(cs Poly &a,cs Poly &b){
	if(a.size()<b.size())return a;
	Poly c=a-(a/b)*b;
	c.resize(b.size()-1);
	return c;
}

cs int PP=2e5+5;

Poly f[PP<<2|1];

inline void build(int k,int l,int r,cs Poly &v){
	if(l==r){f[k].resize(2);f[k][0]=dec(0,v[l]),f[k][1]=1;return ;}
	int mid=l+r>>1;build(k<<1,l,mid,v);build(k<<1|1,mid+1,r,v);
	f[k]=f[k<<1]*f[k<<1|1];
}

inline void solve(int k,int l,int r,cs Poly &a,Poly &y){
	if(l==r){y[l]=a[0];return ;}int mid=l+r>>1;
	solve(k<<1,l,mid,a%f[k<<1],y);
	solve(k<<1|1,mid+1,r,a%f[k<<1|1],y);
}

inline void multi_point_evaluation(cs Poly &A,cs Poly &x,Poly &y){
	build(1,0,x.size()-1,x);y.resize(x.size());
	solve(1,0,x.size()-1,A%f[1],y);
}

cs int M=1e6,N=M+7,P=1e4+7;

int ifac[PP];

inline std::pair<Poly,Poly> interpolation(int l,int r,cs Poly &A,cs Poly &coef){
	if(l==r){
		Poly h(2,0);h[0]=dec(0,l),h[1]=1;
		Poly g(1,mul(A[l],coef[l]));
		return std::pair<Poly,Poly>(h,g);
	}
	int mid=l+r>>1;
	std::pair<Poly,Poly> a=interpolation(l,mid,A,coef);
	std::pair<Poly,Poly> b=interpolation(mid+1,r,A,coef);
	return std::pair<Poly,Poly>(a.fi*b.fi,a.fi*b.se+b.fi*a.se);
}

inline Poly interpolation(cs Poly &G){
	int n=G.size();Poly coef;coef.resize(n);ifac[0]=1;
	for(int re i=1;i<n;++i)ifac[i]=mul(ifac[i-1],power(i,mod-2));
	for(int re i=0;i<n;++i)coef[i]=mul(ifac[i],((n-i)&1)?ifac[n-i-1]:mod-ifac[n-i-1]);
	return interpolation(0,n-1,G,coef).se;
}

int n,m;

int p[N],pct;bool mrk[N];
int mu[N],g[N],s[N];

inline void linear_sieves(int lim=M){
	mu[1]=g[1]=s[1]=1;
	for(int re i=2;i<=lim;++i){
		if(!mrk[i])p[++pct]=i,g[i]=power(i,m),mu[i]=mod-1;
		s[i]=add(s[i-1],mul(g[i],mu[i]));
		for(int re j=1;i*p[j]<=lim;++j){
			mrk[i*p[j]]=true;
			g[i*p[j]]=mul(g[i],g[p[j]]);
			if(i%p[j]==0)break;
			mu[i*p[j]]=mu[i]?mod-mu[i]:0;
		}
	}
	for(int re i=2;i<=lim;++i)Inc(g[i],g[i-1]);
}

int ss[P],gg[P],vis[P];

inline void work(int n){
	if(n<=M)return ;Poly x,tx,y;
	for(int i=0;i<=m+1;++i)y.push_back(g[i]);
	Poly F=interpolation(y);
	for(int i=1,j;i<=n;i=j+1){
		j=n/(n/i);if(j>M)x.push_back(j%mod),tx.push_back(j);
	}
	multi_point_evaluation(F,x,y);
	for(int re i=0;i<x.size();++i)gg[::n/tx[i]]=y[i];
}

inline int G(int n){
	if(n<=M)return g[n];
	else return gg[::n/n];
}

inline int S(int n){
	if(n<=M)return s[n];
	if(vis[::n/n])return ss[::n/n];
	vis[::n/n]=true;int ans=1;
	for(int re i=2,j;i<=n;i=j+1){
		j=n/(n/i);
		Dec(ans,mul(S(n/i),dec(G(j),G(i-1))));
	}
	return ss[::n/n]=ans;
}

#ifdef zxyoi
#undef zxyoi
#endif

signed main(){
#ifdef zxyoi
	freopen("calc.in","r",stdin);
#else
#ifndef ONLINE_JUDGE
	freopen("calc.in","r",stdin);freopen("calc.out","w",stdout);
#endif
#endif
	scanf("%d%d",&n,&m);
	linear_sieves();init_NTT();
	work(n);cout<<S(n);
	return 0;
}