模板 数学（2）

发表于 | 分类于 |
字数统计: 2,234 字

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

求 $\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;
}

Powered By Valine
v1.5.2