Lagrange 插值 & 高斯消元

Lagrange 插值 \(\And\) 高斯消元

~~你说得对，高斯消元是后面补得，因为没啥可说的，放在一起吧。~~

~~你说得对，我是 sb，先学拉插再学高斯消元~~

拉插是用来进行多项式插值的有力工具。

统一设最后的函数为 \(g\)，点分别为 \((x_1,y_1)...(x_n,y_n)\)。

基础

普通拉插：

考虑构造函数 \(f(i)\)，使其满足：

\[\forall x=x_j(j \in [1,n], j\not=i),f(x)=0 \land f(x_i)=y_i \]
有构造：

\[f(i)=y_i\prod_{j\not=i}\frac{x-x_i}{x_j-x_i} \]
比较显然。

然后可以让 \(g(i)=\sum f(i)=\sum\limits_{i=1}^n y_i\prod\limits_{j\not=i}\frac{x-x_i}{x_j-x_i}\)

暴力实现是 \(O(n^2)\) 的。

可以优化到 \(O(n\log^2 n)\)，挺麻烦，我不会，link
连续拉插：

对于 \(x=[1,n]\) 有 \(O(n)\) 拉插。

将原式替换一下

\[g(i)=\sum_{i=1}^n y_i\prod_{j\not=i}\frac{x-i}{j-i} \]
我们预处理出 \(x-i\) 的前后缀积，分子就干完了。发现分母可以拆成 \((i-1)!\times (n-i)! \times (-1)^{n-i}\)，预处理阶乘逆元就行了。
重心拉插：

可以让 \(O(n)\) 加新点~~感觉可以平替牛顿插值了~~。

将原式变换一下：

\[\begin{aligned} g(i)&=\sum_{i=1}^n y_i\prod_{j\not=i}\frac{x-x_i}{x_j-x_i}\\&=\prod_{i=1}^n(x-x_i)\sum_{i=1}^n (\frac{1}{x-x_i} \times \prod_{i!=j}\frac{y_i}{x_i-x_j})\end{aligned} \]
每次加点直接求 \(\prod_{i!=j}\frac{y_i}{x_i-x_j}\) 即可。

但是做题没用过。
求系数：

直接模拟是 \(O(n^3)\) 的，先吧 \(\prod\limits_{i=1}^n (x-x_i)\) 处理出来在模拟多项式除 \(x-t\) 即可 \(O(n^2)\)

有用定理：

若 \(\Delta f(i)=f(i)-f(i-1)\) 是一个 \(k\) 次多项式，则 \(f\) 是 \(k\) 次多项式。

用法：

一个 \(\sum\) 会让 \(f\) 次数 \(+1\)
对于天然查分的 \(dp\) 可以设一个次数，用这个解出来，见最后一个题。

题目

一个人的数论：

和莫反书接上回是吧。

莫反，发现有 \(\sum_{i=1}^n i^k\)，这是一个 \(k+1\) 次多项式，形如 \(\sum_{i=1}^{d+1}c_in^i\) （\(c\) 是系数），将其带入在推即可，最后拉插求系数（也可以高斯消元）。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int W=1003,D=103,MOD=1e9+7;
int d,w,n,ans,cp[W],ca[W];
vector<int> ck;

namespace MT{
	Il int Fpw(int a,int b){
		int ans=1;
		while(b){
			if(b&1) ans=1ll*ans*a%MOD;
			a=1ll*a*a%MOD,b>>=1;
		}
		return ans;
	}
	Il int Inv(Ct int &a){return Fpw(a,MOD-2);}
	class LGG{
	private:
		int cy[D];
		vector<int> tmp,mul;
	public:
		LGG(){tmp.reserve(D),mul.reserve(D);}
		Il void Gety(){For(i,1,d+2,1) cy[i]=(cy[i-1]+Fpw(i,d))%MOD;}
		Il void operator()(vector<int> &f){
			int len=d+2;
			f.clear(),f.resize(len),tmp.clear(),mul.clear(),mul.emplace_back(1);
			For(i,1,len,1){
				tmp.emplace_back(0); for(Ct int &k:mul) tmp.emplace_back(k);
				For(j,0,i-1,1) tmp[j]=(tmp[j]-1ll*mul[j]*i%MOD+MOD)%MOD;
				swap(tmp,mul),tmp.clear();
			}
			For(i,1,len,1){
				int fd=1; tmp.clear(),tmp.resize(len);
				For_(j,len,2,1) tmp[j-1]=(tmp[j-1]+mul[j])%MOD,tmp[j-2]=(tmp[j-2]+1ll*tmp[j-1]*i%MOD)%MOD;
				tmp[0]=(tmp[0]+mul[1])%MOD;
				For(j,1,len,1) if(i^j) fd=1ll*fd*(i-j)%MOD;
				fd=1ll*Inv((fd+MOD)%MOD)*cy[i]%MOD;
				For(j,0,len-1,1) f[j]=(f[j]+1ll*tmp[j]*fd%MOD)%MOD;
			}
		}
	}Lgg;
}

int main(){
	read(d,w); n=1;
	For(i,1,w,1) read(cp[i],ca[i]),n=1ll*n*MT::Fpw(cp[i],ca[i])%MOD;
	MT::Lgg.Gety(); MT::Lgg(ck);
	For(i,1,d+1,1){
		int tmp=1ll*ck[i]*MT::Fpw(n,i)%MOD;
		For(j,1,w,1) tmp=1ll*tmp*(1-MT::Fpw(cp[j],(d-i+MOD-1)%(MOD-1)))%MOD;
		ans=(ans+tmp)%MOD;
	}
	write((ans+MOD)%MOD);
}

成绩比较

拉插优化 dp。

设 \(dp_{i,j}\) 表示前 i 门课，有 \(j\) 个被碾压方案数。

\[dp_{i,j}=\sum_{l=j}^kdp_{i-1,l}\binom{j}{l}\binom{n-r_i-j}{n-1-l}\sum_{p=1}^{u_i}p^{n-r_i}(u_i-p)^{r_i-1} \]

发现后面统计分数是一个 \(n\) 次多项式，直接拉插求值就可以避免枚举到 \(u_i\)。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

const int N=103,MOD=1e9+7;
int n,m,ck,u[N],cr[N],dp[N][N];

namespace MT{
	int fac[N],ivf[N];
	Il int Fpw(int a,int b){
		int ans=1;
		while(b){
			if(b&1) ans=1ll*ans*a%MOD;
			b>>=1,a=1ll*a*a%MOD;
		}
		return ans;
	}
	Il int Inv(Ct int &a){return Fpw(a,MOD-2);}
	struct INIT{INIT(){
		fac[0]=1; For(i,1,N-3,1) fac[i]=1ll*fac[i-1]*i%MOD;
		ivf[N-3]=Inv(fac[N-3]);
		For_(i,N-3,1,1) ivf[i-1]=1ll*ivf[i]*i%MOD;
	}}Initer;
	Il int C(int a,int b){return (a<b||b<0)?0:1ll*fac[a]*ivf[b]%MOD*ivf[a-b]%MOD;}
	class LGG{
	private: int cy[N],lmul[N],rmul[N];
	public:
		Il void Gety(Ct int &r,Ct int &u){
			For(i,1,n+1,1) cy[i]=(cy[i-1]+1ll*Fpw(i,n-r)*Fpw(u-i,r-1)%MOD)%MOD;
		}
		Il int operator()(Ct int &x){
			int len=n+1,y=0;
			lmul[0]=1; For(i,1,len,1) lmul[i]=1ll*lmul[i-1]*(x-i)%MOD;
			rmul[len+1]=1; For_(i,len,1,1) rmul[i]=1ll*rmul[i+1]*(x-i)%MOD;
			For(i,1,len,1) y=(y+1ll*cy[i]*lmul[i-1]%MOD*rmul[i+1]%MOD*ivf[i-1]%MOD*ivf[len-i]%MOD*(len-i&1?-1:1)+MOD)%MOD;
			return y;
		}
	}Lgg;
}

int main(){
	read(n,m,ck);
	For(i,1,m,1) read(u[i]);
	For(i,1,m,1) read(cr[i]);
	dp[0][n-1]=1;
	For(i,1,m,1){
		MT::Lgg.Gety(cr[i],u[i]);
		int tmp=MT::Lgg(u[i]);
		For(j,ck,n-1,1) For(k,j,n-1,1)
			dp[i][j]=(dp[i][j]+1ll*dp[i-1][k]*MT::C(k,j)%MOD*MT::C(n-1-k,n-cr[i]-j)%MOD*tmp%MOD)%MOD;
	}
	write(dp[m][ck],'\n');
}

教科书般的亵渎

sb 题面描述，k 始终不变。

式子显然，用拉插求 \(\sum_{i=1}^n i^k\) 即可。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int M=55,MOD=1e9+7;
int t,n,m,ca[M];

namespace MT{
	Il int Fpw(int a,int b){
		int ans=1;
		while(b){
			if(b&1) ans=1ll*ans*a%MOD;
			b>>=1,a=1ll*a*a%MOD;
		}
		return ans;
	}
	Il int Inv(Ct int &a){return Fpw(a,MOD-2);}
	class LGG{
	private:
		int cy[M],lmul[M],rmul[M],fac[M],ivf[M];
	public:
		LGG(){
			fac[0]=1; For(i,1,M-3,1) fac[i]=1ll*fac[i-1]*i%MOD;
			ivf[M-3]=Inv(fac[M-3]);
			For_(i,M-3,1,1) ivf[i-1]=1ll*ivf[i]*i%MOD;
		}
		Il void Gety(){int k=m+1; For(i,1,k+2,1) cy[i]=(cy[i-1]+Fpw(i,k))%MOD;}
		Il int operator()(Ct int &x){
			int len=m+3,y=0;
			lmul[0]=1; For(i,1,len,1) lmul[i]=1ll*lmul[i-1]*(x-i)%MOD;
			rmul[len+1]=1; For_(i,len,1,1) rmul[i]=1ll*rmul[i+1]*(x-i)%MOD;
			For(i,1,len,1) y=(y+1ll*cy[i]*lmul[i-1]%MOD*rmul[i+1]%MOD*ivf[i-1]%MOD*ivf[len-i]%MOD*(len-i&1?-1:1)+MOD)%MOD;
			return y;
		}
	}Lgg;
}

int main(){
	read(t);
	while(t--){
		read(n,m); MT::Lgg.Gety();
		For(i,1,m,1) read(ca[i]); ca[0]=0;
		sort(ca,ca+m+1);
		int ans=0;
		For(i,0,m,1){
			ans=(ans+MT::Lgg(n-ca[i]))%MOD;
			For(j,i+1,m,1) ans=(ans-MT::Fpw(ca[j]-ca[i],m+1)+MOD)%MOD;
		}
		write(ans,'\n');
	}
}

tyvj 1858 XLkxc

奇妙题目。

拉插套拉插，由有用定理（见上）可知，其是一个 \(k+3\) 次多项式，因为其中有 \(a+id\) 太大，无法直接求，需要提前插出后面函数。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)
#define int long long

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int K=150,MOD=1234567891;
int t,k,ca,n,d;

namespace MT{
	int fac[K],ivf[K];
	Il int Fpw(int a,int b){
		int ans=1;
		while(b){
			if(b&1) ans=1ll*ans*a%MOD;
			b>>=1,a=1ll*a*a%MOD;
		}
		return ans;
	}
	Il int Inv(Ct int &a){return Fpw(a,MOD-2);}
	struct INIT{INIT(){
		fac[0]=1; For(i,1,K-3,1) fac[i]=1ll*fac[i-1]*i%MOD;
		ivf[K-3]=Inv(fac[K-3]);
		For_(i,K-3,1,1) ivf[i-1]=1ll*ivf[i]*i%MOD;
	}}Initer;
	class LGG{
	private:
		int cy[K],lmul[K],rmul[K];
	public:
		Il int &Y(Ct int &x){return cy[x];}
		Il int operator()(Ct int &x,Ct int &len){
			int y=0;
			lmul[0]=1; For(i,1,len,1) lmul[i]=1ll*lmul[i-1]*(x-i)%MOD;
			rmul[len+1]=1; For_(i,len,1,1) rmul[i]=1ll*rmul[i+1]*(x-i)%MOD;
			For(i,1,len,1) y=(y+1ll*cy[i]*lmul[i-1]%MOD*rmul[i+1]%MOD*ivf[i-1]%MOD*ivf[len-i]%MOD*(len-i&1?-1:1)+MOD)%MOD;
			return y;
		}
	}La,Lb;
	Il void Init(){
		int pw[K]={0}; For(i,1,k+3,1) pw[i]=Fpw(i,k);
		For(i,1,k+3,1){
			Lb.Y(i)=0;
			For(j,1,i,1) Lb.Y(i)=(Lb.Y(i)+1ll*(i-j+1)*pw[j]%MOD)%MOD;
		}
		La.Y(0)=Lb(ca,k+3); For(i,1,k+4,1) La.Y(i)=(La.Y(i-1)+Lb((ca+i*d)%MOD,k+3))%MOD;
	}
}

signed main(){
	read(t);
	while(t--){
		read(k,ca,n,d);
		MT::Init();
		write(MT::La(n,k+4),'\n');
	}
}

calc

另一种拉插优化 dp。

设 \(dp_{i,j}\) 表示值域到 \(i\)，选 \(j\) 个方案数，枚举第 \(i\) 个选不选，有。

\[dp_{i,j}=dp_{i-1,j-1}\times i\times j+dp_{i-1,j} \]

发现它天然差分，有：

\[dp_{i,j}-dp_{i-1,j}=dp_{i-1,j-1}\times i\times j \]

设 \(dp_{i,j}\) 是关于 \(i\) 的 \(t_j\) 次多项式，有：

\[t_j-1=t_{j-1}+1 \]

显然 \(t\) 是等差序列，有 \(t_n=2n\)。

所以只求 \(2*n+1\) 项即可。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int N=1010;
int MOD,dp[N][N],n,ca;

namespace MT{
	int fac[N],ivf[N];
	Il int Fpw(int a,int b){
		int ans=1;
		while(b){
			if(b&1) ans=1ll*ans*a%MOD;
			b>>=1,a=1ll*a*a%MOD;
		}
		return ans;
	}
	Il int Inv(Ct int &a){return Fpw(a,MOD-2);}
	Il void Init(){
		fac[0]=1; For(i,1,N-3,1) fac[i]=1ll*fac[i-1]*i%MOD;
		ivf[N-3]=Inv(fac[N-3]);
		For_(i,N-3,1,1) ivf[i-1]=1ll*ivf[i]*i%MOD;
	}
	class LGG{
	private: int cy[N],lmul[N],rmul[N];
	public:
		Il int &Y(Ct int &i){return cy[i];}
		Il int operator()(Ct int &x,Ct int &len){
			int y=0;
			lmul[0]=1; For(i,1,len,1) lmul[i]=1ll*lmul[i-1]*(x-i)%MOD;
			rmul[len+1]=1; For_(i,len,1,1) rmul[i]=1ll*rmul[i+1]*(x-i)%MOD;
			For(i,1,len,1) y=(y+1ll*cy[i]*lmul[i-1]%MOD*rmul[i+1]%MOD*ivf[i-1]%MOD*ivf[len-i]%MOD*(len-i&1?-1:1)+MOD)%MOD;
			return y;
		}
	}Lgg;
}

int main(){
	read(ca,n,MOD);
	MT::Init();
	int m=n*2+1;
	For(i,0,m,1) dp[i][0]=1;
	For(i,1,m,1) For(j,1,n,1) dp[i][j]=(1ll*dp[i-1][j-1]*i*j%MOD+dp[i-1][j])%MOD;
	For(i,1,m,1) MT::Lgg.Y(i)=dp[i][n];
	write(MT::Lgg(ca,m));
}

高斯消元

其实真正的高斯消元没啥用，真正用的多的都是高斯约旦了。

高斯约旦分别判无解和自由元：

每次新消元时从一开始选，因为可以选之前因为是 \(0\) 而跳过的，注意一下优先级，具体见代码~~十分压行~~：

CODE

Il int Gss(double a[N][N]){
	For(i,1,n,1){
		int nw=i; For(j,1/*not i*/,n,1) if((i<=j||Isz(a[j][j])/*judge it not used*/)&&fabs(a[j][i])>fabs(a[nw][i])) nw=j; 
		if(Isz(a[nw][i])) continue; swap(a[i],a[nw]);
		For(j,1,n,1) if(i^j){double tmp=a[j][i]/a[i][i]; For(k,i,n+1,1) a[j][k]-=tmp*a[i][k];}
	}
	int fg=1; For(i,1,n,1) Isz(a[i][i])?fg=-(!Isz(a[i][n+1])||fg==-1/*priority!!*/):a[i][n+1]/=a[i][i]; return fg;
}

题都是板子。

线性方程组

真·板子。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int N=103;
constexpr double eps=1e-7;
int n; double c[N][N];
Il bool Isz(Ct double f){return fabs(f)<eps;}

Il int Gss(double a[N][N]){
	For(i,1,n,1){
		int nw=i; For(j,1,n,1) if((i<=j||Isz(a[j][j]))&&fabs(a[j][i])>fabs(a[nw][i])) nw=j; 
		if(Isz(a[nw][i])) continue; swap(a[i],a[nw]);
		For(j,1,n,1) if(i^j){double tmp=a[j][i]/a[i][i]; For(k,i,n+1,1) a[j][k]-=tmp*a[i][k];}
	}
	int fg=1; For(i,1,n,1) Isz(a[i][i])?fg=-(!Isz(a[i][n+1])||fg==-1):a[i][n+1]/=a[i][i]; return fg;
}

int main(){
	read(n);
	For(i,1,n,1) For(j,1,n+1,1) c[i][j]=read<int>();
	int ans=Gss(c);
	if(ans<=0) write(ans,'\n');
	else For(i,1,n,1) Isz(c[i][n+1])?printf("x%d=0\n",i):printf("x%d=%.2lf\n",i,c[i][n+1]);
}

球形空间产生器sphere

题意模拟有 \(n+1\) 个二次方程。

发现给了 \(n+1\) 个，直接差分消去二次项，就变成 \(n\) 个一次方程。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

#define UN_FAST

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int N=13;
constexpr double eps=1e-7;
int n;
double cpos[N][N],c[N][N];

namespace MT{
	Il bool Isz(Ct double f){return fabs(f)<eps;}
	Il void Gss(double a[N][N]){
		For(i,1,n,1){
			int nw=i; For(j,i+1,n,1) if(fabs(a[j][i])>fabs(a[nw][i])) nw=j; swap(a[i],a[nw]); 
			For(j,1,n,1) if(i^j){double tmp=a[j][i]/a[i][i]; For(k,i,n+1,1) a[j][k]-=tmp*a[i][k];}
		} For(i,1,n,1) a[i][n+1]/=a[i][i];
	}
}

int main(){
	read(n);
	For(i,1,n+1,1) For(j,1,n,1) scanf("%lf",&cpos[i][j]);
	For(i,1,n,1){
		double tmp=0;
		For(j,1,n,1) c[i][j]=2*(cpos[i+1][j]-cpos[i][j]),tmp+=cpos[i+1][j]*cpos[i+1][j]-cpos[i][j]*cpos[i][j];
		c[i][n+1]=tmp;
	}
	MT::Gss(c);
	For(i,1,n,1) printf("%.3f ",c[i][n+1]);
}

臭气弹

考虑类似 dp 的转移，考虑从其他点转移来的概率，设概率，解方程，注意起点加一。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int N=303;
constexpr double eps=1e-7;
int n,m,gph[N][N],dgr[N];
double cpos[N][N],c[N][N],p,q;

namespace MT{
	Il bool Isz(Ct double f){return fabs(f)<eps;}
	Il void Gss(double a[N][N]){
		For(i,1,n,1){
			int nw=i; For(j,i+1,n,1) if(fabs(a[j][i])>fabs(a[nw][i])) nw=j; swap(a[i],a[nw]); 
			For(j,1,n,1) if(i^j){double tmp=a[j][i]/a[i][i]; For(k,i,n+1,1) a[j][k]-=tmp*a[i][k];}
		} For(i,1,n,1) a[i][n+1]/=a[i][i];
	}
}

int main(){
	read(n,m); p=read<int>(),q=read<int>();
	For(i,1,m,1){
		int u,v; read(u,v);
		gph[u][v]=gph[v][u]=1,++dgr[u],++dgr[v];
	}
	For(i,1,n,1){
		c[i][i]=-1;
		For(j,1,n,1) if(gph[i][j]) c[i][j]=(q-p)/q/dgr[j];
	} c[1][n+1]=-1;
	MT::Gss(c);
	For(i,1,n,1) printf("%.6lf\n",c[i][n+1]*p/q);
}

开关问题

设每个开关状态，容易列出异或方程。

统计个数就是 \(2^\text{自由元个数}\) 因为每个自由元有两个状态（这个卡了我好久，我是 sb）。

CODE

#include<bits/stdc++.h>
#include<sys/mman.h>
#include<fcntl.h>
using namespace std;
using llt=long long;
using llf=long double;
using ull=unsigned long long;
#define Ct const
#define Il __always_inline
#define For(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
#define For_(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
#define For_it(i,a,b) for(auto i=(a);i!=(b);++i)

namespace IO{
#ifdef ONLINE_JUDGE
	int Fin=fileno(stdin); FILE* Fout=stdout;
#elif defined(UN_FAST)
	FILE *Fin=freopen("in_out/in.in","r",stdin),*Fout=freopen("in_out/out.out","w",stdout);
#else
	int Fin=open("in_out/in.in",0); FILE *Fout=freopen("in_out/out.out","w",stdout);
#endif // file
#ifdef UN_FAST
	char cc;
	#define G (cc=getchar())
	#define C cc
#else
	const char *I=(char*)mmap(0,1<<28,1,2,Fin,0)-1; 
	#define G (*++I)
	#define C (*I)
#endif // fast (mmap)
#define P(x) putc_unlocked(x,Fout)
	template<class T> Il void read(T &x){x=0;bool f=0; while(f|=G==45,C<48); while(x=x*10+(C&15),G>47); f?x=-x:0;}
	template<class T> void write(T x){if(x<0) P('-'),x=-x; if(x/10) write(x/10); P('0'+x%10);}
	Il void read(char &c){while(G<33); c=C;} void read(char *s){char *p=s-1; while(G<33); while(*++p=C,G>32); *++p='\0';}
	Il void read(string &s){s.clear(); while(G<33); while(s.push_back(C),G>32);}
	template<class T=int> Il T read(){T a; return read(a),a;}
	template<class T,class ...Argc> Il void read(T &x,Argc &...argc){read(x),read(argc...);}
	Il void write(const char &c){P(c);} Il void write(const char *s){for(const char *c=s;*c!='\0';++c) P(*c);}
	Il void write(const string &s){for(const char &c:s) P(c);}
	template<class T,class ...Argc> Il void write(const T &x,const Argc &...argc){write(x),write(argc...);}
	template<class T> Il void Write(const T &x){write(x),P(' ');} Il void Write(const char &c){P(c); if(c>32) P(' ');}
	template<class T,class ...Argc> Il void Write(const T &x,const Argc &...argc){Write(x),Write(argc...);}
#undef G
#undef C
#undef P
}using IO::read; using IO::write; using IO::Write;

constexpr int N=303;
int t,n;
bool c[N][N],bg[N];

Il int Gss(bool a[N][N]){
	For(i,1,n,1){
		int nw; For(j,1,n,1) if((i<=j||!a[j][j])&&a[j][i]) nw=j;
		if(!a[nw][i]) continue; swap(a[i],a[nw]); 
		For(j,1,n,1) if(i^j&&a[j][i]) For(k,i,n+1,1) a[j][k]^=a[i][k];
	}
	int tmp=0; For(i,1,n,1) if(!a[i][i]){if(a[i][n+1]) tmp=-1; else if(tmp!=-1) ++tmp;} return tmp;
}

int main(){
	read(t);
	while(t--){
		read(n); memset(c,0,sizeof(c));
		For(i,1,n,1) c[i][i]=1,read(bg[i]);
		For(i,1,n,1) c[i][n+1]=read()^bg[i];
		int u,v; while(read(v,u),u||v) c[u][v]=1;
		int ans=Gss(c);
		if(ans==-1) write("Oh,it's impossible~!!\n");
		else write((llt)pow(2,ans),'\n');
	}
}

带状矩阵

这就不能约旦了，有交换行和交换列两种做法。

交换列的要维护当前元的编号，交换行的有二倍常数（因为最多会让这一行长度乘二），并且在只有一半带状时没法做。

全部带状 $nd^2$，一半带状 $n^2d$。

posted @ 2024-05-16 16:47 xrlong 阅读(7) 评论(0) 编辑收藏举报

会员力量，点亮园子希望

刷新页面返回顶部

xrlong

Lagrange 插值 & 高斯消元

Lagrange 插值 \(\And\) 高斯消元

基础

题目

高斯消元

带状矩阵

公告