模板 - 数学

快速幂

递归

long long quick_pow(long long x,long long y,long long p)
{
	long long tmp=quick_pow(x,y>>1,p);
	tmp=tmp*tmp%p;
	if(y&1) return tmp*x%p;
	else return tmp;
}

迭代

long long quick_pow(long long x,long long y,long long p)
{
	long long res=1;
	while(y)
	{
		if(y&1) res=res*x%p;
		x=x*x%p;
		y>>=1;
	}
	return res;
}

矩阵快速幂

配合【矩阵】部分模板使用。

Matrix quick_pow(Matrix x,long long y)
{
	Matrix res=x; y--;
	while(y)
	{
		if(y&1) res=res*x;
		x=x*x;
		y>>=1;
	}
	return res;
}

欧几里得算法

gcd

int gcd(int x,int y){return y?gcd(y,x%y):x;}

lcm

int lcm(int x,int y){return x*y/gcd(x,y);}

扩展欧几里得算法 exgcd

LL exgcd(LL x,LL y,LL &u,LL &v)
{
	if(y)
	{
		LL u_,v_;
		LL res=exgcd(y,x%y,u_,v_);
		u=v_,v=u_-v_*(x/y);
		return res;
	}
	else
	{
		u=1,v=0;
		return x;
	}
}

乘法逆元

快速幂（费马小定理）

inline long long inv(long long x){return quick_pow(x,P-2);}

数列所有数的逆元

\[a_i^{-1} = \prod_{j=1}^{i-1}a_j \times \left( \prod_{j=1}^{i}a_j \right)^{-1} \]

fac_a[0]=1;
for(int i=1;i<=n;i++)
	fac_a[i]=fac_a[i-1]*a[i]%p;
LL tmp=inv(fac_a[n]);
for(int i=n;i>=1;i--)
{
	inv_a[i]=tmp*fac_a[i-1]%p;
	tmp=tmp*a[i]%p;
}

排列组合

乘法逆元求组合数

long long fac[N],inv_fac[N];
void Init_fac()
{
	fac[0]=1; inv_fac[0]=inv(fac[0]);
	for(int i=1;i<=n;i++)
	{
		fac[i]=fac[i-1]*i%P;
		inv_fac[i]=inv(fac[i]);
	}
}
long long C(int n,int m)
{
	return fac[n]*inv_fac[n-m]%P*inv_fac[m]%P;
}

高斯消元

Gauss 消元法

【占位】暂无

Gauss-Jordan 消元法

const double EPS=1e-6;
bool Gauss_Jordan()
{
	for(int i=1;i<=n;i++)
	{
		int p=i;
		while(abs(a[p][i])<EPS && p<=n) p++;
		if(p>n) return false;
		for(int j=1;j<=n+1;j++)
			swap(a[i][j],a[p][j]);
		double si=a[i][i];
		for(int j=i;j<=n+1;j++)
			a[i][j]/=si;
		for(int j=1;j<=n;j++)
		{
			if(j==i) continue;
			double sj=a[j][i];
			for(int k=i;k<=n+1;k++)
				a[j][k]-=sj*a[i][k];
		}
	}
	return true;
}

const double EPS = 1e-6;
int Gauss_Jordan()
{
	int row = 1; // 当前还未处理的最早的行数
	for (int col = 1; col <= n; col++) // 枚举列数
	{
		/* 第 1 步 */
		int p = row; // 在没处理过的行里找 p
		while (p <= n && abs(a[p][col]) < EPS) p++;
		if (p > n) continue; // 找不到合法 p，跳过这一行
		/* 第 2 步 */
		for (int i = 1; i <= n + 1; i++)
			swap(a[p][i], a[row][i]); // 交换两行
		/* 第 3 步 */
		double tmp = a[row][col]; // 预存，因为后面会修改其值
		for (int i = 1; i <= n + 1; i++)
			a[row][i] /= tmp; // 把 a[row][col] 调成 1
		/* 第 4 步 */
		for (int i = 1; i <= n; i++)
		{
			if (i == row) continue; // 不能把自己消了
			tmp = a[i][col]; // （同上）预存，因为后面会修改其值
			for (int j = 1; j <= n + 1; j++)
				a[i][j] -= a[row][j] * tmp; // 目的是消掉 a[i][row]
		}
		row++;
	}
	if (row == n + 1) return 1; // 有唯一合法解
	for (int i = row; i <= n; i++)
		if (abs(a[i][n + 1]) >= EPS) return -1; // 无合法解
	return 0; // 有无穷多合法解
}

矩阵

struct Matrix{
	int h,w;
	long long a[3][3];
	Matrix(int h_=0,int w_=0)
	:h(h_),w(w_){memset(a,0,sizeof(a));}
};

Matrix operator * (const Matrix &x,const Matrix &y)
{
	Matrix res(x.h,y.w);
	for(int i=1;i<=res.h;i++)
		for(int j=1;j<=res.w;j++)
		{
			for(int k=1;k<=x.w;k++)
				res.a[i][j]+=x.a[i][k]*y.a[k][j]%m;
			res.a[i][j]%=m;
		}
	return res;
}

扩展中国剩余定理 exCRT

I128 A=a[1],B=b[1];
bool uns=false;
for(int i=2;i<=n;i++)
{
	I128 p,q,d=exgcd(B,-b[i],p,q);
	if((a[i]-A)%d){uns=true;break;}
	I128 k=(a[i]-A)/d;
	B=lcm(B,b[i]);
	A=(a[i]+q*k*b[i])%B;
}
if(uns) printf("-1\n");
else write((A+B)%B,'\n');