简单复习一下数学的几个模板。
总感觉数学还是要多复习，不然很容易忘记。
简单数论
(O(sqrt{n}))判定质数
不用多讲。
Miller-Rabin算法
bool mr(ll x,ll b)
{
if(qp(b,x-1,x)!=1)
return false;
ll k=x-1;
while((k&1)==0)
{
k>>=1;
ll d=qp(b,k,x);//快速幂,计算b^k mod x
if(d!=1 && d!=x-1)
return false;
if(d==x-1)
return true;
}
return true;
}
bool mr(ll x)
{
if(x==46856248255981)
return false;
if(x==2||x==3||x==7||x==61||x==24251)
return true;
return mr(x,2)&&mr(x,3)&&mr(x,7)&&mr(x,61)&&mr(x,24251);
}
对于一个数(x)，我们选一个质数(b)作为基底，依次计算(b^{x-1},b^{frac{x-1}{2}},b^{frac{x-1}{2^2}},cdots)。注意：
- 在第一次计算时不考虑(b^{x-1} equiv x-1)的情况。
- 遇到(b^{frac{x-1}{2^k}} equiv x - 1)或者(frac{x-1}{2^k})为奇数的情况就直接判定正确。
注意(5)个素因子积：(2,3,5,7,24251)。那个很强的伪素数(46856248255981)我倒觉得真心没必要背。
约数(O(sqrt{N}))筛法
又叫试除法。每次围绕(sqrt{N})对称地找到两个因子(p,q)，把这两个因子加入约数列表中。注意当(p^2=N)时要特判一下，不然会加多了。
倍数法
和线性筛的代码比较类似，适合一次筛多个质数的问题。
gcd
埃式筛
比较好理解。
线性筛
for(rg int i = 1; i <= N; ++ i)
{
if(v[i] < 0)
{
v[i] = i;
prime[++ pn] = i;
}
for(rg int j = 1; j <= pn; ++ j)
{
if(prime[j] * i > N || prime[j] > v[i])
break;
v[i * prime[j]] = prime[j];
}
}
暂时还没有试过这个代码。从理论上来说应该是没错的。
进阶数论
冷知识
设一个数的分解式：
我们可以把指数提取出来：
它可以帮助我们从全新的角度观察数论。其实说白了就是把数论中的运算变成对数形式。
举个例子，求两个数的(gcd):
这样我们就可以把方程(gcd(x,a)=b)变成若干个整数不等式了。
再举个例子，整除关系：
这样就也把整除关系变成若干个不等式了。从本质上来说，(<mathbb{Z}^{*},|>)本身就是一个偏序集，因此才会有这样的特性。
可以结合这道题看一下：CSP-S 2009 Hankson的趣味题
重要函数的计算与预处理
欧拉函数(phi(n))。
int phi(int n)
{
int ret = n;
for(rg int i = 2; i * i <= n; ++ i)
if(n % i == 0)
{
ret = ret / i * (i - 1);
while(n % i == 0)
n /= i;
}
if(n > 1)
ret = ret / n * (n - 1);//n是质数
return ret;
}
这个是直接根据定义式计算的。
它的预处理版本：
void init()
{
for(rg int i = 2; i <= n; ++ i)
phi[i] = i;
for(rg int i = 2; i <= n; ++ i)
if(phi[i] == i)
for(int j = i; j <= n; j += i)
phi[j] = phi[j] / i * (i - 1);
}
当然，还有基于线性筛的版本：
void init()
{
for(rg int i = 2; i <= N; ++ i)
{
if(v[i] == 0)
{v[i] = i, prime[++ pn] = i, phi[i] = i - 1;}
for(rg int j = 1; j <= pn; ++ j)
{
if(prime[j] * i > N || prime[j] > v[i])
break;
v[i * prime[j]] = prime[j];
phi[i * prime[j]] = phi[i] * (i % prime[j] == 0 ? prime[j] : prime[j] - 1);
}
}
}
这里用到了(phi(i))的两个性质，这里也不多赘述。
莫比乌斯函数(mu(n))。
for(rg int i = 1; i <= N; ++ i)
mu[i] = 1, v[i] = 0;
for(rg int i = 2; i <= N; ++ i)
{
if(v[i])
continue;
mu[i] = -1;
for(rg int j = (i << 1); j <= N; j += i)
{
v[j] = 1;
if((j / i) % i == 0)
mu[j] = 0;
else
mu[j] *= -1;
}
}
顺带一提，学长zsy讲过一个非常有趣的理解方式。这里不妨粘贴之。
我们都知道莫比乌斯反演：（简便起见，这里直接写成卷积形式）
从之前的素因子式(N=(alpha_i))的角度理解，((alpha_i))其实是一个高维向量，而(g(alpha_i))就是这个向量的点权，(f(alpha_i))就是(g)的前缀和。联想到二维前缀和(a_i = S_i - S_{i-1})，莫比乌斯反演事实上就是对于高维前缀和(f(alpha_i))的差分变换。
实际上，我之前补充过一个冷知识：广义的莫比乌斯反演。对于偏序关系(preceq)，定义它的莫比乌斯函数：
则有：
这实际上说明了莫比乌斯反演实质上就是一个局部有限偏序集上的高维前缀和差分变换。这也进一步说明了莫比乌斯函数在容斥原理和数论上的链接作用。
模算术
正如你所见，(a mod b)这个函数在设计时出现了一些问题，使得两个字母之间间隔太大。因此，下文中均用C++写法(a\% b)代表(amod b)。同理，为了方便书写，我们统一用(a/b)代表(lfloorfrac{a}{b} floor)。
欧拉定理
若正整数(a,p)互质，则：
对于任意正整数(b)，还有：
如果(a)和(p)不一定互质，则应该使用下面这个公式：
证明均略去。
裴蜀定理和exgcd
对于整数方程(ax+by=c)，存在整数解((x,y))的充要条件是(gcd(a,b)|c)。它的另一个版本是(ax+by = gcd(a,b))一定有整数解。
它的证明是基于(gcd(a,b)=gcd(b,a\%b))的。在递归终点(b=0)时，显然有(x=1,y=0)。否则，
根据数学归纳法，假设下面这个方程成立：
根据模数的另一种定义：
我们可以得到：
整理得到：
因此，我们令(x'=y, y'=(x-(b/a)y))，就可以得到 (ax'+by'=gcd(a,b))了。
在程序中可以这样简写：
inline ll exgcd(ll a, ll b, ll &x, ll &y)
{
if(b == 0)
{
x = 1, y = 0;
return a;
}
ll d = exgcd(b, a % b, y, x);
y -= (a / b) * x;
return d;
}
一次同余方程FCE
一般形式为(ax equiv b pmod{p})。可以转换成(ax + py = b)的形式，然后用exgcd解出((x,y))。
假设(gcd(a,p)
ot mid b)，那么方程无解。
注意一下得到特解(x_0)后，方程的通解(x equiv x_0 pmod{frac{p}{gcd(a,p)}})。
逆元
(a)在模(p)意义下的逆元为(a^{-1})，则：(acdot a^{-1} equiv 1 pmod{p})。这是一个FCE，直接求解就可以了。
值得注意的是逆元的几种非常特殊的求法：
-
如果(p)是质数，那么(a)的逆元就是(a^{p-2}pmod{p})。
-
线性求逆元：(i^{-1} equiv (p-p/i) imes (p \% i)^{-1} pmod{p})
-
线性求阶乘逆元：(frac{1}{(i+1)!} imes (i+1) equiv frac{1}{i!} pmod{p})
稍微注意一下，在用阶乘求逆元时可能会出现(p | i!)的情况，此时阶乘算出来都是(0)。
中国剩余定理
设方程组：
当(p_1,p_2,cdots,p_k)两两互质时，设(P=prod_{i=1}^{k}p_i)，则方程一定有整数解：
其中((frac{P}{p_i})^{-1}(frac{P}{p_i})equiv 1 pmod{p_i})。
假设(p_1,p_2,cdots,p_k)不满足两两互质，那么应该用扩展版本。
设(x_{k-1})是通过EXCRT计算出的前(k-1)个方程的特解，而(P_{k-1} = prod_{i=1}^{k-1}p_i)。为了让它满足第(k)个方程，有：
已知(x_{k-1},P_{k-1},a_k,p_k)，这个方程就是一个关于(lambda)的FCE，可以直接计算。当然，如果它无解，那整个方程组就都无解了。新的特解(x_k=x_{k-1}+lambda P_{k-1})通过加上或减去(P_k)就是前(k)个方程的特解。
另外有个小技巧：求乘积(P_k)可以改为求前(k)个模数的(operatorname{lcm})。粘贴核心代码：
int main()
{
n=qr(1);
RP(i,1,n)
m[i]=qr(1ll),a[i]=qr(1ll);//x = a[i] (mod m[i])
ans=0,M=1;//M = prod
RP(i,1,n)
{
ll p=FCE(M,a[i]-ans,m[i]);
ans=ans+p*M;
M*=m[i]/gcd(M,m[i]);
ans=((ans%M)+M)%M;
}
printf("%lld",((ans%M)+M)%M);
return 0;
}
高次同余方程（指数型）
求解形如(a^{x}equiv b pmod{p})类型的，关于(x)的方程。其中有(a,p)互质，(x)非负。
这里没有一般的解析方法。但是由于(a,p)互质，我们可以自由计算逆元，进行和(a)有关的乘除法。
设(x = epsilon t - eta)，其中(t=lceilsqrt{p}
ceil)，(0leq eta < t)。事先声明，这样的记法其实有点“混用”的意味，但事实上相比普通的拉丁字母，它更有助于记忆理解。
方程变成了：
即：
对于所有的(bcdot a^{eta},0 leq eta < t)，我们将它插入一个Hash表中。枚举(epsilon)的所有取值(0,1,cdots,t)，依次判断在Hash表中是否有对应的(bcdot a^{eta})，然后拼凑得到答案。
这种方法叫做Baby step，Giant step。是一种类分块的优化枚举。
简单的线性代数
矩阵可以看成一个二维数表，可以用两个紧挨的下标(ij)来标识位置((i,j))处的元素。
矩阵加法就是对应的元素相加，即：(C_{ij}=A_{ij}+B_{ij})。
矩阵乘法相对复杂。只有形如(A_{(n imes k)})和(B_{(k imes m)})才能相乘得到一个形如(C_{(n imes m)})的矩阵。标准的公式：
可以理解为：(C_{ij}=A)的第(i)行行向量(\,cdot\,B)的第(j)列列向量。正所谓“行列式”，它说明行在前，列在后。
矩阵乘法的应用
对于一个(k)阶的齐次线性递推关系：
由于是线性递推关系，我们可以把所有的状态(x_n sim x_{n-k})记录下来，封装到一个向量中去。即：
此时这个(k)阶方矩阵：
又叫作转移矩阵。
如果知道这个递推关系，想要从(x_0,x_1,cdots,x_k)开始递推(c)次，这个过程就等价于对这个列向量左乘(c)次，即(T^c)。由于矩阵乘法满足结合律，这个过程可以用矩阵快速幂优化到(O(k^3log c))。
如果上面不是齐次递推关系，而是带了一个常数(p)（如(a_n=2a_n+3)），该怎么办呢？可以直接把常数接在状态向量的最下面，转移矩阵赋值(T_{k+1,k+1}=1)，原封不动地转移到新状态去就可以了。
当然，直接代入矩阵乘法的公式，对着公式依次赋值，才是最稳妥的做法。
当然，矩阵乘法的应用不止于此。考虑这样一个问题：
给定一个带边权无向连通图，求过(k)个点（不算起点），从(1)走到(n)的最短路？
设邻接矩阵(E_{ij})表示点(i)和(j)初始时的距离。如果两点无边相连，则赋值(infty)。由于(+)和( imes)都具有交换律和结合律，而(min)和(+)同样具有交换律和结合律，我们可以把矩阵乘法(langle+, imes angle)重定义成(langle min, + angle)。也就是说，每一次“乘法”：
因此最终的答案就是(E^{k}_{1,n})。可以用矩阵快速幂优化到(O(n^3log k))。
高斯消元
为节省空间，这里直接接一个链接过去。
线性基
狭义的线性基仅指异或空间下的线性基。其将大小为(N)的二进制集合，转换为(log)值域大小级别的集合，大大降低了时间和空间复杂度。但是线性基往往会破坏顺序结构：你可以得到集合中若干数异或后的最大值，最小值，但却无法具体得知你到底异或了哪些数。
线性基可以看成数表({p_0,p_1,cdots,p_k})，其中(p_i)表示为(1)的最高位为(i)的数。它可能是由原来若干个数异或得到的。
当你想新插入一个数(a)到线性基中，你需要从高到低枚举每一个已有的基底(p_i)。如果这个基底还没有数插入，你当然可以选择将这个位置加入(a)。但如果这一位上有数了，为了提高存储效率，你需要令(a = a operatorname{xor} p_i)从而消去这一位的值，然后继续检查低位的基底。
举个例子。这个线性基：
就是一个不错的线性基。而：
则非常糟糕。最高位为(2,4)的基底发生了重复，这两对基底应该分别异或起来，再插入这个线性基中。
粘贴一下核心代码：
void bits(ll x)
{
DRP(i,62,0)
{
if(!(x >> (ll)i))
continue;
if(!p[i])
{p[i] = x; break;}
x ^= p[i];
}
}
int main()
{
scanf("%d", &n);
RP(i,1,n)
scanf("%lld", &a[i]), bits(a[i]);
DRP(i,62,0)
if((ans ^ p[i]) > ans)
ans ^= p[i];
cout<<ans;
return 0;
}