模板 数学(2)

一些简单的数学相关算法模板。
本文主要集中于模相关的问题。

求 $\gcd$

欧几里得算法(辗转相除法)

辗转相除法基于 GCD 递归定理:

该定理可以用反证法进行证明,也可以通过证明左右两数分别被对方整除证明。此处不列出证法。

1
2
3
int gcd(int a, int b){
return (!b) ? a : gcd(b, a % b);
}

Stein算法

该算法相对欧几里得算法而言,在计算机上实现的效率更高,尤其是对于一些大数的计算。

该算法的步骤也非常简单,即为以下式子:

以上式子的正确性可以用同余式的性质证明,此处不给出。

1
2
3
4
5
6
7
8
9
10
11
int gcd(int a, int b){
if(!b) return a;
if(!a) return b;
if(a & 1){
if(b & 1) return gcd(b, abs(a - b) >> 1);
else return gcd(a, b >> 1);
}else{
if(b & 1) return gcd(a >> 1, b);
else return gcd(a >> 1, b >> 1) << 1;
}
}

迭代版本:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ll gcd(ll a, ll b){
ll c = 0;
while(a && b){
if(a & 1){
if(b & 1){
if(a > b) a = (a - b) >> 1;
else b = (b - a) >> 1;
}else b >>= 1;
}else{
if(b & 1) a >>= 1;
else a >>= 1, b >>= 1, c++;
}
}
return (a + b) << c;
}


扩展欧几里得

扩展欧几里得可以用来求解如下形式的线性同余方程组(由贝祖定理,解一定存在):

它的思想如下所示:假设我们已经求出了 $bx’+(a \bmod b)y’ = \gcd (b, a\bmod b)$ 的解 $x’, y’$,由 GCD 递归定理可以知道 $\gcd (a, b) = \gcd (b, a\bmod b)$,则:

比较两侧系数可以求出 $y=x’-\lfloor \frac{a}{b} \rfloor y’, x=y’$。当 $b = 0$ 时很容易求出 $x=1, y=0$,将这个结果向前迭代就可以求出答案。

1
2
3
4
5
6
7
8
9
int extgcd(int a, int b, int &x, int &y){
if(!b){
x = 1, y = 0;
return a;
}
int d = extgcd(b, a % b, y, x);
y -= x * (a / b);
return d;
}

需要说明的是,通过以上算法求出的 $x, y$ 满足:在所有可能的解中,$|x|+|y|$ 最小。


欧拉函数

欧拉定理

欧拉定理叙述如下:

欧拉定理降幂公式


解线性同余方程组

中国剩余定理

中国剩余定理(CRT)可以用来求解线性同余方程组。它的描述如下:
有$n$个线性同余方程组:

其中 两两互质,那么存在唯一的一个 $x$ 在模 使得上述方程组成立。
证明如下:首先构造一个这样的 $x$ 证明存在性。令 为除 外所有 的积,那么 。解 。然后令 ,则这个 $x$ 就是答案。(容易验证这一点)
然后再证明唯一性。假设存在另一个不同的解 $y$ ,那么 $x-y$ 同时是 的倍数。因为 两两互质。所以 $x-y$ 是 $M$ 的倍数。所以 $x\equiv y \pmod M$。矛盾。
这样就可以利用这个构造性算法构造出这个问题的解了:

1
2
3
4
5
6
7
8
9
10
11
12
int CRT(int m[], int a[], int n){
int M = 1, x = 0;
for(int i = 0; i < n; ++i)
M *= m[i];
for(int i = 0; i < n; ++i){
int Mi = M / m[i], b, q;
extgcd(Mi, m[i], b, q);
Mi = ((Mi * a[i]) % M) * b % M;
x = (x + Mi) % M;
}
return (x + M) % M;
}

扩展中国剩余定理

当 $m_i$ 不两两互质怎么办?依然可以求出线性同余方程组的解。
先考虑只有两个方程的情况。假如只有这样两个方程:

写 $x=b_1+tm_1$,带入第二个方程得 $m_1 a t\equiv b_2 - ab_1 \pmod {m_2}$。解出 $t$,如 $t\equiv c \pmod {m_2}$,就可以得到 $t=c+dm_2$,带入,从而 $x\equiv b_1+cm_1 \pmod {m_1m_2}$。这样就可以依次解开第三,四,……个方程。
实际上为了防止爆数据类型,完全可以把 $\pmod {m_1m_2}$ 改写为 $\pmod {lcm(m_1, m_2)}$。下面的代码就是这么做的。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
ll extCRT(ll b[], ll m[], int n){
ll M = m[1], ans = b[1], x, y;
for(int i = 2; i <= n; ++i){
ll d = extgcd(M, m[i], x, y), a_ = (b[i] - ans) % m[i];
if(a_ < 0) a_ += m[i];
if(a_ % d) return -1;
x = fstmul(x, a_ / d, m[i] / d);
if(x < 0) x += m[i] / d;
ll MM = M;
M /= d, M *= m[i];
ans = (ans + fstmul(x, MM, M)) % M;
if(ans < 0) ans += M;
}
return ans;
}


组合数取模

卢卡斯定理

一句话描述卢卡斯定理,就是:

其中 $p$ 是质数。
容易发现卢卡斯定理的表达式里面有一个递归结构,所以可以把它做成递归的。

1
2
3
ll Lucas(ll n, ll m){
return (!m) ? 1 : C(n % p, m % p) * Lucas(n / p, m / p) % p;
}

扩展卢卡斯定理

如果要取模的数 $p$ 不是质数怎么办?可以考虑将 $p$ 进行质因数分解,分解成 的形式,算出组合数对每一个 取模的值,然后用中国剩余定理合并即可。方便起见,记 $P=p_i^e$,即下面讨论的都是对某个质数的幂取模的情况。
那么问题就是怎么求 $\tbinom{n}{m} \bmod P$。先考虑 $n!, m!, (n-m)!$ 各自含因子 $p$ 的次数。设三者分别含 $p$ 的次数是 ,那么只有当 的时候才有 $\tbinom{n}{m} \bmod P \neq 0$。 可以在 $O(\log_p n)$ 的时间内求出来。
再考虑除去 $p$ 后的阶乘如何计数。可以发现阶乘中,乘的数在大于 $P$ 时会形成循环节,如 $1\times 2\times 4\times 5\times 7\times 8 \equiv 10\times 11\times 13\times 14\times 16\times 17 \pmod 9$。我们预处理出

这一循环节会形成总共 $\lfloor \frac{n}{P}\rfloor$ 次,带来的贡献就是 $fac[p^e]^{\lfloor\frac{n}{P}\rfloor}$。此外,还有不完整的循环节,它带来的贡献是 $fac[n\bmod P]$。而不在完整循环节或者不完整循环节中的项就都是 $p$ 的倍数,它们的贡献是 $p\times 2p\times \cdots \times p\lfloor \frac{n}{p}\rfloor=p^{\lfloor \frac{n}{p}\rfloor} \times \lfloor \frac{n}{p}\rfloor !$。看这一项,因为 $p$ 要算指数所以要单独拿开,而后面那个阶乘可以递归处理。所以这一部分的计算也能在 $O(\log_p n)$ 时间内完成。
于是答案就是

至此,扩展Lucas定理(虽然到这里已经和Lucas定理没什么关系了)算法完成。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
/* pair的first记录的是阶乘模P的值,second记录该阶乘中p的指数 */
pair<ll, ll> mod_fac(ll n, ll p, ll p_){
if(n == 0) return make_pair(1, 0);
pair<ll, ll> res = mod_fac(n / p, p, p_);
ll fir = res.first * poww(fac[p_], n / p_, p_) % p_;
fir = fir * fac[n % p_] % p_;
return make_pair(fir, res.second + n / p);
}
ll exLucas(ll n, ll m, ll p, ll k, ll p_){
ll cnt = 1;
fac[0] = 1;
for(ll i = 1; i <= p_; ++i){
if(cnt == p) cnt = 1, fac[i] = fac[i - 1];
else cnt++, fac[i] = fac[i - 1] * i % p_;
}
pair<ll, ll> fz = mod_fac(n, p, p_), fm1 = mod_fac(m, p, p_), fm2 = mod_fac(n - m, p, p_);
if(fz.second - fm1.second - fm2.second >= k) return 0;
ll ps = poww(p, fz.second - fm1.second - fm2.second, p_);
ps = ps * fz.first % p_;
ps = ps * inv(fm1.first, p_) % p_;
ps = ps * inv(fm2.first, p_) % p_;
return ps;
}


求逆元

单个逆元

$x$ 在模 $p$ 意义下有逆元的充要条件:$(x, p) = 1$。
如果 $p$ 是质数,根据费马小定理,可知 $x^{-1}\equiv x^{p-2} \pmod p$。如果 $p$ 不是质数,根据欧拉定理理论上也可以解,但求欧拉函数要花费 $O(\sqrt p)$ 的时间,不划算。因此可以用扩展欧几里得,按照求线性同余方程的形式求逆元。
两种情况下的时间复杂度都是 $O(\log n)$。

1
2
3
4
5
6
int inv(int x, int p){
int u, v;
if(extgcd(x, p, u, v) != 1)
return -1;
return u < 0 ? u + p: u;
}

求 $[1, n]$ 逆元

因为逆元是积性函数,所以理论上能用线性筛筛出。但是这么做效率不高。
考虑除法定理:对于 $i<p$,记 $p=ki+j(1\le j<p)$。那么 $ki+j\equiv 0 \pmod p$。移项有 $-ki\equiv j \pmod p\implies i^{-1}\equiv (-k)\times j^{-1} \pmod p$。
时间复杂度:从小到大处理,$O(n)$。

1
2
3
inv[1] = 1;
for(int i = 2; i < p; ++i)
inv[i] = 1ll * (p - p / i) * inv[p % i] % p;

求多个数逆元

如果数据并不符合 $[1,n]$ 的形式,应当怎么办?此时,我们设有 $n$ 个数,$a[1\cdots n]$,被模的是一个质数 $p$。

。因为逆元是完全积性的,从而 。那么,我们可以先求出 $pre[n]^{-1}$,然后根据下面的递推式:

可以求出给定的所有数的逆元。时间复杂度:$O(n+\log p)$。

1
2
3
4
5
6
7
8
9
void quickinv(int *a, int *pre, int n, int p, int *inv){
pre[0] = 1;
for (int i = 1; i <= n; ++i)
pre[i] = 1ll * pre[i - 1] * a[i] % p;
int preinv = poww(pre[n], p - 2);
for (int i = n; i >= 1; --i)
inv[i] = 1ll * pre[i - 1] * preinv % p,
preinv = 1ll * a[i] * preinv % p;
}