近期算法竞赛中的数学问题

CodeJam 2021 qualified

Reversort Engineering
本例涉及到摊还分析，动态规划中有一个很常见的相关技巧就是未来费用提前计算

• 首先很容易想到对 $[1\cdots n]$ 的排序代价 $C$ 满足

$$n-1 \leqslant C \leqslant \frac{n(n+1)}{2} - 1$$

• 可以反着来构造，最开始令序列有序 $a = \{1, 2, \cdots n\}$
  $d(i)$ 表示从 $i$ 出发，找到最小元素下标为 $j$，需要走的距离，即 $\text{len}[i \cdots j]$
  我们是反着构造的，所以实际上 $d(i)=\text{len}[i \cdots j]$，$j$ 为从 $i \to n$ 最大元素的下标
  初始阶段，令 $C \leftarrow C-(n-1) \quad d(i) = 1$

  • 如果 $C > 0$，就说明 $\text{len}[i \cdots j] > 1$，产生了额外代价

• 于是可以贪心构造，如果在排序的某个阶段，$\forall i \in [n-2 \to 0]$
  $[i \cdots j]$，让 $a_j$ 为最小，$a_i$ 为第 $2$ 小
  反着构造就是 $[i \cdots j]$，$a_i$ 为当前最小，$a_j$ 为第 $2$ 小，反转 $[i \cdots j]$
  此时代价就是每一段 $[i \cdots j]$ 能够取到的最坏代价 $C'$

• $\textbf{for} \ \forall i \in [n-2 \to 0]$，$i$ 能往前走的最长距离是 $\max d(i) = n-1-i$
  只要 $C > 0$，执行

  $$\textbf{for} [1 \cdots n-1-i]: \ d(i) += 1, \ C -= 1$$

  即让 $i$ 尽可能往前走，贪心地让 $[i\cdots n-1]$ 中最小的元素尽可能一直位于 $a_{n-1}$ 处
  $d(i)$ 初始时为 $1$，$\max d(i) = n-i$
  此时需要反转的区间是 $[i, j] = [i \cdots i+d(i)-1]$

  • 反转 $[i, j] = [i, i+d(i)-1]$，在第 $i-1$ 轮迭代中

    $$\begin{gathered} a_{i-1} < a[i \cdots j], \quad \text{对于区间} [i-1\cdots j] \\ a_{i-1} \text{为最小,} \quad a_{j} \text{第 2 小} \end{gathered}$$

    $C > 0$，第 $i-1$ 阶段，$[i\cdots n-1]$ 中最小的元素尽可能一直位于 $a_{n-1}$ 处
    又因为 $a_{i-1} < \forall \ a[i \cdots n-1]$
    $a_{i-1}$ 最小，$a_{n-1}$ 第 $2$ 小，反转 $a[i-1, n-1]$ 即可

• 反着迭代 $n-1$ 次，即可得到原序列
  最后对于不足以填充 $[i \cdots n-1]$ 的 $d(i)$，也不会影响结果
  因为此时 $[i \cdots i+d(i)-1]$ 一定是最后一次迭代的区间，反转之后对应原序列第 $i$ 次迭代
  恰好满足 $a_{i+d(i)-1}$ 是 $[i \cdots i+d(i)-1]$ 区间最小元素，对 $i-1$ 及之前的迭代没有影响
  因为前面的元素 $a[0 \cdots i-1] < a[i \cdots n]$，$\forall j \in [0, i-1], d(j) = 1$

#define _for(i, l, r) for (int i = (l); i < (r); i++)
using namespace std;
const int maxn = 100 + 10;
int n, C;

void solve() {
    vector<int> d(n-1, 1);
    C -= (n-1);
    for (int i = n-2; i >= 0; i--) {
        for (int _ = 1; _ <= n-1-i; _ ++) {
            if (C > 0) {
                C -= 1, d[i] += 1;
            }
        }
    }

    vector<int> a(n);
    for (int i = 0; i < n; i++) a[i] = i+1;
    for (int i = n-2; i >= 0; i--) {
        reverse(a.begin()+i, a.begin()+i+d[i]);
    }
    for (auto x : a) printf("%d ", x);
    printf("\n");
}

int main() {
    freopen("input.txt", "r", stdin);
    int T, kase = 0;
    cin >> T;
    while (T--) {
        printf("Case #%d: ", ++kase);
        cin >> n >> C;

        // impossible
        int lo = n-1, hi = n*(n+1) / 2 - 1;
        if (C < lo || C > hi) {
            puts("IMPOSSIBLE");
            continue;
        }

        // solve
        solve();
    }
}

数学基础

筛法求素数

朴素筛法

vector<int> primes;
bool st[maxn]; // init st[...] = 0
// st = true 表示这个数是合数

void get_primes(int n) {
    for (int i = 2; i <= n; i++) {
        if (st[i]) continue;
        primes.push_back(i);
        for (int j = 2*i; j <= n; j += i) st[j] = true;
    }
}

线性筛法

vector<int> primes;
bool st[maxn];

void get_primes() {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) primes.push_back(i);
        for (int j = 0; primes[j] <= n/i; j++) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
}

约数

试除法分解素因子

vector<PII> div;
void divide(int x) {
    for (int i = 2; i <= x/i; i++) {
        if (x % i == 0) {
            int s = 0;
            while (x % i == 0) s++, x /= i;
            // i^s
            div.push_back({i, s});
        }
    }
    if (x > 1) div.push_back({x, 1});
    // x^1
}

阶乘分解

inline int get(int n, int p) {
    int s = 0;
    while (n) s += n/p, n /= p;
    return s;
}

组合数的几种求法

递推求组合数

$$\dbinom{n}{k} = \binom{n-1}{k} + \binom{n-1}{k-1}$$

const int maxn = 1e4 + 10;
const int mod = 100000007;
int c[maxn][maxn], n;
void calc() {
    c[0][0] = 1;
    for (int i = 1; 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;
        }
}

乘法逆元求组合数

乘法逆元

若 $(b, m) = 1$，并且 $b \mid a$，那么 $a / b = a \cdot (b^{-1} \bmod m)$
$(b^{-1} \bmod m)$ 即为 $b$ 的 $\bmod \ m$ 的乘法逆元

注意前提，$(b, m) = 1$，即 $b, m$ 互素

$$\begin{aligned} a/b &= a \cdot b^{-1} \equiv a/b \cdot (b \cdot b^{-1}) (\bmod \ m) \\ &\Rightarrow b \cdot b^{-1} \equiv 1 (\bmod \ m) \end{aligned}$$

$\textbf{algorithm}$
于是在求组合数，求 $a / b$ 并且对一个比较大的素数 $p$ 取模的时候，注意一定要 $(b, p) = 1$

• $a / b \ (\bmod \ p) \Longrightarrow a \cdot b^{-1} \ (\bmod \ p)$，以下省略 $(\bmod \ p)$
  除以 $b$ 等价于乘以 $b$ 的乘法逆元 $b^{-1}$
• 求$b$ 的乘法逆元 $b^{-1}$，实际上就是求 $bx \equiv 1 (\bmod \ p)$ 的解
  该式的唯一解为

  $$\begin{aligned} x &\equiv b^{\phi (p) - 1} \ (\bmod \ p) \equiv b^{p-2} \ (\bmod \ p) \\ b^{-1} &= \text{qpow}(b, p-2) \quad \text{可以用快速幂迅速求出} \end{aligned}$$

证明

1. Euler-Fermat定理 设 $(a, m) = 1$，则有

   $$a^{\phi(m)} \equiv 1 (\bmod \ m)$$

   设 $A = \{b_1, b_2 \cdots b_{\phi (m)}\}$ 为模 $m$ 的一个简化剩余系，因为 $(a, m) = 1$，所以
   $B= \{ab_1, ab_2, \cdots ab_{\phi (m)}\}$ 也是模 $m$ 的一个简化剩余系，于是 $\prod A$ 和 $\prod B$ 同余

   $$b_1b_2 \cdots b_{\phi (m)} \equiv a^{\phi (m)} b_1b_2 \cdots b_{\phi (m)} (\bmod \ m)$$

   每一个 $b_i$ 都与 $m$ 互素，即可得到定理
2. 如果一个素数 $p$ 不能整除 $a$，那么

   $$a^{p-1} \equiv 1 (\bmod \ p)$$

   根据欧拉函数 $\phi (p) = p-1$ 即可得到
3. 对任一整数 $a$ 与任一素数 $p$，我们有

   $$a^{p} \equiv a(\bmod \ p)$$

   这也是很显然的，因为 $p \not | \ a$，那么这个结论就是 $2$，如果 $p \mid a$
   那么 $a^{p}, a$ 都同于 $0 \bmod p$
4. 如果 $(a, m) = 1$，那么一次同余式 $ax \equiv b \ (\bmod \ m)$ 的唯一解

   $$x \equiv ba^{\phi (m) - 1} (\bmod \ m)$$

   根据 Euler-Fermat 定理，上式满足 $ax \equiv b \ (\bmod \ m)$，又因为 $(a, m) = 1$
   所以这个解对于 $\bmod \ m$ 是唯一的

逆元表

当取模数很大的时候，需要利用递推式，线性求出 $[1\cdots n]$ 中每个数在 $\bmod \ p$ 意义下的乘法逆元

$$(p \bmod \ i) + \left\lfloor \frac{p}{i} \right\rfloor \cdot i = p$$

两边同时对 $p$ 取 $\bmod$

$$\left\lfloor \frac{p}{i} \right\rfloor \cdot i \equiv -(p \bmod \ i)$$

从而推出

$$i^{-1} \equiv -\left\lfloor \frac{p}{i} \right\rfloor \cdot (p \bmod \ i)^{-1} \bmod \ p$$

注意 $- \displaystyle \left\lfloor \frac{p}{i} \right\rfloor \bmod p$ 一般写成 $\displaystyle \left(p - \displaystyle \left\lfloor \frac{p}{i} \right\rfloor \right) \bmod p$

const int maxn = 1e6 + 10;
ll inv[maxn];

void init_inv() {
    inv[1] = 1;
    for (int i = 2; i < n; i++) inv[i] = (mod - mod/i) * inv[mod % i] % mod;
}

预处理逆元求组合数

$$\binom{x}{y} = \frac{x!}{y!(x-y)!} \ \bmod \ p$$

• 预处理 $x \in [1\cdots N]$ 所有数的阶乘 $\bmod \ p$，得到 $\bold{fac}[x] = x! \bmod \ p$
• 预处理 $x \in [1\cdots N]$ 阶乘的 $\bmod p$ 乘法逆元，得到 $\bold{infac}[x] = (x!)^{-1} \bmod p$
  这可以通过快速幂 $\bold{infac}[x] = \bold{infac}[x-1] \cdot \textbf{qpow}(x, p-2)$
• 最后的结果为

  $$\bold{fac}[x] \cdot \bold{infac}[y] \cdot \bold{infac}[x-y]$$

const int maxn = 1e5 + 10, N = 1e5;
const int mod = 10000007;
ll fac[maxn], infac[maxn];

void calc() {
    fac[0] = infac[0] = 1;
    for (int i = 1; i <= N; i++) {
        fac[i] = fac[i-1] * i % mod;
        infac[i] = infac[i-1] * ksm(i, mod-2, mod) % mod;
    }
}
int C(int x, int y) {
    ll ans = fac[x] * infac[y] * infac[x-y] % mod;
    return ans;
}

也可以直接利用组合数公式，预处理阶乘，逆元

ll inv[maxn], fac[maxn], infac[maxn];
inline ll mul(ll a, ll b) {
    return a * b % mod;
}

inline ll mul(ll a, ll b, ll c) {
    return a * b % mod * c % mod;
}

void pre() {
    inv[0] = inv[1] = 1;
    for (int i = 2; i < maxn; i++) inv[i] = (mod - mod/i) * inv[mod%i] % mod;

    fac[0] = infac[0] = 1;
    for (int i = 1; i < maxn; i++) fac[i] = mul(fac[i-1], i);
    for (int i = 1; i < maxn; i++) infac[i] = mul(infac[i-1], inv[i]);
}

inline int C(int x, int y) {
    if (x < y) return 0;
    return mul(fac[x], infac[x-y], infac[y]);
}

Lucas 定理求组合数

Lucas 定理

对于素数 $p$

$$\binom{n}{m} \bmod \ p = \binom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor} \binom{n \bmod \ p}{m \bmod p} \bmod p$$

证明
构造生成函数 $g(x) = (1 + x)^p$，有如下式子成立

$$(1 + x)^p \equiv 1 + x^p (\bmod \ p)$$

根据二项式展开，$\forall \ k \in [1, p-1], p \mid \binom{p}{k}$

$$\begin{aligned} (1+x)^p &= 1 + \binom{p}{1}x + \binom{p}{2} x^2 + \cdots \binom{p}{p-1}x^{p-1} + x^p \\ &\equiv 1 + x^p (\bmod \ p) \end{aligned}$$

又因为 $n = \lfloor \frac{n}{p} \rfloor \cdot p + a_0, \quad a_0 = n \bmod p$，定义等式 $(1)$

$$\begin{aligned} (1 + x)^n &= (1+x)^{\lfloor \frac{n}{p} \rfloor \cdot p} \cdot (1+x)^{a_0} \\ &\equiv (1 + x^p)^{\lfloor \frac{n}{p} \rfloor} \cdot (1 + x)^{a_0} (\bmod \ p) \\ &\equiv \left( \sum_{i=0}^{\lfloor n/p \rfloor} \binom{\lfloor n/p \rfloor}{i} x^{i \cdot p} \right) \cdot \left( \sum_{j=0}^{a_0} \binom{a_0}{j} x^j \right) \end{aligned}$$

根据费马小定理，定义函数

$$\begin{aligned} f(x) &\equiv (1 + bx^m)^p (\bmod \ p) \\ &\equiv (1 + bx^m) (\bmod \ p) \end{aligned}$$

又根据以上推论可知

$$\begin{aligned} f(x) &\equiv (1 + bx^m)^p (\bmod \ p) \\ &\equiv (1 + bx^{pm}) (\bmod \ p) \equiv (1 + bx^{pm})^{p} (\bmod \ p) \\ &\equiv f(x^p) \end{aligned}$$

• $n \choose m$ 即为左边第 $m$ 项，$x$ 的幂指数为 $m$，那么 $m = i \cdot p + j$
• 因为 $\bmod \ p$，对于右边，因为 $a_0 = n \bmod p \leqslant p-1$，对于 $\binom{a_0}{j}$ 没有 $p$ 的倍数项
  $(1+x)^n$ 二项展开，根据上面的结论 $f(x) \equiv f(x^p)$，只有 $x$ 的幂指为 $p$ 的倍数时，$(1)$ 才成立

  $$\binom{\lfloor n/p \rfloor}{i} = \binom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor}, \quad j = m \bmod \ p, \ a_0 = n \bmod p$$

• 生成函数系数相乘，即有

  $$\binom{n}{m} \bmod p = \binom{\lfloor n/p \rfloor}{\lfloor m/p \rfloor} \binom{n \bmod p}{m \bmod p} \bmod p$$

Lucas 定理求组合数

$$C(a, b) = \binom{a}{b} = \frac{a(a-1) \cdots (a-b+1)}{b!}$$

const int mod = 1000007;
int ksm(int a, int k) {
    int res = 1;
    for (; k; k >>= 1) {
        if (k & 1) res = (ll)res * a % mod;
        a = (ll)a * a % mod;
    }
    return res;
}
int C(int a, int b, int mod) {
    if (a < b) return 0;

    ll x = 1, y = 1;
    for (int i = a, j = 1; j <= b; i--, j++) {
        x = x * i % mod;
        y = y * j % mod;
    }
    return x * (ll)ksm(y, mod-2, mod) % mod;
}
int lucas(ll a, ll b) {
    if (a < mod && b < mod) return C(a, b);
    return (ll)C(a % mod, b % mod) * lucas(a/mod, b/mod) % mod;
}

分解质因数求组合数

当我们需要求出组合数的真实值，而非对某个数的余数时，分解质因数比较好用

• 筛法预处理范围内的所有质数
• 根据组合数公式

  $$\binom{a}{b} = a! / b! / (a-b)!$$

  可以用阶乘分解来求，阶乘分解函数 $\textbf{get}(n, p)$ 表示 $n!$ 中 $p$ 的次数
• 高精度乘法将所有素因子相乘

const int maxn = 1e5 + 10;
int st[maxn];
vector<int> primes, sum;
int a, b;

void get_primes(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) primes.push_back(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) {
    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 (auto x : a) {
        t += x * b;
        c.push_back(t % 10);
        t /= 10;
    }
    while (t) c.push_back(t % 10), t /= 10;
    return c;
}

void calc() {
    get_primes(a);
    sum.resize(primes.size());
    for (int i = 0; i < primes.size(); 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 < primes.size(); i++)
        _for(_, 0, sum[i]) res = mul(res, primes[i]);
}

数学基础问题

GCD LCM
$G \cdot L = a \cdot b$

• $G > L \ \textbf{or} \ G \bmod L \neq 0$，无解
• 否则如果 $a \leqslant b$，那么 $a$ 最小一定是 $\text{gcd}(a, b)$

Benefit
$A \cdot B = C \cdot \text{gcd}(A, B)$

$$B = \frac{C}{A} \cdot \text{gcd}(A, B) = C' \cdot \text{gcd}(A, B)$$

• 如果 $B$ 最小必然有 $\text{gcd}(A, B) = 1$，否则
  $B \leftarrow B / \text{gcd}(A, B)$ 可以得到一个更小的 $B$

• $\textbf{let} \ C' \leftarrow C/A, \quad B \leftarrow C'$
  $d \leftarrow \textbf{gcd}(A, B)$
  $\textbf{while} \quad d > 1$:
  $\quad C' \times = d, \quad \text{gcd}(A, B) / = d$

• $\text{gcd}(A, B) /= d \Leftrightarrow A /= d$
  最后令 $B = C' \times \prod d$

• 迭代过程如下

  $$\begin{aligned} B &= C' \cdot 1 \\ B &= C'/d \cdot (1 \cdot d) = B/d \cdot (1 \cdot d) \\ B &= B/d^2 \cdot (1 \cdot d^2) \\ \vdots \\ B &= B/\text{gcd} \cdot (1 \cdot \text{gcd}) \end{aligned}$$

• 算法执行过程实际上是倒着执行上述迭代

How do you add
本例实际上是一个隔板法

$$\binom{n+k-1}{k-1}$$

就是本例的答案