CodeJam 2021 qualified

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

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

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

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

    • 如果C>0C > 0,就说明len[ij]>1\text{len}[i \cdots j] > 1,产生了额外代价
  • 于是可以贪心构造,如果在排序的某个阶段,i[n20]\forall i \in [n-2 \to 0]
    [ij][i \cdots j],让aja_j 为最小,aia_i 为第22
    反着构造就是[ij][i \cdots j]aia_i 为当前最小,aja_j 为第22 小,反转[ij][i \cdots j]
    此时代价就是每一段[ij][i \cdots j] 能够取到的最坏代价CC'

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

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

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

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

      ai1<a[ij],对于区间[i1j]ai1为最小,aj第 2 小\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>0C > 0,第i1i-1 阶段,[in1][i\cdots n-1] 中最小的元素尽可能一直位于an1a_{n-1}
      又因为ai1< a[in1]a_{i-1} < \forall \ a[i \cdots n-1]
      ai1a_{i-1} 最小,an1a_{n-1}22 小,反转a[i1,n1]a[i-1, n-1] 即可
  • 反着迭代n1n-1 次,即可得到原序列
    最后对于不足以填充[in1][i \cdots n-1]d(i)d(i),也不会影响结果
    因为此时[ii+d(i)1][i \cdots i+d(i)-1] 一定是最后一次迭代的区间,反转之后对应原序列第ii 次迭代
    恰好满足ai+d(i)1a_{i+d(i)-1}[ii+d(i)1][i \cdots i+d(i)-1] 区间最小元素,对i1i-1 及之前的迭代没有影响
    因为前面的元素a[0i1]<a[in]a[0 \cdots i-1] < a[i \cdots n]j[0,i1],d(j)=1\forall j \in [0, i-1], d(j) = 1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
#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();
}
}

数学基础

筛法求素数

朴素筛法

1
2
3
4
5
6
7
8
9
10
11
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;
}
}

线性筛法
01

1
2
3
4
5
6
7
8
9
10
11
12
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;
}
}
}

约数

试除法分解素因子

1
2
3
4
5
6
7
8
9
10
11
12
13
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
}

阶乘分解

02

1
2
3
4
5
inline int get(int n, int p) {
int s = 0;
while (n) s += n/p, n /= p;
return s;
}

组合数的几种求法

递推求组合数

(nk)=(n1k)+(n1k1)\dbinom{n}{k} = \binom{n-1}{k} + \binom{n-1}{k-1}

1
2
3
4
5
6
7
8
9
10
11
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, m) = 1,并且bab \mid a,那么a/b=a(b1modm)a / b = a \cdot (b^{-1} \bmod m)
(b1modm)(b^{-1} \bmod m) 即为bbmod m\bmod \ m 的乘法逆元

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

a/b=ab1a/b(bb1)(mod m)bb11(mod 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}

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

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

    xbϕ(p)1 (mod p)bp2 (mod p)b1=qpow(b,p2)可以用快速幂迅速求出\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, m) = 1,则有

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

    A={b1,b2bϕ(m)}A = \{b_1, b_2 \cdots b_{\phi (m)}\} 为模mm 的一个简化剩余系,因为(a,m)=1(a, m) = 1,所以
    B={ab1,ab2,abϕ(m)}B= \{ab_1, ab_2, \cdots ab_{\phi (m)}\} 也是模mm 的一个简化剩余系,于是A\prod AB\prod B 同余

    b1b2bϕ(m)aϕ(m)b1b2bϕ(m)(mod m)b_1b_2 \cdots b_{\phi (m)} \equiv a^{\phi (m)} b_1b_2 \cdots b_{\phi (m)} (\bmod \ m)

    每一个bib_i 都与mm 互素,即可得到定理
  2. 如果一个素数pp 不能整除aa,那么

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

    根据欧拉函数ϕ(p)=p1\phi (p) = p-1 即可得到
  3. 对任一整数aa 与任一素数pp,我们有

    apa(mod p)a^{p} \equiv a(\bmod \ p)

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

    xbaϕ(m)1(mod m)x \equiv ba^{\phi (m) - 1} (\bmod \ m)

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

预处理逆元求组合数

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

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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] = fac[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;
}

Lucas 定理求组合数

Lucas 定理

对于素数pp

(nm)mod p=(n/pm/p)(nmod pmmodp)modp\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)pg(x) = (1 + x)^p,有如下式子成立

(1+x)p1+xp(mod p)(1 + x)^p \equiv 1 + x^p (\bmod \ p)

根据二项式展开, k[1,p1],p(pk)\forall \ k \in [1, p-1], p \mid \binom{p}{k}

(1+x)p=1+(p1)x+(p2)x2+(pp1)xp1+xp1+xp(mod p)\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=npp+a0,a0=nmodpn = \lfloor \frac{n}{p} \rfloor \cdot p + a_0, \quad a_0 = n \bmod p,定义等式(1)(1)

(1+x)n=(1+x)npp(1+x)a0(1+xp)np(1+x)a0(mod p)(i=0n/p(n/pi)xip)(j=0a0(a0j)xj)\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}

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

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

又根据以上推论可知

f(x)(1+bxm)p(mod p)(1+bxpm)(mod p)(1+bxpm)p(mod p)f(xp)\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}

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

    (n/pi)=(n/pm/p),j=mmod p, a0=nmodp\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

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

    (nm)modp=(n/pm/p)(nmodpmmodp)modp\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)=(ab)=a(a1)(ab+1)b!C(a, b) = \binom{a}{b} = \frac{a(a-1) \cdots (a-b+1)}{b!}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
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 res = 1;
for (int i = 1, j = a; i <= b; i++, j--) {
res = (ll)res * j % mod;
res = (ll)res * qpow(i, mod-2) % mod;
}
return res;
}
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;
}

分解质因数求组合数

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

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

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

    可以用阶乘分解来求,阶乘分解函数get(n,p)\textbf{get}(n, p) 表示n!n!pp 的次数
  • 高精度乘法将所有素因子相乘
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
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
GL=abG \cdot L = a \cdot b

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

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

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

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

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

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

  • 迭代过程如下

    B=C1B=C/d(1d)=B/d(1d)B=B/d2(1d2)B=B/gcd(1gcd)\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
本例实际上是一个隔板法
choose

(n+k1k1)\binom{n+k-1}{k-1}

就是本例的答案