知识点 - 多项式的前缀和

知识点 - 多项式的前缀和

解决问题类型：

$n^k$ 的前缀和

hdu多校第一场 Sequence 貌似n阶前缀和/矩阵快速幂但其实要用组合数。。。

复杂度：

O(n^2)?

讲义

1. 多项式的前缀和

1.1 $n^k$ 的前缀和

首先来看几个众所周知的公式
$$\sum _{i=1}^{n}1=n$$

$$\sum _{i=1}^{n}i=\frac{n(n+1)}{2}$$

$$\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$

用数学归纳法很容易验证上述公式的正确性，但是对于任何给定的非负整数 $k$ ，如何求出 ${\sum }_{i=1}^n{i}^k$ 呢？下面给出一种利用杨辉三角的计算方法。
$$\begin{array}{ccccc}1& 1& 1& 1& 1\\ 1& 2& 3& 4& 5\\ 1& 3& 6& 10& 15\\ 1& 4& 10& 20& 35\\ 1& 5& 15& 35& 70\end{array}$$

杨辉三角的自然数形式

$$\begin{array}{ccccc}{C}_0^0& C_1^1& C_2^2& C_3^3& C_4^4\\ C_1^0& C_2^1& C_3^2& C_4^3& C_5^4\\ C_2^0& C_3^1& C_4^2& C_5^3& C_6^4\\ C_3^0& C_4^1& C_5^2& C_6^3& C_7^4\\ C_4^0& C_5^1& C_6^2& C_7^3& C_8^4\end{array}$$

杨辉三角的组合数形式

不难发现，第$i$列的前$n$个数字之和刚好等于第$i+1$列第$n$个数字，即
$$C_k^k+C_{k+1}^k+\cdots +C_{k+n}^k=C_{k+n+1}^{k+1}$$
证 根据杨辉恒等式 $C_n^k=C_{n-1}^{k-1}+C_{n-1}^k$ ，$左边=C_k^k+C_{k+1}^k+\cdots +C_{k+n}^k=C_{k+1}^{k+1}+C_{k+1}^k+\cdots +C_{k+n}^k=C_{k+2}^{k+1}+\cdots +C_{k+n}^k=\cdots =C_{k+n+1}^{k+1}=右边$

换言之，除第一列外，每一列都是前一列的“前缀和”，而每一列的数字，都是 $C_n^k(k=列数-1)$ 的形式，其本质就是 $n$ 的多项式，例如：
$$第一列:C_n^0=1$$

$$第二列：C_n^1=\sum _{i=0}^{n-1}{C}_i^0=\sum _{i=0}^{n-1}1=n$$

$$第三列：C_n^2=\sum _{i=1}^{n-1}{C}_i^1=\sum _{i=1}^{n-1}i=\frac{n(n-1)}{2}$$

$$第四列：C_n^3=\sum _{i=2}^{n-1}{C}_i^2=\sum _{i=2}^{n-1}\frac{i(i-1)}{2}=\frac{n(n-1)(n-2)}{6}$$

根据第三列，得到
$$\sum _{i=1}^{n}i=\sum _{i=1}^{n-1}i+n=\frac{n(n+1)}{2}$$
根据第四列，得到
$$\sum _{i=2}^{n-1}\frac{i(i-1)}{2}=\frac{1}{2}\sum _{i=2}^{n-1}{i}^2-\frac{1}{2}\sum _{i=2}^{n-1}i=\frac{1}{2}\left(\sum _{i=1}^n{i}^2-1-n^2\right)-\frac{(n+1)(n-2)}{4}=\frac{n(n-1)(n-2)}{6}$$
进而得到
$$\sum _{i=1}^n{i}^2=\frac{n(n+1)(2n+1)}{6}$$
至于更高次的求和，以此类推即可，上述计算过程比较机械化，因此不难用计算机来实现

#include <bits/stdc++.h>
#define rep(i,a,b) for(ll i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll gcd(ll a, ll b)
{
	return b ? gcd(b, a%b) : a;
}
class frac
{
	public:
	ll x,y;
	frac(){}frac(ll x,ll y):x(x),y(y){}
	bool operator < (const frac &b)const{return x*b.y<y*b.x;}
	bool operator > (const frac &b)const{return x*b.y>y*b.x;}
	bool operator ==(const frac &b)const{return x*b.y==y*b.x;}
	frac operator + (const frac &b)const{ll d=gcd(x*b.y+b.x*y,y*b.y);return frac((x*b.y+b.x*y)/d,(y*b.y)/d);}
	frac operator - (const frac &b)const{ll d=gcd(x*b.y-b.x*y,y*b.y);return frac((x*b.y-b.x*y)/d,(y*b.y)/d);}
	frac operator * (const frac &b)const{ll d=gcd(x*b.x,y*b.y);return frac((x*b.x)/d,(y*b.y)/d);}
	frac operator / (const frac &b)const{ll d=gcd(x*b.y,b.x*y);return frac((x*b.y)/d,(b.x*y)/d);}
	frac operator * (ll b)const{ll d=gcd(x*b,y);return frac((x*b)/d,(y)/d);}
	frac operator / (ll b)const{ll d=gcd(x,y*b);return frac((x)/d,(y*b)/d);}
	frac operator = (ll b){*this=frac(b,1);return *this;}
};
ostream &operator <<(ostream &out,const frac &a)
{
	if(a.y==1)out<<a.x;
	else out<<a.x<<"/"<<a.y;
	return out;
}
typedef frac type;
bool isZero(type x){
	return x.x==0;
}
class Poly{
public:
	vector<type>a={frac(0,1)};
	Poly(){}
	Poly(vector<type> b):a(b){}
	ll n(){
		return a.size()-1;
	}
	Poly operator = (type b){
		this->a.resize(1);
		this->a[0]=b;
		return *this;
	}
	Poly operator = (vector<type> b){
		this->a=b;
		return *this;
	}
	friend ostream &operator << (ostream &o,const Poly &f){
		for(int i=f.a.size()-1;~i;i--){
			if(!i)cout<<"("<<f.a[i]<<")";
			else cout<<"("<<f.a[i]<<")"<<"x^"<<i<<"+";
		}
		cout<<endl;
	}
	type coef(int i){
		if(i>=a.size() || i<0)return frac(0,1);
		return a[i];
	}
	type& operator [] (int i){
		if(i>=a.size() || i<0)cout<<" Warning: Index out of range\n";
		return a[i];
	}
	type operator () (type x){
		type ans;
		ans=0;
		for(int i=n();~i;i--)ans=ans*x+a[i];
		return ans;
	}
	Poly operator () (Poly x){
		Poly ans,t;
		for(int i=n();~i;i--){
			t=Poly((vector<type>){a[i]});
			ans=ans*x+t;
		}
		return ans;
	}
	Poly operator + (Poly &b){
		Poly c;
		c.a.resize(max(a.size(),b.a.size()));
		for(int i=0;i<c.a.size();i++)c.a[i]=coef(i)+b.coef(i);
		while(c.a.size()>1 && isZero(*(c.a.end()-1)))c.a.erase(c.a.end()-1);
		return c;
	}
	Poly operator - (Poly &b){
		Poly c;
		c.a.resize(max(a.size(),b.a.size()));
		for(int i=0;i<c.a.size();i++)c.a[i]=coef(i)-b.coef(i);
		while(c.a.size()>1 && isZero(*(c.a.end()-1)))c.a.erase(c.a.end()-1);
		return c;
	}
	Poly operator * (Poly &b){
		Poly c;
		c.a.resize(a.size()+b.a.size()-1);
		for(int i=0;i<c.a.size();i++)c.a[i]=0;
		for(int i=0;i<a.size();i++)
			for(int j=0;j<b.a.size();j++)
				c.a[i+j]=c.a[i+j]+a[i]*b.a[j];
		while(c.a.size()>1 && isZero(*(c.a.end()-1)))c.a.erase(c.a.end()-1);
		return c;
	}
};
Poly Cn(ll k){
	Poly ans,t;
	ans=frac(1,1);
	for(ll i=0;i<k;i++){
		t=Poly((vector<type>){frac(-i,1),frac(1,1)});
		ans=ans*t;
		t=Poly((vector<type>){frac(1,i+1)});
		ans=ans*t;
	}
	return ans;
}
Poly sum(ll k){
	if(k==0)return Poly((vector<type>){frac(0,1),frac(1,1)});
	Poly ans=Cn(k+1),f=Cn(k),t,p;
	ans=ans+f;
	for(int i=1;i<k;i++){
		t=f(Poly((vector<type>){frac(i,1)}));
		ans=ans+t;
	}
	for(int i=0;i<k;i++){
		p=sum(i);
		t=f[i];
		t=t*p;
		ans=ans-t;
	}
	t=frac(1,1)/f[k];
	ans=ans*t;
	return ans;
}
Poly sum(Poly &f){
	Poly ans,t,p;
	ans=frac(0,1);
	for(int i=0;i<=f.n();i++){
		t=f[i];p=sum(i);t=t*p;
		ans=ans+t;
	}
	return ans;
}
Poly Lagrange(vector<type> x,vector<type> y){
	int n=x.size()-1;
	Poly ans;
	rep(k,0,n){
		Poly t,p;
		t=y[k];
		rep(j,0,n)if(j!=k){
			p=(vector<type>){frac(0,1)-x[j]/(x[k]-x[j]),frac(1,1)/(x[k]-x[j])};
			t=t*p;
		}
		ans=ans+t;
	}
	return ans;
}
int main(){
	cout<<sum(1);
	cout<<sum(2);
	cout<<sum(3);
	cout<<sum(4);
	return 0;
}

运行结果：

(1/2)x^2+(1/2)x^1+(0)
(1/3)x^3+(1/2)x^2+(1/6)x^1+(0)
(1/4)x^4+(1/2)x^3+(1/4)x^2+(0)x^1+(0)
(1/5)x^5+(1/2)x^4+(1/3)x^3+(0)x^2+(1/-30)x^1+(0)

1.2 一般多项式的前缀和

对于一般的多项式
$$f(x)=\sum _{i=0}^k{a}_i{x}^i$$

$$\sum _{i=1}^{n}f(i)=\sum _{i=0}^k\left(a_i\sum _{j=1}^n{j}^i\right)$$

多项式的前缀和依然是多项式，这是多么美妙的结论！

例题

hdu多校第一场 Sequence 貌似n阶前缀和/矩阵快速幂但其实要用组合数。。。

题解：令$b_i={\sum }_{j=i-k\cdot x}{a}_j(0\le x,1\le j\le i)$，则$\sum b_i{x}^i=(\sum a_i{x}^i)(\sum x^{ik})$，因此操作的顺序没有影响，可以统计各类操作次数，然后批量处理，也就是$\sum x^{ik}$的幂。

由于$n,m$比较大，直接做快速幂或者$\ln +\exp$大概会$T$，但是$(\sum x^{ik}{)}^n=\sum (\genfrac{}{}{0ex}{}{n-1+i}{i})x^{ik}$,因此可以省掉求幂，对于一种操作只做一次卷积即可。

代码

#include <bits/stdc++.h>

const int XN=2e6+11, P=998244353;

int Pow(long long base,int v) {
	long long res;
	for(res=1;v;v>>=1,(base*=base)%=P)
		if(v&1)
			(res*=base)%=P;
	return res;
}

int D(int x) {
	((x>=P) && (x-=P)) || ((x<0) && (x+=P));
	return x;
}

void NTT(int a[],int n,int op) {
	for(int i=1,j=n>>1;i<n-1;++i) {
		if(i<j)
			std::swap(a[i],a[j]);
		int k=n>>1;
		while(k<=j) {
			j-=k;
			k>>=1;
		}
		j+=k;
	}
	for(int len=2;len<=n;len<<=1) {
		int rt=Pow(3,(P-1)/len);
		for(int i=0;i<n;i+=len) {
			int w=1;
			for(int j=i;j<i+len/2;++j) {
				int u=a[j],t=1LL*a[j+len/2]*w%P;
				a[j]=D(u+t),a[j+len/2]=D(u-t);
				w=1LL*w*rt%P;
			}
		}
	}
	if(op==-1) {
		std::reverse(a+1,a+n);
		int in=Pow(n,P-2);
		for(int i=0;i<n;++i)
			a[i]=1LL*a[i]*in%P;
	}
}

std::vector<int> Conv(std::vector<int> const &A,std::vector<int> const &B,int N) {
	static int a[XN],b[XN];
	auto Make2=[](int x)->int {
		return 1<<((32-__builtin_clz(x))+((x&(-x))!=x));
	};
	int n=Make2(A.size()+B.size()-1);
	for(int i=0;i<n;++i) {
		a[i]=i<A.size()?A[i]:0;
		b[i]=i<B.size()?B[i]:0;
	}
	NTT(a,n,1);NTT(b,n,1);
	for(int i=0;i<n;++i)
		a[i]=1LL*a[i]*b[i]%P;
	NTT(a,n,-1);
	std::vector<int> C(N);
	for (int i=0;i<N;i++)
		C[i]=a[i];
	return C;
}

const int M=2e6,XM=M+11;

int fact[XM],ifact[XM];

int C(int n,int m) {
	return m<=n?(long long)fact[n]*ifact[m]%P*ifact[n-m]%P:0;
}

int main() {
	std::ios::sync_with_stdio(0);
	std::cin.tie(0);
	fact[0]=1;
	for(int i=1;i<=M;++i)
		fact[i]=(long long)fact[i-1]*i%P;
	ifact[M]=Pow(fact[M],P-2);
	for(int i=M-1;i>=0;--i)
		ifact[i]=(long long)ifact[i+1]*(i+1)%P;
	int T;std::cin>>T;
	while(T--) {
		int n,m;std::cin>>n>>m;
		
		std::vector<int> a(n);
		for(int i=0;i<n;++i) {
			std::cin>>a[i];
		}

		int cnt[]={0,0,0,0};
		for(int i=1;i<=m;++i) {
			int x;std::cin>>x;
			cnt[x]++;
		}

		for(int c=1;c<=3;++c) if(cnt[c]) {
			std::vector<int> h(n);
			for(int i=0;i*c<n;++i)
				h[i*c]=C(cnt[c]-1+i,i);
			a=Conv(h,a,n);
		}
		long long Ans=0;
		for(int i=0;i<n;++i)
			Ans^=1LL*(i+1)*a[i];
		std::cout<<Ans<<'\n';
	}
	return 0;
}