山东游记 - 五一数学专题

CHANGE LOG

• 2023.5.4 更新扩展中国剩余定理代码与解释。
• 2023.5.5 添加线性预处理逆元解释说明。
• 2023.5.7 重构初等矩阵部分。
• 2023.5.11 添加快读模板。
• 2023.5.24 更新高斯消元法代码，添加高斯 - 约旦消元法。
• 2023.6.3 补充逆元部分概念，修改时间复杂度事实性错误。
• 2023.6.9 补充逆元一般性求法，增加扩展卢卡斯定理（exLucas）。
• 2023.7.29 优化求逆元部分文章结构，修改扩展欧几里得事实性错误，对其代码进行优化。

写在前面

本博客中所有 ll 代指 long long，ksm 代指 快速幂 代码，如下：

template <typename _Tp> _Tp ksm(_Tp base,int ind,int mod=mod) {
	_Tp ans=1;
	while(ind) {
		if(ind&1) ans=ans*base%mod;
		base=base*base%mod,ind>>=1;
	}
	return ans%mod;
}

请注意，当前写法是每次 乘法后取模，未拆解大数乘法为加法，这意味着 int 型快速幂计算 不能涉及过程量大于 46341（即 $\lfloor \sqrt{2147483647}\rfloor$）的数字，大数快速幂请自觉送参 long long。

本快速幂模板运用了运算符 * 与 %，并使用了转换构造函数，呼吁广大自编结构体完善构造函数和重载 operator * 与 operator %。

同时，推荐使用如下快读模板：

template <typename _Tp> inline void read(_Tp &x) {
    int f=1; char ch=getchar(); x=0;
    while(!isdigit(ch)) { if(ch=='-') f*=-1; ch=getchar(); }
    while(isdigit(ch)) x=(x<<3)+(x<<1)+(ch^48),ch=getchar();
    x*=f;
}

---

从接触竞赛到现在，不知不觉 8 个月了；如今，春回大地，万物复苏——是时候出去看看了！

于是——

TSOI2022 进军山东！

4.28 启程

上午八点，TSOI2022 师生 5 人从唐一出发，开始了征程。

上午十点半左右正式沿秦滨高速跨入山东，久违地看见了浪漫的浓积云~

下午五点左右安置完毕，晚上还去了海边。在辛安河特大桥边到沙滩上跑了跑。

归程下起了细雨，意外的有趣。

这两天还了解了一些东西，写一下：

• bool 其实是 char，所以理论上是 true 和 false，事实上是 0 或 非 0，是可能出现 0 和 1 之外的数字的。

• 数组越界后操作会在 内存中已有的位置 进行，所以在你数组越界后，其他变量是有可能遭殃的。

• 结构体内的函数可以 随意调用，即使是 顺序在当前函数之后。这意味着你 无需提前声明。

• 结构体和类可以用 *this 来代表当前所在的变量。

---

基础前置

4.29 - Day 1

上午

见到了钟皓曦学长，NOI2013 金牌，清华姚班大佬。
讲了一些比较基础的东西，大部分还都有了解，阿弥陀佛。

那就说几个点吧，一些零碎的知识：

构造函数

在结构体或类中，每次声明变量时执行的函数。声明的 等号 = 和 括号 其实是在调用构造函数。如：

struct high {int len,x[2023];};
high a (0,{0});//显式调用初始化构造函数
high b={0,{0}};//隐式调用初始化构造函数
high c=a;//隐式调用构造函数
high d;//调用默认构造函数
high e=1;//调用转换构造函数
a=1;//赋值
b=a;//赋值

其实都是调用了构造函数。

由上面代码可知，构造函数分为四类：

1. 默认构造函数：声明的时候啥也不说。

2. 初始化构造函数：使用传参构造。

3. 复制构造函数：用同类构造同类。

4. 转换构造函数：使用其他类型的变量构造。

在我们不显式编写构造函数时，系统会 自动生成隐式构造函数。但一旦我们自己编写了构造函数，所有 隐式构造函数全部 失效。

所以上面代码中，就应在结构体中编写：

struct high {
	int len,bit[2023];
	high() {len=1,memset(bit,0,sizeof bit);}//默认构造函数 
	high(int given_len,std::vector <int> given_bit) {
		len=given_len;
		std::copy(given_bit.begin(),given_bit.begin()+given_bit.size(),bit);
	}//初始化构造函数 
	high(int x) {
		int now=0;
		while(x) bit[++now]=x%10,x/=10;
		len=now;
	}//转换构造函数 
};

隐式的复制构造函数始终存在，可以不写。

转换函数

顾名思义，就是用于将当前变量转换为其他类型变量的函数。如：

operator int() const {
		int now=0;
		for(int i=len;i>0;--i)
			now=now*10+bit[i];
		return now;
	}

在刚才的结构体 high 中加入这段代码，就可以实现 high 转 int 了。

取模运算的修正与优化

负数取模后还是负数。正确的做法是，对于 int 型变量 a 和正整数 b，使用 (a%b+b)%b 得出正确结果。

值得一提的是，乘法的速度比除法和取模快。所以，当多次取模时，因为你两个加数 a 和 b 本身是取过模的，所以对于结果 ans 和取模数 mod，可以直接 ans=a+b>mod?a+b-mod:a+b; 总之，面对几乎极限的时间要求时，取模是有可能可以榨出一些油水的。

5.2 补：等比数列快速幂

一种使用类似快速幂求 $\sum _{i=0}^{n-1}{p}^i$ 的方法，我愿称此为 阶乘快速幂（fac_pow）。

出现在 杰杰的女性朋友 一题中，后来好不容易研究明白了。

我们设 $S_n$ 为前 $n$ 项的和，即 $\sum _{i=0}^n{p}_i$，就有

$$S_{\lfloor n/2\rfloor -1}=1+p+p^2+\cdots +p^{\lfloor n/2\rfloor -1},$$

$$S_{n-1}\{\begin{array}{ll}n\equiv 0(\mathrm{mod}2)& =1+p+p^2+\cdots +p^{n-1}\\ & =1+p+\cdots +p^{n/2-1}+p^{n/2}+p^{n/2+1}+\cdots +p^{n-1}\\ & =(1+p^{n/2})S_{n/2-1}\\ n\equiv 1(\mathrm{mod}2)& =1+p+p^2+\cdots +p^{n-1}\\ & =1+p+\cdots +p^{\lfloor n/2\rfloor -1}+p^{\lfloor n/2\rfloor }+p^{\lfloor n/2\rfloor +1}+\cdots +p^{n-1}\\ & =(1+p^{\lfloor n/2\rfloor })S_{\lfloor n/2\rfloor -1}+p^{n-1}\end{array}.$$

至此，我们找到了 $n$ 和 $n/2$ 的关系。然后，我们就可以考虑采用倍增的方法，采用类似快速幂的方式计算。

如果采用右推，即完全类似于快速幂的方法，以对指数 $ind$ 的右移 >>=1 为基本操作，每次操作之间并不是严格的二倍关系，需要记录大量的过程量，这对复杂度并不友好，所以我们考虑采用正推。

左推，即以对 $1$ 对左移操作 <<=1 为基本操作，直至左移到与指数相同为止。由于左移就是严格的二倍关系，只在当前位为 $1$ 时加上 $p^{n-1}$ 即可，而 $p^{n-1}$ 正是 $p^{2\lfloor n/2\rfloor }$，同样满足倍增关系，所以一切就变得非常好办了。

template <typename _Tp> _Tp fac_pow(_Tp base,int ind,int mod=INT_MAX) {
	_Tp ans=0,now=1; int st[100] {0};
	while(ind) {st[++st[0]]=ind&1,ind>>=1;}
	while(st[0]) {
		ans=(ans+(now*ans%mod))%mod,now=now*now%mod;
		if(st[st[0]]) ans=(ans+now)%mod,now=(now*base)%mod;
		st[0]--;
	}
	return ans;
}

其实快速幂也可以同理左推实现，同时，我们也可以通过寻找倍增关系快速求出 $\sum _{i=0}^n{p}^i$ 和 $\sum _{i=1}^n{p}^i$， 他们都是异曲同工的。

---

下午

其实第一天下午就讲了很多线性代数啊……

讲了素数筛法，矩阵等等。大部分还都算是见过的，比较好理解。

一起上课的很多同学年龄都比我们小很多，和我们听同样的课程，山东的竞赛确实比河北要强啊。

矩阵快速幂优化

矩阵运算都比较简单，需要注意的是，矩阵满足结合律但 不满足交换律。（abel 群笑话 hhhh）

矩阵乘法在 OI 中有很多用处：

1. 加速递推：数列的前几项用向量表示，使用矩阵快速幂加速计算过程，$O(n)\to O(\log n)$。

2. 加速一维相同的二维数组计算：比如下面这个问题：

   给定一个有向图，从 $A$ 点走 $k$ 步恰好走到 $B$ 点的方案数？

   这道题用 $E(x,y)\leftarrow 1$ 表示 $x$ 和 $y$ 之间有一条有向边，可得对于 $C=A\times A$ 有 $C(i,j)=\sum A(i,k)\times A(k,j)$，把 $k$ 作为中继点，不难看出其实这表达的就是从 $i$ 走两步恰好到 $j$ 的方案数。所以，矩阵快速幂，直接解决。

   可以看一把 P4159 [SCOI2009] 迷路 ，这道题用前缀优化建边和矩阵快速幂来解决。

数组遍历顺序优化

在计算机中，所有临时数据将优先存储在缓存中。遇到数组时，缓存中会将当前访问的数组下标的后几个数同时写入缓存，所以在遍历时，按顺序逐个遍历确实是最快的。

比如经典的矩阵乘法：

for(int i=1;i<=ans.n;++i)
	for(int k=1;k<=m;++k)
		for(int j=1;j<=ans.m;++j)
			ans.x[i][j]+=x[i][k]*a.x[k][j];

i,k,j 的遍历顺序就是速度最快的。

---

初等数论

4.30 - Day 2

上午

一觉醒来啊…7：05，爽。

逆元

众所周知，加减乘法满足 过程中多次取模 和 对结果取模 等效。我愿称此为 模运算的盖斯定律。然而，除法并不满足模运算的盖斯定律。我们可以计算 $12÷4mod5$，却对 $2÷4mod5$ 无计可施。所以，我们提出了 逆元 解决这一问题：

如果 $gcd(a,m)=1$ 且存在唯一的 $b$ 使得 $a\times b\equiv 1(\mathrm{mod}m)$ 且 $1\le b<m$，则 $b$ 为 $a$ 在 $modm$ 意义下的 逆元。

事实上，逆元是一个极广阔的概念，他是指一个可以取消另一给定元素运算的元素。在数学里，逆元素广义化了加法中的加法逆元和乘法中的倒数。模意义下的逆元又称为 数论倒数，这也很容易理解。

求解逆元，我们有如下几种方法：

1. 费马小定理

   对于任意一个 质数 $p$，若 $a,p$ 互质，有 $a^{p-1}\equiv 1(\mathrm{mod}p).$

   由费马小定理，我们可以推出 $a^{p-2}\equiv {\displaystyle \frac{1}{a}}(\mathrm{mod}p)$，快速幂求逆元即可。

2. 欧拉定理

   对于任意一个 自然数 $m$，若 $a,m$ 互质，有 $a^{\phi (m)}\equiv 1(\mathrm{mod}m).$ 其中，$\phi (m)$ 为欧拉函数。

   由于对于任意质数 $p$ 有 $\phi (p)=p-1$，费马小定理其实是欧拉定理在 $m$ 为质数下的 特殊情况。

   接下来，我们就可以得出 ${\displaystyle \frac{1}{a}}\equiv a^{\phi (m)-1}(\mathrm{mod}m)$，于是将 $÷a$ 转换为 $\times {\displaystyle \frac{1}{a}}$，最后转换成 $a^{\phi (m)-1}$，就得到了 $a$ 的逆元。

   Q: 怎么求欧拉函数？
   A: 你往下看。

3. 线性预处理序列逆元

   如果我们需要且仅需要求一组数的逆元，那么这还算是个不错的方法。

   我们可以线性处理前缀积 $fac$，（用这个名称是因为，当这一组数是从 $1$ 开始的公差为 $1$ 的等差数列时，前缀积就是 阶乘 了。）费马小 / 欧拉 / 扩欧求前缀积末项逆元 $ifa{c}_n$ ，向前逆推求出所有前缀积逆元 $ifac$ ，再用 $in{v}_i=fa{c}_{i-1}\cdot ifa{c}_i$ 的方法求出从 $1$ 到 $n$ 的所有逆元。本质是不断 约分。

   以从 $1$ 开始的公差为 $1$ 的等差数列为例，$Code$ 如下：

   void cal_inv() {
   	fac[0]=1;
   	for(int i=1;i<=n;++i) fac[i]=1ll*fac[i-1]*i%p
   		ifac[n]=ksm(fac[n],p-2,p);
   	for(int i=n;i>=0;--i) ifac[i]=1ll*ifac[i+1]*(i+1)%p;
   	for(int i=1;i<=n;++i) inv[i]=1ll*fac[i-1]*ifac[i]%p;
   }
   

   这个方法非常巧妙，我们可以具体的板演一下。

   对于需要求逆元的序列 $a_1,a_2,a_3,\cdots ,a_n$，设

   $$f_x=\prod _{i=1}^x{a}_i,$$

   再设

   $$g_x={\displaystyle \frac{1}{f_x}},$$

   则有

   $$g_i={\displaystyle \frac{1}{a_1{a}_2{a}_3\cdots a_i}}={\displaystyle \frac{1}{a_1{a}_2{a}_3\cdots a_i{a}_{i+1}}}\times a_{i+1}=g_{i+1}\times a_{i+1},$$

   故求 $g_n$ 逆元，后可一直逆推至 $g_1$ 逆元。后可知

   $${\displaystyle \frac{1}{a_i}}={\displaystyle \frac{1}{a_1{a}_2\cdots a_{i-1}{a}_i}}\times a_1{a}_2\cdots a_{i-1}=f_{i-1}\cdot g_i$$

   就求出逆元了。

   建议练习 洛谷P5431 【模板】乘法逆元 2。

4. 另一种线性预处理逆元方法（不常用）

   我们将 $pmodi$ 做如下变换：

   设 $k=\lfloor p÷i\rfloor ,r=pmodi,$

   得：

   $$p=ki+r$$
   $$0\equiv ki+r(\mathrm{mod}p)$$
   $$-r\equiv ki(\mathrm{mod}p)$$
   $${\displaystyle \frac{r}{i}}\equiv -k(\mathrm{mod}p)$$
   $${\displaystyle \frac{1}{i}}\equiv {\displaystyle \frac{p-k}{r}}(\mathrm{mod}p)$$

   逆元可求。

   $Code$：

   inv[1]=1;
   void cal_inv() {
   	for(int i=2;i<=n;++i) {
   		int k=p/i,r=p%i;
   		inv[i]=1ll*(p-k)*inv[r]%p;
   	}
   }
   

   UPD. 2023.5.5 关于预处理逆元，应用场景有限，特殊时候特殊方法对待。在 $p$ 为质数且不连续的时候，首选 费马小 + 快速幂，$p$ 不为质数……祝你好运。

5. 扩展欧几里得

   UPD. 2023.6.9 隆重介绍！应用场景更广泛且复杂度比肩费马小的 扩欧求逆元！

   根据逆元的定义，易得

   $$a\cdot x\equiv 1(\mathrm{mod}p),$$

   其中 $x$ 为 a 逆元。
   那么对于任意整数 $k,$ 就有

   $$k\cdot p+1\equiv 1(\mathrm{mod}p)$$
   $$k\cdot p+a\cdot x\equiv 1(\mathrm{mod}p),$$

   由于逆元定义写到 $gcd(a,p)=1,$ 那么我们就可以用下文中提到的 exgcd 求解方程

   $$p\cdot k+a\cdot x=1.$$

   $x$ 即为所求。

   这可真是太强了，即使 $p$ 不是素数仍然可以这么求解，所以对于 $p$ 不是素数的情况 认定他就对了。（稍微会比费马小慢，所以可以是质数的时候可以考虑费马小。）

欧拉函数

众所周知对于质数 $p$ 有 $\phi (p)=p-1$。

后可证：

$$\phi (p^2)=p(p-1)$$

$$\phi (p^n)={\displaystyle \frac{p-1}{p}}\times p^n$$

$$\phi (n)=n-{\displaystyle \frac{n}{p_1}}-{\displaystyle \frac{n}{p_2}}+{\displaystyle \frac{n}{p_1{p}_2}}$$

$$\phi (n)=n\times {\displaystyle \frac{p_1-1}{p_1}}\times \cdots \times {\displaystyle \frac{p_k-1}{p_k}}$$

其中第三式 $n$ 为质数 $p_1$ 与 $p_2$ 的积，利用第二式结合容斥原理推出；第四式在第三式基础上因式分解，然后合理外推，$n$ 表示任意自然数，$p_i$ 表示 $n$ 的质因数。

至此，我们可以 $O(\sqrt{n})$ 计算欧拉函数了。

int getphi(int x) {
	int ret=x,tmp=x;
	for(int i=2;i*i<=x;++i) {
		if(tmp%i==0) {
			ret=ret/i;
			ret=ret*(i-1);
		}
		while(!tmp%i) tmp=tmp/i;
	}
	if(tmp>1) ret=ret/tmp*(tmp-1);
	return ret;
}

Miller-Rabin 素性测试

如果 $n$ 是质数，则以下两条至少一条成立：

取 $a<n$，设 $n-1=d\times 2^r$，

$$a^d\equiv 1(\mathrm{mod}n)$$

$$\text{or}$$

$$\mathrm{\exists }0\le i<r,\text{s.t. }{a}^{d\times 2^i}\equiv 1(\mathrm{mod}n)$$

但是，有一些合数也是有记录通过 Miller-Rabin 测试的，所以我们可以多次测试，提高概率。

利用 Miller-Rabin，我们可以 $O(\log n)$ 判断质数。

$Code$：

bool miller_rabin(int n,int a) {
	int d=n-1,r=0;
	while(!d%2) d=d/2,r++;
	int x=ksm(a,d,n);
	if(x==1) return 1;
	for(int i=0;i<r;++i) {
		if(x==n-1) return 1;
		x=1ll*x*x%n;
	}
	return 0;
}
bool is_prime(int n) {
	if(n<2) return 0;
	for(int i=1;i<=23;++i) {
		int a=rand()%(n-1)+1;
		if(!miller_rabin(n,a)) return 0;
	}
	return 1;
}

或者，为提高成功几率，我们可以制定一个 $prime\mathrm{\_}list$，不用太大，几个数即可，然后用这些数进行检验，只要能保证他在 int 范围内正确就行。

int prime_list[]={2,3,5,7,13,23,37,73};
bool miller_rabin(int n,int a) {
	int d=n-1,r=0;
	while(!d%2) d=d/2,r++;
	int x=ksm(a,d,n);
	if(x==1) return 1;
	for(int i=0;i<r;++i) {
		if(x==n-1) return 1;
		x=1ll*x*x%n;
	}
	return 0;
}
bool is_prime(int n) {
	if(n<2) return 0;
	for(int i=0;i<8;++i) {
		if(n==prime_list[i]) return 1;
		if(!n%prime_list[i]) return 0;
		if(!miller_rabin(n,prime_list[i]%n)) return 0;
	}
	return 1;
}

扩展欧几里得 - exgcd

给定 $a,b$，已知 $gcd(a,b)=g$，求 $x,y,\text{ s.t. }$

$$xa+yb=g.$$

如何解决？

之前学到的辗转相除：

$$gcd(a,b)=gcd(b,amodb)=g$$

$$xa+yb=x^{\prime }b+y^{\prime }(amodb)$$

$$xa+yb=x^{\prime }b+y^{\prime }(a-b\cdot \lfloor {\displaystyle \frac{a}{b}}\rfloor )$$

$$xa+yb=x^{\prime }b+y^{\prime }(a-b\cdot \lfloor {\displaystyle \frac{a}{b}}\rfloor )$$

$$xa+yb=y^{\prime }a+(x^{\prime }-y^{\prime }\cdot \lfloor {\displaystyle \frac{a}{b}}\rfloor )b$$

$$x=y^{\prime },y=(x^{\prime }-y^{\prime }\cdot \lfloor {\displaystyle \frac{a}{b}}\rfloor )$$

问题解决。

$Code$：

int exgcd(int a,int b,int &x,int &y) {
	if(!b) {x=1,y=0; return a;}
	int g=exgcd(b,a%b,x,y),tmp=x;
	x=y,y=tmp-(a/b)*y;
	return g;
}

注意，这样递归的边界条件是 $x\cdot (a,b)+y\cdot 0=(a,b)$，所以在判断 !b 后 $x\leftarrow 1,y\leftarrow 0$，否则，$y$ 不影响结果（废话，因为是 $y\cdot 0$），$x$ 赋的值是几 $ax+by$ 就是几倍 $gcd(a,b)$ 。

UPD. 2023.7.29 谢邀，人在 ZROI。以前对 $x,y$ 我们都是手动 swap，其实可以直接动个小手脚：

int exgcd(int a,int b,int &x,int &y) {
	if(!b) { x=1,y=0; return a; 	}
	int g=exgcd(b,a%b,y,x);
	return y-=a/b*x,g;
}

一个意思。

裴蜀定理

害，就是你可以一步走 $x$ 米，也可以一步走 $y$ 米，问你离原点最近几米，就是 $gcd(x,y)$。

我要提出 碘酊定理 ！

$$\mathrm{\forall }x\in N^{\ast },1\equiv 1(\mathrm{mod}x)$$

伟大！伟大啊！！！

成功解决中国剩余定理。

下午

跟大家说一声，这个 $LATEX$ 不念雷泰克斯，

补：关于 exgcd

我们使用扩欧解决 P1082 [NOIP2012 提高组] 同余方程，找到的解不一定是最小的解。我们用找到的解除以他们的 $gcd$，寻找到一组 互质 的解 $x,y$，就可以找到当前方程的所有解了：

$$\{\begin{array}{l}{x}^{\prime }=x+kb\\ y^{\prime }=y+ka\end{array},k\in R$$

（扩展）中国剩余定理 - CRT & exCRT

解决

$$\{\begin{array}{l}x\equiv a_1(\mathrm{mod}{p}_1)\\ x\equiv a_2(\mathrm{mod}{p}_2)\\ \cdots \\ x\equiv a_k(\mathrm{mod}{p}_k)\end{array}$$

的问题。

不难得出，事实上同余方程组的解是一个同余方程。所以，其实中国剩余定理解决的是 同余方程的合并问题。

我们先以解

$$\{\begin{array}{l}x\equiv a_1(\mathrm{mod}{p}_1)\\ x\equiv a_2(\mathrm{mod}{p}_2)\end{array}$$

为例，介绍两种解法：

1. 大数翻倍法（暴力）

   《小学奥数题常用骗分方法》

   枚举第一个同余方程的解 $a_1+k\cdot p_1$，直至找到一个解满足第二个同余方程，或枚举至 大于 $\mathrm{lcm}(p_1,p_2)$，枚举至此意味着该方程组 无解。

   就是优化点的暴力呗。时间复杂度 $O(p_2)$。

   之所以叫 大数 翻倍法，就是枚举 $p$ 大的那个方程的解。这样到达无解的极端情况的时间更快。即 $O(min(p_1,p_2))$。

   有人说“卡常小妙招”？ZHX：事实上，他还没被成功卡过。整体来说，他的复杂度是 $\mathrm{\Theta }((\sum _{i=1}^{k}p)-\underset{i=1}{\overset{k}{max}}\{p_i\})$。之所以没被卡，是因为大部分时候都会出 $p$ 是质数的情况，而如果要保证 $\prod p$ 在 long long 范围内，那么其实最大的 $p$ 对于这个积来说还是蛮小的，不是很好卡。

2. 扩展中国剩余定理（正解）

   正解好闪，拜谢正解。

   将方程写成：

   $$x=k_1\cdot p_1+a_1=k_2\cdot p_2+a_2$$
   $$k_1\cdot p_1-k_2\cdot p_2=a_2-a_1$$

   扩欧！

   没错！扩展欧几里得求 $k_1,k_2$。

   ZHX 调整后代码：

   #include<bits/stdc++.h>
   #define int long long
   const int N=1e5+5;
   int n;
   int p[N],a[N];
   int exgcd(int a,int b,int &x,int &y) {
   	if(!b) {x=1,y=0; return a;}
   	int g=exgcd(b,a%b,x,y),tmp=x;
   	x=y,y=tmp-(a/b)*y;
   	return g;
   }
   void combine(int x,int y) {
   	int k1,k2;
   	int g=exgcd(p[x],p[y],k1,k2);
   	int k=(a[y]-a[x])/g; k1*=k;
   	int z=k1*p[x]+a[x];
   	p[y]=p[x]/g*p[y]; a[y]=(z%p[y]+p[y])%p[y];
   }
   int solve() {
   	for(int i=2;i<=n;++i)
   		combine(i-1,i);
   	return a[n];
   }
   signed main() {
   	std::ios::sync_with_stdio(0);
   	std::cin.tie(0);std::cout.tie(0);
   	std::cin>>n;
   	for(int i=1;i<=n;++i)
   		std::cin>>p[i]>>a[i];
   	int ans=solve();
   	std::cout<<ans;
   	return 0;
   }
   

   这种做法被称为 扩展中国剩余定理（exCRT）。扩展中国剩余定理不要求 $p$ 互质，nice。

   UPD. 2023.5.4 关于 ZHX 的代码，他 死 了。原因是过程量爆 long long 了。所以我们需要一个很新的乘法，就是类似于快速幂的方法，然后处处取模。

   #include<bits/stdc++.h>
   #define ll long long
   const int N=1e5+5;
   int n;
   ll p[N],a[N];
   template <typename T> void read(T &x)
   {
   	int f=1;x=0;char ch=getchar();
   	while(ch>'9'||ch<'0'){if(ch=='-')f=-1;ch=getchar();}
   	while(ch>='0'&&ch<='9'){x=(x<<3)+(x<<1)+(ch^48);ch=getchar();}
   	x*=f;
   }
   ll mul (ll a,ll b,ll mod) {
   	ll ans=0;
   	while(b) {
   		if(b&1) ans=(ans+a)%mod;
   		a=(a+a)%mod; b>>=1;
   	}
   	return ans;
   }
   ll exgcd(ll a,ll b,ll &x,ll &y) {
   	if(!b) {x=1,y=0; return a;}
   	ll g=exgcd(b,a%b,x,y),tmp=x;
   	x=y,y=tmp-(a/b)*y;
   	return g;
   }
   ll excrt() {
   	ll x,y,g,k,tmp,ans=a[1],P=p[1];
   	for(int i=2;i<=n;++i) {
   		ll B=p[i],C=(a[i]-ans%B+B)%B;
   		g=exgcd(P,B,x,y),tmp=B/g;
   		if(C%g) return -1;
   		x=mul(x,C/g,tmp),ans+=x*P,P*=tmp,ans=(ans%P+P)%P;
   	}
   	return (ans%P+P)%P;
   }
   int main() {
   	std::ios::sync_with_stdio(0);
   	std::cin.tie(0);std::cout.tie(0);
   	read(n);
   	for(int i=1;i<=n;++i)
   		read(p[i]),read(a[i]);
   	std::cout<<excrt();
   	return 0;
   }
   

质数筛法

嗯，TSOI 暑假入门内容。
先讲的埃氏筛：

for(int i=2;i<=n;++i)
	if(!not_prime[i])
		for(int b=a+a;b<=n;b+=a)
			 not_prime[b]=1;
for(int i=2;i<=n;++i)
	if(!not_prime[i]) prime[++cnt]=i;

根据调和级数，这个玩意复杂度 $O(n\log n)$，加第二行优化是 $O(\log (\log n))$，Vergil 讲了但我不会证。

然后是线性筛。

for(int i=2;i<=n;++i) {
	if(!not_prime[i]) cnt++,prime[++cnt]=i;
	for(int j=1;j<=cnt;++j) {
		int x=prime[j]*i;
			if(x>n) break;
		not_prime[x]=1;
		if(!i%prime[j]) break;//保证每个数只被最小的质因数筛掉
	}
}

就是把埃氏筛的内层拉出来放外层，然后加上关键性的加注释的那一行优化。

积性函数

$\mathrm{\forall }$ 互质的 $a,b$，若有 $f(a)\cdot f(b)=f(ab)$，则 $f(x)$ 是 积性函数。

$\mathrm{\forall }a,b$，若有 $f(a)\cdot f(b)=f(ab)$，则 $f(x)$ 是 完全积性函数。

易证欧拉函数是 积性函数，所以可以在线性筛的同时计算欧拉函数。

for(int i=2;i<=n;++i) {
	if(!not_prime[i]) cnt++,prime[++cnt]=i,phi[i]=i-1;
	for(int j=1;j<=cnt;++j) {
		int x=prime[j]*i;
		if(x>n) break;
		not_prime[x]=1;
		phi[x]=phi[prime[j]]*phi[i];//如果x与prime[j]互素
		if(!i%prime[j]) {
			phi[x]=prime[j]*phi[i];//如果x与prime[j]不互素
			break;
		}
	}
}

今天，ZHX说他 $\text{\#临时}$ 出了一道 $\text{\#凑数题}$，是$\text{\#T3}$，结果 $\text{\#得分率最低}$，是 $\text{\#斐波那契}$，靠。

5.1 - Day 3

上午

大步小步算法 - Baby Step Giant Step

AKA 北上广深 (BSGS)。

给定质数 $a,b,p$，求解方程 $a^x\equiv b(\mathrm{mod}p)$，保证 $a,b,p\le {10}^9$。

考虑使用暴力，枚举至已经出现过的模数，即枚举一个循环内是否有解。由于 $p$ 是素数，$a^{p-1}\equiv 1(\mathrm{mod}p)$，即一个循环的开始必然为 $1$，所以枚举至取模后为 $1$ 即可。

int solve(int a,int b,int p) {
	int v=1;
	for(int x=0;x<p-1;++x) {
		if(x==b) return x;
		v=1ll*v*a%p;
	}
	return -1;
}

srds，他复杂度 $O(p)$，算了。

所以，介绍 Baby Step Giant Step 算法。

将一个循环（$0\sim p-2$）分为一组 $s$ 个，如果

$$\mathrm{\exists }x\in [(j-1)s,js)\text{, s.t. }{a}^x\equiv b(\mathrm{mod}p),$$

那么令 $x=(j-1)s+k(k\in [0,s)\bigcap \mathbf{Z}),$ 就有

$$\begin{array}{rl}{a}^{(j-1)s+k}& \equiv b\\ a^k& \equiv b\cdot a^{-(j-1)s}.\end{array}$$

即如果第 $j$ 组有解可以满足该同余方程，则 $b\cdot a^{-(j-1)s}$ 必然在第 $1$ 组出现。最小非负整数解即为当前组的基数加上在第 $1$ 组出现的位置，即 $(j-1)s+k.$

int BSGS(int a,int b,ll p) {
	int s=sqrt(p),v=1;
	std::map<int,int> mp;
	for(int i=0;i<s;++i) mp[v]=i,v=(ll)v*a%p;
	for(int i=0;(ll)i*s<=p;++i) {
		int chk=(ll)b*ksm(ksm((ll)a,(ll)i*s,p),p-2,p)%p;
		if(mp.count(chk)) return mp[chk]+i*s;
	}
	return -1;
}

复杂度 $O(max({\displaystyle \frac{p}{s}},s))$，所以 $s=\sqrt{p}$ 时均摊复杂度最低。

$$\mathsf{\text{初等数论结课！}}$$

---

组合数学

好在主播有所了解。

排列数与组合数

排列数：在一个含 $n$ 个元素的集合中选取 $m$ 的元素，考虑顺序（同样的元素不同顺序视作不同集合），组成的子集的数量，记作 $P(n,m)$ 或 $A(n,m)$。

$$P(n,m)={\displaystyle \frac{n!}{(n-m)!}}$$

组合数：在一个含 $n$ 个元素的集合中选取 $m$ 的元素，不考虑顺序（同样的元素不同顺序视作相同集合），组成的子集的数量，记作 $C(n,m)$。

$$C(n,m)={\displaystyle \frac{n!}{m!(n-m)!}}$$

下文中对排列数和组合数用下标表示 $n,m$，即 $P_n^m$ 和 $C_n^m$。

一些有意思的结论：

$$C_n^0=C_n^n=1$$

$$C_n^m=C_n^{n-m}$$

当然，还有就是组合数的 递推式：

$$C_n^m=C_{n-1}^{m-1}+C_{n-1}^m$$

应该都不难理解吧。

根据递推式，我们不难发现 杨辉三角 就是 组合数下三角阵，也易得可以通过 动态规划 计算组合数。也可以用 预处理阶乘 和 预处理组合数 计算排列数，完美解决 除法不满足取模盖斯定律 的问题。

还有一些性质：

$$C_n^0+C_n^1+C_n^2+\cdots +C_n^n=2^n$$

对于 $n\ge 1$，有

$$C_n^1-C_n^2+C_n^3-\cdots C_n^n=0$$

证明如下：

$$C_n^1=C_{n-1}^0+C_{n-1}^1$$

$$C_n^2=C_{n-1}^1+C_{n-1}^2$$

$$\cdots$$

$$C_n^n=C_{n-1}^{n-1}$$

一直枚举，刚好可以抵消。

二项式定理

将 $(x+y{)}^n$ 展开可以发现，其各项系数刚好就是杨辉三角的第 $n+1$ 行。所以其实组合数也正好是 二项式系数。方便书写，我们也将 $C_n^m$ 写作 ${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}$，代表二项式系数。

于是有：

$$(x+y{)}^n=\sum _{i=0}^n{\displaystyle (\genfrac{}{}{0ex}{}{n}{i})}{x}^{n-i}{y}^i$$

同时，如果将 ${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}$ 展开 $k$ 次，也有：

$${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}=\sum _{i=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{i})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n-k}{m-k+i})}$$

令 $j=k-i$，得：

$$\begin{array}{rl}\sum _{i=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{i})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n-k}{m-k+i})}& =\sum _{i=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{k-j})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n-k}{m-j})}\\ & =\sum _{i=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n-k}{m-j})}\end{array}$$

使用换元法，将 $j$ 换成 $i$，得到最终结果：

$${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}=\sum _{i=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{i})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n-k}{m-k+i})}=\sum _{i=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{i})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n-k}{m-i})}$$

$6.$

组合数的 6 种求法

$ZHX$ 显出极高兴的样子，将两个指头的长指甲敲着讲台，点头说，“对呀对呀！……组合数有六样求法，你知道么？”

接下来我们一个一个说：

$\mathbf{\text{Question. }}{\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}modk$

1. 特殊的模数

$\mathbf{\text{Sol. }}\alpha -n,m\le {10}^{18},k=1$

$k=1$ 啊兄弟！就是 $0$！

2. 递推式

$\mathbf{\text{Sol. }}\beta -n,m\le {10}^3,k\le {10}^9$

使用杨辉三角，$O(n^2)$ 的复杂度。

for(int i=0;i<=n;++i) {
	C[i][0]=1;
	for(int j=1;j<=i;++j)
		C[i][j]=(C[i-1][j-1]+C[i-1][j])%p;
}

3. 预处理阶乘逆元

$\mathbf{\text{Sol. }}\gamma -n,m\le {10}^6,k\le {10}^9\text{ and }k\text{ is prime.}$

预处理阶乘，然后除法使用逆元（$k$ 是质数，费马小），定义式求解组合数。

fac[0]=1;
for(int i=1;i<=1000000;++i)
	fac[i]=(ll)fac[i-1]*i%p;
C(n,m) =(ll)fac[n]*ksm(fac[m],p-2,p)%p*ksm(fac[n-m],p-2,p)%p;

4. 定义式约分

$\mathbf{\text{Sol. }}\delta -n\le {10}^9,m\le {10}^3,k\le {10}^9$

ZHX OI 社会学定律：没有思路看数据范围，一切尽在数据范围最诡异的数中。

本范围的 1e3 就够诡异。于是，我们尝试利用一下这个 1e3。

我们可以将组合数约分，消去 $(n-m)!$ ，就有

$$\begin{array}{rl}{\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}& ={\displaystyle \frac{n!}{m!(n-m)!}}\\ & ={\displaystyle \frac{n\cdot (n-1)\cdot (n-2)\cdot ...\cdot (n-m+1)}{1\cdot 2\cdot ...\cdot m}},\end{array}$$

这串式子的上下都只有 $m$ 项，复杂度是可以压住的；考虑组合数必定是整数，我们就可以每次寻找分子与分母的 $gcd$ 进行约分，暴力但有效地解决这一问题。

for(int i=1;i<=m;i++)
	deno[i]=i,nume[i]=n-i+1;
for(int i=1;i<=m;++i)
	for (int j=1;j<=m;++j) {
		int g=gcd(nume[i],deno[j]);
		nume[i]/=g;
		deno[j]/=g;
	}
int ans=1;
for(int i=1;i<=m;++i)
ans=(ll)ans*nume[i]%p;

5. Lucas 定理

$\mathbf{\text{Sol. }}ϵ-n,m\le {10}^9,k\le 100\text{ and k is prime.}$

隆重介绍 卢卡斯定理：

对于质数 $p$，有：

$${\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}modp={\displaystyle (\genfrac{}{}{0ex}{}{\lfloor n/p\rfloor }{\lfloor m/p\rfloor })}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{nmodp}{mmodp})}modp$$

所以，我们可以将 $n,m$ 转换为 $p$ 进制，然后对于 $n_{(p)}$ 的每位 $a_i$ 和 $m_{(p)}$ 的对应为 $b_i$，计算

$$\prod {\displaystyle (\genfrac{}{}{0ex}{}{a_i}{b_i})}modp$$

即可。

int solve() {
	while(n) x[++x[0]]=n%p,n/=p;
	while(m) y[++y[0]]=m%p,m/=p;
	int ans=1;
	for(int i=1;i<=std::max(x[0],y[0]);++i) {
		if(i>x[0]) x[i]=0;
		if(i>y[0]) y[i]=0;
		if(x[i]<y[i]) { ans=0; break; }
		ans=(ll)ans*fac[x[i]]%p*ksm((ll)fac[y[i]]*fac[x[i]-y[i]]%p,p-2,p)%p;
	}
	return ans;
}

6. 扩展卢卡斯定理 - exLucas

$\mathbf{\text{Sol. }}\zeta -n,m\le {10}^9,k\le 100$

$k$ 不是质数，怎么办呢？

把 $k$ 拆分成质因数 $k=p_1^{a_1}\cdot p_2^{a_2}\cdot \cdots \cdot p_k^{a_k}$，然后就会有同余方程组：

$$\{\begin{array}{l}{\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}mod{p}_1^{a_1}=x_1\\ {\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}mod{p}_2^{a_2}=x_2\\ \cdots \\ {\displaystyle (\genfrac{}{}{0ex}{}{n}{m})}mod{p}_k^{a_k}=x_k\end{array}$$

解 中国剩余定理 即可。

具体过程详见 $dalao$ 题解。

$\mathbf{\text{AC Code:}}$

#include<bits/stdc++.h>
#define ll long long
#define ld long double
#define fsp(x) std::fixed<<std::setprecision(x)
#define forE(u) for(ll p=head[u],v=E[p].to;p;p=E[p].next,v=E[p].to)
const int N=1e5+5;
ll cnt,P[N],A[N];
template <typename _Tp> _Tp ksm(_Tp base,ll ind,ll mod) {
	_Tp ans=1;
	while(ind) {
		if(ind&1) ans=ans*base%mod;
		base=base*base%mod; ind>>=1;
	}
	return ans%mod;
}
ll exgcd(ll a,ll b,ll &x,ll &y) {
	if(!b) return x=1,y=0,a;
	ll g=exgcd(b,a%b,x,y),tmp=x;
	x=y,y=tmp-(a/b)*y;
	return g;
}
inline ll inv(ll a,ll p) { ll x,y; exgcd(a,p,x,y); return (x+p)%p; }
ll f_pk(ll n,ll p,ll pk) {
	if(!n) return 1;
	ll ans=1;
	for(ll i=1;i<=pk;++i)
		if(i%p) ans=ans*i%pk;
	ans=ksm(ans,n/pk,pk);
	for(ll i=pk*(n/pk);i<=n;++i)
		if(i%p) ans=ans*(i%pk)%pk;
	return f_pk(n/p,p,pk)*ans%pk;
}
ll g_p(ll n,ll p) { return n<p?0:n/p+g_p(n/p,p); }
inline ll C_pk(ll n,ll m,ll p,ll pk) { return f_pk(n,p,pk)*inv(f_pk(m,p,pk)*f_pk(n-m,p,pk),pk)%pk*ksm(p,g_p(n,p)-g_p(m,p)-g_p(n-m,p),pk)%pk; }
inline ll excrt() {
	ll x,y,g,tmp,ans=A[1],p=P[1];
	for(ll i=2;i<=cnt;++i) {
		ll B=P[i],C=(A[i]-ans%B+B)%B;
		g=exgcd(p,B,x,y),tmp=B/g;
		x=C/g*x%tmp,ans+=x*p,p*=tmp,ans=(ans%p+p)%p;
	}
	return (ans%p+p)%p;
}
inline ll exLucas(ll n,ll m,ll p) {
	for(ll i=2;i*i<=p;++i) {
		ll tmp=1;
		while(!(p%i)) tmp*=i,p/=i;
		if(tmp==1) continue;
		P[++cnt]=tmp,A[cnt]=C_pk(n,m,i,tmp);
	}
	if(p!=1) P[++cnt]=p,A[cnt]=C_pk(n,m,p,p);
	return excrt();
}
int main() {
	std::ios::sync_with_stdio(false);
	std::cin.tie(nullptr); std::cout.tie(nullptr);
	ll n,m,p;
	std::cin>>n>>m>>p;
	ll ans=exLucas(n,m,p);
	std::cout<<ans<<'\n';
	return 0;
}

tips：关于非常大的 ${\displaystyle (\genfrac{}{}{0ex}{}{n_1}{m_1})},{\displaystyle (\genfrac{}{}{0ex}{}{n_2}{m_2})}$ 比大小，可以转换成 $\log {\displaystyle (\genfrac{}{}{0ex}{}{n_1}{m_1})},\log {\displaystyle (\genfrac{}{}{0ex}{}{n_2}{m_2})}$ 进行比较，开 double 就行。

求和变形

遇到多层求和 $\sum \sum$ 的时候，复杂度压不住啊。

于是有了求和变形。

1. 增加枚举量

2. 交换枚举顺序

3. 分离无关变量

4. 换元法

比如 P3746 [六省联考 2017] 组合数问题，他想求：

$$\left(\sum _{i=0}^{\mathrm{\infty }}{C}_{nk}^{ik+r}\right)modp,$$

我就可以

$$\begin{array}{rl}& \sum _{i=0}^{\mathrm{\infty }}{C}_{nk}^{ik+r}\\ & =\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{nk}{ik+r})}\\ & =\sum _{i=0}^{\mathrm{\infty }}\sum _{j=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{nk-k}{ik+r-j})}& \mathit{1}\\ & =\sum _{j=0}^k\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{nk-k}{ik+r-j})}& \mathit{2}\\ & =\sum _{j=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{(n-1)k}{ik+r-j})}& \mathit{3}& ,\end{array}$$

其中 $\mathit{1}$ 式使用了增加枚举量，$\mathit{2}$ 式使用了交换枚举顺序，$\mathit{3}$ 式使用了分离无关变量。那么，到了这里如何解呢？

我们使用第 4 招——换元法，设

$$f(n,r)=\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{nk}{ik+r})},$$

就有

$$f(n-1,r-j)=\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{(n-1)k}{ik+r-j})},$$

接下来，添加辅助纬度，将原式表达为

$$f_n(1,r)=\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{nk}{ik+r})},$$

再使用多项式定理，就有

$$\begin{array}{rl}{f}_n(1,r)& =\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{nk}{ik+r})}\\ & =\sum _{i=0}^{\mathrm{\infty }}\sum _{j=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}{\displaystyle (\genfrac{}{}{0ex}{}{(n-1)k}{ik+r-j})}\\ & =\sum _{j=0}^k\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{(n-1)k}{ik+r-j})}{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}\\ & =\sum _{j=0}^k{f}_{n-1}(1,r-j){\displaystyle (\genfrac{}{}{0ex}{}{k}{j})},\end{array}$$

怎么样，是不是有点意思？

到了这儿，其实复杂度并没有降，所以是时候想个方法了——矩阵快速幂。但是现在还不满足条件，所以我们开一个新的矩阵 $D$，并使

$$D(r-j,r)={\displaystyle (\genfrac{}{}{0ex}{}{k}{j})},$$

于是就有了

$$f_n(1,r)=\sum _{j=0}^k{f}_{n-1}(1,r-j)\cdot D(r-j,r),$$

完全满足矩阵乘法。

所以，我们就可以使用

$$f_n=f_{n-1}\times D$$

结合矩阵快速幂完美解决这道题，也就是说我们使用矩阵快速幂，然后

$$\begin{array}{rl}\sum _{i=0}^{\mathrm{\infty }}{C}_{nk}^{ik+r}& =\sum _{j=0}^k{\displaystyle (\genfrac{}{}{0ex}{}{k}{j})}\sum _{i=0}^{\mathrm{\infty }}{\displaystyle (\genfrac{}{}{0ex}{}{(n-1)k}{ik+r-j})}\\ & =f_n(1,r).\end{array}$$

泰裤辣！

抽屉原理

有 $n+1$ 个物品和 $n$ 个抽屉，把这些物品放进抽屉，至少有 一个抽屉 里面有 多于一个 物品。

没了。

还在？

真没了。

就把其他的题抽象成这个模型就行了，没了。

容斥原理

最基本的，

$$|A_1\cup A_2|=|A_1|+|A_2|-|A_1\cap A_2|,$$

$$\begin{array}{rl}|A_1\cup A_2\cup A_3|& =|A_1|+|A_2|+|A_3|\\ & -|A_1\cap A_2|-|A_1\cap A_3|-|A_2\cap A_3|\\ & +|A_1\cap A_2\cap A_3|,\end{array}$$

反正就是，一层加，一层减，然后就有一个你看不懂但是可以意会其精神的式子：

$$\left|\bigcup _{i=1}^n{S}_i\right|=\sum _{m=1}^n(-1{)}^{m-1}\sum _{a_i<a_{i+1}}\left|\bigcap _{i=1}^m{S}_{a_i}\right|,$$

理论就结束了。

全靠思考，抽象模型。

在容斥原理的运用中，我们可以将较强的条件（苛刻的）转换为较弱的条件（如不满足其中一个条件），然后使用容斥原理求解。

5.2 - Day 4

上午

容斥定理解决春季赛 T2

没错！其实春季赛 T2 也是一道 容斥原理 的题。

[春季测试 2023] 幂次

题目描述

小 Ω 在小学数学课上学到了“幂次”的概念：$\mathrm{\forall }a,b\in {\mathbb{N}}^{+}$，定义 $a^b$ 为 $b$ 个 $a$ 相乘。

她很好奇有多少正整数可以被表示为上述 $a^b$ 的形式？由于所有正整数 $m\in N^{+}$ 总是可以被表示为 $m^1$ 的形式，因此她要求上述的表示中，必须有 $b\ge k$，其中 $k$ 是她事先选取好的一个正整数。

因此她想知道在 $1$ 到 $n$ 中，有多少正整数 $x$ 可以被表示为 $x=a^b$ 的形式，其中 $a,b$ 都是正整数，且 $b\ge k$？

输入格式

第一行包含两个正整数 $n,k$，意义如上所述。

输出格式

输出一行包含一个非负整数表示对应的答案。

考虑使用枚举法，依次枚举 $b=k,k+1,k+2\cdots$ 时所有的数，可以推出在 $n$ 以内这样的数有 $\lfloor \sqrt[b]{n}\rfloor$ 个。当然，其中会有重复计数。比如，$a^6$ 就会被 $a^2$ 和 $a^3$ 同时计数。所以，我们需要扣去一次。

那么，如何确定当前幂次是应该加上还是减去呢？理论上，这需要用到 莫比乌斯反演。但是显然这对目前阶段的我们来说可能还是有点超纲了。所以我们采用一种更加暴力的方法。

使用 $nu{m}^b$ 表示 $a^b$ 被计数了几次，根据 $num$ 决定当前幂次是加是减。核心部分如下：

for(int a=2;a<=64;a++)
	num[a] = 0;
for(int a=2;a<=64;a++) {
//	long long v=pow(n,1.0/a)-1; 可能导致精度问题
	long long v=(long long)floor(exp(log(n)/a))-1; //精度优化
	int d=1-num[a];
	ans+=v*d;
	for(int b=a;b<=64;b+=a)
		num[b]+=d;
}

冷知识:

$$2^{\frac{{\log }_{2}n}{a}}=\sqrt[a]{n}$$

可以解决精度问题。

---

线性代数

高斯消元

现在我们要解决

$$\{\begin{array}{l}{a}_1^1{x}_1+a_1^2{x}_2+\cdots +a_1^n{x}_n=b_1\\ a_2^1{x}_1+a_2^2{x}_2+\cdots +a_2^n{x}_n=b_2\\ \cdots \\ a_n^1{x}_1+a_n^2{x}_2+\cdots +a_n^n{x}_n=b_n\end{array}$$

这么一个 $n$ 元一次方程组，how？

高斯消元。

其实就是把咱们手解方程的过程转换成代码了。

void gauss() {
	for(int i=1;i<=n;++i) {
		for(int j=i;j<=n;++j)
			if(a[j][i]!=0) {
				for(int k=1;k<=n+1;k++) swap(a[i][k],a[j][k]);
				break;
			}
		for(int j=i+1;j<=n;++j) {
			if(a[j][i]==0) continue;
			int l=a[i][i]/gcd(abs(a[i][i]),abs(a[j][i]))*a[j][i];
			int ratioi=l/a[i][i];
			int ratioj=l/a[j][i];
			for (int k=1;k<=n+1;++k)
				a[j][k] = a[j][k] * ratioj - a[i][k] * ratioi;
		}	
	} 
	for(int i=n;i>=1;--i) {
		for(int j=i+1;j<=n;++j)
			a[i][n+1]-=a[i][j]*x[j];
		x[i]=a[i][n+1]/a[i][i];
	}
}

还有一个优化一点：

bool gauss() {
	for(int i=1;i<=n;++i) {
		int r=i;
		for(int j=i;j<=n;++j)
			if(fabs(a[j][i])>fabs(a[r][i]))
				r=j;
		if(fabs(a[r][i])<eps) return 0;
		if(i!=r) std::swap(a[i],a[r]);
		for(int j=i+1;j<=n;++j) {
			double ratio=a[j][i]/a[i][i];
			for(int k=1;k<=n+1;++k)
				a[j][k]-=a[i][k]*ratio;
		}
	}
	for(int i=n;i>=1;--i) {
		for(int j=i+1;j<=n;++j)
			a[i][n+1]-=a[i][j]*x[j];
		x[i]=a[i][n+1]/a[i][i];
	}
	return 1;
}

最后是升级版的 高斯 - 约旦消元法：

bool jordan() {
	for(int i=1;i<=n;++i) {
		int r=i;
		for(int j=i;j<=n;++j)
			if(fabs(a[j][i])>fabs(a[r][i])) r=j;
		if(fabs(a[r][i])<eps) return 0;
		if(i!=r) std::swap(a[i],a[r]);
		for(int j=1;j<=n;++j) {
			if(i!=j) {
				double ratio=a[j][i]/a[i][i];
				for(int k=i+1;k<=n+1;++k)
					a[j][k]-=a[i][k]*ratio;
			}
		}
	}
	for(int i=1;i<=n;++i) x[i]=a[i][n+1]/a[i][i];
	return 1;
}

年轻人！珍惜约旦优美的"$\setminus$"型吧。

单位矩阵

就是这样一个矩阵：

$$\left[\begin{array}{ccc}1& 0& \cdots \\ 0& 1& \ddots \\ ⋮& \ddots & \ddots \end{array}\right]$$

有这么一个性质：

$$\left[\begin{array}{ccc}x& 0& \cdots \\ 0& x& \ddots \\ ⋮& \ddots & \ddots \end{array}\right]\times A=xA$$

然后就可以把数乘转化成矩阵乘法了，可以运用这么一个构造函数对矩阵赋 int 型的初值。

matrix(int z=0,int a=0,int b=0) {
	n=a,m=b,memset(x,0,sizeof x);
	for(int i=0;i<=24;++i)
	x[i][i]=z;
}

那个 24 是我开的矩阵的数组的长度，不同情况不同对待。

单位矩阵只能是 $n\times n$ 的，即 行数与列数相同。

初等矩阵

一个奇妙的矩阵：

$$\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]$$

你会发现他可以这样：

$$\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\times \left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]=\left[\begin{array}{cc}3& 4\\ 1& 2\end{array}\right]$$

$$\left[\begin{array}{cc}1& 2\\ 3& 4\end{array}\right]\times \left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}2& 1\\ 4& 3\end{array}\right]$$

这个矩阵是将单位矩阵的第 $1$ 行和第 $2$ 行交换，比如还有这个矩阵：

$$\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\end{array}\right]$$

他相比单位矩阵，第 $2$ 行和第 $3$ 行做了交换。他乘矩阵 $X$ 后，结果就是 $X$ 第 $2$ 行和第 $3$ 行 交换。

而如果单位矩阵的第 $i$ 行第 $j$ 列增加 $x$，将该矩阵乘矩阵 $X$，就会将 $X$ 的第 $j$ 行的 $x$ 倍加到第 $i$ 行。

这些矩阵都是 初等矩阵，他们是由单位矩阵经过 初等变换 得来的。

初等变换共有三种：

1. 交换两行或两列；

2. 用一个数 $k$ 乘某一行；

3. 用一个数 $k$ 乘某一行后加到另一行去（乘行不变）。

哎，不行就看看 这篇 博吧，反正对于 OIer 来说这些应该已经够了。

逆矩阵

对一个 $n$ 阶矩阵 $A$ ，如果存在另一个 $n$ 阶矩阵 $B$，它们满足：$AB=BA=E$（$E$ 为单位矩阵），那么 $A$ 和 $B$ 互为逆矩阵。

如果你把

$$\{\begin{array}{l}{a}_1^1{x}_1+a_1^2{x}_2+...+a_1^n{x}_n=b_1\\ a_2^1{x}_1+a_2^2{x}_2+...+a_2^n{x}_n=b_2\\ ...\\ a_n^1{x}_1+a_n^2{x}_2+...+a_n^n{x}_n=b_n\end{array}$$

写成

$$\left[\begin{array}{cccc}{a}_1^1& a_1^2& \cdots & a_1^n\\ a_2^1& a_2^2& \cdots & a_2^n\\ ⋮& ⋮& ⋮& ⋮\\ a_n^1& a_n^2& \cdots & a_n^n\end{array}\right]\times \left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right],$$

然后就不难发现，解方程组其实可以看作求逆矩阵的过程。

同时，也有一些题会问到你，让你去求逆矩阵，可以把原矩阵 $A$ 和单位矩阵 $E$ 放在一起进行高斯消元，就会有

$$A\times B=E$$

$$E\times B=Y$$

然后这个矩阵 $Y$ 就是你要求的逆矩阵了。

---

初等概率

下午

嗨害，直接就是概率啊。

$$P(A)\cdot P(B)=P(AB)$$

• 随机试验：

  1. 不能预先确知结果；

  2. 试验之前可以预测所有可能结果或范围；

  3. 可以在相同条件下重复实验。

• 样本空间：随机试验所有可能结果组成的集合。

• 事件：样本空间的任意一个子集。

• 事件发生：在一次试验中，事件的一个样本点发生。

• 必然事件：样本空间全集。

• 不可能事件：空集。

• 条件概率：$B$ 发生时，$A$ 发生的概率。记作 $P(A|B).$

  $$P(A|B)={\displaystyle \frac{P(AB)}{P(B)}}$$
  $$P(A|B)P(B)=P(AB)$$

• 独立事件：$A$ 是否发生与 $B$ 无关。

• 期望：我懂了，但我不好说。大概是是每种结果贡献与概率的乘积的和。

  期望有一个非常重要的性质，就是

  $$E[x_1+x_2]=E[x_1]+E[x_2],$$

  即 期望的和等于和的期望。

最后的最后，我们来一道题，实战一下：

小胡站在原点，手里拿着两枚硬币。抛第一枚硬币正面向上的概率为 $p$，第二枚正面向上的概率为 $q$。
小胡开始抛第一枚硬币，每次抛到反面小胡就向 $x$ 轴正方向走一步，直到抛到正面。
接下来小胡继续抛第一枚硬币，每次抛到反面小胡就向 $y$ 轴正方向走一步，直到抛到正面。
现在小胡想回来了，于是他开始抛第二枚硬币，如果小胡抛到正面小胡就向 $x$ 轴的负方向走一步，否则小胡就向 $y$ 轴的负方向走一步。
现在小胡想知道他在往回走的时候经过原点的概率是多少呢？

OK，我们可以把计算概率的式子列出来：

$$\sum _{x=0}^{\mathrm{\infty }}\sum _{y=0}^{\mathrm{\infty }}(1-p{)}^{x}p\cdot (1-p{)}^{y}p\cdot q^x(1-q{)}^y{\displaystyle (\genfrac{}{}{0ex}{}{x+y}{x})}$$

因为有可能走到无限远，所以求和从 $0$ 到 $\mathrm{\infty }.$ $(1-p{)}^{x}p$ 是走到 $(x,0)$ 的概率，$(1-p{)}^{y}p$ 是在此基础上走到 $(x,y)$ 的概率，$q^x(1-q{)}^y{\displaystyle (\genfrac{}{}{0ex}{}{x+y}{x})}$ 是再走回 $(0,0)$ 的概率。接下来，进行史诗级变形：

$$\begin{array}{rl}& \sum _{x=0}^{\mathrm{\infty }}\sum _{y=0}^{\mathrm{\infty }}(1-p{)}^{x}p\cdot (1-p{)}^{y}p\cdot q^x(1-q{)}^y{\displaystyle (\genfrac{}{}{0ex}{}{x+y}{x})}\\ =& \sum _{x=0}^{\mathrm{\infty }}\sum _{y=0}^{\mathrm{\infty }}{p}^2(1-p{)}^{x+y}\cdot q^x(1-q{)}^y{\displaystyle (\genfrac{}{}{0ex}{}{x+y}{x})}\\ =& \sum _{t=0}^{\mathrm{\infty }}\sum _{x=0}^t{p}^2(1-p{)}^t\cdot q^x(1-q{)}^{t-x}{\displaystyle (\genfrac{}{}{0ex}{}{t}{x})}\\ =& p^2\sum _{t=0}^{\mathrm{\infty }}(1-p{)}^t\sum _{x=0}^t\cdot q^x(1-q{)}^{t-x}{\displaystyle (\genfrac{}{}{0ex}{}{t}{x})}\\ =& p^2\sum _{t=0}^{\mathrm{\infty }}(1-p{)}^t(q+1-q{)}^t\\ =& p^2\sum _{t=0}^{\mathrm{\infty }}(1-p{)}^t\\ =& p^2\cdot {\displaystyle \frac{1-(1-p{)}^{\mathrm{\infty }+1}}{1-(1-p)}}\\ =& p^2\cdot {\displaystyle \frac{1}{p}}\\ =& p.\end{array}$$

我知道你现在的心情，一方面是像疯了一样，震惊于这个变形绝妙的美感；一方面是像傻了一样，尝试用毕生所学去理解这个式子。接下来，我们一步一步来说。

第一步，合并同类项，即 $p^2$ 和 $(1-p{)}^{x+y}$；

第二步，换元法，设 $t=x+y$，更换枚举变量；

第三步，分离无关变量，改变了 $p^2$ 和 $(1-p{)}^{x+y}$ 的位置；

第四步，运用 二项式定理，换 $\sum _{x=0}^t\cdot q^x(1-q{)}^{t-x}{\displaystyle (\genfrac{}{}{0ex}{}{t}{x})}$ 为 $(q+1-q{)}^t$；

$$\mathbf{\text{二项式定理：}}(a+b{)}^n=\sum _{x=0}^n{a}^x\cdot b^{n-x}\cdot {\displaystyle (\genfrac{}{}{0ex}{}{n}{x})}.$$

第五步，化简该式；

第六步，运用 等差数列求和公式，化 $\sum _{t=0}^{\mathrm{\infty }}(1-p{)}^t$ 为 ${\displaystyle \frac{1-(1-p{)}^{\mathrm{\infty }+1}}{1-(1-p)}}$;

$$\begin{array}{rl}\mathbf{\text{等差数列求和公}}& \mathbf{\text{式推导：}}\\ x& =1+a+a^2+\cdots +a^n\\ ax& =a^1+a^2+a^3+\cdots +a^{n+1}\\ (1-a)x& =1-a^{n+1}\\ x& ={\displaystyle \frac{1-a^{n+1}}{1-a}}.\end{array}$$

第七步，运用极限思想，认为 $(1-p{)}^{\mathrm{\infty }+1}$ 无限趋近于 $0$，将 ${\displaystyle \frac{1-(1-p{)}^{\mathrm{\infty }+1}}{1-(1-p)}}$ 化简为 ${\displaystyle \frac{1}{p}}$;

最后一步，化简得出结果。

没想到吧！如此繁杂的式子，变形化简之后竟只剩一个 $p$！

---

结语

五一数学专题集训就到此为止了，但是探索的旅途永远不会结束。不论是在数学上，还是在信竞中，人类正以昂扬的姿态谱写着新的历史——而 $TSOI$ 的新的历史，将由我们书写。

王重阳，2023.5.4 署

$$\mathbf{\text{- The End. -}}$$