清北学堂五一数学营笔记

矩阵

https://www.luogu.com.cn/blog/531930/shou-xie-yi-ge-ju-zhen-lei。

P4159 [SCOI2009] 迷路

如果边权是 $1$ 的话，答案就是邻接矩阵的 $k$ 次方的第 $1$ 行第 $n$ 列。

这个题有边权且很小（$0\le w\le 10$），考虑拆成 $1$ 边权。拆边运算量太大，考虑拆点。

我们把一个点拆成 $9$ 个。记 $id(i,j)$ 表示点 $i$ 拆出来的第 $j+1$ 个点，其中 $0\le j\le 8$，满足 $id(i,j)$ 到 $i$ 的距离为 $j$。$id(i,0)$ 就是点 $i$ 本身。

我们把 $id$ 看做一个 $n\times 9$ 的二维数组，这样我们就可以给每一个位置一个编号：$id(i,j)=9(i-1)+j$。容易发现点的范围是 $[0,9n-1]$。

建边我们“向内”建，由 $id(i,j)$ 向 $id(i,j-1)$ ；连权值为 $1$ 的边，这样从 $id(i,j)$ 往前走 $j$ 条边就可以到达 $id(i,0)$。

然后是原图的建边。如果 $w(i,j)\not=0$，那么由 $id(i,0)$ 向 $id(j,w(i,j)-1)$ 连边。也就是先从 $id(i,0)$ 经过一条边到达 $id(j,w(i,j)-1)$，然后再经过 $w(i,j)-1$ 条边到达 $id(j,0)$，一共 $w(i,j)$ 条边，和题目相符。

最后我们就得到了一个大小为 $9n\times9n$ 的邻接矩阵，直接用矩阵快速幂计算该矩阵的 $t$ 次方即可。答案就是该矩阵的第 $id(1,0)$ 行第 $id(n,0)$ 列的数。

由于点的范围是 $[0,9n-1]$ 的，所以矩阵乘法的时候直接从 $0$ 枚举到 $9n-1$ 就可以了。最后不要忘记取模。

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int mod=2009;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
   	int x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
int n,tot,t;
inline int id(int x,int y){
	return (x-1)*9+y;
}
struct Matrix{
	int val[105][105];
	Matrix(){
		memset(val,0,sizeof(val));
	}
	inline void reset(){
		for(int i=0;i<tot;++i) val[i][i]=1;
	}
}tmp;
inline Matrix operator*(const Matrix &a,const Matrix &b){
	Matrix c;
	for(int i=0;i<tot;++i){
		for(int k=0;k<tot;++k){
			for(int j=0;j<tot;++j){
				c.val[i][j]=(c.val[i][j]+a.val[i][k]*b.val[k][j]%mod)%mod; 
			}
		}
	}
	return c;
}
inline Matrix ksm(Matrix a,int b){
	Matrix c;
	c.reset();
	while(b){
		if(b&1)c=c*a;
		a=a*a,b>>=1;
	}
	return c;
}
signed main(){
	n=read(),t=read(),tot=9*n;
	for(int i=1;i<=n;++i){
		for(int j=1;j<=8;++j) tmp.val[id(i,j)][id(i,j-1)]=1;
		for(int j=1;j<=n;++j){
			int x;
			scanf("%1d",&x);
			if(x) tmp.val[id(i,0)][id(j,x-1)]=1;
		}
	}
	tmp=ksm(tmp,t);
	println(tmp.val[id(1,0)][id(n,0)]);
	return 0;
}

一个小语法

如果你在主函数中定义了一个全局变量 $k$，然后在某个循环中有定义了一次，如：

int k=1;
signed main(){
	for(int k=2;k<=3;++k) println(k);
	return 0;
}

这时候就会输出内部循环中的 $k$。但是只要改成 ::k 就可以访问到外部的 $k$。

数论

判断质数

一、试除法

先判断 $1,0$，然后直接枚举 $[2,x]$ 的数试除即可。$O(x)$

显而易见的优化是，质因子总是成对出现的，因此只需要枚举较小的那个质因子（不超过 $\sqrt{x}$）。试除 $[2,\sqrt{x}]$ 的数，$O(\sqrt{x})$。

同样的，我们只需要试除质数。可以先 $O(\sqrt{x})$ 筛出 $[2,\sqrt{x}]$ 的质数然后试除即可。根据素数密度就知道复杂度约为 $O(\dfrac{\sqrt{x}}{\ln x})$。

代码太简单，不写了。

二、Miller - Rabin

https://www.luogu.com.cn/paste/0ffg5egd

分解质因数

一、朴素做法

直接枚举 $[2,\sqrt{x}]$ 的所有数 $a$ 然后除就可以了。如果 $a$ 是合数，那么 $x$ 已经被 $a$ 的某个更小的质因子除过了，所以 $x$ 必然不是 $a$ 的倍数。因此这样找到的 $a$ 必然是质数。最后，由于 $x$ 最多只有一个 $>\sqrt{x}$ 的质因子，所以特判一下即可。时间复杂度 $O(\sqrt{x})$ 同样可以通过先筛质数来优化复杂度。

二、Pollard's rho 算法

还没过。

UPD：没讲

CF45G - Prime Problem

纯纯的找规律 & 打表题。

首先我们需要知道哥德巴赫猜想，并且知道它在本题的数据范围内已经得到了验证（也可以自己手动验证）。

首先计算出数列的和 $sum=\dfrac{n(n+1)}{2}$。初始都归在同一组。

如果 $sum$ 为质数（试除法），那么直接输出，结束；
如果 $sum$ 为偶数，根据哥猜，可以拆成两个质数之和。通过打表可以发现必然有一个不超过 $n$。考虑每一个 $i$，验证即可；
如果 $sum$ 为奇数，由于 $n\ge2$，所以必然能拆出一个 $3$ 来，然后 $sum$ 就变成偶数，此时再考虑上一种情况即可。注意不要考虑上 $3$。

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=6005;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
   	int x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
int n,sum,ans[N];
bool vis[N];
inline bool is_prime(int x){
	if(x<2)return 0;
	for(int i=2;i*i<=x;++i) if(x%i==0)return 0;
	return 1;
}
inline void output(){
	for(int i=1;i<=n;++i) printsp(ans[i]);
	printf("\n");
}
signed main(){
	n=read(),sum=n*(n+1)/2;
	for(int i=1;i<=n;++i) ans[i]=1;
	if(is_prime(sum)){
		output();
		return 0;
	}
	if(sum%2==1&&!is_prime(sum-2))sum-=3,ans[3]=3,vis[3]=1;
	for(int i=2;i<=n;++i){
		if(vis[i])continue;
		if(is_prime(i)&&is_prime(sum-i)){
			ans[i]=2;
			break;
		}
	}
	for(int i=1;i<=n;++i) printsp(ans[i]);
	return 0;
}

逆元

若 $b\times c \equiv 1 \pmod p$，那么就称 $c$ 为 $b$ 在模 $p$ 意义下的乘法逆元，类似于倒数。借助逆元，我们就可以完成模意义下的除法操作了。

那么怎么求逆元呢？

费马小定理：$a^{p-1}\equiv1\pmod p$，满足 $p$ 为质数，$\gcd(a,p)=1$。

我们发现 $a\times a^{p-2}=a^{p-1}$，所以：

\[a\times a^{p-2}\equiv1\pmod p \]

然后这个 $a^{p-2}$ 不就是 $a$ 的逆元了吗！

快速幂计算即可，时间复杂度 $O(\log p)$。

所以，在 $p$ 为质数，$a$ 于 $p$ 互质的情况下，$a\div b \bmod c=a\times b^{p-2} \bmod c$。

欧拉函数：$\phi(p)$ 表示 $[1,p]$ 中有多少个数与 $p$ 互质。

欧拉定理：$a^{\phi(p)}\equiv1\pmod p$，满足 $\gcd(a,p)=1$

这个时候我们会发现如果一个数 $p$ 是素数，那么 $\phi(p)=p-1$（除 $p$ 外都与它互质），这样就是费马小定理了。所以费马小定理实际上是欧拉定理的特殊情况。

这样 $a$ 的逆元就是 $a^{\phi(p)-1}$。

注意：逆元存在的充分条件是 $\gcd(a,p)=1$，其他情况下不一定存在逆元。

最后的问题就是如何算出 $\phi(n)$ 了。

一、暴力

枚举 $[1,n-1]$ 的所有数并和 $n$ 取 $\gcd$，时间复杂度 $O(n\log n)$。

二、优化

暴力的复杂度太高了，我们考虑优化。

当 $n=$ 一个质数 $p$ 时，显然 $\phi(n)=p-1$。
当 $n=p^2$ 时，考虑用总数减去不互质的。显然只有 $p$ 的倍数不互质，$[1,n]$ 中一共有 $p$ 个数是 $p$ 的倍数，所以 $\phi(n)=n-p=p^2-p=p(p-1)$。
当 $n=p^k$ 时，考虑将 $[1,n]$ 分为若干长度为 $p$ 的段，一共有 $\dfrac{n}{p}$ 段。每段只有最后一个数（$p$ 的倍数）与 $n$ 不互质，其他 $p-1$ 个数都互质，所以 $\phi(p)=\dfrac{n}{p}(p-1)$。稍微整理一下（把 $p$ 移到 $p-1$ 下面）得 $\phi(n)=n(1-\dfrac{1}{p})$。
加一个质因子：当 $n=p_1^kp_2^k$ 时，还是按照上面的方法分段，容斥一下得 $\phi(n)=n-\dfrac{n}{p_1}-\dfrac{n}{p_2}+\dfrac{n}{p_1p_2}$。因式分解一下：先提 $n$，得 $\phi(n)=n(1-\dfrac{1}{p_1}-\dfrac{1}{p_2}+\dfrac{1}{p_1p_2})$；分组分解，得 $\phi(n)=n(1-\dfrac{1}{p_1})(1-\dfrac{1}{p_2})$。
然后就可以得出规律了：将 $n$ 质因数分解为 $\sum_{i=1}^kp_i^{a_i}$，那么答案就是 $n\prod_{i=1}^k(1-\dfrac{1}{p_i})$。

所以直接分解就可以了，一个小技巧：先除后乘防止溢出。在求 $lcm$ 中也有用。

inline int get_phi(int n){
	int phi=n;
	for(int i=2;i*i<=n;++i){
		if(n%i==0){
			phi=phi/i*(i-1);
			while(n%i==0) n/=i;
		}
	}
	if(n!=1) phi=phi/n*(n-1);
	return phi;
}

这样时间复杂度就是 $O(\sqrt{n})$ 了。你想先筛一下也行。

求 $[1,n]$ 逆元

$1.$ 暴力：枚举每个数，快速幂计算逆元即可。时间复杂度 $O(n\log n)$。

$2.$ 顺推：假设我们求出了 $[1,i-1]$ 所有数的逆元，现在考虑求出 $i$ 的。$p$ 明显可以写成 $ki+r$ 的形式。而且由于 $i<p$ 且 $p$ 为质数，所以 $0<r<i$。现在把等式两边同时除 $p$：

\[ki+r\equiv0\pmod p \]

移项，得：

\[ki\equiv-r\pmod p \]

除 $i$ ，得:

\[k\equiv\dfrac{-r}{i}\pmod p \]

最后再把 $-r$ 移过去：

\[\dfrac{1}{i}\equiv-\dfrac{k}{r}\pmod p \]

最后，把 $k$ 和 $r$ 重新表示一下：

\[\dfrac{1}{i}\equiv-\dfrac{\lfloor\dfrac{p}{i}\rfloor}{p \bmod i}\pmod p \]

写成递推式就是 $inv_i=(-\lfloor\dfrac{p}{i}\rfloor)\times inv_{p \bmod i}$。

最后，由于 $-\lfloor\dfrac{p}{i}\rfloor$ 是负的，所以加个 $p$ 即可：

$inv_i=(p-\lfloor\dfrac{p}{i}\rfloor)\times inv_{p \bmod i}$

$3.$ 倒推：首先把每个数的阶乘算出来，然后调用一次快速幂算出 $n!$ 的逆元。由于 $\dfrac{1}{i!}=\dfrac{i+1}{(i+1)!}$，因此我们就可以倒推求出所有 $\dfrac{1}{i!}$ 了。最后 $\dfrac{(i-1)!}{i!}=\dfrac{1}{i}$，于是就做完了。时间复杂度 $O(n)$。

这种做法还可以拓展到任意 $n$ 个数的情况。把阶乘换成前缀积即可。

exgcd

给定 $a,b,c$，求任意一组整数解 $x,y$，满足 $ax+by=\gcd(a,b)$

假设我们现在要求 $ax+by=\gcd(a,b)$。令 $d=\gcd(a,b)$，由于 $\gcd(a,b)=\gcd(b,a \bmod b)$，所以我们考虑把 $b$ 和 $a\bmod b$ 也写成类似的形式：

\[bx'+(a\bmod b)y'=d \]

现在我们要从这个式子推出 $x,y$。把 $a\bmod b$ 拆开，得：

\[bx'+(a-\lfloor\dfrac{a}{b}\rfloor\times b)y'=d \]

拆项，得：

\[bx'+ay'-\lfloor\dfrac{a}{b}\rfloor\times b\times y'=d \]

把带 $a$ 的凑到一起，带 $b$ 的凑到一起，得：

\[ay'+b(x'-\lfloor\dfrac{a}{b}\rfloor\times y')=d \]

这个时候，我们发现这个式子和原来的式子形式一样了！因此：

\[x=y',y=x'-\lfloor\dfrac{a}{b}\rfloor\times y' \]

因此只要递归算出 $bx+(a\bmod b)\times y = d$ 的解即可。

那递归的终止条件呢？当 $b=0$ 时，$\gcd(a,b)=a$，因此令 $x=1,y=0$ 即可。

代码实现时可以用一个三个元素的结构体维护 $x$，$y$ 和 $\gcd$，这样就可以顺便求出 $\gcd(a,b)$ 了。

代码。

用途？

现在有一个问题：$x,y,a,b$ 为整数，那么 $ax+by$ 能表示出的最小正整数是多少？

答案是 $\gcd(a,b)$，现在考虑证明：

设 $\gcd(a,b)=d$，那么 $a$ 可以表示为 $a'd$，$b$ 可以表示为 $b'd$，代回原式：

\[a'dx+b'dy \]

提 $d$，得：

\[d(a'x+b'y) \]

所以得数一定是 $\gcd(a,b)$ 的倍数，所以 $\gcd(a,b)$ 就是最小正整数。这也说明了，只有得数是 $\gcd(a,b)$ 的倍数的时候，方程才有解。这就是裴蜀定理。

现在要求 $ax+by=c$ ，满足 $c$ 是 $\gcd(a,b)$ 的倍数。怎么做？

把得到的 $x,y$ 同乘 $\dfrac{c}{\gcd(a,b)}$ 即可。

现在只求出了一组解，求通解怎么办？

假设 $x_0,y_0$ 是一开始 exgcd 的出来的一组解，那么通解可以写成 $x=x_0+k\times\dfrac{b}{\gcd(a,b)},y=y_0-k\times\dfrac{a}{\gcd(a,b)}$ 的形式，其中 $k$ 为任意整数。

而且有结论：上述算法取得的一组解 $(x,y)$ 满足 $|x|\le|\frac{b}{2\gcd(a,b)}|,|y|\le|\frac{a}{2\gcd(a,b)}|$，证明可以考虑归纳。具体见《深进》P298-P299。

P1082 [NOIP2012 提高组] 同余方程

这个题既然保证有解，那么 $\gcd(a,b)=1$。

做法一：容易发现这个 $x$ 就是 $a$ 的逆元。由于 $b$ 不一定是质数，因此用欧拉定理即可。

做法二：我们改写一下它的形式。既然 $ax\bmod b=1$，那么 $ax$ 一定可以写成 $by+1$ 的形式。移个项，变为 $ax-by=1$，那么直接上 exgcd 就行了。

但是这个题要求最小正整数解。考虑通解公式，容易发现 $x\equiv x_0\pmod \dfrac{b}{\gcd(a,b)}$，由于 $\gcd=1$，所以 $x\equiv x_0\pmod b$，容易发现答案就是 $x_0\bmod b$，所以 (%mod+mod)%mod 防止负数即可。

CRT/EXCRT

用途是解如下同余方程组

\[\begin{cases} x \equiv b_1\ ({\rm mod}\ a_1) \\ x\equiv b_2\ ({\rm mod}\ a_2) \\ ... \\ x \equiv b_n\ ({\rm mod}\ a_n)\end{cases} \]

如果不要求 $a$ 两两互质就是 EXCRT，否则就是普通 CRT。

对于普通 CRT，有一个小构造：记 $M=\prod_{i=1}^{n}a_i$，$m_i=\dfrac{M}{a_i},t_i\equiv\dfrac{1}{m_i}\pmod {a_i}$，那么答案为 $\sum_{i=1}^{n}b_it_im_i$。

证明：对于第 $i$ 个方程组，对于 $j\not=i$，$m_j$ 都是 $a_i$ 的倍数，所以 $b_jt_jm_j=0$；而对于 $i$，这三项都不是 $a_i$ 的倍数，所以 $b_it_im_i=b_i\times1=b_i$。

对于 $a$ 不互质的情况，$m_i$ 与 $a_i$ 也不一定互质，也就不一定存在 $t_i$（逆元），这时就需要用 EXCRT 了。

EXCRT 在这里。

筛法

一、暴力

暴力判断每个数是否是质数即可。试除法可以做到 $O(n\sqrt{n})$，Miller - Rabin 可以做到 $O(n\log n)$。

二、埃氏筛

把每个数的倍数筛掉。

not_prime[1]=1;
for(int i=2;i<=n;++i){
	for(int j=i+i;j<=n;j+=i){
		not_prime[j]=1;
	}
}
for(int i=2;i<=n;++i) if(!not_prime[i]) prime[++cnt]=i;

复杂度是多少？$\sum_{i=2}^{n}\dfrac{n}{i}$。补一个 $\dfrac{n}{1}$，提出 $n$，得 $n\sum_{i=1}^{n}\dfrac{1}{i}$。有调和级数可得约为 $O(n\log n)$。

优化：只需要筛质数的倍数。

not_prime[1]=1;
for(int i=2;i<=n;++i){
	if(!not_prime[i]){
		for(int j=i+i;j<=n;j+=i){
			not_prime[j]=1;
		}
	}
}
for(int i=2;i<=n;++i) if(!not_prime[i]) prime[++cnt]=i;

不加证明地给出算法的复杂度：约为 $O(n\log\log n)$。（是 $\log$ 套 $\log$，不是 $\log$ 乘 $\log$）。

其实还可以优化：设 $j=ki$，当 $k<i$ 时，$j$ 已经被 $k$ 筛过了，因此从 $j=i^2$ 开始筛即可。

not_prime[1]=1;
for(int i=2;i<=n;++i){
	if(!not_prime[i]){
		for(int j=i*i;j<=n;j+=i){
			not_prime[j]=1;
		}
	}
}
for(int i=2;i<=n;++i) if(!not_prime[i]) prime[++cnt]=i;

三、线性筛

埃氏筛的缺点在于，每个数可能被它的质因子筛了多次。所以它的时间复杂度必然不可能达到 $O(n)$。

线性筛的思想就是，让每个数被它的最小质因子筛掉。

过程：维护 bool 数组，判断某个数是否是质数；再维护一张质数表，里面存放 $[2,n]$ 的所有质数（省略 $1$）。从小到大扫描 $i$，如果 $i$ 不是质数，存放到素数表里。然后，取出素数表里的数，把该质数 $\times i$ 筛掉（超过 $n$ 就 break）。

const int N=1e6+5;
int n,prime[N],cnt,T;
bool not_prime[N];
inline void do_prime(){
	not_prime[0]=not_prime[1]=1,cnt=0;
	for(int i=2;i<=N-5;++i){
		if(!not_prime[i])prime[++cnt]=i;
		for(int j=1;j<=cnt&&i*prime[j]<=N-5;++j){
			not_prime[i*prime[j]]=1;
		}
	}
}

可惜，这样还是 $O(n\log\log n)$ 的。

但是我们只需要在第 $9$ 行下面添加如下代码：

if(i%prime[j]==0)break;

这样就 $O(n)$ 了。

为什么？

当 $i\bmod prime_j=0$ 时，那么 $i$ 的最小质因子就是 $prime_j$（否则前面已经 break 掉了）。而如果 $j+1$，那么被筛掉的数是 $prime_{j+1}\times i$，由于 $prime_{j+1}>prime_j$，所以 $prime_j$ 也是 $prime_{j+1}\times i$ 的最小质因子。这样它就不是由它的最小质因子筛掉了，需要 break 掉。这样就可以保证每个数都是由它的最小质因子筛掉，时间复杂度 $O(n)$。

代码：

const int N=1e6+5;
int n,prime[N],cnt,T;
bool not_prime[N];
inline void do_prime(){
	not_prime[0]=not_prime[1]=1,cnt=0;
	for(int i=2;i<=N-5;++i){
		if(!not_prime[i])prime[++cnt]=i;
		for(int j=1;j<=cnt&&i*prime[j]<=N-5;++j){
			not_prime[i*prime[j]]=1;
			if(i%prime[j]==0)break;
		}
	}
}

这玩意还能用来筛积性函数。

积性函数

定义：对于一个函数 $f$，若 $f(a)f(b)=f(ab)$，且 $\gcd(a,b)=1$，那么称 $f$ 为积性函数。如果没有 $\gcd(a,b)=1$ 的限制，那么 $f$ 就是完全积性函数。

举个例子：$\phi$ （欧拉函数）是积性函数。

回忆一下公式：设

\[a=\sum_{i=1}^{n}p_i^{a_i},b=\sum_{i=1}^{m}q_i^{b_i} \]

那么

\[\phi(a)=a\prod_{i=1}^{n}(1-\dfrac{1}{p_i}),\phi(b)=b\prod_{i=1}^{m}(1-\dfrac{1}{q_i}) \]

假设 $a,b$ 互质，所以 $a,b$ 没有公共的质因子，所以

\[ab=\sum_{i=1}^{n}p_i^{a_i}\sum_{i=1}^{m}q_i^{b_i} \]

所以

\[\phi(ab)=ab\prod_{i=1}^{n}(1-\dfrac{1}{p_i})\prod_{i=1}^{m}(1-\dfrac{1}{q_i}) \]

而

\[\phi(a)\phi(b)=ab\prod_{i=1}^{n}(1-\dfrac{1}{p_i})\prod_{i=1}^{m}(1-\dfrac{1}{q_i}) \]

所以 $\phi(ab)=\phi(a)\phi(b)$，所以 $\phi$ 是积性函数。

现在我们要线性求出 $[1,n]$ 中所有数的 $\phi$ 值。

首先，$\phi(1)=1$。并且当 $p$ 为质数时，$\phi(p)=p-1$。这样我们就只需要考虑计算合数的欧拉函数值。

由于我们会筛掉 $i\times prime_j$，所以直接更新 $\phi(i\times prime_j)=\phi(i)\times\phi(prime_j)$ 就可以了……吗？

并不行，因为两者不一定互质。什么时候不互质？$i\bmod prime_j=0$ 时。所以我们重新算一下 $\phi$ 就好了。等于多少呢？等于 $prime_j\times\phi(i)$。

考虑证明，令 $i=\sum_{j=1}^{n}p_j^{a_j}$，那么

\[\phi(i)=i\prod_{j=1}^{n}(1-\dfrac{1}{p_j}) \]

由于 $prime_j$ 是 $i$ 的质因子，所以乘上它之后只会让它的指数 $+1$，而不会多出一个新的质因子。所以：

\[\phi(i\times prime_j)=prime_j\times i\prod_{k=1}^{n}(1-\dfrac{1}{p_k}) \]

容易发现后面的部分就是 $\phi(i)$，所以：

\[\phi(i\times prime_j)=prime_j\times\phi(i) \]

这样我们就可以筛出 $[1,n]$ 所有数的欧拉函数值了。

const int N=1e6+5;
int n,prime[N],cnt,T,phi[N];
bool not_prime[N];
inline void do_prime(){
	not_prime[0]=not_prime[1]=phi[1]=1,cnt=0;
	for(int i=2;i<=N-5;++i){
		if(!not_prime[i])prime[++cnt]=i,phi[i]=i-1;
		for(int j=1;j<=cnt&&i*prime[j]<=N-5;++j){
			not_prime[i*prime[j]]=1;
			phi[i*prime[j]]=phi[i]*phi[prime[j]];
			if(i%prime[j]==0){
				phi[i*prime[j]]=prime[j]*phi[i];
				break;
			}
		}
	}
}

BSGS 算法

~~Baidu Search Google Search 算法~~

原题链接。

一、暴力

直接枚举 $x$，并迭代计算 $a^x$ 并检验即可。

时间复杂度为 $O(\text{答案})$，明显不靠谱。

代码还是放一下吧：

inline int solve(int a,int b,int p){
	// 求解 a^x % p = b
	int val=1;
	for(int x=0;x;++x){
		if(val==b)return x;
		val=1ll*val*a%p;
		if(val==1)return -1;// 出现循环，无解 
	} 
}

为什么循环一定从 $1$ 开始呢？因为费马小定理（$p$ 是质数）。同理，我们还可以知道循环节长度为 $p-1$。因此，加上无解判断的代码时间复杂度实际上是 $O(p)$ 的。

所以我们改一下：

inline int solve(int a,int b,int p){
	// 求解 a^x % p = b
	int val=1;
	for(int x=0;x<p-1;++x){
		if(val==b)return x;
		val=1ll*val*a%p;
	}
	return -1;
}

这样我们就可以喜提 $0$ 分了。

二、正解

由于循环节长度只有 $p-1$，所以我们只需要考虑 $p-1$ 个数。

我们将这些数分组，假设每组有 $s$ 个数。

暴力扫描每一个组，复杂度还是一样的。

现在第一组还是暴力，但是我们需要考虑用更快的方法求出后面的组。

很明显，第二组的每个数都是第一组对应数的 $a^s$ 倍。换句话说，把第二组中的数全部除 $a^s$，就得到了第一组。

我们可以发现，如果第二组存在 $b$，那么第一组一定存在 $\dfrac{b}{a^s}$！

换句话说，如果第 $i$ 组存在 $b$，那么第一组一定存在 $\dfrac{b}{a^{s(i-1)}}$。

这样我们只需要在第一组里面查找了。怎么查找？哈希或者二分。

容易发现这是一个类似分块的操作。我们需要扫描第一组，用哈希维护，时间复杂度 $O(s)$；然后需要扫描每一组找答案，时间复杂度 $O(\dfrac{p}{s})$。根据均值不等式很容易知道 $s=\sqrt{p}$ 最优，这样时间复杂度就是 $O(\sqrt{p})$。

我代码里偷懒用了 set，还写了快速幂，时间复杂度是 $O(\sqrt{p}\log p)$ 的，很逊。

代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline ll read(){
   	ll x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
ll a,b,p;
inline ll ksm(ll a,ll b,ll p){
	ll res=1;
	a%=p;
	while(b){
		if(b&1)res=res*a%p;
		a=a*a%p,b>>=1;
	}
	return res;
}
inline ll solve(ll a,ll b,ll p){
	// 求解 a^x % p = b
	ll s=sqrt(p),val=1;
	set<ll>S;
	for(int i=0;i<s;++i){
		S.insert(val);
		val=val*a%p;
	}
	for(int i=0;1ll*i*s<=p;++i){
		// 查找第 i 组（从 0 开始，这样就查找 b/a^si 即可） 
		ll inv=ksm(a,i*s,p);
		ll tmp=b*ksm(inv,p-2,p)%p;
		if(S.count(tmp)){
			// inv 就是第 i 行的第一个数 
			for(int j=i*s;;++j){
				// 暴力查找第 i 组中满足条件的数 
				if(inv==b)return j;
				inv=inv*a%p;
			}
		}
	}
	return -1;
}
signed main(){
	p=read(),a=read(),b=read();
	int ans=solve(a,b,p);
	if(ans==-1)printf("no solution\n");
	else println(ans);
	return 0;
}

容易发现，每次的循环的 $tmp$ 只是在上一次的基础上乘上 $\dfrac{1}{a^s}$。所以快速幂预处理出 $\dfrac{1}{a^s}$ 就可以去掉那只 $\log$。

当 $i=0$ 的时候情况稍微有点复杂，特判一下即可。

提供一个哈希表板子。

由于 Hash 表的时间复杂度是期望线性的，所以时间复杂度 $O(\sqrt{p}+\log p)=O(\sqrt{p})$。

代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=1e5+5;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline ll read(){
    ll x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
ll a,b,p;
inline ll ksm(ll a,ll b,ll p){
    ll res=1;
    a%=p;
    while(b){
        if(b&1)res=res*a%p;
        a=a*a%p,b>>=1;
    }
    return res;
}
template<class P,class Q>
struct Hash_Table{
    static const int mod=5e5+9;
    struct Node{
        P key;Q val;
        int nxt;
    }node[N];
    int head[mod+5],cnt;
    inline int get_hash(P x){return (x%mod+mod)%mod;}
    inline void clear(){
        for(int i=1;i<=cnt;++i){
            head[get_hash(node[i].key)]=0,node[i].key=node[i].nxt=0;
            Q tmp;
            swap(node[i].val,tmp);
        }
        cnt=0;
    }
    inline void insert(P key){
        int val=get_hash(key);
        node[++cnt].key=key;
        node[cnt].nxt=head[val];
        head[val]=cnt;
    }
    inline Q &operator[](const P key){
        int u=get_hash(key);
        for(int i=head[u];i;i=node[i].nxt){
            P v=node[i].key;
            if(v==key)return node[i].val;
        }
        insert(key);
        return node[cnt].val;
    }
    inline bool count(const P key){
        int u=get_hash(key);
        for(int i=head[u];i;i=node[i].nxt){
            P v=node[i].key;
            if(v==key)return 1;
        }
        return 0;
    } 
};
Hash_Table<ll,int>mp;
inline ll BSGS(ll a,ll b,ll p){
    // 求解 a^x % p = b
    ll s=sqrt(p),val=1;
    for(int i=0;i<s;++i) mp[val]=i,val=val*a%p;
    if(mp.count(b))return mp[b];
    ll inv=ksm(ksm(a,s,p),p-2,p),tmp=b%p;
    for(int i=1;1ll*i*s<=p;++i){
        // 查找第 i 组（从 0 开始，这样就查找 b/a^si 即可） 
        tmp=(tmp*inv)%p;
        if(mp.count(tmp))return mp[tmp]+i*s;
    }
    return -1;
}
signed main(){
    p=read(),a=read(),b=read();
    int ans=BSGS(a,b,p);
    if(ans==-1)printf("no solution\n");
    else println(ans);
    return 0;
}

也可以换成 unorderep_map 或者 gp/cc_hash_table，但是性能会差一点。

组合数学

基础知识

加法原理：两种选择是有关系/矛盾的，不能同时选择（同一阶段）。

乘法原理：两种选择是没有关系的，可以同时选择（不同阶段）。

排列：从 $n$ 个人里选出 $m$ 个人站成一列（考虑顺序, $12$ 和 $21$ 是两种情况），有几种选法？

对于第一个人，有 $n$ 种选法（都可以选），第二个人可以选 $n-1$ 种（第一个人已经选了）……一直到第 $m$ 个人，有 $n-m+1$ 种选法。

由于每个人的选择是独立的，所以乘法原理：$P(n,m)=n\times n-1\times n-2\times ... \times (n-m+1)=\dfrac{n!}{m!}$（排列数）

ZHX 语录：小学奥数里用 $A$ 是因为怕我们不认识 $P$。

组合：从 $n$ 个人里选出 $m$ 个人，不考虑顺序（$12 $ 和 $21$ 是一种情况）。要排除重复的方案数，有 $\binom{n}{m}=\dfrac{n!}{m!(n-m)!}$

组合

定义式：$\binom{n}{m}=\dfrac{n!}{m!(n-m)!}$

现在来考虑一些性质：

$\binom{n}{0}=\binom{n}{n}=1$（只有一种方案，都不选和都选）。
$\binom{n}{m}=\binom{n}{n-m}$（选 $m$ 个相当于不选 $n-m$ 个，也可以直接代进公式求出来）。
递推式：$\binom{n}{m}=\binom{n-1}{m}+\binom{n-1}{m}$（DP 思想，考虑第 $n$ 个选不选，同一阶段，加法原理）。
$\sum_{i=0}^{n}\binom{n}{i}=2^n$（选 $[0,n]$ 任意多个物品的方案数，每个物品可选可不选）。
当 $n\ge1$ 时，$\sum_{i=0}^{n}(-1)^i\binom{n}{i}=0$。

相当于证明选偶数个物品的方案数等于选奇数个物品的方案数的方案数。

考虑递推式：$\binom{n}{m}=\binom{n-1}{m}+\binom{n-1}{m}$。那么，选奇数个物品的方案数实际上就是
$ \sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\binom{n-1}{2i+1}+\binom{n-1}{2i+1-1}=\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\binom{n-1}{2i+1}+\binom{n-1}{2i} $，和选偶数个物品的方案数一样。

二项式定理：$(x+y)^n=\sum_{i=0}^{n}\binom{n}{i}x^iy^{n-i}$。证明考虑乘法分配率的过程。
把 $\binom{n}{m}$ 按照递推式进行 $k$ 次展开，得：$\sum_{i=0}^{k}\binom{k}{i}\times\binom{n-k}{m-k+i}=\sum_{i=0}^{k}\binom{k}{i}\times\binom{n-k}{m-i}$，这个就是组合数的展开式，好像也是一种卷积。

例题

一、$n$ 个物品选 $m$ 个，一个数可以被选任意多次，求方案数（不考虑顺序）。

首先考虑每个数只能选一次的情况。将选的数排序，实际上就是求不等式 $1\le a_1< a_2<a_3<....<a_m\le n$ 的解的数量。明显答案就是 $\binom{n}{m}$。而如果可以重复，那么就相当于把 $a$ 之间的 $<$ 换成了 $\le$。

那么我们怎么把 $\le$ 变回去呢？令 $b_i=a_i+i-1$，那么原式就变成了 $1\le b_1<b_2<...<b_m\le n+m-1$。这个时候我们就会惊喜地发现，该不等式的解就是 $\binom{n+m-1}{m}$！然后把 $b_i$ 表示为 $a_i$，所以 $a$ 数组的数量（也就是本题的答案）也是 $\binom{n+m-1}{m}$！

二、求 $\binom{n}{m}\bmod p$ 的值

部分分一： $n,m\le 10^{18},p=1$

弱智档，输出 $0$ 即可。

部分分二：$n,m\le 1000$，$p$ 没有限制

递推求组合数即可。

部分分三：$n,m\le10^6$，$p$ 为质数

递推求出阶乘和阶乘逆元，然后套公式即可。

部分分四：$n\le 10^9,m\le 10^3$，$p$ 没有限制

$m$ 这么小，突破口一定在 $m$ 上。并且可能是一个 $O(m^2)$ 级别的算法。

先约分，得：

\[\binom{n}{m}=\dfrac{\prod_{i=0}^{m-1}n-i}{m!} \]

容易发现，把乘法拆开，分子和分母的项数都是 $O(m)$ 级别的。

由于最后的答案一定是整数，所以我们 $O(m^2)$ 两两约分，直到分母全为 $1$ 为止。最后把分子乘起来就是答案。由于约分要求 gcd，所以时间复杂度 $O(m^2\log n)$。

还是写一下代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
   	int x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
int n,m,fz[1005],fm[1005],mod;
ll ans; 
int gcd(int x,int y){
	return y?gcd(y,x%y):x;
}
signed main(){
	n=read(),m=read(),mod=read();
	for(int i=1;i<=m;++i) fz[i]=n-i+1,fm[i]=i;
	for(int i=1;i<=m;++i){
		for(int j=1;j<=m;++j){
			int g=gcd(fz[i],fm[j]);
			fz[i]/=g,fm[i]/=g;
		}
	}
	ans=1;
	for(int i=1;i<=m;++i) ans=(ans*fz[i])%mod;
	println(ans);
	return 0;
}

部分分五：$n,m\le10^9,p\le100$ 且为质数

这里引入一种新的定理：Lucas 定理。

内容：

\[\binom{n}{m}\equiv\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\times\binom{n\bmod p}{m\bmod p}\pmod p \]

其中因为 $p$ 比较小，$\binom{n\bmod p}{m\bmod p}$ 可以直接 $O(p^2)$ 递推算出来，而剩下的部分继续递归。

终止条件：$m=0$ 时，返回 $1$。

实际上，观察上述式子，发现答案实际上就是将 $n$ 和 $m$ 的 $p$ 进制表示按位求组合数，然后相乘，比如：

\[n=221,m=110(\text{三进制}),p=3 \]

那么

\[\binom{n}{m}\equiv\binom{2}{1}\times\binom{2}{1}\times\binom{1}{0}\pmod p \]

这样我们就得到了一个迭代做法：先暴力预处理出组合数，用短除法将 $n,m$ 分解为 $p$ 进制，然后逐位计算即可。

时间复杂度 $O(p^2+\log_p n)$。

代码：

#include<bits/stdc++.h>
using namespace std;
int n,m,p,C[105][105];
int num_a[105],tot_a,num_b[105],tot_b;
inline int Lucas(int a,int b,int p){
	memset(num_a,0,sizeof(num_a));
	memset(num_b,0,sizeof(num_b));
	tot_a=tot_b=0;
	while(a) num_a[++tot_a]=a%p,a/=p;
	while(b) num_b[++tot_b]=b%p,b/=p;
	int ans=1;
	for(int i=1;i<=tot_a;++i) ans=1ll*ans*C(num_a[i],num_b[i])%p;
	return ans;
}
signed main(){
	scanf("%d%d%d",&n,&m,&p);
	C[0][0]=1;
	for(int i=1;i<p;++i){
		C[i][0]=1;
		for(int j=1;j<p;++j) C[i][j]=(C[i-1][j]+C[i-1][j-1])%p;
	}
	printf("%d\n",Lucas(n,m,p));
	return 0;
}

实际上直接处理阶乘和阶乘逆元就可以达到 $O(p+\log_p n)$。下面给出 Luogu P3807 的代码：

#include<bits/stdc++.h>
const int N=1e5+5;
using namespace std;
int T,n,m,p,fac[N],invfac[N];
int num_a[105],tot_a,num_b[105],tot_b;
inline int ksm(int a,int b){
	int res=1;a%=p;
	while(b){
		if(b&1)res=1ll*res*a%p;
		a=1ll*a*a%p,b>>=1;
	}
	return res;
}
inline int C(int n,int m){
	if(n<m||n<0||m<0)return 0;
	return 1ll*fac[n]*invfac[m]%p*invfac[n-m]%p;	
}
inline int Lucas(int a,int b,int p){
	memset(num_a,0,sizeof(num_a));
	memset(num_b,0,sizeof(num_b));
	tot_a=tot_b=0;
	while(a) num_a[++tot_a]=a%p,a/=p;
	while(b) num_b[++tot_b]=b%p,b/=p;
	int ans=1;
	for(int i=1;i<=tot_a;++i) ans=1ll*ans*C(num_a[i],num_b[i])%p;
	return ans;
}
signed main(){
	scanf("%d",&T);
	while(T--){
		scanf("%d%d%d",&n,&m,&p);
		fac[0]=1,n+=m;
		for(int i=1;i<p;++i) fac[i]=1ll*fac[i-1]*i%p;
		invfac[p-1]=ksm(fac[p-1],p-2);
		for(int i=p-2;i>=0;--i) invfac[i]=1ll*invfac[i+1]*(i+1)%p;
		printf("%d\n",Lucas(n,m,p));
	}
	return 0;
}

假设要求 $\binom{n}{m}\bmod 30$ 的结果，怎么办？

质因数分解，$30=2\times3\times5$，然后分别计算出 $\binom{n}{m}\bmod 2,3,5$ 的结果 $x_1,x_2,x_3$，那么就有：

\[\begin{cases} \binom{n}{m} \equiv x_1\ ({\rm mod}\ 2) \\ \binom{n}{m}\equiv x_2\ ({\rm mod}\ 3) \\ \binom{n}{m} \equiv x_3\ ({\rm mod}\ 5)\end{cases} \]

然后直接上 CRT 就行。

实际上，当 $p$ 可以分解为若干个不相同的质数相乘时，都可以这么做。

但是如果 $p$ 的某个质因子的指数 $>1$，那么模数就不是质数，就不行了。

习题

一、给定 $n\le10^9,k\le10^3$，让你构造一组解，使得 $n$ 可以表示为 $k$ 个不同的组合数之和（上面和下面有一个不同就行）。

$n<k$，无解；因为组合数最小都是 $1$；
$n\ge k$，构造 $k-1$ 个 $1$，然后 $n-k+1$。$1$ 直接狂填 $\binom{x}{0}$，$n-k+1$ 直接 $\binom{n-k+1}{1}$ 即可。

原题是 Luogu P4369。

代码：

#include<bits/stdc++.h>
using namespace std;
int x,k;
signed main(){
	scanf("%d%d",&x,&k);
	for(int i=1;i<k;++i) printf("%d %d\n",i,0);
	printf("%d %d\n",x-k+1,1); 
	return 0;
}

二、比较 $\binom{n_1}{m_1}$ 和 $\binom{n_2}{m_2}$ 的大小（$n_1,n_2,m_1,m_2\le10^6$）。

对数函数的性质：$\log_x(ab)=\log_x a+\log_x b$。

同理，有 $\log_x(\frac{a}{b})=\log_x a-\log_x b$。

如果 $a<b$，那么 $\log_x a<\log_x b$，反过来同理。

这样，我们只需要求 $\log_x\binom{n_1}{m_1}$ 和 $\log_x\binom{n_2}{m_2}$ 并比较大小即可。

我们算一下：

\[\log_x\binom{n}{m} \]

\[=\log_xn!-\log_xm!-\log_x(n-m)! \]

现在怎么算阶乘的 $\log$？

\[n!=(n-1)!\times n \]

所以

\[\log_xn!=\log_x(n-1)!+\log_xn \]

然后我们就可以递推算出阶乘的 $\log$ 了。

这样我们就可以比较大小了。

至于底数 $x$ 选多少，随便。

实际上类似于放缩法，缩小答案的范围然后计算。

代码（把缺省源砍了）：

#include<bits/stdc++.h>
//typedef __uint128_t ulll;
const int N=1e6+5;
using namespace std;
int n1,m1,n2,m2;
double fac[N];
inline double logC(int n,int m){
	return fac[n]-fac[m]-fac[n-m];
}
signed main(){
	cin>>n1>>m1>>n2>>m2;
	fac[0]=0;
	for(int i=1;i<=N-5;++i) fac[i]=fac[i-1]+log(i);
	double _1=logC(n1,m1),_2=logC(n2,m2);
	if(_1>_2) printf("A > B\n");
	else if(_1==_2) printf("A = B\n");
	else printf("A < B\n");
	return 0;
}

时间复杂度 $O(n)$。

还有个问题：这里用的 double，精度怎么办？

精度误差影响的是小数点后好几位，和这个题没关系。

三、选 $k$ 个不同的组合数 $\binom{a}{b}$，使得它们的和最大。$1\le k\le 10^5,0\le b\le a\le 10^6$。

首先，最大的组合数是什么？

考虑杨辉三角的性质，明显是最后一行最中间的那个数（有偶数列的话就除二即可）。

那第二大的呢？

明显是在最大的上下左右。现在不知道谁更大，比较一下即可（上一题刚说过）。

然后从第二大的数又可以往四周扩充几个点。维护一个候选集，每次从候选集里面挑选一个最大的，取出来之后更新候选集，重复 $k$ 次即可。用结构体维护行和列，重载运算符然后优先队列维护即可。

很明显，这是一个类似 BFS 的过程。

注意判重。

Luogu 原题是 P4370。

代码（我每次加入的是当前点周围八联通的点）：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=1e6+5;
const int dx[8]={-1,-1,0,1,1,1,0,-1};
const int dy[8]={0,1,1,1,0,-1,-1,-1};
const int mod=1e9+7;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
   	int x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
int n,k,facc[N],ifac[N];
double fac[N];
map<int,map<int,bool> >mp;
inline ll ksm(ll a,ll b){
	ll res=1;
	a%=mod;
	while(b){
		if(b&1)res=res*a%mod;
		a=a*a%mod,b>>=1;
	}
	return res;
}
inline double logC(int n,int m){
    return fac[n]-fac[m]-fac[n-m];
}
inline bool outmp(int x,int y){
	return x<0||x>n||y<0||y>n;
}
struct Binom{
	int x,y;
	// C(x,y)
	Binom(int _x=0,int _y=0):x(_x),y(_y){}
	inline bool operator<(const Binom &b)const{return logC(x,y)<logC(b.x,b.y);}
};
priority_queue<Binom,vector<Binom>,less<Binom> >pq;
inline ll C(int n,int m){
	if(n<m||n<0||m<0)return 0;
	return 1ll*facc[n]*ifac[m]%mod*ifac[n-m]%mod;
}
signed main(){
	n=read(),k=read(),fac[0]=0,facc[0]=1;
    for(int i=1;i<=n;++i) fac[i]=fac[i-1]+log(i),facc[i]=1ll*facc[i-1]*i%mod;
    ifac[n]=ksm(facc[n],mod-2);
    for(int i=n-1;i>=0;--i) ifac[i]=1ll*ifac[i+1]*(i+1)%mod;
    pq.push(Binom(n,(n+1)/2));
    int cnt=0;
    ll ans=0;
    while(cnt<k){
    	int x=pq.top().x,y=pq.top().y;
    	pq.pop();
    	if(mp[x][y])continue;
    	mp[x][y]=1;
    	ans=(ans+C(x,y))%mod,cnt++;
    	for(int i=0;i<8;++i){
    		int xx=x+dx[i],yy=y+dy[i];
    		if(outmp(xx,yy))continue;
    		pq.push(Binom(xx,yy));
		}
	}
	println(ans);
	return 0;
}

四、P3746 [六省联考 2017] 组合数问题

原题链接。

这个题 $k$ 很小，我们考虑从 $k$ 入手。

我们记

\[f(n,r)=\sum_{i=0}^{\infty}\binom{nk}{ik+r}\bmod p \]

现在我们来介绍几个求和变形的技巧：

一、增加枚举量

把 $f(n,r)$ 展开 $k$ 次，得：

\[f(n,r)=\sum_{i=0}^{\infty}\sum_{j=0}^{k}\binom{k}{j}\times\binom{nk-k}{ik+r-j}\bmod p \]

二、交换枚举顺序

明显交换两个 $\sum$ 不影响答案，所以：

\[f(n,r)=\sum_{j=0}^{k}\sum_{i=0}^{\infty}\binom{k}{j}\times\binom{nk-k}{ik+r-j}\bmod p \]

三、分离无关变量

容易发现此时 $\binom{k}{j}$ 和内层的 $\sum$ 没有关系了，提出来：

\[f(n,r)=\sum_{j=0}^{k}\binom{k}{j}\sum_{i=0}^{\infty}\binom{nk-k}{ik+r-j}\bmod p \]

然后我们把式子整理一下：

\[f(n,r)=\sum_{j=0}^{k}\binom{k}{j}\sum_{i=0}^{\infty}\binom{(n-1)k}{ik+(r-j)}\bmod p \]

我们会发现，第二个 $\sum$ 正好是 $f(n-1,r-j)$！所以：

\[f(n,r)=\sum_{j=0}^{k}\binom{k}{j}f(n-1,r-j)\bmod p \]

我们使用一个经典套路：给他扩充一维，凑成矩阵乘法：

\[f(n,0,r)=\sum_{j=0}^{k}f(n-1,1,r-j)\binom{k}{j}\bmod p \]

这个时候我们把 $f(n)$ 看成一个 $1\times (k+1)$ 的矩阵，那么我们实际上需要凑出一个矩阵 $D$，使其满足：

\[f_n[0][r]=\sum_{j=0}^{k}f_{n-1}[0][r-j]\times D[r-j][r]\bmod p \]

$D$ 是一个 $(k+1)\times (k+1)$ 的矩阵。

对比原式，只需要令 $D[r-j][j]=\binom{k}{j}$ 即可。

这样就能用矩阵乘法了。

有 $f_n=f_0\times D^n$。$D$ 已经构造出来了，现在的问题就是构造 $f_0$。

我们有

\[f(0,0,r)=\sum_{i=0}^{\infty}\binom{0k}{nk+r}=\sum_{i=0}^{\infty}\binom{0}{nk+r} \]

很明显，只有当 $r=0$ 的时候式子是 $1$，其余为 $0$。

所以，$f_0$ 是一个 $1\times(k+1)$ 的矩阵，且只有 $f_0[0][0]=1$，其余都是 $0$。

但还有一个问题：$r-j$ 可能 $<0$。

我们考虑证明：$j>0$ 时，$f(i,0,-j)=f(i,0,k-j)$。

把 $f(i,0,k-j)$ 展开，得：

\[\sum_{l=0}^{\infty}\binom{ik}{lk+k-j} \]

\[=\sum_{l=0}^{\infty}\binom{ik}{(l+1)k-j} \]

\[=\sum_{l=1}^{\infty}\binom{ik}{lk-j} \]

实际上，因为 $\binom{ik}{0k-j}=\binom{ik}{-j}=0$，所以原式

\[=\sum_{l=0}^{\infty}\binom{ik}{lk-j} \]

\[=f(i,0,-j) \]

这样，如果 $r-j<0$，我们就可以直接给 $D[r-j+k][j]$ 加上 $\binom{k}{j}$ 了（不能直接赋值，因为一个位置可能会被重复算多次）。此时，$D$ 的两维就都是 $[0,k]$ 之间的数了。

然后我们就可以直接进行矩阵乘法了（由于预处理时已经解决了负数的情况，所以直接从 $0$ 到 $k$ 枚举下标即可）。

这样我们就做完了。时间复杂度 $O(k^3\log n)$。

代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=55;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
   	int x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
int n,p,k,r,C[N][N];
struct Matrix{
	int n,m;
	ll val[N][N];
	Matrix(int _n=0,int _m=0){
		n=_n,m=_m;
		memset(val,0,sizeof(val));
	}
}f0,D;
inline Matrix operator*(const Matrix &a,const Matrix &b){
	Matrix c(a.n,b.m);
	for(int i=0;i<c.n;++i){
		for(int k=0;k<a.m;++k){
			for(int j=0;j<c.m;++j){
				c.val[i][j]=(c.val[i][j]+a.val[i][k]*b.val[k][j]%p)%p;
			}
		}
	}
	return c;
}
inline Matrix ksm(Matrix a,int b){
	Matrix c=a;
	b--;
	while(b){
		if(b&1)c=c*a;
		a=a*a,b>>=1;
	}
	return c;
}
signed main(){
	n=read(),p=read(),k=read(),r=read();
	C[0][0]=1;
	for(int i=1;i<=k;++i){
		C[i][0]=1;
		for(int j=1;j<=k;++j) C[i][j]=(C[i-1][j]+C[i-1][j-1])%p;
	}
	f0.n=1,f0.m=k+1;
	f0.val[0][0]=1;
	for(int j=1;j<=k;++j) f0.val[0][j]=0;
	D.n=D.m=k+1;
	for(int r=0;r<=k;++r){
		for(int j=0;j<=k;++j){
			if(r>=j) D.val[r][r-j]+=C[k][j];
			else D.val[r][r-j+k]+=C[k][j];
		}
	}
	Matrix F=f0*ksm(D,n);
	println(F.val[0][r]);
	return 0;
}

抽屉原理

把 $n+1$ 个东西放到 $n$ 个抽屉里，必然有一个抽屉有 $2$ 个东西（反证）。

推广：把 $kn+1$ 个东西放到 $n$ 个抽屉里，必然有一个抽屉有 $k+1$ 个东西（同理）。

习题

一、POJ 3370

给定 $n$ 个数，要求从中选出任意多个数（至少一个），使得它们的和是 $c$ 的倍数。$c\le n\le10^5$。

奇怪的地方是 $c\le n$，我们考虑把 $c$ 作抽屉，$n$ 作物品。

构造：强制令选的数连续。做一遍前缀和，那么需要满足 $sum_r-sum_{l-1}\equiv0\pmod c$。这样会产生 $n+1$ 个前缀和（$0$ 到 $n$），考虑把它们 $\bmod c$ 之后的结果分类。由于 $n\ge c$，所以 $n+1>c$，有抽屉原理可知必然有两个前缀和会存在同一组里。这两个前缀和对应的区间就是一组解。那么把每一组开一个 vector，把前缀和的编号插入里面，最后找到 size()>1 的一组，取出里面的两个编号即可。

代码：

#include<cstdio>
#include<vector>
const int N=1e5+5;
using namespace std;
int c,n,a[N],sum[N];
vector<int>box[N];
signed main(){
	while(1){
		scanf("%d%d",&c,&n);
		if(!c&&!n)break;
		for(int i=0;i<c;++i) box[i].clear();
		sum[0]=0;
		for(int i=1;i<=n;++i){
			scanf("%d",a+i);
			sum[i]=(sum[i-1]+a[i]%c)%c;
		}
		for(int i=0;i<=n;++i) box[sum[i]%c].push_back(i);
		for(int i=0;i<c;++i){
			if(box[i].size()>1){
				int l=box[i][0],r=box[i][1]; 
				for(int i=l+1;i<=r;++i) printf("%d ",i);
				break;
			}
		}
		printf("\n");
	}
	return 0;
}

二、P2218 [HAOI2007]覆盖问题

原题链接。

平面上有 $n$ 个点 $(x_i,y_i)$，要用三个大小为 $L\times L$ 的正方形覆盖所有点（可以不在格点上，可重叠，平行于坐标轴），求最小的 $L$（$L$ 不一定是整数）。

首先发现答案满足单调性，考虑二分。

现在的问题就是如何 check：找到最左、最右、最上、最下的四个点，然后画一个大长方形；根据抽屉原理，必然有一个正方形要盖住两个。很明显这个正方形要放在大长方形的角上。由于点数很小，可以直接暴力枚举盖哪个角。使用 dfs 递归三层即可，具体来说：

第一层递归：找到四个极限点，枚举盖住哪个角，进入下一层；

第二层递归：找到没被覆盖的四个极限点，枚举盖住哪个角；

第三层递归：找到没被覆盖的四个极限点，判断横纵坐标的最大差值是否 $\le L$。

注意使用下一个正方形的时候要清空当前被覆盖的点。

代码：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=2e4+5;
const int inf=1e9+1; 
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline int read(){
   	int x=0;bool sgn=0;char ch=gt();
   	while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
   	while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
	return sgn?-x:x;
}
template<class T>
inline void print(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
	static char st[70];short top=0;
	if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
	int siz=s.size();
	for(int i=0;i<siz;++i) pt(s[i]);
	printf("\n");
}
int n,flag[N];
struct Point{
	int x,y;
	Point(int _x=0,int _y=0):x(_x),y(_y){}
	inline void reset(){x=y=0;}
}a[N];
struct Square{
	Point up_left,down_right; 
	Square(){
		up_left.reset();
		down_right.reset();
	}
};
inline bool in(Point a,Square b){
	// 判断 a 是否在 b 内
	return a.x>=b.up_left.x&&a.x<=b.down_right.x&&a.y>=b.down_right.y&&a.y<=b.up_left.y;
}
bool dfs(int k,int dep){
	// k * k，当前是第 dep 个正方形
	Point L(inf,0),R(-inf,0),U(0,-inf),D(0,inf);
	for(int i=1;i<=n;++i){
		if(a[i].x<L.x&&!flag[i]) L=a[i];
		if(a[i].x>R.x&&!flag[i]) R=a[i];
		if(a[i].y>U.y&&!flag[i]) U=a[i];
		if(a[i].y<D.y&&!flag[i]) D=a[i];
	}
	if(max(1ll*R.x-L.x,1ll*U.y-D.y)<=k) return 1;
	if(dep==3) return 0;
	Square r[4];
	r[0].up_left=Point(L.x,U.y),r[0].down_right=Point(L.x+k,U.y-k);// 左上角
	r[1].up_left=Point(R.x-k,U.y),r[1].down_right=Point(R.x,U.y-k);// 右上角
	r[2].up_left=Point(L.x,D.y+k),r[2].down_right=Point(L.x+k,D.y);// 左下角
	r[3].up_left=Point(R.x-k,D.y+k),r[3].down_right=Point(R.x,D.y);// 右下角
	for(int i=0;i<4;++i){
		// 使用第 i 个正方形
		for(int j=1;j<=n;++j) if(!flag[j]&&in(a[j],r[i])) flag[j]=dep;// 覆盖 
		if(dfs(k,dep+1)) return 1;
		for(int j=1;j<=n;++j) if(flag[j]==dep) flag[j]=0;
	}
	return 0;
}
inline bool check(int k){
	for(int i=1;i<=n;++i) flag[i]=0;
	return dfs(k,1);
}
signed main(){
	n=read();
	for(int i=1;i<=n;++i) a[i].x=read(),a[i].y=read();
	int l=0,r=2*inf;
	while(l<r){
		int mid=l+((r-l)>>1);
		if(check(mid)) r=mid;
		else l=mid+1;
	}
	println(l);
	return 0;
}

容斥原理

$A\cap B$：$A$ 和 $B$ 的交集。

$A\cup B$：$A$ 和 $B$ 的并集。

$|A|$：集合 $A$ 的大小。

那么我们有：

\[|A\cup B|=|A|+|B|-|A\cap B| \]

同理可得：

\[|A\cup B\cup C|=|A|+|B|+|C|-|A \cap B|-|A\cap C|-|B\cap C|+|A\cap B\cap C| \]

推广到 $n$ 个集合 $A_1,A_2,...,A_n$：

\[\sum_{B\in(A_1,A_2,...,A_n)}(-1)^{|B|+1}\times|\cap_{A_i\in B}A_i| \]

$n$ 个人排成一排的方案数：$n!$

$n$ 个人排成一圈，且旋转后相同算一种方案的方案数：$(n-1)!$（每一种方案有 $n$ 种旋转结果；或定死一个人的位置，就是剩下 $n-1$ 个人的排列）。

习题

一、N 对夫妻问题

$n$ 对夫妻，总共 $2n$ 个人。把他们排成一圈，满足每对夫妻都不能相邻的位置且旋转后相同的方案算一种，求方案数。

$n$ 对夫妻都不相邻是一个很强的条件，那么考虑容斥原理：首先，$2n$ 个人的圆排列数量是 $(2n-1)!$，然后，考虑强制让一对夫妻相邻，有：

\[ans=(2n-1)!-n\times(2n-2)!\times2 \]

$n$ 是从 $n$ 对夫妻挑一对，$(2n-2)!$ 是因为可以把这对夫妻看成一个人，圆排列数量是 $(2n-2)!$，$2$ 是因为一对夫妻有两种排法。

做完了吗？并没有，因为情况会算重。因此，加回来（使两对夫妻强制相邻）：

\[ans=(2n-1)!-n\times(2n-2)!\times2+\binom{n}{2}\times (2n-3)!\times2^2 \]

同理，最终可得：

\[ans=\sum_{i=0}^{n}(-1)^i\times\binom{n}{i}\times(2n-i-1)!\times2^i \]

这个 $i$ 的含义是强制让 $i$ 对夫妻相邻。

二、P9118 [春季测试 2023] 幂次/ HDU 2204

Luogu P9118 链接。

可以啊，春测出原题。

首先我们发现 $1$ 这个数情况比较复杂，所以我们直接拿出来。

记 $ans_i$ 表示指数 $i$ 的答案。对于指数 $i$，有 $ans_i=\lfloor\sqrt[i]{n}\rfloor-1$（$-1$ 是把 $1$ 去掉），这个用 powl(1.0*n,1.0/i) 就可以求出来（是 long long，一定要开 powl！！！）。

然后发现 $i$ 的倍数次幂会被算多次，所以减去，有：

\[ans_i=\lfloor\sqrt[i]{n}\rfloor-1-\sum_{j=2}^{ij\le100}dp_{i\times j} \]

不超过一百是因为 $k$ 最大只有 $100$。

我们发现编号较小的 $dp$ 值是由编号较大的 $dp$ 值推出来的，所以从大到小枚举即可。

最后再加上 $1$，答案为：

\[(\sum_{i=k}^{100}ans_i)+1 \]

代码：

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
ll n,ans,dp[105];
int k;
signed main(){
	scanf("%lld%d",&n,&k);
	for(int i=100;i>=k;--i){
		dp[i]=powl(1.0*n,1.0/i)-1;
		for(int j=2;j*i<=100;++j) dp[i]-=dp[i*j];
	}
	ans=1;
	for(int i=k;i<=64;++i) ans+=dp[i];
	printf("%lld\n",ans);
	return 0;
}

HDU 2204 就把 $k$ 改成 $2$ 即可。

矩阵 - 二

高斯消元

原理已经写过了，见本文最开头的链接。

用处：

一个题，输入一个数，输出一个数，通常会怎么做？

打表！

$n\le10^9$，怎么办！

~~分块打表！~~

~~OEIS！~~

找规律！

这个过程就可以用高消。

比如，一个数列在 $n=1,2,3,4,5$ 时候的值是 $0,5,20,51,104$，现在我们要找出这个数列的通项公式（禁用 OEIS）。

我们猜测：答案是形如 $an+b$ 的形式。把 $n=1,2$，带进去：

\[\begin{cases} a+b=0 \\ 2a+b=5 \end{cases} \]

解出来，得：

\[\begin{cases} a=5 \\ b=-5 \end{cases} \]

试一下，好像不对。那就加一次，加个未知数。最终可以试出来 $a_n=n^3-n^2+n-1$。

这个解方程的过程就可以使用高斯消元。

限制条件：答案必须是形如若干个 $n$ 的次方相加的形式。

这实际上是一种骗分技巧。

矩阵求逆

定义：见模板题。

考虑如下矩阵：

\[\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \]

容易发现，它和单位矩阵的唯一区别就是第 $2$ 行的 $1$ 和第 $3$ 行的 $1$ 的列号互换。

如果用它左乘一个矩阵，那么将会交换这个矩阵的第二、第三行。

再考虑如下矩阵：

\[\begin{bmatrix} 1 & 0\\ k & 1 \\ \end{bmatrix} \]

如果它左乘一个矩阵，那么将会把那个矩阵的第一行 $\times k$ 再加到第二行上。

我们会发现，这正好是高斯消元的操作！

这样，高斯消元的过程就可以转化为矩阵乘法了。

我们考虑把 $A$ 消成 $I$。由逆矩阵的性质，这相当于乘了一个 $A^{-1}$。对 $I$ 做一遍同样的操作，我们会发现它最终就是 $A^{-1}$！（$I\times A^{-1}=A^{-1}$）

怎么一起消呢？把两个矩阵放一块即可。

代码还是在文章最开头的博客里。

概率与期望

概率

假设我们要扔骰子，那么扔出来的结果 $1,2,3,4,5,6$ 就被称为样本空间。而其中的每个情况被称为样本点。事件是若干样本点的集合。

既然事件是集合，那么它们之间也可以进行运算。设 $A=\{1,2,3\},B=\{2,3,4\}$，那么

\[A\cap B=A\cdot B=\{2,3\} \]

\[A\cup B=A+B=\{1,2,3,4\} \]

那减法呢？

$A-B$ 就是把 $B$ 里的东西在 $A$ 中扣掉之后剩下来的东西，即：

\[A-B=A-AB=\{1\} \]

样本空间内每个样本点都有一个概率，每个样本点的概率不一定相等。一个事件的概率就是事件内样本点的概率之和。

接下来介绍概率的一些性质：

$P(A)$ 表示事件 $A$ 发生的概率。很明显，有 $0\le P(A)\le 1$。
假设样本空间内共有 $n$ 个样本点 $B_1,B_2,...,B_n$，那么有 $\sum_{i=1}^{n}P(B_i)=1$.
如果 $P(A)=1$，称 $A$ 为必然事件；若 $P(A)=0$，那么 $A$ 为不可能事件。

条件概率：$P(A|B)$ 表示 $B$ 发生时 $A$ 发生的概率。比如，$P(\text{扔出 1}|\text{扔出的点数在 1 到 3 之间})=\dfrac{1}{3}$。

再举个例子，如果 $A=\{1,2,3\},B=\{1,3,5\}$，那么 $P(A|B)=\dfrac{2}{3}$（$3$ 个里面有 $2$ 个满足条件）。

实际上，条件概率有一个更简单的公式：

\[P(A|B)=\dfrac{P(AB)}{P(B)} \]

解释：如图，

如果 $B$ 要发生，必然要落在绿色区域里面；而 $A$ 也要发生，所以最终必然要落在蓝色区域里面。由面积比可得，$P(A|B)=\dfrac{P(AB)}{P(B)}$

所以，还可以得出

\[P(A|B)P(B)=P(AB) \]

独立事件：如果 $A$ 事件不影响 $B$ 事件发生，那么就称它们为独立事件。实际上，如果 $A,B$ 独立，有：

\[P(A)P(B)=P(AB) \]

实际上，对比条件概率的公式，有：

\[P(A)P(B)=P(A|B)P(B)=P(AB) \]

\[P(A)=P(A|B) \]

这也可以解释：由于这两个事件独立，所以 $B$ 事件的发生对 $A$ 事件的发生没有影响。

比如，扔两次骰子，这两个事件就是互相独立的。

期望

还是扔骰子，每个数的概率都是 $\dfrac{1}{6}$，求扔出来的数的期望。

期望，实际上就是每个事件的权值（比如扔骰子的期望） $\times$ 概率再加起来，所以扔骰子的期望值就是：

\[(1+2+3+4+5+6)\times\dfrac{1}{6}=3.5 \]

那如果是扔出来的数的平方的期望呢？同理，为：

\[(1^2+2^2+3^2+4^2+5^2+6^2)\times\dfrac{1}{6}=\dfrac{91}{6} \]

接下来，介绍一下期望的性质：

期望的和 $=$ 和的期望。比如我们有两颗骰子：第一颗是正常的，每种情况概率都是 $\dfrac{1}{6}$；第二颗灌了铅，扔出 $6$ 的概率是 $\dfrac{95}{100}$，其他都是 $\dfrac{1}{100}$，那么：

\[E(x_1+x_2)=E(x_1)+E(x_2) \]

实际上，这在任何情况下都适用：

\[E(x_1+x_1^2)=E(x_1)+E(x_1^2) \]

这也是对的。这样，就可以大大降低运算量。

严谨证明？不会

习题

一、课件 P111

某小区有男 $51$ 人女 $49$ 人，男人色盲的概率为 $\dfrac{2}{100}$，女人色盲的概率为 $\dfrac{0.25}{100}$，现在抽到一个色盲，问这个人是男人的概率。

很明显，这是求 $P(\text{这个人是男人}|\text{这个人是色盲})$，根据条件概率的公式，有：

\[P(\text{这个人是男人}|\text{这个人是色盲})=\dfrac{P(\text{这个人是男人}\cdot \text{这个人是色盲})}{P(\text{这个人是色盲})} \]

\[=\dfrac{\dfrac{51}{100}\times\dfrac{2}{100}}{P(\text{这个人是男色盲})+P(\text{这个人是女色盲})} \]

\[=\dfrac{\dfrac{51}{100}\times\dfrac{2}{100}}{\dfrac{51}{100}\times\dfrac{2}{100}+\dfrac{49}{100}\times\dfrac{0.25}{100}} \]

剩下懒得约分了。

二、课件 P117 - 小葱过河

问题在于比较 $(\dfrac{99}{100})^{100}$ 和 $(\dfrac{999}{1000})^{1000}$ 的大小。

发现指数可以约掉，问题变为比较 $\dfrac{99}{100}$ 和 $(\dfrac{999}{1000})^{10}$ 的大小。

考虑放缩法，由于 $\dfrac{999}{1000}>\dfrac{998}{999}>...>\dfrac{990}{991}$，所以

\[(\dfrac{999}{1000})^{10}>\prod_{i=990}^{999}\dfrac{i}{i+1} \]

约分，得：

\[(\dfrac{999}{1000})^{10}>\dfrac{99}{100} \]

所以，答案是~~不过河~~走那 $1000$ 个石头。

三、课件 Question 11 - 小胡抛硬币

操作二和操作三之间的那个点非常关键，我们考虑枚举它的坐标。

显然，最坏情况下会走无限步，所以 $x,y$ 都有可能达到无限。

首先，我们需要走到 $(x,0)$，这意味着我们需要连续抛 $x$ 次反面，并再抛一次正面停下来。这部分的概率为 $(1-p)^xp$
然后，我们需要走到 $(x,y)$，这意味着我们需要连续抛 $y$ 次反面，并再抛一次正面停下来。这部分的概率为 $(1-p)^yp$
最后，我们需要原路返回。很明显，我们一共需要抛 $x+y$ 次，且需要抛 $x$ 次正面， $y$ 次反面，方案数为 $\binom{x+y}{x}$。这部分的概率为 $q^x(1-q)^y\binom{x+y}{x}$。这实际上就是二项式分布。

把上面三个部分乘起来（不是一个阶段），然后枚举 $x,y$，有：

\[ans=\sum_{x=0}^{\infty}\sum_{y=0}^{\infty}(1-p)^xp (1-p)^ypq^x(1-q)^y\binom{x+y}{x} \]

我们考虑用求和变形技巧化简式子。

首先，把能合并的合并，有：

\[ans=\sum_{x=0}^{\infty}\sum_{y=0}^{\infty}(1-p)^{x+y}p^2q^x(1-q)^y\binom{x+y}{x} \]

然后我们介绍第四个求和变形的技巧：

换元法：我们令 $t=x+y$，然后枚举 $t$ 和 $x$。有：

\[ans=\sum_{t=0}^{\infty}\sum_{x=0}^{t}(1-p)^tp^2(1-q)^{t-x}q^x\binom{t}{x} \]

分离无关变量：

\[ans=p^2\sum_{t=0}^{\infty}(1-p)^t\sum_{x=0}^{t}(1-q)^{t-x}q^x\binom{t}{x} \]

这是什么？二项式定理！

所以：

\[ans=p^2\sum_{t=0}^{\infty}(1-p)^t[q+(1-q)]^t \]

\[ans=p^2\sum_{t=0}^{\infty}(1-p)^t1^t \]

\[ans=p^2\sum_{t=0}^{\infty}(1-p)^t \]

剩下的是什么？是等比数列。

我们可以使用等比数列求和公式：

\[\sum_{i=0}^{n-1}a^i=\begin{cases}n & a=1\\\dfrac{a^n-1}{a-1}& a\not=1 \end{cases} \]

所以：

\[ans=p^2\dfrac{(1-p)^{\infty}-1}{1-p-1} \]

分子分母取相反数，得：

\[ans=p^2\dfrac{1-(1-p)^{\infty}}{1-(1-p)} =p^2\dfrac{1-(1-p)^{\infty}}{p} \]

在极限意义下，$(1-p)^{\infty}$ 可以看做 $0$，所以：

\[ans=p^2\dfrac{1-0}{p}=p^2\dfrac{1}{p}=p \]

问题来了：当 $p=0$ 时，$\dfrac{1}{p}$ 无意义，怎么办？

实际上，当 $p=0$ 时，在使用等比数列求和之前就可以变为 $0$，所以还是对的。

四、BZOJ 2396

可以在 Dark BZOJ 上交。网址。

梯子在这里。

我们可以构造一个 $n\times 1$ 的矩阵（向量） $D$，很明显，如果 $A\times B=C$，那么：

\[A\times B\times D=C\times D \]

首先，$C\times D$ 是 $O(n^2)$ 的，对于左边的式子，由于矩阵有结合律，所以我们可以先 $O(n^2)$ 算出 $B\times D$，再 $O(n^2)$ 算出 $A\times B\times D$，可以接受。

但是反过来对吗？

明显不对，比如 $D$ 是零矩阵。

但是我们多随几次，正确率就很高了。

代码（我随了十次）：

#include<bits/stdc++.h>
//#pragma GCC optimize("Ofast")
#define gt getchar
#define pt putchar
#define y1 y233
//typedef unsigned int uint;
typedef long long ll;
typedef unsigned long long ull;
//typedef __int128 lll;
//typedef __uint128_t ulll;
const int N=1005;
using namespace std;
inline bool __(char ch){return ch>=48&&ch<=57;}
inline ll read(){
    ll x=0;bool sgn=0;char ch=gt();
    while(!__(ch)&&ch!=EOF){sgn|=(ch=='-');ch=gt();}
    while(__(ch)){x=(x<<1)+(x<<3)+(ch-48);ch=gt();}
    return sgn?-x:x;
}
template<class T>
inline void print(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);
}
template<class T>
inline void printsp(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(32);
}
template<class T>
inline void println(T x){
    static char st[70];short top=0;
    if(x<0)pt('-');
    do{st[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top)pt(st[top--]);pt(10);
}
inline void put_str(string s){
    int siz=s.size();
    for(int i=0;i<siz;++i) pt(s[i]);
    printf("\n");
}
int n;
struct Matrix{
    int n,m;
    ll val[N][N];
    Matrix(int _n=0,int _m=0){
        n=_n,m=_m;
        memset(val,0,sizeof(val));
    }
    inline void Reader(){
        for(int i=0;i<n;++i) for(int j=0;j<m;++j) val[i][j]=read();
    }
    inline void Output(){
        for(int i=0;i<n;++i){
            for(int j=0;j<m;++j) printsp(val[i][j]);
            printf("\n");
        }
    }
    inline bool operator==(const Matrix &b)const{
        for(int i=0;i<n;++i) for(int j=0;j<m;++j) if(val[i][j]!=b.val[i][j]) return 0;
        return 1;
    }
}A,B,C;
inline Matrix operator*(const Matrix &a,const Matrix &b){
    Matrix c(a.n,b.m);
    for(int i=0;i<c.n;++i){
        for(int k=0;k<a.m;++k){
            for(int j=0;j<c.m;++j){
                c.val[i][j]+=a.val[i][k]*b.val[k][j];
            }
        }
    }
    return c;
}
mt19937 Rand(time(0));
inline bool check(Matrix a,Matrix b,Matrix c){
    // 检查 a * b 是否 = c
    Matrix Vector(c.m,1); 
    for(int i=0;i<Vector.n;++i) Vector.val[i][0]=Rand();
    return a*(b*Vector)==c*Vector;
}
inline void solve(){
	A.n=C.n=A.m=B.n=B.m=C.m=n;
   	A.Reader(),B.Reader(),C.Reader();
    for(int i=1;i<=10;++i){
       	if(!check(A,B,C)){
           	printf("No\n");
           	return;
       	}
    }
    printf("Yes\n");
}
signed main(){
    while(~scanf("%d",&n)) solve();
    return 0;
}

注意 check 函数的最后，要加括号，不然复杂度就伪了。

这题在 luogu 上有个原：P10102。

不过这题只需要随一次就行了，题解说错误率是 $\frac{1}{\bmod}$，没看懂，好深刻。

贴个码：

#include<bits/stdc++.h>
#include<ext/pb_ds/assoc_container.hpp>
#include<ext/pb_ds/tree_policy.hpp>
#include<ext/pb_ds/hash_policy.hpp>
#define gt getchar
#define pt putchar
#define fst first
#define scd second
#define SZ(s) ((int)s.size())
#define all(s) s.begin(),s.end()
#define pb push_back
#define eb emplace_back
typedef long long ll;
typedef double db;
typedef long double ld;
typedef unsigned long long ull;
typedef unsigned int uint;
const int N=3005;
const int mod=998244353;
using namespace std;
using namespace __gnu_pbds;
typedef pair<int,int> pii;
template<class T,class I> inline void chkmax(T &a,I b){a=max(a,(T)b);}
template<class T,class I> inline void chkmin(T &a,I b){a=min(a,(T)b);}
inline bool __(char ch){return ch>=48&&ch<=57;}
template<class T> inline void read(T &x){
	x=0;bool sgn=0;static char ch=gt();
	while(!__(ch)&&ch!=EOF) sgn|=(ch=='-'),ch=gt();
	while(__(ch)) x=(x<<1)+(x<<3)+(ch&15),ch=gt();
	if(sgn) x=-x;
}
template<class T,class ...I> inline void read(T &x,I &...x1){
	read(x);
	read(x1...);
}
template<class T> inline void print(T x){
	static char stk[70];short top=0;
	if(x<0) pt('-');
	do{stk[++top]=x>=0?(x%10+48):(-(x%10)+48),x/=10;}while(x);
    while(top) pt(stk[top--]);
}
template<class T> inline void printsp(T x){
	print(x);
	putchar(' ');
}
template<class T> inline void println(T x){
	print(x);
	putchar('\n');
}
int T,n,a[N][N],b[N][N],c[N][N],d[N];
int res1[N],res2[N],res3[N];
inline void add(int &a,int b){
	a+=b;
	if(a>=mod) a-=mod; 
}
mt19937 Rand(time(0));
inline void solve(){
	read(n);
	for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) read(a[i][j]);
	for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) read(b[i][j]);
	for(int i=1;i<=n;++i) for(int j=1;j<=n;++j) read(c[i][j]);
	for(int i=1;i<=n;++i) d[i]=Rand()%mod,res1[i]=res2[i]=res3[i]=0;
	for(int i=1;i<=n;++i){
		for(int j=1;j<=n;++j){
			add(res1[i],1ll*d[j]*a[j][i]%mod);
		}
	}
	for(int i=1;i<=n;++i){
		for(int j=1;j<=n;++j){
			add(res2[i],1ll*res1[j]*b[j][i]%mod);
		}
	}
	for(int i=1;i<=n;++i){
		for(int j=1;j<=n;++j){
			add(res3[i],1ll*d[j]*c[j][i]%mod);
		}
	}
	for(int i=1;i<=n;++i) if(res2[i]!=res3[i]) return printf("No\n"),void();
	printf("Yes\n");
}
signed main(){
	read(T);
	while(T--) solve();
	return 0;
}

显然 $O(n^2)$。

但是原来的代码过不了，说明不用封装的时候别 xjb 封装净增常数。

posted @ 2024-02-28 15:22 Southern_Dynasty 阅读(13) 评论(0) 编辑收藏举报

刷新页面返回顶部

syzqwq

清北学堂五一数学营笔记

矩阵

P4159 [SCOI2009] 迷路

一个小语法

数论

判断质数

一、试除法

二、Miller - Rabin

分解质因数

一、朴素做法

二、Pollard's rho 算法

CF45G - Prime Problem

逆元

一、暴力

二、优化

求 \([1,n]\) 逆元

exgcd

CRT/EXCRT

筛法

一、暴力

二、埃氏筛

三、线性筛

积性函数

BSGS 算法

一、暴力

二、正解

组合数学

基础知识

组合

例题

习题

抽屉原理

习题

容斥原理

习题

矩阵 - 二

高斯消元

矩阵求逆

概率与期望

概率

期望

习题

公告