矩阵乘法

二维 乘 二维

void mul(vector<vector<ll>>& c, vector<vector<ll>> a, vector<vector<ll>> b, int p) {

    int n = a.size();

    for (int i = 0; i < n; ++ i) {
        for (int j = 0; j < n; ++ j) {
            c[i][j] = 0;
            for (int k = 0; k < n; ++ k) {
                c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % p;
            }
        }
    } 

}

一维 乘 二维

void mul(vector<ll>& c, vector<ll> a, vector<vector<ll>> b, int p) {
    int n = a.size();

    for (int i = 0; i < n; ++ i) {
        c[i] = 0;
        for (int k = 0; k < n; ++ k) {
            c[i] = (c[i] + a[k] * b[k][i]) % p;
        }
    }
}

高斯消元

// a[N][N]是增广矩阵
int gauss()
{
    int c, r;
    for (c = 0, r = 0; c < n; c ++ )
    {
        int t = r;
        for (int i = r; i < n; i ++ )   // 找到绝对值最大的行
            if (fabs(a[i][c]) > fabs(a[t][c]))
                t = i;

        if (fabs(a[t][c]) < eps) continue;

        for (int i = c; i <= n; i ++ ) swap(a[t][i], a[r][i]);      // 将绝对值最大的行换到最顶端
        for (int i = n; i >= c; i -- ) a[r][i] /= a[r][c];      // 将当前行的首位变成1
        for (int i = r + 1; i < n; i ++ )       // 用当前行将下面所有的列消成0
            if (fabs(a[i][c]) > eps)
                for (int j = n; j >= c; j -- )
                    a[i][j] -= a[r][j] * a[i][c];

        r ++ ;
    }

    if (r < n)
    {
        for (int i = r; i < n; i ++ )
            if (fabs(a[i][n]) > eps)
                return 2; // 无解
        return 1; // 有无穷多组解
    }

    for (int i = n - 1; i >= 0; i -- )
        for (int j = i + 1; j < n; j ++ )
            a[i][n] -= a[i][j] * a[j][n];

    return 0; // 有唯一解
}

高精度

高精度加法

// 倒着存放
vector<int> add(vector<int> a, vector<int> b) {
    vector<int> ret;
    int len = min((int)a.size(), (int)b.size());
    int at = 0, i = 0;
    for (; i < len; ++ i) {
        int t = a[i] + b[i] + at;
        at = t / 10;
        ret.push_back(t % 10);
    } 
    while (i < (int)a.size()) {
        int t = a[i] + at;
        at = t / 10;
        ret.push_back(t % 10);
        ++ i;
    }
    while (i < (int)b.size()) {
        int t = b[i] + at;
        at = t / 10;
        ret.push_back(t % 10);
        ++ i;
    }
    while (at) {
        ret.push_back(at % 10);
        at /= 10;
    }
    return ret; 
}

高精度减法

vector<int> sub(vector<int>& a, vector<int>& b) {

    vector<int> ret;
    int t = 0;
    for (int i = 0; i < a.size(); ++ i) {
        t = a[i] - t;
        if (i < b.size()) t -= b[i];
        ret.push_back((t + 10) % 10);
        if (t < 0) t = 1;
        else t = 0;
    }
    while (ret.size() > 1 && ret.back() == 0) ret.pop_back();
    return ret;
}

高精度乘法

// 高精度 乘 低精度
vector<int> mul(vector<int>& a, int b) {
    vector<int> ret;
    int t = 0;
    for (int i = 0; i < a.size() || t; ++ i) {
        if (i < a.size()) t += a[i] * b;
        ret.push_back(t % 10);
        t /= 10;
    }

    while (ret.size() > 1 && ret.back() == 0) ret.pop_back();
    return ret; 
}
// 高精度 乘 高精度
vector<int> mul(vector<int> a, vector<int> b) {
    vector<int> ret(a.size() + b.size() + 5, 0);
    for (int i = 0; i < a.size(); ++ i) {
        for (int j = 0; j < b.size(); ++ j) {
            ret[i + j] += a[i] * b[j];
        }
    }   

    for (int i = 0; i + 1 < ret.size(); ++ i) {
        ret[i + 1] += ret[i] / 10;
        ret[i] %= 10;
    } 

    while (ret.size() > 1 && ret.back() == 0) ret.pop_back();

    return ret;
}

高精度除法

int yu = 0; // 除不尽后的余数
vector<int> divv(vector<int>& a, int b) {

    vector<int> ret;
    int t = 0;
    for (int i = a.size() - 1; i >= 0; -- i ) {
        t = t * 10 + a[i];
        ret.push_back(t / b);
        t %= b;
    }
    yu = t;
    reverse(ret.begin(), ret.end());
    while (ret.size() > 1 && ret.back() == 0) ret.pop_back();
    return ret;
}

欧几里得算法

详细理论

常用

求 a, b 最大公约数即 $gcd(a, b)$
有 $gcd(a, b) = gcd(b, a\;mod\; b)$

int gcd(int a, int b)
{
    return b ? gcd(b, a % b) : a;
}

扩展

---

（裴蜀定理）

设 a, b 为正整数，则关于$ax + by = c$ 的方程有整数解当且仅当 c 是 gcd(a, b) 的倍数。

---

求$ax + by = c$ 可以通过求 $b{x_0} + (a\;mod\;b){y_0} = c$
即可得$\begin{cases}x = {y_0} \\ y = {x_0} - \lfloor \frac{a}{b} \rfloor {y_0} \end{cases}$

该方程通解
$\begin{cases}{x_k} = x + k \cdot \frac{b}{gcd(a,b)} \\ {y_k} = y - k \cdot \frac{a}{gcd(a, b)} \end{cases}(k \in z)$

int exgcd(int a, int b, int &x, int &y)
{
    if (b == 0)
    {
        x = 1;
        y = 0;
        return a;
    }
    int d = exgcd(b, a % b, y, x); //这里交换了x和y
    y -= (a / b) * x;
    return d;
}

逆元

详细理论

欧几里德算法

同余方程

$ax \equiv 1\; (mod\; b)$

有解条件：$gcd(a, b) = 1$
即解 $ax - by = 1$
通解：${x_k} = x + kb$
即：求 $x \; mod \;b$

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;
}
ll inv(ll a, ll p)
{
    ll x, y;
    if (exgcd(a, p, x, y) != 1) // 无解的情形
        return -1;
    return (x % p + p) % p;
}

费马小定理

若 p 为质数且 $gcd(a, p) = 1$ 则有 ${a^{p-1} \equiv 1 \;(mod\;p)}$
可知$inv(a) = {a^{p-2}}\;(mod\;p)$

#define ll long long
ll qpow(ll a, ll n, ll p)// 快速幂
{
    ll ans = 1;
    while (n)
    {
        if (n & 1)
            ans = ans % p * a % p;
        a = a % p * a % p;
        n >>= 1;
    }
    return ans;
}
ll inv(ll a, ll p)
{
    return qpow(a, p - 2, p);
}

线性递推

$inv(a) = - \lfloor \frac{p}{a} \rfloor \cdot inv(p\;mod\;a)\;(mod\;p)$

#define ll long long
// 多次对不同的p使用需要清空Inv数组
ll Inv[N] = {0, 1};
ll mod(ll a, ll p)
{
    return (a % p + p) % p;
}
ll inv(ll a, ll p)
{
    if (Inv[a])
        return Inv[a];
    Inv[a] = mod(-p / a * inv(p % a, p), p);
    return Inv[a];
}

中国剩余定理

解同余方程组$\begin{cases}x\equiv 2 \;(mod\;3) \\ x\equiv 3 \;(mod\;5) \\ x\equiv 2 \;(mod\;7)\end{cases}$

余数数组为 b, 模数数组为 a
设

$p = \prod_{i = 1} ^ n {a_i}$

和

${r_i} = \frac{p}{{a_i}}$

则有

$x \equiv \sum_{i=1}^n{b_i}\cdot {r_i} \cdot inv({r_i})\mid \;(mod\;p)$

#define ll long long
ll CRT(ll a[], ll b[], ll n)  // a是模数数组，b是余数数组，n是数组长度
{
    ll p = 1, x = 0;
    for (int i = 0; i < n; ++i)
        p *= a[i];
    for (int i = 0; i < n; ++i)
    {
        ll r = p / a[i];
        x += (b[i] * r * inv(r, a[i])) % p;   // 逆元的求法参见上篇文章，或者下面有完整代码
    }
    return x % p;
}

质数筛

bool isnp[MAXN];
vector<int> primes; // 质数表
void init(int n)
{
    for (int i = 2; i <= n; i++)
    {
        if (!isnp[i])
            primes.push_back(i);
        for (int p : primes)
        {
            if (p * i > n)
                break;
            isnp[p * i] = 1;
            if (i % p == 0)
                break;
        }
    }
}

欧拉函数

详细理论

int phi(int n)
{
    int res = n;
    for (int i = 2; i * i <= n; i++)
    {
        if (n % i == 0)
            res = res / i * (i - 1); // 先除再乘防止溢出
        while (n % i == 0) // 每个质因数只处理一次，可以把已经找到的质因数除干净
            n /= i;
    }
    if (n > 1)
        res = res / n * (n - 1); // 最后剩下的部分也是原来的n的质因数
    return res;
}

int phi[N];
bool isnp[N];
vector<int> primes;
void init(int n)
{
    phi[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!isnp[i])
            primes.push_back(i), phi[i] = i - 1; // 性质一，指数为1的情形
        for (int p : primes)
        {
            if (p * i > n)
                break;
            isnp[p * i] = 1;
            if (i % p == 0)
            {
                phi[p * i] = phi[i] * p; // 性质二
                break;
            }
            else
                phi[p * i] = phi[p] * phi[i]; // 这时肯定互质，用性质三
        }
    }
}

Min25 筛

#define ll long long

inline ll V2IDX(ll v, ll N, ll Ndr, ll nv) {
    return v >= Ndr ? (N/v - 1) : (nv - v);
}

ll primesum(ll N) {     //求取1~N的所有质数和
    ll *S;
    ll *V;

    ll r = (ll)sqrt(N);
    ll Ndr = N/r;

    assert(r*r <= N and (r+1)*(r+1) > N);

    ll nv = r + Ndr - 1;

    V = new ll[nv];
    S = new ll[nv];

    for (ll i=0; i<r; i++) {
        V[i] = N/(i+1);
    }
    for (ll i=r; i<nv; i++) {
        V[i] = V[i-1] - 1;
    }

    for (ll i=0; i<nv; i++) {
        S[i] = V[i] * (V[i] + 1) / 2 - 1;
    }

    for (ll p=2; p<=r; p++) {
        if (S[nv-p] > S[nv-p+1]) {
            ll sp = S[nv-p+1];
            ll p2 = p*p;
            for (ll i=0; i<nv; i++) {
                if (V[i] >= p2) {
                    S[i] -= p * (S[V2IDX(V[i]/p, N, Ndr, nv)] - sp);
                } else {
                    break;
                }
            }
        }
    }

    return S[0];
}

扩展欧拉定理

详细理论

$a^p \equiv a^{b\;mod\;\varphi(m)}\;(mod\;m), 其中 gcd(a, m) = 1$

求组合数

递推

C_a^b = C_{a - 1} ^ {b} + C_{a - 1} ^ {b - 1}

// c[a][b] 表示从a个苹果中选b个的方案数
for (int i = 0; i < N; i ++ )
    for (int j = 0; j <= i; j ++ )
        if (!j) c[i][j] = 1;
        else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;

预处理逆元

$C_a^b = \frac{a!}{(a - b)!b!}$

首先预处理出所有阶乘取模的余数fact[N]，以及所有阶乘取模的逆元infact[N]
如果取模的数是质数，可以用费马小定理求逆元
int qmi(int a, int k, int p)    // 快速幂模板
{
    int res = 1;
    while (k)
    {
        if (k & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        k >>= 1;
    }
    return res;
}

// 预处理阶乘的余数和阶乘逆元的余数
fact[0] = infact[0] = 1;
for (int i = 1; i < N; i ++ )
{
    fact[i] = (LL)fact[i - 1] * i % mod;
    infact[i] = (LL)infact[i - 1] * qmi(i, mod - 2, mod) % mod;
}

分解质因数

$C_a^b = \frac{a!}{(a - b)!b!}$

当我们需要求出组合数的真实值，而非对某个数的余数时，分解质因数的方式比较好用：
    1. 筛法求出范围内的所有质数
    2. 通过 C(a, b) = a! / b! / (a - b)! 这个公式求出每个质因子的次数。 n! 中p的次数是 n / p + n / p^2 + n / p^3 + ...
    3. 用高精度乘法将所有质因子相乘

int primes[N], cnt;     // 存储所有质数
int sum[N];     // 存储每个质数的次数
bool st[N];     // 存储每个数是否已被筛掉

void get_primes(int n)      // 线性筛法求素数
{
    for (int i = 2; i <= n; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i;
        for (int j = 0; primes[j] <= n / i; j ++ )
        {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

int get(int n, int p)       // 求n！中的次数
{
    int res = 0;
    while (n)
    {
        res += n / p;
        n /= p;
    }
    return res;
}

vector<int> mul(vector<int> a, int b)       // 高精度乘低精度模板
{
    vector<int> c;
    int t = 0;
    for (int i = 0; i < a.size(); i ++ )
    {
        t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }

    while (t)
    {
        c.push_back(t % 10);
        t /= 10;
    }

    return c;
}

get_primes(a);  // 预处理范围内的所有质数

for (int i = 0; i < cnt; i ++ )     // 求每个质因数的次数
{
    int p = primes[i];
    sum[i] = get(a, p) - get(b, p) - get(a - b, p);
}

vector<int> res;
res.push_back(1);

for (int i = 0; i < cnt; i ++ )     // 用高精度乘法将所有质因子相乘
    for (int j = 0; j < sum[i]; j ++ )
        res = mul(res, primes[i]);

Lucas 定理

详细理论

C_m^n \equiv C_{m\;mod\;p}^{n\;mod\;p}\cdot C_{\lfloor\frac{m}{p} \rfloor}^{\lfloor\frac{n}{p}\rfloor}\;(mod\;p)

// 需要先预处理出fact[]，即阶乘
inline ll C(ll m, ll n, ll p)
{
    return m < n ? 0 : fact[m] * inv(fact[n], p) % p * inv(fact[m - n], p) % p;
}
inline ll lucas(ll m, ll n, ll p)
{
    return n == 0 ? 1 % p : lucas(m / p, n / p, p) * C(m % p, n % p, p) % p;
}

大步小步算法

详细理论

$a^x\equiv b\;(mod\;m), 其中 gcd(a, m) = 1$

ll BSGS(ll a, ll b, ll m)
{
    static unordered_map<ll, ll> hs;
    hs.clear();
    ll cur = 1, t = sqrt(m) + 1;
    for (int B = 1; B <= t; ++B)
    {
        (cur *= a) %= m;
        hs[b * cur % m] = B; // 哈希表中存B的值
    }
    ll now = cur; // 此时cur = a^t
    for (int A = 1; A <= t; ++A)
    {
        auto it = hs.find(now);
        if (it != hs.end())
            return A * t - it->second;
        (now *= cur) %= m;
    }
    return -1; // 没有找到，无解
}

扩展

// 修改版的BSGS，额外带一个系数
ll BSGS(ll a, ll b, ll m, ll k = 1)
{
    static unordered_map<ll, ll> hs;
    hs.clear();
    ll cur = 1, t = sqrt(m) + 1;
    for (int B = 1; B <= t; ++B)
    {
        (cur *= a) %= m;
        hs[b * cur % m] = B; // 哈希表中存B的值
    }
    ll now = cur * k % m;
    for (int A = 1; A <= t; ++A)
    {
        auto it = hs.find(now);
        if (it != hs.end()) return A * t - it->second;
        (now *= cur) %= m;
    }
    return -INF; // 这里因为要多次加1，要返回更小的负数
}
ll exBSGS(ll a, ll b, ll m, ll k = 1)
{
    ll A = a %= m, B = b %= m, M = m;
    if (b == 1) return 0;
    ll cur = 1 % m;
    for (int i = 0;; i++)
    {
        if (cur == B) return i;
        cur = cur * A % M;
        ll d = gcd(a, m);
        if (b % d) return -INF;
        if (d == 1) return BSGS(a, b, m, k * a % m) + i + 1;
        k = k * a / d % m, b /= d, m /= d; // 相当于在递归求解exBSGS(a, b / d, m / d, k * a / d % m)
    }
}