数学 / Math（ACM-Template 2.0）

数学

数论

整数环相关 and 扩展欧几里得

ll mul(ll a,ll b,ll m=mod){return a*b%m;} // 模乘
ll qpow(ll a,ll b,ll m=mod){ // 快速幂
	ll ans=1;
	for(;b;a=mul(a,a,m),b>>=1)
		if(b&1)ans=mul(ans,a,m);
	return ans;
}
void exgcd(ll a,ll b,ll &d,ll &x,ll &y){ // ax + by = gcd(a, b) = d
	if(!b)d=a,x=1,y=0;
	else exgcd(b,a%b,d,y,x),y-=x*(a/b);
}
ll gcdinv(ll v,ll m=mod){ // 扩欧版逆元
	ll d,x,y;
	exgcd(v,m,d,x,y);
	return (x%m+m)%m;
}
ll getinv(ll v,ll m=mod){ // 快速幂版逆元，m 必须是质数!!
	return qpow(v,m-2,m);
}
ll qpows(ll a,ll b,ll m=mod){
	if(b>=0)return qpow(a,b,m);
	else return getinv(qpow(a,-b,m),m);
}

阶乘 and 组合数

  • \(O(n)\) 初始化，\(O(1)\) 查询
struct Comb {
	static const int N = 100010;
	int fac[N], inv[N];
	Comb() {
		fac[0] = 1;
		repeat (i, 1, N)
			fac[i] = 1ll * fac[i - 1] * i % mod;
		inv[N - 1] = qpow(fac[N - 1], mod - 2, mod);
		repeat_back (i, 1, N)
			inv[i - 1] = 1ll * inv[i] * i % mod;
	}
	int operator()(int a, int b) { // a >= b
		if (a < b || b < 0) return 0;
		return 1ll * fac[a] * inv[a - b] % mod * inv[b] % mod;
	}
	int A(int a, int b) { // a >= b
		if (a < b || b < 0) return 0;
		return 1ll * fac[a] * inv[a - b] % mod;
	}
} C;

防爆模乘

// int128版本
ll mul(ll a,ll b,ll m=mod){return (__int128)a*b%m;}
// long double版本（欲防爆，先自爆）
ll mul(ll a,ll b,ll m){
	ll c=a*b-(ll)((long double)a*b/m+0.5)*m;
	return c<0?c+m:c;
}
// 每位运算一次版本，注意这是真·龟速乘，O(logn)
ll mul(ll a,ll b,ll m=mod){
	ll ans=0;
	while(b){
		if(b&1)ans=(ans+a)%m;
		a=(a+a)%m;
		b>>=1;
	}
	return ans;
}
// 把 b 分成两部分版本，要保证 m 小于 1<<42（约等于 4e12），a, b < m
ll mul(ll a,ll b,ll m=mod){
	a%=m,b%=m;
	ll l=a*(b>>21)%m*(1ll<<21)%m;
	ll r=a*(b&(1ll<<21)-1)%m;
	return (l+r)%m;
}

最大公约数

__gcd(a, b) // 内置 gcd，推荐
ll gcd(ll a, ll b) { return b == 0 ? a : gcd(b, a % b); } // 不推荐，比内置 gcd 慢
ll gcd(ll a, ll b) { // 卡常 gcd！！
	#define tz __builtin_ctzll
	if (!a || !b) return a | b;
	int t = tz(a | b);
	a >>= tz(a);
	while (b) {
		b >>= tz(b);
		if (a > b) swap(a, b);
		b -= a;
	}
	return a << t;
	#undef tz
}
  • 实数 gcd
lf fgcd(lf a,lf b){return abs(b)<1e-5?a:fgcd(b,fmod(a,b));}

CRT + extra

// CRT，m[i]两两互质
ll crt(ll a[],ll m[],int n){ // ans%m[i]==a[i]
	repeat(i,0,n)a[i]%=m[i];
	ll M=1,ans=0;
	repeat(i,0,n)
		M*=m[i];
	repeat(i,0,n){
		ll k=M/m[i],t=gcdinv(k%m[i],m[i]); // 扩欧!!
		ans=(ans+a[i]*k*t)%M; // 两个乘号可能都要mul
	}
	return (ans+M)%M;
}
// exCRT，m[i]不需要两两互质，基于扩欧exgcd和龟速乘mul
ll excrt(ll a[],ll m[],int n){ // ans%m[i]==a[i]
	repeat(i,0,n)a[i]%=m[i]; // 根据情况做适当修改
	ll M=m[0],ans=a[0],g,x,y; // M是m[0..i]的最小公倍数
	repeat(i,1,n){
		ll c=((a[i]-ans)%m[i]+m[i])%m[i];
		exgcd(M,m[i],g,x,y); // Ax=c(mod B)
		if(c%g)return -1;
		ans+=mul(x,c/g,m[i]/g)*M; // 龟速乘
		M*=m[i]/g;
		ans=(ans%M+M)%M;
	}
	return (ans+M)%M;
}

离散对数 using BSGS + extra

  • 求 \(a^x \equiv b \pmod m\) ，\(O(\sqrt m)\)
// BSGS，a 和 mod 互质
ll bsgs(ll a,ll b,ll mod){ // a^ans%mod==b
	a%=mod,b%=mod;
	static unordered_map<ll,ll> m; m.clear();
	ll t=(ll)sqrt(mod)+1,p=1;
	repeat(i,0,t){
		m[mul(b,p,mod)]=i; // p==a^i
		p=mul(p,a,mod);
	}
	a=p; p=1;
	repeat(i,0,t+1){
		if(m.count(p)){ // p==a^i
			ll ans=t*i-m[p];
			if(ans>0)return ans;
		}
		p=mul(p,a,mod);
	}
	return -1;
}
// exBSGS，a 和 mod不需要互质，基于 BSGS
ll exbsgs(ll a,ll b,ll mod){ // a^ans%mod==b
	a%=mod,b%=mod;
	if(b==1)return 0;
	ll ans=0,c=1,g;
	while((g=__gcd(a,mod))!=1){
		if(b%g!=0)return -1;
		b/=g,mod/=g;
		c=mul(c,a/g,mod);
		ans++;
		if(b==c)return ans;
	}
	ll t=bsgs(a,mul(b,getinv(c,mod),mod),mod); // 必须扩欧逆元!!
	if(t==-1)return -1;
	return t+ans;
}

阶与原根

  • 一些原根
if (m == 167772161) return 3;
if (m == 469762049) return 3;
if (m == 754974721) return 11;
if (m == 998244353) return 3;
  • 判断是否有原根：若 m 有原根，则 m 一定是下列形式：\(2,4,p^a,2p^a\)（ p 是奇素数， a 是正整数）
  • 求所有原根：若 g 为 m 的一个原根，则 \(g^s\space(1\le s\le\varphi(m),\gcd(s,\varphi(m))=1)\) 给出了 m 的所有原根。因此若 m 有原根，则 m 有 \(\varphi(\varphi(m))\) 个原根
  • 求一个原根，\(O(n\log\log n)\) 实际远远不到
ll getG(ll n){ // 求 n 最小的原根
	static vector<ll> a; a.clear();
	ll k=n-1;
	repeat(i,2,sqrt(k+1)+1)
	if(k%i==0){
		a.push_back(i); // a 存放 (n-1) 的质因数
		while(k%i==0)k/=i;
	}
	if(k!=1)a.push_back(k);
	repeat(i,2,n){ // 枚举答案
		bool f=1;
		for(auto j:a)
		if(qpow(i,(n-1)/j,n)==1){
			f=0;
			break;
		}
		if(f)return i;
	}
	return -1;
}

N 次剩余

  • 求 \(x^a \equiv b \pmod m\) ，基于 BSGS、原根
// 只求一个
ll residue(ll a,ll b,ll mod){ // ans^a%mod==b
	ll g=getG(mod),c=bsgs(qpow(g,a,mod),b,mod);
	if(c==-1)return -1;
	return qpow(g,c,mod);
}
// 求所有N次剩余
vector<ll> ans;
void allresidue(ll a,ll b,ll mod){ // ans^a%mod==b
	ll g=getG(mod),c=bsgs(qpow(g,a,mod),b,mod);
	ans.clear();
	if(b==0){ans.push_back(0);return;}
	if(c==-1)return;
	ll now=qpow(g,c,mod);
	ll step=(mod-1)/__gcd(a,mod-1);
	ll ps=qpow(g,step,mod);
	for(ll i=c%step;i<mod-1;i+=step,now=mul(now,ps,mod))
		ans.push_back(now);
	sort(ans.begin(),ans.end());
}

数论函数

单个欧拉函数

  • \(\varphi(n)=\) 小于 n 且与 n 互质的正整数个数
  • 令 n 的唯一分解式 \(n=\prod({p_k}^{a_k})\)，则 \(\varphi(n)=n\cdot \prod(1-\dfrac 1 {p_k})\)
  • \(O(\sqrt n)\)
int getphi(int n){
	int ans=n;
	repeat(i,2,sqrt(n)+2)
	if(n%i==0){
		while(n%i==0)n/=i;
		ans=ans/i*(i-1);
	}
	if(n>1)ans=ans/n*(n-1);
	return ans;
}

离线乘法逆元

  • 求 \(1..(n-1)\) 的逆元，\(O(n)\)。
void get_inv(int n,int m=mod){
	inv[1]=1;
	repeat(i,2,n)inv[i]=m-m/i*inv[m%i]%m;
}
  • 求 \(a_{1..n}\) 的逆元，离线，\(O(n)\)。
void get_inv(int a[],int n){ // 求a[1..n]的逆元，存在inv[1..n]中
	static int pre[N];
	pre[0]=1;
	repeat(i,1,n+1)
		pre[i]=(ll)pre[i-1]*a[i]%mod;
	int inv_pre=qpow(pre[n],mod-2,mod);
	repeat_back(i,1,n+1){
		inv[i]=(ll)pre[i-1]*inv_pre%mod;
		inv_pre=(ll)inv_pre*a[i]%mod;
	}
}

线性筛

  • 定理：求出 \(f(p)\)（p 为质数）的复杂度不超过 \(O(\log p)\) 的积性函数可以被线性筛。

筛素数和最小质因数

  • prime[i] 表示第 \(i+1\) 个质数，vis[i] == 0 表示 i 是素数，lpf[i] 为 i 的最小质因数。
  • \(O(n)\)。
struct Sieve {
	static const int N = 1000010;
	bool vis[N]; int lpf[N]; vector<int> prime;
	Sieve() {
		vis[1] = 1;
		repeat (i, 2, N) {
			if (!vis[i]) prime.push_back(i), lpf[i] = i;
			for (auto j : prime) {
				if (i * j >= N) break;
				vis[i * j] = 1; lpf[i * j] = j;
				if (i % j == 0) break;
			}
		}
	}
} sieve;

筛欧拉函数

  • 线性版，\(O(n)\)
bool vis[N]; int phi[N]; vector<int> a;
void get_phi(){
	vis[1]=1; phi[1]=1;
	repeat(i,2,N){
		if(!vis[i])a.push_back(i),phi[i]=i-1;
		for(auto j:a){
			if(i*j>=N)break;
			vis[i*j]=1;
			if(i%j==0){phi[i*j]=phi[i]*j; break;}
			phi[i*j]=phi[i]*(j-1);
		}
	}
}
  • 不是线性但节省力气和空间版，\(O(n\log\log n)\)
void get_phi(){
	phi[1]=1; // 其他的值初始化为0
	repeat(i,2,N)if(!phi[i])
	for(int j=i;j<N;j+=i){
		if(!phi[j])phi[j]=j;
		phi[j]=phi[j]/i*(i-1);
	}
}

筛莫比乌斯函数

  • \(O(n)\)
bool vis[N]; int mu[N]; vector<int> a;
void get_mu(){
	vis[1]=1; mu[1]=1;
	repeat(i,2,N){
		if(!vis[i])a.push_back(i),mu[i]=-1;
		for(auto j:a){
			if(i*j>=N)break;
			vis[i*j]=1;
			if(i%j==0){mu[i*j]=0; break;}
			mu[i*j]=-mu[i];
		}
	}
}

筛约数个数

bool vis[N]; int d[N]; vector<int> a;
void get_d(){
	vector<int> c(N); vis[1]=1; d[1]=1,c[1]=0;
	repeat(i,2,N){
		if(!vis[i])a.push_back(i),d[i]=2,c[i]=1;
		for(auto j:a){
			if(i*j>=N)break;
			vis[i*j]=1;
			if(i%j==0){
				d[i*j]=d[i]/(c[i]+1)*(c[i]+2);
				c[i*j]=c[i]+1;
				break;
			}
			d[i*j]=d[i]*2,c[i*j]=1;
		}
	}
}

筛约数之和

bool vis[N]; int d[N]; vector<int> a;
void get_d(){
	vector<int> c(N),sum(N);
	vis[1]=1; d[1]=1,c[1]=0,sum[1]=0;
	repeat(i,2,N){
		if(!vis[i])a.push_back(i),d[i]=i+1,c[i]=i,sum[i]=1+i;
		for(auto j:a){
			if(i*j>=N)break; vis[i*j]=1;
			if(i%j==0){
				c[i*j]=c[i]*j;
				sum[i*j]=sum[i]+c[i*j];
				d[i*j]=d[i]/sum[i]*sum[i*j];
				break;
			}
			d[i*j]=d[i]*d[j],c[i*j]=j,sum[i*j]=1+j;
		}
	}
}

筛 gcd

int gcd[N][N];
void get_gcd(int n,int m){
	repeat(i,1,n+1)
	repeat(j,1,m+1)
	if(!gcd[i][j])
	repeat(k,1,min(n/i,m/j)+1)
		gcd[k*i][k*j]=k;
}

杜教筛

  • 杜教筛只能筛部分积性函数。

\[g(1)S(n)=\sum_{i=1}^n(f*g)(i)-\sum_{i=2}^n g(i)S(\lfloor\dfrac n i \rfloor),S(n)=\sum_{i=1}^nf(i) \]

  • 如果能找到合适的 \(g(n)\)，能快速计算 \(\displaystyle\sum_{i=1}^n(f*g)(i)\)，就能快速计算 \(S(n)\)。

\[\begin{aligned}f(n)&=\mu(n)&g(n)&=1&(f*g)(n)&=[n=1]\\ f(n)&=\varphi(n)&g(n)&=1&(f*g)(n)&=n\\ f(n)&=n\cdot\varphi(n)&g(n)&=n&(f*g)(n)&=n^2\\ f(n)&=d(n)&g(n)&=\mu(n)&(f*g)(n)&=1\\ f(n)&=\sigma(n)&g(n)&=\mu(n)&(f*g)(n)&=n\end{aligned} \]

（但是用公式 \(\sum_{i=1}^{n}\sigma(n)=\sum_{i=1}^{n}i\cdot\lfloor\dfrac n i\rfloor\) 更好）

  • \(O(n^{\tfrac 2 3})\)，注意有递归的操作就要记忆化
struct DU{
	static const int N=2000010;
	int sum[N];
	DU(){
		vector<int> a,mu(N,1),vis(N,0);
		repeat(i,2,N){
			if(!vis[i])a.push_back(i),mu[i]=-1;
			for(auto j:a){
				if(i*j>=N)break;
				vis[i*j]=1;
				if(i%j==0){mu[i*j]=0; break;}
				mu[i*j]=-mu[i];
			}
		}
		repeat(i,1,N)sum[i]=sum[i-1]+mu[i];
	}
	ll sum_mu(ll n){
		if(n<N)return sum[n];
		static map<ll,ll> rec; if(rec.count(n))return rec[n];
		ll ans=1;
		for(ll l=2,r;l<=n;l=r+1){
			r=n/(n/l);
			ans-=sum_mu(n/l)*(r-l+1);
		}
		return rec[n]=ans;
	}
	ll sum_phi(ll n){
		ll ans=0;
		for(ll l=1,r;l<=n;l=r+1){
			r=n/(n/l);
			ans+=(sum_mu(r)-sum_mu(l-1))*(n/l)*(n/l);
		}
		return ((ans-1)>>1)+1;
	}
	ll sum_d(ll n){
		static map<ll,ll> rec; if(rec.count(n))return rec[n];
		ll ans=n;
		for(ll l=2,r;l<=n;l=r+1){
			r=n/(n/l);
			ans-=(sum_mu(r)-sum_mu(l-1))*sum_d(n/l);
		}
		return rec[n]=ans;
	}
	ll sum_sigma(ll n){
		ll ans=0;
		for(ll l=1,r;l<=n;l=r+1){
			r=n/(n/l);
			ans+=(l+r)*(r-l+1)/2*(n/l);
		}
		return ans;
	}
}du;

min_25 筛

  • 学不会。
  • 求 \([1,n]\) 内的素数个数。
#include<cstdio>
#include<math.h>
#define ll long long
const int N = 316300;
ll n, g[N<<1], a[N<<1];
int id, cnt, sn, prime[N];
inline int Id(ll x){return x<=sn?x:id-n/x+1;}
int main() {
	scanf("%lld", &n), sn=sqrt(n);
	for(ll i=1; i<=n; i=a[id]+1) a[++id]=n/(n/i), g[id]=a[id]-1;
	for(int i=2; i<=sn; ++i) if(g[i]!=g[i-1]){
		// 这里 i 必然是质数，因为 g[] 是前缀质数个数
		// 当 <i 的质数的倍数都被筛去，让 g[] 发生改变的位置只能是下一个质数
		// 别忘了 i<=sn 时，ID(i) 就是 i。
		prime[++cnt]=i;
		ll sq=(ll)i*i;
		for(int j=id; a[j]>=sq; --j) g[j]-=g[Id(a[j]/i)]-(cnt-1);
	}
	return printf("%lld\n", g[id]), 0;
}
  • 求 \([1,n]\) 内的素数之和。
namespace Min25 {
	int prime[N],id1[N],id2[N],flag[N],cnt,m;
	ll g[N],sum[N],a[N],T,n;
	int ID(ll x){return x<=T?id1[x]:id2[n/x];}
	ll getsum(ll x){return x*(x+1)/2-1;}
	ll f(ll x){return x;}
	void work(){
		T=sqrt(n+0.5);cnt=0;fill(flag,flag+T+1,0);m=0;
		for(int i=2;i<=T;i++){
			if(!flag[i]) prime[++cnt]=i,sum[cnt]=sum[cnt-1]+i;
			for(int j=1;j<=cnt && i*prime[j]<=T;j++){
				flag[i*prime[j]]=1;
				if(i%prime[j]==0) break;
			}
		}
		for(ll l=1;l<=n;l=n/(n/l)+1){
			a[++m]=n/l;
			if(a[m]<=T)id1[a[m]]=m;else id2[n/a[m]]=m;
			g[m]=getsum(a[m]);
		}
		for(int i=1;i<=cnt;i++)
			for(int j=1;j <= m && 1ll*prime[i]*prime[i]<=a[j];j++)
				g[j]=g[j]-1ll*prime[i]*(g[ID(a[j]/prime[i])]-sum[i-1]);
	}
	ll solve(ll x){
		if(x<=1) return x;
		return n=x,work(),g[ID(n)];
	}
}

素数约数相关

唯一分解 / 质因数分解

  • 用数组表示数字唯一分解式的素数的指数，如 \(50=[0,0,1,0,0,2,0,\ldots]\)。
  • 可以用来计算阶乘和乘除操作。
void fac(int a[],ll n){
	repeat(i,2,(int)sqrt(n)+2)
	while(n%i==0)a[i]++,n/=i;
	if(n>1)a[n]++;
}
  • set 维护版
struct fac{
	#define facN 1010
	ll a[facN]; set<ll> s; // 乘法就是multiset
	fac(){mst(a,0); s.clear();}
	void lcm(ll n){ // self=lcm(self,n)
		repeat(i,2,facN)
		if(n%i==0){
			ll cnt=0;
			while(n%i==0)cnt++,n/=i;
			a[i]=max(a[i],cnt); // 改成a[i]+=cnt就变成了乘法
		}
		if(n>1)s.insert(n);
	}
	ll value(){ // return self%mod
		ll ans=1;
		repeat(i,2,facN)
			if(a[i])ans=ans*qpow(i,a[i],mod)%mod;
		for(auto i:s)ans=ans*i%mod;
		return ans;
	}
}f;

素数判定 using Miller-Rabin

  • \(O(\log^3 n)\)
bool mr(ll x,ll b){
	ll k=x-1;
	while(k){
		ll cur=qpow(b,k,x);
		if(cur!=1 && cur!=x-1)return 0;
		if(k%2==1 || cur==x-1)return 1;
		k>>=1;
	}
	return 1;
}
bool isprime(ll x){
	if(x<2 || x==46856248255981ll)
		return 0;
	if(x<4 || x==61) // raw: if(x==2 || x==3 || x==7 || x==61 || x==24251)
		return 1;
	return mr(x,2) && mr(x,61); // 2,3,5,7,13,29,37,89
}

大数分解 using Pollard-rho

  • \(O(n^{\tfrac 1 4})\)，基于 Miller-Rabin 素性测试
ll pollard_rho(ll x){
	ll s=0,t=0,c=rnd()%(x-1)+1;
	int stp=0,goal=1; ll val=1;
	for(goal=1;;goal<<=1,s=t,val=1){
		for(stp=1;stp<=goal;++stp){
			t=((__int128)t*t+c)%x;
			val=(__int128)val*abs(t-s)%x;
			if(stp%127==0){
				ll d=__gcd(val,x);
				if(d>1)return d;
			}
		}
		ll d=__gcd(val,x);
		if(d>1)return d;
	}
}
vector<ll> ans; // result
void rho(ll n){
	if(isprime(n)){
		ans.push_back(n);
		return;
	}
	ll t;
	do{t=pollard_rho(n);}while(t>=n);
	rho(t);
	rho(n/t);
}

求约数 / 因数

  • \(O(\sqrt n)\)
int get_divisor(int n){
	int ans=0;
	for(int i=1;i<n;i=n/(n/(i+1)))
	if(n%i==0)
		ans++; // v.push_back(i);
	return ans+1; // v.push_back(n);
}
  • 小常数版（要求 \(n\le 10^7\)），基于线性筛
vector<pii> pd; vector<ll> v; // pd: <k, p>; v: divisors
void dfs(int x,int y,ll s){
	if(x==(int)pd.size()){v.push_back(s); return;}
	dfs(x+1,0,s);
	if(y<pd[x].se)dfs(x,y+1,s*pd[x].fi);
}
void get_divisor(ll n){
	pd.clear(); v.clear();
	while(n!=1){
		if(!pd.empty() && pd.back().fi==lpf[n])pd.back().se++;
		else pd.push_back({lpf[n],1});
		n/=lpf[n]; // needs initialized
	}
	dfs(0,0,1);
}

组合数学

组合数取模 using Lucas+extra

  • Lucas定理用来求模意义下的组合数。
  • 真·Lucas，p 是质数。
ll lucas(ll a,ll b,ll p){ // a>=b
	if(b==0)return 1;
	return mul(C(a%p,b%p,p),lucas(a/p,b/p,p),p);
}
  • 特例：如果p=2，可能lucas失效。（？）
ll C(ll a,ll b){ // a>=b，p=2的情况
	return (a&b)==b;
}
  • 快速阶乘和 exLucas
  • qfac.A(x),qfac.B(x) 满足 \(A\equiv \dfrac{x!}{p^B}\pmod {p^k}\)。
  • qfac.C(a,b) \(\equiv C_a^b \pmod {p^k}\)。
  • \(\text{exlucas}(a,b,m)\equiv C_a^b \pmod m\)，函数内嵌中国剩余定理。
struct Qfac{
	ll s[2000010];
	ll p,m;
	ll A(ll x){ // 快速阶乘的A值
		if(x==0)return 1;
		ll c=A(x/p);
		return s[x%m]*qpow(s[m],x/m,m)%m*c%m;
	}
	ll B(ll x){ // 快速阶乘的B值
		int ans=0;
		for(ll i=x;i;i/=p)ans+=i/p;
		return ans;
	}
	ll C(ll a,ll b){ // 组合数，a>=b
		ll k=B(a)-B(b)-B(a-b);
		return A(a)*gcdinv(A(b),m)%m
			*gcdinv(A(a-b),m)%m
			*qpow(p,k,m)%m;
	}
	void init(ll _p,ll _m){ // 一定要满足m=p^k
		p=_p,m=_m;
		s[0]=1;
		repeat(i,1,m+1)
			if(i%p)s[i]=s[i-1]*i%m;
			else s[i]=s[i-1];
	}
}qfac;
ll exlucas(ll a,ll b,ll mod){
	ll ans=0,m=mod;
	for(ll i=2;i<=m;i++) // 不能repeat
	if(m%i==0){
		ll p=i,k=1;
		while(m%i==0)m/=i,k*=i;
		qfac.init(p,k);
		ans=(ans+qfac.C(a,b)*(mod/k)%mod*gcdinv(mod/k,k)%mod)%mod;
	}
	return (ans+mod)%mod;
}

博弈论

SG 定理

  • 有向无环图中，两个玩家轮流推多颗棋子，不能走的判负。
  • 假设 x 的后继状态为 \(y_1,y_2,...,y_k\)。
  • 则 \(SG[x]=mex\{SG[y_i]\}\)，\(mex(S)\) 表示不属于集合 S 的最小自然数。
  • 当且仅当所有起点SG值的异或和为 0 时先手必败。
  • （如果只有一个起点，SG的值可以只考虑01）。
  • 例题：拿 n 堆石子，每次只能拿一堆中的斐波那契数颗石子。
void getSG(int n){
	mst(SG,0);
	repeat(i,1,n+1){
		mst(S,0);
		for(int j=0;f[j]<=i && j<=N;j++)
			S[SG[i-f[j]]]=1;
		for(int j=0;;j++)
		if(!S[j]){
			SG[i]=j;
			break;
		}
	}
}

置换群

  • 定理：排列逆序对数奇偶性和 \(n-\) 环数奇偶性相同。
  • 求 \(A^x\)，编号从 0 开始，\(O(n)\)。
void qpow(int a[],int n,int x){
	static int rec[N],c[N];
	static bool vis[N];
	fill(vis,vis+n,0);
	repeat(i,0,n)if(!vis[i]){
		int cnt=0; rec[cnt++]=i;
		for(int p=a[i];p!=i;p=a[p])
			rec[cnt++]=p,vis[p]=1;
		repeat(J,0,cnt)
			c[rec[J]]=a[rec[(J+x-1)%cnt]];
		repeat(J,0,cnt)
			a[rec[J]]=c[rec[J]];
	}
}
  • \(A^k=B\) 求任一 A，编号从 0 开始，\(O(n)\)。（暂无判断有解操作）
repeat(i,0,n){
	a[read()-1]=i;
	vis[i]=0;
}
repeat(i,0,n)if(!vis[i]){
	int cnt=0; rec[cnt++]=i;
	for(int p=a[i];p!=i;p=a[p])
		rec[cnt++]=p,vis[p]=1;
	repeat(J,0,cnt)
		c[1ll*J*k%cnt]=rec[J];
	repeat(J,0,cnt)
		ans[c[(J+1)%cnt]]=c[J];
}

多项式

技能树：拉格朗日反演，快速插值和多点求值。

拉格朗日插值

  • 函数曲线通过n个点 \((x_i,y_i)\)，求 \(f(k)\)。
  • 拉格朗日插值：\(\displaystyle f(x)=\sum_{i=1}^n\left[y_i\prod_{j!=i}\dfrac{x-x_j}{x_i-x_j}\right]\)。
  • \(O(n^2)\)。
ll solve(int n,int x0){
	ll ans=0; x0%=mod;
	repeat(i,0,n)x[i]%=mod,y[i]%=mod;
	repeat(i,0,n){
		int s1=y[i],s2=1;
		repeat(j,0,n)
		if(i!=j){
			s1=s1*(x0-x[j])%mod;
			s2=s2*(x[i]-x[j])%mod;
		}
		ans=(ans+s1*qpow(s2,mod-2)%mod+mod)%mod;
	}
	return ans;
}

ll solve(int n,int x0){ // (i,y[i]),i=1..n的优化
	ll ans=0,up=1; x0%=mod;
	if(x0>=1 && x0<=n)return y[x0];
	repeat(i,1,n+1)
		up=up*(x0-i)%mod;
	repeat(i,1,n+1){
		ans+=y[i]*up%mod*qpow((x0-i)*((n+i)%2?-1:1)*C.fac[i-1]%mod*C.fac[n-i]%mod,mod-2)%mod;
	}
	return ans%mod;
}

快速傅里叶变换 / FFT + 任意模数

  • 离散傅里叶变换(DFT)即求 \((\omega_n^k,f(\omega_n^k))\)，多项式 \(\displaystyle d_k=\sum_{i=0}^{n-1}a_i(\omega_n^k)^i\)
  • 离散傅里叶反变换(IDFT)即求多项式 \((\omega_n^{-k},g(\omega_n^{-k}))\)，多项式 \(\displaystyle c_k=\sum_{i=0}^{n-1}d_i(\omega_n^{-k})^i\)，最后 \(a_i=\dfrac {c_i}{n}\)
  • 求两个多项式的卷积，\(O(n\log n)\)
struct FFT{
	static const int N=1<<20;
	struct cp{
		long double a,b;
		cp(){}
		cp(const long double &a,const long double &b):a(a),b(b){}
		cp operator+(const cp &t)const{return cp(a+t.a,b+t.b);}
		cp operator-(const cp &t)const{return cp(a-t.a,b-t.b);}
		cp operator*(const cp &t)const{return cp(a*t.a-b*t.b,a*t.b+b*t.a);}
		cp conj()const{return cp(a,-b);}
	};
	cp wn(int n,int f){
		static const long double pi=acos(-1.0);
		return cp(cos(pi/n),f*sin(pi/n));
	}
	int g[N];
	void dft(cp a[],int n,int f){
		repeat(i,0,n)if(i>g[i])swap(a[i],a[g[i]]);
		for(int i=1;i<n;i<<=1){
			cp w=wn(i,f);
			for(int j=0;j<n;j+=i<<1){
				cp e(1,0);
				for(int k=0;k<i;e=e*w,k++){
					cp x=a[j+k],y=a[j+k+i]*e;
					a[j+k]=x+y,a[j+k+i]=x-y;
				}
			}
		}
		if(f==-1){
			cp Inv(1.0/n,0);
			repeat(i,0,n)a[i]=a[i]*Inv;
		}
	}
	#ifdef CONV
	cp a[N],b[N];
	vector<ll> conv(const vector<ll> &u,const vector<ll> &v){ // 一般fft
		const int n=(int)u.size()-1,m=(int)v.size()-1;
		const int k=__lg(n+m+1)+1,s=1<<k;
		g[0]=0; repeat(i,1,s)g[i]=(g[i/2]/2)|((i&1)<<(k-1));
		repeat(i,0,s){
			a[i]=cp(i<=n?u[i]:0,0);
			b[i]=cp(i<=m?v[i]:0,0);
		}
		dft(a,s,1); dft(b,s,1);
		repeat(i,0,s)a[i]=a[i]*b[i];
		dft(a,s,-1);
		vector<ll> ans;
		repeat(i,0,n+m+1)ans<<llround(a[i].a);
		return ans;
	}
	#endif
	#ifdef CONV_MOD
	cp a[N],b[N],Aa[N],Ab[N],Ba[N],Bb[N];
	vector<ll> conv_mod(const vector<ll> &u,const vector<ll> &v,ll mod){ // 任意模数fft
		const int n=(int)u.size()-1,m=(int)v.size()-1,M=sqrt(mod)+1;
		const int k=__lg(n+m+1)+1,s=1<<k;
		g[0]=0; repeat(i,1,s)g[i]=(g[i/2]/2)|((i&1)<<(k-1));
		repeat(i,0,s){
			a[i]=i<=n?cp(u[i]%mod%M,u[i]%mod/M):cp();
			b[i]=i<=m?cp(v[i]%mod%M,v[i]%mod/M):cp();
		}
		dft(a,s,1); dft(b,s,1);
		repeat(i,0,s){
			int j=(s-i)%s;
			cp t1=(a[i]+a[j].conj())*cp(0.5,0);
			cp t2=(a[i]-a[j].conj())*cp(0,-0.5);
			cp t3=(b[i]+b[j].conj())*cp(0.5,0);
			cp t4=(b[i]-b[j].conj())*cp(0,-0.5);
			Aa[i]=t1*t3,Ab[i]=t1*t4,Ba[i]=t2*t3,Bb[i]=t2*t4;
		}
		repeat(i,0,s){
			a[i]=Aa[i]+Ab[i]*cp(0,1);
			b[i]=Ba[i]+Bb[i]*cp(0,1);
		}
		dft(a,s,-1); dft(b,s,-1);
		vector<ll> ans;
		repeat(i,0,n+m+1){
			ll t1=llround(a[i].a)%mod;
			ll t2=llround(a[i].b)%mod;
			ll t3=llround(b[i].a)%mod;
			ll t4=llround(b[i].b)%mod;
			ans+=(t1+(t2+t3)*M%mod+t4*M*M)%mod;
		}
		return ans;
	}
	#endif
}fft;

多项式全家桶 数组版

  • 多项式开根 \(g[0]=1\)，若 \(g[0]\) 是其他二次剩余则另需二次剩余板子。
  • 若 \(\ln f(x)\) 存在，则 \([x^0]f(x)=1\)；若 \(\exp f(x)\) 存在，则 \([x^0]f(x)=0\)。
void eachfac(ll a[],int n){ // ans[i]=a[i]*i!
	ll p=1;
	repeat(i,1,n)p=p*i%mod,a[i]=a[i]*p%mod;
}
void eachfacinv(ll a[],int n){ // ans[i]=a[i]/i!
	ll p=1;
	repeat(i,1,n)p=p*i%mod,a[i]=a[i]*qpow(p,mod-2)%mod;
}
void shift(ll a[],int n,int p){ // ans=a*x^p
	if(p>0){
		repeat_back(i,p,n)a[i]=a[i-p];
		repeat(i,0,p)a[i]=0;
	}
	if(p<0){
		repeat(i,0,n+p)a[i]=a[i-p];
		repeat(i,n+p,n)a[i]=0;
	}
}

const int mod=998244353;
inline ll D(ll x){return x>=mod?x-mod:x;}
#define ex(a) fill(a+n,a+n*2,0)
int polyinit(ll a[],int n1){ // a[0..n1-1], n>=n1, n=2^k
	int n=1; while(n<n1)n<<=1;
	fill(a+n1,a+n,0);
	return n;
}
void read(ll a[],int n){
	repeat(i,0,n)a[i]=read();
}
void print(const ll a[],int n){
	repeat(i,0,n)print(a[i],i==n-1);
}
void der(const ll a[],int n,ll b[]){ // n=2^k, b = d a(x) / dx
	repeat(i,1,n){b[i-1]=i*a[i]%mod;} b[n-1]=0;
}
void cal(const ll a[],int n,ll b[]){ // n=2^k, b = integral a(x) dx
	repeat_back(i,1,n){b[i]=qpow(i,mod-2,mod)*a[i-1]%mod;} b[0]=0;
}
void ntt(ll a[],ll n,ll op){ // n=2^k
	for(int i=1,j=n>>1;i<n-1;++i){
		if(i<j)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){
		ll rt=qpow(3,(mod-1)/len,mod);
		for(int i=0;i<n;i+=len){
			ll w=1;
			repeat(j,i,i+len/2){
				ll u=a[j],t=1ll*a[j+len/2]*w%mod;
				a[j]=D(u+t),a[j+len/2]=D(u-t+mod);
				w=1ll*w*rt%mod;
			}
		}
	}
	if(op==-1){
		reverse(a+1,a+n);
		ll in=qpow(n,mod-2,mod);
		repeat(i,0,n)a[i]=1ll*a[i]*in%mod;
	}
}
void conv(ll a[],ll b[],int n,ll c[],const function<ll(ll,ll)> &f=[](ll a,ll b){return a*b%mod;}){ // n=2^k, c=a*b
	ntt(a,n,1); if(b!=a)ntt(b,n,1);
	repeat(i,0,n)c[i]=f(a[i],b[i]);
	ntt(c,n,-1);
}
void inv(const ll a[],int n,ll b[]){ // n=2^k, a!=b, a*b=1
	static ll t[N];
	if(n==1){b[0]=qpow(a[0],mod-2); b[1]=0; return;}
	inv(a,n/2,b); copy(a,a+n,t);
	ex(b); ex(t);
	conv(b,t,n*2,b,[](ll a,ll b){
		return a*(2-a*b%mod+mod)%mod;
	});
	ex(b);
}
const int inv2=qpow(2,mod-2);
void sqrt(const ll a[],int n,ll b[]){ // n=2^k, a!=b, b*b=a
	static ll f[N],g[N];
	if(n==1){b[0]=1; /*sqrtmod(a[0])*/ return;}
	sqrt(a,n/2,b);
	inv(b,n,f);
	copy(a,a+n,g); ex(g);
	conv(g,f,n*2,g);
	repeat(i,0,n)b[i]=inv2*(g[i]+b[i])%mod;
}
void ln(const ll a[],int n,ll b[]){ // n=2^k, a!=b
	static ll t[N];
	der(a,n,t); inv(a,n,b); ex(t);
	conv(t,b,n*2,t);
	cal(t,n,b); ex(b);
}
void exp(const ll a[],int n,ll b[]){ // n=2^k, a!=b
	static ll t[N];
	if(n==1){b[0]=1; return;}
	exp(a,n/2,b); ln(b,n,t);
	repeat(i,0,n){t[i]=D(a[i]-t[i]+mod);} t[0]++;
	ex(b);
	conv(b,t,n*2,b);
}
void divmod(const ll a[],const ll b[],int n,int m,ll d[],ll r[]){ // n=2^k, m<n, |a|=n, |b|=m, |d|=n-m+1, |r|=m-1, a=d*b+r
	#define er(a,m) fill(a+m,a+n*2,0)
	static ll f[N],g[N];
	reverse_copy(b,b+m,f); er(f,m);
	inv(f,n,g);
	reverse_copy(a,a+n,f);
	conv(f,g,n*2,d);
	reverse(d,d+n-m+1); er(d,n-m+1);
	copy(d,d+n*2,f);
	copy(b,b+n,g); er(g,n);
	conv(f,g,n*2,r);
	repeat(i,0,n*2)r[i]=D(a[i]-r[i]+mod);
}
const ll im=911660635; // im = sqrtmod(-1), imaginary number
namespace tri{
	ll f[N],g[N];
	void getexp(const ll a[],int n){ // n=2^k, f=exp(ia), g=exp(-ia)
		repeat(i,0,n)g[i]=a[i]*im%mod;
		exp(g,n,f);
		inv(f,n,g);
	}
	void sin(const ll a[],int n,ll b[]){ // n=2^k
		getexp(a,n);
		repeat(i,0,n)b[i]=D(f[i]-g[i]+mod)*inv2%mod*(mod-im)%mod;
	}
	void cos(const ll a[],int n,ll b[]){ // n=2^k
		getexp(a,n);
		repeat(i,0,n)b[i]=D(f[i]+g[i])*inv2%mod;
	}
	void tan(const ll a[],int n,ll b[]){ // n=2^k
		getexp(a,n);
		repeat(i,0,n)
			tie(f[i],g[i])=make_pair(
				D(f[i]-g[i]+mod)*inv2%mod*(mod-im)%mod,
				D(f[i]+g[i])*inv2%mod
			);
		inv(g,n,b); ex(f);
		conv(f,b,n*2,b);
	}
	void asin(const ll a[],int n,ll b[]){ // n=2^k
		der(a,n,f); ex(f);
		copy(a,a+n,g); ex(g);
		conv(g,g,n*2,g,[](ll a,ll b){
			return D(1-a*b%mod+mod);
		});
		sqrt(g,n,b); inv(b,n,g);
		conv(f,g,n*2,g); cal(g,n,b);
	}
	void acos(const ll a[],int n,ll b[]){ // n=2^k
		asin(a,n,b);
		repeat(i,0,n)b[i]=D(mod-b[i]);
	}
	void atan(const ll a[],int n,ll b[]){ // n=2^k
		der(a,n,f); ex(f);
		copy(a,a+n,g); ex(g);
		conv(g,g,n*2,g,[](ll a,ll b){
			return D(1+a*b%mod);
		});
		inv(g,n,b); conv(f,b,n*2,g);
		cal(g,n,b);
	}
}
ll getmod(const char s[],int mod){ // ans=s%mod
	ll ans=0;
	repeat(i,0,strlen(s))ans=(ans*10+s[i]-'0')%mod;
	return ans;
}
void qpow_trivial(const ll a[],ll p,int n,ll b[]){ // n=2^k, b=a^p
	static ll t[N];
	ln(a,n,t);
	repeat(i,0,n)(t[i]*=p)%=mod;
	exp(t,n,b);
}
void qpow(const ll a[],char s[],int n,ll b[]){ // n=2^k, b=a^s
	static ll f[N],g[N];
	ll m=getmod(s,mod),m1=getmod(s,mod-1);
	ll d=0; while(d<n && a[d]==0)d++;
	if(d*m>=n || (d && strlen(s)>=8)){
		fill(b,b+n,0); return;
	}
	int in=qpow(a[d],mod-2,mod),owe=qpow(a[d],m1,mod);
	repeat(i,0,n-d)f[i]=a[i+d]*in%mod; fill(f+n-d,f+n,0);
	qpow_trivial(f,m,n,g); d*=m;
	repeat(i,0,d)b[i]=0;
	repeat_back(i,d,n)b[i]=g[i-d]*owe%mod;
}

快速沃尔什变换 / FWT

  • 计算 \(\displaystyle c_i=\sum_{i=f(j,k)}a_jb_k\)，\(O(n\log n)\)。
void fwt(ll a[],int n,int flag,char c){ // flag = -1 / 1
	if(c=='|'){
		for(int w=1;w<n;w<<=1)
		for(int i=0;i<n;i+=w*2)
		repeat(j,0,w)
			ad(a[i+j+w]+=a[i+j]*flag);
	}
	else if(c=='&'){
		for(int w=1;w<n;w<<=1)
		for(int i=0;i<n;i+=w*2)
		repeat(j,0,w)
			ad(a[i+j]+=a[i+j+w]*flag);
	}
	else if(c=='^'){
		if(flag==-1)flag=qpow(2,mod-2,mod);
		for(int w=1;w<n;w<<=1)
		for(int i=0;i<n;i+=w*2)
		repeat(j,0,w){
			ad(a[i+j]+=a[i+j+w]);
			ad(ad(a[i+j+w]=a[i+j]-a[i+j+w]*2));
			(a[i+j]*=flag)%=mod;
			(a[i+j+w]*=flag)%=mod;
		}
	}
}
void bitmul(ll a[],ll b[],int n,char c){ // n=2^k
	fwt(a,n,1,c); fwt(b,n,1,c);
	repeat(i,0,n)a[i]=(a[i]*b[i])%mod;
	fwt(a,n,-1,c);
}

多项式复合

  • 求 \(H(x)\equiv F(G(x)) (\bmod x^{n+1})\)，\(n=20000\)，有点卡常
ll RT[17][N];
struct INIT{INIT(){ // 预处理原根
	for(int i=0;i<17;++i){
		ll *G=RT[i]; G[0]=1;
		const int gi=G[1]=qpow(3,(mod-1)/(1<<i+1));
		for(int j=2;j<1<<i;++j)G[j]=G[j-1]*gi%mod;
	}
}}Init; 
void ntt(ll a[],ll n,ll op){
	for(int i=1,j=n>>1;i<n-1;++i){
		if(i<j)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){
		const ll *G=RT[__builtin_ctz(len)-1];
		for(int i=0;i<n;i+=len){
			ll w=1;
			repeat(j,i,i+len/2){
				ll u=a[j],t=1ll*a[j+len/2]*G[j-i]%mod;
				a[j]=D(u+t),a[j+len/2]=D(u-t+mod);
			}
		}
	}
	if(op==-1){
		reverse(a+1,a+n);
		ll in=qpow(n,mod-2,mod);
		repeat(i,0,n)a[i]=1ll*a[i]*in%mod;
	}
}
void conv(ll a[],ll b[],int n,ll c[],const function<ll(ll,ll)> &f=[](ll a,ll b){return a*b%mod;}){ // n=2^k, c=a*b
	ntt(a,n,1); if(b!=a)ntt(b,n,1);
	repeat(i,0,n)c[i]=f(a[i],b[i]);
	ntt(c,n,-1);
}
const int L=142; // sqrt(n1)
ll f[N],g[L+1][N],ng[L+1][N],G[L+1][N],nG[L+1][N];
void prework(ll g[][N],ll ng[][N],int n){
	#define cpy(a,b) copy(a,a+n,b)
	n*=2; g[0][0]=1;
	static ll e[N]; cpy(g[1],e); ntt(e,n,1);
	repeat(i,1,L+1){
		cpy(g[i-1],ng[i-1]); ntt(ng[i-1],n,1);
		repeat(j,0,n)g[i][j]=e[j]*ng[i-1][j]%mod;
		ntt(g[i],n,-1); fill(g[i]+n/2,g[i]+n,0);
	}
}
void Solve(){
	int n1=read()+1,m1=read()+1;
	read(f,n1); int n=polyinit(f,max(n1,m1));
	read(g[1],m1); prework(g,ng,n);
	cpy(g[L],G[1]); prework(G,nG,n);
	static ll ans[N];
	repeat(i,0,L){
		static ll s[N]; fill(s,s+n*2,0);
		repeat(j,0,L){
			int x=i*L+j;
			if(x<n1){
				repeat(k,0,n1)
					(s[k]+=f[x]*g[j][k])%=mod;
			}
		}
		ntt(s,n*2,1);
		repeat(j,0,n*2)s[j]=s[j]*nG[i][j]%mod;
		ntt(s,n*2,-1);
		repeat(k,0,n1)ad(ans[k]+=s[k]);
	}
	print(ans,n1);
}

多项式多点求值

  • 已知多项式 f 和序列 a，求 \(f(a_1),f(a_2),\ldots,f(a_m)\)
  • 线性算法指输入 n 维向量 x，经过 \(m\times n\) 矩阵 A 变换后输出 m 维向量 \(y=Ax\) 的算法
  • 转置原理指出，如果存在 \(x'=A^Ty'\) 的算法，那么就有存在相同复杂度的 \(y=Ax\) 的算法。将 \(A^T\) 分解为三种指令 x[i]+=x[j],x[i]*=c,swap(x[i],x[i])，那么 A 即倒着执行这些指令，并且将第一种指令变为 x[j]+=x[i]。（\(A=E_1E_2\ldots E_k\rightarrow A^T=E_k^TE_{k-1}^T\ldots E_1^T\)）
  • \(O(n\log^2n)\)，常数极大（1.8s 跑 64000）
ll D(ll x){return x>=mod?x-mod:x<0?x+mod:x;}
ll &ad(ll &x){return x=D(x);}
typedef vector<ll> vi;
#define rs(a) [&]{if((int)a.size()<n)a.resize(n,0);}()
#define cut(a) fill(a.begin()+n/2,a.begin()+n,0)
int polyn(int n1){ // return 2^k >= n1
	return 1<<(__lg(n1-1)+1);
}
void ntt(vi &a,ll n,ll op){ // n=2^k
	rs(a);
	for(int i=1,j=n>>1;i<n-1;++i){
		if(i<j)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){
		ll rt=qpow(3,(mod-1)/len,mod);
		for(int i=0;i<n;i+=len){
			ll w=1;
			repeat(j,i,i+len/2){
				ll u=a[j],t=1ll*a[j+len/2]*w%mod;
				a[j]=D(u+t),a[j+len/2]=D(u-t+mod);
				w=1ll*w*rt%mod;
			}
		}
	}
	if(op==-1){
		reverse(a.begin()+1,a.begin()+n);
		ll in=qpow(n,mod-2,mod);
		repeat(i,0,n)a[i]=1ll*a[i]*in%mod;
	}
}
vi conv(vi a,vi b,int n,const function<ll(ll,ll)> &f=[](ll a,ll b){return a*b%mod;}){ // n=2^k, ans=a*b
	n*=2; rs(a),rs(b); cut(a),cut(b);
	ntt(a,n,1); ntt(b,n,1);
	repeat(i,0,n)a[i]=f(a[i],b[i]);
	ntt(a,n,-1); cut(a);
	return a;
}
vi inv(const vi &a,int n){ // n=2^k, ans=1/a
	if(n==1)return vi(1,qpow(a[0],mod-2,mod));
	return conv(inv(a,n/2),a,n,[](ll a,ll b){
		return a*(2-a*b%mod+mod)%mod;
	});
}
vi convauto(vi a,vi b,const function<ll(ll,ll)> &f=[](ll a,ll b){return a*b%mod;}){
	int n1=a.size()+b.size()-1,n=polyn(n1);
	rs(a),rs(b);
	ntt(a,n,1); ntt(b,n,1);
	repeat(i,0,n)a[i]=f(a[i],b[i]);
	ntt(a,n,-1);
	a.resize(n1);
	return a;
}
vi convtr(vi a,vi b,const function<ll(ll,ll)> &f=[](ll a,ll b){return a*b%mod;}){
	int n1=a.size()+b.size()-1,n=polyn(n1);
	rs(a),rs(b);
	reverse(a.begin(),a.end());
	ntt(a,n,1); ntt(b,n,1);
	repeat(i,0,n)a[i]=f(a[i],b[i]);
	ntt(a,n,-1);
	reverse(a.begin(),a.end());
	return a;
}
vi prod[N]; ll ans[N],a[N];
#define lc x*2
#define rc x*2+1
void getprod(int x,int l,int r){
	if(l==r){
		prod[x]={1,D(mod-a[l])};
		return;
	}
	int mid=(l+r)/2;
	getprod(lc,l,mid);
	getprod(rc,mid+1,r);
	prod[x]=convauto(prod[lc],prod[rc]);
}
void dfs(int x,int l,int r,vi G){
	G.resize(r-l+1);
	if(l==r){
		ans[l]=G[0];
		return;
	}
	int mid=(l+r)/2;
	dfs(lc,l,mid,convtr(G,prod[rc]));
	dfs(rc,mid+1,r,convtr(G,prod[lc]));
}
void Solve(){
	int n=read()+1,m=read();
	vi f; repeat(i,0,n)f<<read(); // poly
	repeat(i,0,m)
		a[i]=read(); // query
	getprod(1,0,m-1);
	vi v=inv(prod[1],polyn(m+1));
	dfs(1,0,m-1,convtr(f,v));
	repeat(i,0,m)
		print(ans[i],1);
}

多项式快速插值

  • \(O(n\log^2 n)\)，常数极大（2.85s 跑 100000）
const int N=(1<<18)+5;
int D(int x){return x>=mod?x-mod:x;}
int r[19][N],w[2][N],lg[N],inv[19];
void Pre(){
	repeat(d,1,18+1){
		repeat(i,1,(1<<d))r[d][i]=(r[d][i>>1]>>1)|((i&1)<<(d-1));
		lg[1<<d]=d,inv[d]=qpow(1<<d,mod-2);
	}
	for(int t=(mod-1)>>1,i=1,x,y; i<262144; i<<=1,t>>=1){
		x=qpow(3,t),y=qpow(332748118,t),w[0][i]=w[1][i]=1;
		repeat(k,1,i){
			w[1][k+i]=mul(w[1][k+i-1],x);
		    w[0][k+i]=mul(w[0][k+i-1],y);
		}
	}
}
int lim,e,n,m;
void init(int len){
	lim=1,e=0;
	while(lim<len)lim<<=1,++e;
}
void NTT(int *a,int ty){
	repeat(i,0,lim)if(i<r[e][i])swap(a[i],a[r[e][i]]);
	for(int mid=1; mid<lim; mid<<=1)
	for(int j=0,t; j<lim; j+=(mid<<1))
	repeat(k,0,mid){
		t=mul(w[ty][mid+k],a[j+k+mid]);
		a[j+k+mid]=D(a[j+k]-t+mod),a[j+k]=D(a[j+k]+t);
	}
	if(!ty)repeat(i,0,lim)a[i]=mul(a[i],inv[e]);
}
void Inv(int *a,int *b,int len){
	if(len==1)return b[0]=qpow(a[0],mod-2),void();
	Inv(a,b,len>>1),lim=(len<<1),e=lg[lim];
	static int c[N],d[N];
	repeat(i,0,len)c[i]=a[i],d[i]=b[i];
	repeat(i,len,lim)c[i]=d[i]=0;
	NTT(c,1),NTT(d,1);
	repeat(i,0,lim)c[i]=mul(c[i],mul(d[i],d[i]));
	NTT(c,0);
	repeat(i,0,len)b[i]=D(D(b[i]*2)-c[i]+mod);
	repeat(i,len,lim)b[i]=0;
}
struct node {
	node *lc,*rc;
	vector<int> vec;
	int deg;
	void Mod(const int *a,int *r,int n){
		static int b[N],c[N],d[N];
		int len=1;
		while(len<=n-deg)len<<=1;
		repeat(i,0,n+1)b[i]=a[n-i];
		repeat(i,0,deg+1)c[i]=vec[deg-i];
		repeat(i,n-deg+1,len)c[i]=0;
		Inv(c,d,len);
		lim=(len<<1),e=lg[lim];
		repeat(i,n-deg+1,lim)b[i]=d[i]=0;
		NTT(b,1),NTT(d,1);
		repeat(i,0,lim)b[i]=mul(b[i],d[i]);
		NTT(b,0);
		reverse(b,b+n-deg+1);
		init(n+1);
		repeat(i,n-deg+1,lim)b[i]=0;
		repeat(i,0,deg+1)c[i]=vec[i];
		repeat(i,deg+1,lim)c[i]=0;
		NTT(b,1),NTT(c,1);
		repeat(i,0,lim)b[i]=mul(b[i],c[i]);
		NTT(b,0);
		repeat(i,0,deg)r[i]=D(a[i]-b[i]+mod);
	}
	void Mul(){
		static int a[N],b[N];
		deg=lc->deg+rc->deg,vec.resize(deg+1),init(deg+1);
		repeat(i,0,lc->deg+1)a[i]=lc->vec[i];
		repeat(i,lc->deg+1,lim)a[i]=0;
		repeat(i,0,rc->deg+1)b[i]=rc->vec[i];
		repeat(i,rc->deg+1,lim)b[i]=0;
		NTT(a,1),NTT(b,1);
		repeat(i,0,lim)a[i]=mul(a[i],b[i]);
		NTT(a,0);
		repeat(i,0,deg+1)vec[i]=a[i];
	}
} pool[N],*rt;
struct bnode {
	bnode *lc,*rc;
	vector<int> vec;
	int l,r;
	void Mul(node* p){
		static int a[N],b[N],c[N],d[N];
		int mid=(l+r)>>1;
		init(r-l+1+1);
		repeat(i,0,mid-l+1+1)a[i]=lc->vec[i],c[i]=p->lc->vec[i];
		repeat(i,mid-l+2,lim)a[i]=c[i]=0;
		repeat(i,0,r-mid+1)b[i]=rc->vec[i],d[i]=p->rc->vec[i];
		repeat(i,r-mid+1,lim)b[i]=d[i]=0;
		NTT(a,1),NTT(b,1),NTT(c,1),NTT(d,1);
		repeat(i,0,lim)a[i]=D(mul(a[i],d[i])+mul(b[i],c[i]));
		NTT(a,0);
		vec.resize(r-l+2);
		repeat(i,0,r-l+1+1)vec[i]=a[i];
	}
} o[N],*brt;
int a[N],tot,cnt;
node *newnode(){
	return &pool[tot++];
}
bnode *newbnode(){
	return &o[cnt++];
}
void solve(node *&p,int l,int r){
	p=newnode();
	if(l==r)return p->deg=1,p->vec.resize(2),p->vec[0]=mod-a[l],p->vec[1]=1,void();
	int mid=(l+r)>>1;
	solve(p->lc,l,mid),solve(p->rc,mid+1,r);
	p->Mul();
}
int b[25],f[N];
void calc(node *p,int l,int r,const int *d){
	if(r-l<=512){
		repeat(i,l,r+1){
			int x=a[i],c1,c2,c3,c4,now=d[r-l];
			b[0]=1;
			repeat(j,1,16+1)b[j]=mul(b[j-1],x);
			for(int j=r-l-1; j-15>=0; j-=16){
				c1=(1ll*now*b[16]+1ll*d[j]*b[15]+1ll*d[j-1]*b[14]+1ll*d[j-2]*b[13])%mod,
				c2=(1ll*d[j-3]*b[12]+1ll*d[j-4]*b[11]+1ll*d[j-5]*b[10]+1ll*d[j-6]*b[9])%mod,
				c3=(1ll*d[j-7]*b[8]+1ll*d[j-8]*b[7]+1ll*d[j-9]*b[6]+1ll*d[j-10]*b[5])%mod,
				c4=(1ll*d[j-11]*b[4]+1ll*d[j-12]*b[3]+1ll*d[j-13]*b[2]+1ll*d[j-14]*b[1])%mod,
				now=(0ll+c1+c2+c3+c4+d[j-15])%mod;
			}
			repeat_back(j,0,(r-l)%16)now=(1ll*now*x+d[j])%mod;
			f[i]=now;
		}
		return;
	}
	int mid=(l+r)>>1,b[p->deg+1];
	p->lc->Mod(d,b,p->deg-1),calc(p->lc,l,mid,b);
	p->rc->Mod(d,b,p->deg-1),calc(p->rc,mid+1,r,b);
}
int x[N],y[N],A[N];
void loli(bnode *&u,node *p,int l,int r){
	u=newbnode(),u->l=l,u->r=r;
	if(l==r)return u->vec.resize(2),u->vec[0]=mul(y[l],qpow(f[l],mod-2)),u->vec[1]=0,void();
	int mid=(l+r)>>1;
	loli(u->lc,p->lc,l,mid),loli(u->rc,p->rc,mid+1,r);
	u->Mul(p);
}
int main(){
	n=read(),Pre();
	repeat(i,1,n+1)x[i]=a[i]=read(),y[i]=read();
	solve(rt,1,n);
	repeat(i,1,n+1)A[i-1]=mul(rt->vec[i],i);
	A[n]=0;
	calc(rt,1,n,A);
	loli(brt,rt,1,n);
	repeat(i,0,n)print(brt->vec[i]);
	return 0;
}

k 进制异或卷积 / k 进制 FWT

  • 令 \(a \oplus_k b\) 为 k 进制异或（k 进制无进位加法）
  • 求 \(\displaystyle c_m=\sum_{i\oplus_k j=m}a_ib_j\)
const int k=7;
const ll rt=qpow(3,(mod-1)/k);
void kfwt(ll a[],int n,int op){ // n = k^(int), op = -1 / 1
	static ll t[k],w[k];
	w[0]=1; repeat(i,1,k)w[i]=w[i-1]*rt%mod;
	for(int len=1;len<n;len*=k)
	for(int i=0;i<n;i+=len*k)
	repeat(j,i,i+len){
		repeat(x,0,k){
			t[x]=0;
			repeat(y,0,k)
				(t[x]+=a[j+y*len]*w[y*(k+op*x)%k])%=mod;
		}
		repeat(x,0,k)a[j+x*len]=t[x];
	}
	if(op==-1){
		ll inv=qpow(n,mod-2);
		repeat(i,0,n)a[i]=a[i]*inv%mod;
	}
}

任意长度 FFT using Bluestein 算法

  • 编号从 0 开始，\(O(n\log n)\)。
typedef complex<lf> cp;
void DFT(cp a[],int n,int op){ // op=-1 时不是真正的 IDFT
	static int r[N]{};
	repeat(i,1,n){
		r[i]=(r[i>>1]>>1)|(i&1?n>>1:0);
		if(i<r[i])swap(a[i],a[r[i]]);
	}
	for(int len=2;len<=n;len<<=1){
		cp w=cp(cos(pi*2/len),op*sin(pi*2/len));
		for(int i=0;i<n;i+=len){
			cp d(1,0);
			repeat(j,i,i+len/2){
				cp x=a[j],y=a[j+len/2]*d;
				a[j]=x+y; a[j+len/2]=x-y;
				d=d*w;
			}
		}
	}
}
void bluestein(cp a[],int n,int op){ // n is arbitrary, op in {-1,1}
	static cp t[N],u[N];
	int k=1;
	while(k<4*n)k<<=1;
	repeat(i,0,k)t[i]=u[i]=0;
	repeat(i,0,n)
		t[i]=a[i]*cp(cos(pi*i*i/n),op*sin(pi*i*i/n));
	repeat(i,0,n*2)
		u[i]=cp(cos(pi*(i-n)*(i-n)/n),-op*sin(pi*(i-n)*(i-n)/n));
	DFT(t,k,1); DFT(u,k,1);
	repeat(i,0,k)t[i]=t[i]*u[i];
	DFT(t,k,-1);
	repeat(i,0,n)
		a[i]=t[i+n]*cp(cos(pi*i*i/n),op*sin(pi*i*i/n))/lf(k);
}

递推式插值 using BM 算法

  • 已知数列前几项，求递推式系数 \(C_0a_i+C_1a_{i+1}+...+C_ka_{i+k}=0,C_k=-1\)
  • 用来找规律
using vtr = vector<ll>;
vtr mul(const vtr &a, const vtr &b, const vtr &m, int k) {
	vtr r(2 * k - 1);
	repeat(i, 0, k) repeat(j, 0, k) (r[i + j] += a[i] * b[j]) %= mod;
	repeat_back(i, 0, k - 1) {
		repeat(j, 0, k) (r[i + j] += r[i + k] * m[j]) %= mod;
		r.pop_back();
	}
	return r;
}
vtr qpow(const vtr &m, ll n) {
	int k = (int)m.size() - 1;
	vtr r(k), x(k);
	r[0] = x[1] = 1;
	for (; n; n >>= 1, x = mul(x, x, m, k))
		if (n & 1) r = mul(x, r, m, k);
	return r;
}
ll go(const vtr &a, const vtr &x, ll n) {
	int k = (int)a.size() - 1;
	if (n <= k) return x[n - 1];
	if (a.size() == 2) return x[0] * qpow(a[0], n - 1, mod) % mod;
	vtr r = qpow(a, n - 1);
	ll ans = 0;
	repeat(i, 0, k) (ans += r[i] * x[i]) %= mod;
	return (ans + mod) % mod;
}
vtr BM(const vtr &x) {
	vtr C = {-1}, B = {-1};
	ll L = 0, m = 1, b = 1;
	repeat(n, 0, x.size()) {
		ll d = 0;
		repeat(i, 0, L + 1) (d += C[i] * x[n - i]) %= mod;
		if (d == 0) { m++; continue; }
		vtr T = C;
		ll c = mod - d * qpow(b, mod - 2, mod) % mod;
		repeat(i, C.size(), B.size() + m) C.push_back(0);
		repeat(i, 0, B.size()) (C[i + m] += c * B[i]) %= mod;
		if (2 * L > n) { m++; continue; }
		L = n + 1 - L; B.swap(T); b = d; m = 1;
	}
	reverse(C.begin(), C.end());
	return C;
}

分治 FFT

  • 比如 \(\displaystyle f[i]=\sum_{j=1}^if[i-j]g[i]\)，把要求的多项式分成两边，先算 \(f[0..n-1]\) 对自己的贡献（此时 \(f[0..n-1]\) 已确定），然后算 \(f[0..n-1]\) 对 \(f[n..2n-1]\) 的贡献，再算 \(f[n..2n-1]\) 对自己的贡献，\(O(n\log^2n)\)
int g[N],f[N],GG[N],FF[N];
void work(int l,int r){
	if(r-l==1)return;
	int m=(l+r)/2;
	work(l,m);
	copy(g,g+r-l,GG);
	copy(f+l,f+m,FF+l); fill(FF+m,FF+r,0);
	conv(GG,FF+l,r-l,FF+l);
	repeat(i,m,r)ad(f[i]+=FF[i]);
	work(m,r);
}
// int n=polyinit(g,n1); fill(f,f+n,0); f[0]=1; work(0,n);
  • 卡特兰数 \(\displaystyle S_n=\sum_{k=0}^{n-1}S_kS_{n-1-k}\)
int f[N],GG[N],FF[N];
void work(int l,int r){
	if(r-l==1)return;
	int m=(l+r)/2;
	work(l,m);
	copy(f,f+r-l,GG+1); GG[0]=0;
	copy(f+l,f+m,FF+l); fill(FF+m,FF+r,0);
	conv(GG,FF+l,r-l,FF+l);
	int t=1+(l!=0);
	repeat(i,m,r)(f[i]+=FF[i]*t)%=mod;
	work(m,r);
}
// fill(f,f+n,0); f[0]=1; work(0,n);
  • 超级卡特兰数 \(\displaystyle S_n=S_{n-1}+\sum_{k=0}^{n-1}S_kS_{n-1-k}\)
ll f[N],GG[N],FF[N],f0[N];
void work(int l,int r){
	if(r-l==1)return;
	int m=(l+r)/2;
	work(l,m);
	copy(f,f+r-l,GG+1); GG[0]=0;
	copy(f+l,f+m,FF+l); fill(FF+m,FF+r,0);
	conv(GG,FF+l,r-l,FF+l);
	int t=1+(l!=0);
	repeat(i,m,r)(f0[i]+=f0[i-1]+FF[i]*t)%=mod;
	ad(f0[r]+=f0[r-1]);
	repeat(i,m,r)ad(f[i]+=f0[i]),f0[i]=0;
	work(m,r);
}
// fill(f,f+n,0); fill(f0,f0+n,0); f[0]=f0[0]=1; work(0,n);

第二类斯特林数·行

  • \(\displaystyle S(n,r)=[x^r](\sum_{i=0}^n\dfrac{(-1)^i}{i!}x^i)(\sum_{i=0}^{n}\dfrac{i^n}{i!}x^i)\)
repeat(i,0,n1+1){
	a[i]=C.inv[i]; if(i%2==1)a[i]=-a[i];
	b[i]=qpow(i,n1)*C.inv[i]%mod;
}
int n=polyinit(a,n1+1); polyinit(b,n1+1);
conv(a,b,n,a);

矩阵

矩阵乘法 and 矩阵快速幂

  • 已并行优化，矩乘 \(O(n^3)\)，矩快 \(O(n^3\log m)\)
struct mat{
	static const int N=110;
	ll a[N][N];
	explicit mat(ll e=0){
		repeat(i,0,n)
		repeat(j,0,n)
			a[i][j]=e*(i==j);
	}
	mat operator*(const mat &b)const{
		mat ans(0);
		repeat(i,0,n)
		repeat(k,0,n){
			ll t=a[i][k];
			repeat(j,0,n)
				(ans.a[i][j]+=t*b.a[k][j])%=mod;
		}
		return ans;
	}
	vector<ll> operator*(const vector<ll> &b)const{
		vector<ll> ans(n);
		repeat(i,0,n)
		repeat(j,0,n)
			(ans[i]+=a[i][j]*b[j])%=mod;
		return ans;
	}
	void input(){
		repeat(i,0,n)repeat(j,0,n)a[i][j]=read();
	}
	void print(){
		cout<<"mat size="<<n<<endl;
		repeat(i,0,n){
			repeat(j,0,n)
				cout<<a[i][j]<<' ';
			cout<<endl;
		}
	}
	friend mat qpow(mat a,ll b){
		mat ans(1); // mat ans; repeat(i,0,n)ans[i][i]=1;
		while(b){
			if(b&1)ans=ans*a;
			a=a*a; b>>=1;
		}
		return ans;
	}
	ll *operator[](int x){return a[x];}
	const ll *operator[](int x)const{return a[x];}
};

高斯消元

struct mat{
	static const int N=410;
	vector<ll> a[N];
	mat(){for(auto &i:a)i.assign(N*2,0);} // if get_inv is needed, N*2
	ll det;
	void r_div(int x,int m,ll k){ // a[x][]/=k
		ll r=qpow(k,mod-2);
		repeat(i,0,m)
			a[x][i]=a[x][i]*r%mod;
		det=det*k%mod;
	}
	void r_plus(int x,int y,int m,ll k){ // a[x][]+=a[y][]*k
		repeat(i,0,m)
			a[x][i]=(a[x][i]+a[y][i]*k)%mod;
	}
	bool gauss(int n,int m){ // n<=m, return whether succuss
		det=1;
		repeat(i,0,n){
			int t=-1;
			repeat(j,i,n)if(a[j][i]){t=j; break;}
			if(t==-1){det=0; return 0;}
			if(t!=i){a[i].swap(a[t]); det=-det;}
			r_div(i,m,a[i][i]);
			repeat(j,0,n)
			if(j!=i && a[j][i])
				r_plus(j,i,m,mod-a[j][i]);
		}
		return 1;
	}
	ll get_det(int n){gauss(n,n); return (det+mod)%mod;} // return det
	bool get_inv(int n){ // self=inv(self), return whether success
		repeat(i,0,n)
		repeat(j,0,n)
			a[i][j+n]=(i==j);
		bool t=gauss(n,n*2);
		repeat(i,0,n)
		repeat(j,0,n)
			a[i][j]=a[i][j+n];
		return t;
	}
	vector<ll> &operator[](int x){return a[x];}
	const vector<ll> &operator[](int x)const{return a[x];}
}a;
  • 浮点版
struct mat{
	static const int N=110;
	vector<lf> a[N];
	mat(){for(auto &i:a)i.assign(N*2,0);} // if get_inv is needed, N*2
	lf det;
	void r_div(int x,int m,lf k){ // a[x][]/=k
		lf r=1/k;
		repeat(i,0,m)a[x][i]*=r;
		det*=k;
	}
	void r_plus(int x,int y,int m,lf k){ // a[x][]+=a[y][]*k
		repeat(i,0,m)a[x][i]+=a[y][i]*k;
	}
	bool gauss(int n,int m){ // n<=m, return whether succuss
		det=1;
		repeat(i,0,n){
			int t=-1;
			repeat(j,i,n)if(abs(a[j][i])>eps){t=j; break;}
			if(t==-1){det=0; return 0;}
			if(t!=i){a[i].swap(a[t]); det=-det;}
			r_div(i,m,a[i][i]);
			repeat(j,0,n)
			if(j!=i && abs(a[j][i])>eps)
				r_plus(j,i,m,-a[j][i]);
		}
		return 1;
	}
	lf get_det(int n){gauss(n,n); return det;} // return det
	bool get_inv(int n){ // self=inv(self), return whether success
		repeat(i,0,n)
		repeat(j,0,n)
			a[i][j+n]=(i==j);
		bool t=gauss(n,n*2);
		repeat(i,0,n)
		repeat(j,0,n)
			a[i][j]=a[i][j+n];
		return t;
	}
	vector<lf> &operator[](int x){return a[x];}
	const vector<lf> &operator[](int x)const{return a[x];}
}a;

int n;
struct mat{
	static const int N=110;
	vector< vector<ll> > a;
	mat():a(N,vector<ll>(N)){}
	ll det(int n){
		ll ans=1;
		repeat(i,0,n){
			repeat(j,i+1,n)
			while(a[j][i]){
				ll t=a[i][i]/a[j][i];
				repeat(k,i,n)a[i][k]=(a[i][k]-a[j][k]*t)%mod;
				swap(a[i],a[j]);
				ans=-ans;
			}
			ans=ans*a[i][i]%mod;
			if(!ans)return 0;
		}
		return (ans+mod)%mod;
	}
}a;

异或方程组

  • 编号从 0 开始，高斯消元部分 \(O(n^3)\)（luogu P2962）
bitset<N> a[N]; bool l[N];
int n,ans;
int gauss(int n){ // -1 : no solution, 0 : multi, 1 : single
	repeat(i,0,n){
		int t=-1;
		repeat(j,i,n)if(a[j][i]){t=j; break;}
		if(t==-1)continue;
		if(t!=i)swap(a[i],a[t]);
		repeat(j,0,n)
		if(i!=j && a[j][i])
			a[j]^=a[i];
	}
	repeat(i,0,n)if(!a[i][i] && a[i][n])return -1;
	repeat(i,0,n)if(!a[i][i])return 0;
	return 1;
}
void dfs(int x=n-1,int num=0){
	if(num>ans)return;
	if(x==-1){ans=num; return;}
	if(a[x][x]){
		bool v=a[x][n];
		repeat(i,x+1,n)
		if(a[x][i])
			v^=l[i];
		dfs(x-1,num+v);
	}
	else{
		dfs(x-1,num);
		l[x]=1;
		dfs(x-1,num+1);
		l[x]=0;
	}
}
int solve(){ // 返回满足方程组的sum(xi)最小值
	ans=inf; gauss(n); dfs(n-1,0);
	return ans;
}

线性基

  • 线性基是一系列线性无关的基向量组成的集合

异或线性基

  • 结论：\(basis.exist(a^b)\) 等价于 a, b 在 \(basis\) 里消去关键位后相等（要求是最简线性基，即第一个板子）
  • 插入、查询 \(O(\log M)\)
struct basis{
	static const int n=63;
	#define B(x,i) ((x>>i)&1)
	ll a[n],sz;
	bool failpush; // 是否线性相关
	void init(){mst(a,0); sz=failpush=0;}
	void push(ll x){ // 插入元素
		repeat(i,0,n)if(B(x,i))x^=a[i];
		if(x!=0){
			int p=__lg(x); sz++;
			repeat(i,p+1,n)if(B(a[i],p))a[i]^=x;
			a[p]=x;
		}
		else failpush=1;
	}
	ll top(){ // 最大值
		ll ans=0;
		repeat(i,0,n)ans^=a[i];
		return ans;
	}
	bool exist(ll x){ // 是否存在
		repeat_back(i,0,n)
		if((x>>i)&1){
			if(a[i]==0)return 0;
			else x^=a[i];
		}
		return 1;
	}
	ll kth(ll k){ // 第k小，不存在返回-1
		if(failpush)k--; // 如果认为0是可能的答案就加这句话
		if(k>=(1ll<<sz))return -1;
		ll ans=0;
		repeat(i,0,n)
		if(a[i]!=0){
			if(k&1)ans^=a[i];
			k>>=1;
		}
		return ans;
	}
}b;
basis operator+(basis a,const basis &b){ // 将b并入a
	repeat(i,0,a.n)
	if(b.a[i])a.push(b.a[i]);
	a.failpush|=b.failpush;
	return a;
}
  • 这个版本中求 kth 需要 rebuild \(O(\log^2 n)\)
struct basis{
	// ...
	void push(ll x){ // 插入元素
		repeat_back(i,0,n)
		if((x>>i)&1){
			if(a[i]==0){a[i]=x; sz++; return;}
			else x^=a[i];
		}
		failpush=1;
	}
	ll top(){ // 最大值
		ll ans=0;
		repeat_back(i,0,n)
			ans=max(ans,ans^a[i]);
		return ans;
	}
	void rebuild(){ // 求第k小的前置操作
		repeat_back(i,0,n)
		repeat_back(j,0,i)
		if((a[i]>>j)&1)
			a[i]^=a[j];
	}
}b;

实数线性基

  • 编号从 0 开始，插入、查询 \(O(n^2)\)
struct basis{
	lf a[N][N]; bool f[N]; int n; // f[i]表示向量a[i]是否被占
	void init(int _n){
		n=_n;
		fill(f,f+n,0);
	}
	bool push(lf x[]){ // 返回0表示可以被线性表示，不需要插入
		repeat(i,0,n)
		if(abs(x[i])>1e-5){ // 这个值要大一些
			if(f[i]){
				lf t=x[i]/a[i][i];
				repeat(j,0,n)x[j]-=t*a[i][j];
			}
			else{
				f[i]=1;
				repeat(j,0,n)a[i][j]=x[j];
				return 1;
			}
		}
		return 0;
	}
}b;

线性基求交

  • 未测试，\(O(\log^2 W)\)。
typedef array<int,64> Base;
Base merge(Base a,Base b){
	Base tmp(a),ans{};
	int cur,d;
	repeat(i,0,64)if(b[i]){
		cur=0,d=b[i];
		repeat_back(j,0,i+1)if(d>>j&1){
			if(tmp[j]){
				d^=tmp[j],cur^=a[j];
				if(d)continue;
				ans[i]=cur;
			}
			else tmp[j]=d,a[j]=cur;
			break;
		}
	}
	return ans;
}

线性规划 using 单纯形法

  • 声明：还没学会
  • \(\left[\begin{array}{ccccccc} a & a & a & a & a & a & b \\ a & a & a & a & a & a & b \\ a & a & a & a & a & a & b \\ c & c & c & c & c & c & v \end{array}\right]\)
  • 每行表示一个约束，\(\sum ax\le b\)，并且所有 \(x\ge 0\)，求 \(\sum cx\) 的最大值
  • 对偶问题：每列表示一个约束，\(\sum ax\ge c\)，并且所有 \(x\ge 0\)，求 \(\sum bx\) 的最小值
  • 先找 \(c[y]>0\) 的 y，再找 \(b[x]>0\) 且 \(\dfrac {b[x]}{a[x][y]}\) 最小的x（找不到 y 则 v，找不到 x 则 INF），用行变换将 \(a[x][y]\) 置 1，将其他 \(a[i][y]\) 和 \(c[y]\) 置 0
  • 编号从 1 开始，\(O(n^3)\)，缺init
const int M=1010; const lf eps=1e-6;
int n,m;
lf a[N][M],b[N],c[M],v; // a[1..n][1..m],b[1..n],c[1..m]
void pivot(int x,int y){
	b[x]/=a[x][y];
	repeat(j,1,m+1)if(j!=y)
		a[x][j]/=a[x][y];
	a[x][y]=1/a[x][y];
	repeat(i,1,n+1)
	if(i!=x && abs(a[i][y])>eps){
		b[i]-=a[i][y]*b[x];
		repeat(j,1,m+1)if(j!=y)
			a[i][j]-=a[i][y]*a[x][j];
		a[i][y]=-a[i][y]*a[x][y];
	}
	v+=c[y]*b[x];
	repeat(j,1,m+1)if(j!=y)
		c[j]-=c[y]*a[x][j];
	c[y]=-c[y]*a[x][y];
}
lf simplex(){ // 返回INF表示无限制，否则返回答案
	while(1){
		int x,y;
		for(y=1;y<=m;y++)if(c[y]>eps)break;
		if(y==m+1)return v;
		lf mn=INF;
		repeat(i,1,n+1)
		if(a[i][y]>eps && mn>b[i]/a[i][y])
			mn=b[i]/a[i][y],x=i;
		if(mn==INF)return INF; // unbounded
		pivot(x,y);
	}
}
void init(){v=0;}

数学杂项

struct of 自动取模

  • 未测，不好用，别用了。
struct mint{
	ll v;
	mint(ll _v){v=_v%mod;}
	mint operator+(const mint &b)const{return v+b.v;}
	mint operator-(const mint &b)const{return v-b.v;}
	mint operator*(const mint &b)const{return v*b.v;}
	explicit operator ll(){return (v+mod)%mod;}
};

struct of 区间

  • 未测，不好用，别用了。
bool isinter(pii a, pii b) {
	int l = max(a.fi, b.fi), r = min(a.se, b.se);
	return l <= r;
}
pii inter(pii a, pii b) {
	int l = max(a.fi, b.fi), r = min(a.se, b.se);
	if (l > r) return {0, -1};
	return {l, r};
}
pii unioned(pii a, pii b) { // 首先得 isinter
	return {min(a.fi, b.fi), max(a.se, b.se)};
}
void unioned(vector<pii> &a) {
	sort(a.begin(), a.end());
	int pre = 0;
	for (int i = 1; i < (int)a.size(); i++) {
		if (isinter(a[pre], a[i]))
			a[pre] = unioned(a[pre], a[i]);
		else
			a[++pre] = a[i];
	}
	a.erase(a.begin() + pre + 1, a.end());
}

struct of 高精度

  • 乘除 \(O(n^2)\)，且除法常数巨大。
  • 一般建议 Java。
struct big{
	vector<ll> a;
	static const ll k=1000000000,w=9;
	int size()const{return a.size();} // memory size
	explicit big(const ll &x=0){ // from ll
		*this=big(to_string(x));
	}
	explicit big(const string &s){ // from string
		static ll p10[9]={1};
		repeat(i,1,w)p10[i]=p10[i-1]*10;
		int len=s.size();
		int f=(s[0]=='-')?-1:1;
		a.resize(len/w+1);
		repeat(i,0,len-(f==-1))
			a[i/w]+=f*(s[len-1-i]-48)*p10[i%w];
		adjust();
	}
	int sgn(){return a.back()>=0?1:-1;} // sign (please used after adjust())
	void shrink(){ // pop zeros (will not release memory)
		while(size()>1 && a.back()==0)a.pop_back();
	}
	void adjust(){ // weak adjust, a[i] in (-k, k)
		repeat(i,0,3)a.push_back(0);
		repeat(i,0,size()-1){
			a[i+1]+=a[i]/k;
			a[i]%=k;
		}
		shrink();
	}
	void final_adjust(){ // strong adjust, a[i] have same sign
		adjust();
		int f=sgn();
		repeat(i,0,size()-1){
			ll t=(a[i]+k*f)%k;
			a[i+1]+=(a[i]-t)/k;
			a[i]=t;
		}
		shrink();
	}
	explicit operator string(){ // to string
		static char s[N]; char *p=s;
		final_adjust();
		if(sgn()==-1)*p++='-';
		repeat_back(i,0,size()){
			sprintf(p,i==size()-1?"%lld":"%09lld",abs(a[i]));
			p+=strlen(p);
		}
		return s;
	}
	const ll &operator[](int n)const{ // visit
		return a[n];
	}
	ll &operator[](int n){ // flexible visit
		repeat(i,0,n-size()+1)a.push_back(0);
		return a[n];
	}
	friend big operator+(big a,const big &b){ // <big + big>
		repeat(i,0,b.size())a[i]+=b[i];
		a.adjust();
		return a;
	}
	friend big operator-(big a,const big &b){ // <big - big>
		repeat(i,0,b.size())a[i]-=b[i];
		a.adjust();
		return a;
	}
	friend big operator*(const big &a,const big &b){ // <big * big>
		big ans;
		repeat(i,0,a.size()){
			repeat(j,0,b.size())
				ans[i+j]+=a[i]*b[j];
			ans.adjust();
		}
		return ans;
	}
	friend big operator*(big a,int b){ // <big * int>
		repeat(i,0,a.size())a[i]*=b;
		a.adjust();
		return a;
	}
	friend int to_abs(big &a){ // used in divide and mod
		a.final_adjust();
		int t=a.sgn(); a=a*t;
		return t;
	}
	friend big operator/(big a,int b){ // <big / int>
		int f=to_abs(a) * ((b>0)-(b<0));
		b=abs(b);
		ll rem=0;
		repeat_back(i,0,a.size()){
			ll s=rem*a.k+a[i];
			a[i]=s/b;
			rem=s%b;
		}
		return a*f;
	}
	friend big div(big &a,const big &b){ // used in <big / big>
		if(a<b)return big();
		big ans=div(a,b*2)*2;
		if(!(a<b)){
			a=a-b;
			ans=ans+big(1);
		}
		return ans;
	}
	friend big operator/(big a,big b){ // <big / big>
		int f=to_abs(a)*to_abs(b);
		return div(a,b)*f;
	}
	friend ll operator%(const big &a,ll mod){ // <big % ll>
		ll ans=0,p=1; mod=abs(mod);
		repeat(i,0,a.size()){
			(ans+=p*a[i])%=mod;
			(p*=a.k)%=mod;
		}
		return (ans+mod)%mod;
	}
	friend big operator%(big a,big b){ // <big % big>
		to_abs(b);
		int f=to_abs(a);
		return (a-a/b*b)*f;
	}
	friend bool operator<(big a,big b){ // <big less than big>
		a.final_adjust();
		b.final_adjust();
		repeat_back(i,0,max(a.size(),b.size()))
			if(a[i]!=b[i])return a[i]<b[i];
		return 0;
	}
	friend bool operator==(big a,big b){ // <big == big>
		a.final_adjust();
		b.final_adjust();
		repeat_back(i,0,max(a.size(),b.size()))
			if(a[i]!=b[i])return 0;
		return 1;
	}
};

struct of 分数

  • （可以直接哈希）（避免0/0，会当成0/1处理）
struct frac{
	ll x,y;
	explicit frac(ll x=0,ll y=1):x(x),y(y){init();}
	void init(){
		if(y<0)x=-x,y=-y;
		if(x==0)y=1;
		else if(y==0)x=1;
		else{
			ll d=abs(__gcd(x,y));
			x/=d; y/=d;
		}
	}
	frac operator-()const{
		return frac(-x,y);
	}
	friend frac operator+(const frac &a,const frac &b){
		return frac(a.x*b.y+a.y*b.x,a.y*b.y);
	}
	friend frac operator-(const frac &a,const frac &b){
		return a+-b;
	}
	friend frac operator*(const frac &a,const frac &b){
		return frac(a.x*b.x,a.y*b.y);
	}
	friend frac operator/(const frac &a,const frac &b){
		return frac(a.x*b.y,a.y*b.x);
	}
	friend ostream &operator<<(ostream &cout,const frac &f){
		return cout<<f.x<<'/'<<f.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;
	}
};

表达式求值

inline int lvl(const string &c){ // 运算优先级，小括号要排最后
	if(c=="*")return 2;
	if(c=="(" || c==")")return 0;
	return 1;
}
string convert(const string &in) { // 中缀转后缀
	stringstream ss;
	stack<string> op;
	string ans,s;
	repeat(i,0,in.size()-1){
		ss<<in[i];
		if(!isdigit(in[i]) || !isdigit(in[i+1])) // 插入空格
			ss<<" ";
	}
	ss<<in.back();
	while(ss>>s){
		if(isdigit(s[0]))ans+=s+" ";
		else if(s=="(")op.push(s);
		else if(s==")"){
			while(!op.empty() && op.top()!="(")
				ans+=op.top()+" ",op.pop();
			op.pop();
		}
		else{
			while(!op.empty() && lvl(op.top())>=lvl(s))
				ans+=op.top()+" ",op.pop();
			op.push(s);
		}
	}
	while(!op.empty())ans+=op.top()+" ",op.pop();
	return ans;
}
ll calc(const string &in){ // 后缀求值
	stack<ll> num;
	stringstream ss;
	ss<<in;
	string s;
	while(ss>>s){
		char c=s[0];
		if(isdigit(c))
			num.push((stoll(s))%mod);
		else{
			ll b=num.top(); num.pop();
			ll a=num.top(); num.pop();
			if(c=='+')num.push((a+b)%mod);
			if(c=='-')num.push((a-b)%mod);
			if(c=='*')num.push((a*b)%mod);
			// if(c=='^')num.push(qpow(a,b));
		}
	}
	return num.top();
}

数值积分 using 自适应辛普森法

  • 求 \(\displaystyle \int_{l}^{r}f(x)\mathrm{d}x\) 的近似值
lf raw(lf l,lf r){ // 辛普森公式
	return (f(l)+f(r)+4*f((l+r)/2))*(r-l)/6;
}
lf asr(lf l,lf r,lf eps,lf ans){
	lf m=(l+r)/2;
	lf x=raw(l,m),y=raw(m,r);
	if(abs(x+y-ans)<=15*eps)
		return x+y-(x+y-ans)/15;
	return asr(l,m,eps/2,x)+asr(m,r,eps/2,y);
}
// 调用方法：asr(l,r,eps,raw(l,r))
posted @ 2021-07-07 01:14  axiomofchoice  阅读(245)  评论(0)  编辑  收藏  举报