离散微积分

差分

设 `x_n` (`n ge 0`) 是一个数列. 定义 `x_n` 的差分为 `Delta x_n = x_(n+1) - x_n`, 其中 `Delta` 称为差分算子. 定义 `x_n` 是自己的零阶差分: `x_n = Delta^0 x_n`. 一般地, 对任意正整数 `k`, `x_n` 的 `k` 阶差分递归定义为 `Delta^k x_n = Delta( Delta^(k-1) x_n)`. 差分可以类比于微分运算.

上面定义的差分又叫向前差分. 还可以定义向后差分与中心差分 `grad x_n = x_n - x_(n-1)`, `quad delta x_n = x_(n+1/2) - x_(n-1/2)`. 下文我们只讨论向前差分.

利用差分表计算 `x_n = 2 n^2 + 3 n + 1` 的各阶差分.

1   6   15  28  45  66  91  ..
  5   9   13  17  21  25  ..
    4   4   4   4   4   ..
	  0   0   0   0   ..
	

使用组合数进行符号计算. 将 `x_n` 改写为 `x_n = 4(n;2) + 5(n;1) + (n;0)`, 从而 `Delta x_n = 4(n;1) + 5(n;0)`,
`Delta^2 x_n = 4(n;0)`,
`Delta^3 x_n = 0`.

如同微积分中我们有一张常用函数的导数表一样, 我们来推导一些常见数列的差分.

指数函数的差分公式 直接计算, `Delta 2^n = 2^n`, `quad Delta a^n = a^n (a-1)`, `quad Delta^k a^n = a^n (a-1)^k`. `2^n` 的差分是它自己, 这一性质对应于微积分中的 `"e"^x`. 如果记 `H_n = sum_(k=1)^n 1/k`, 则 `Delta H_n = 1/(n+1)`, 这一性质对应于微积分中的 `ln x`.

下降幂与上升幂 设 `k` 为非负整数, `n^(ul k) = prod_(0 le i lt k) (n-i) = n(n-1) cdots (n-k+1)`,
`n^(bar k) = prod_(0 le i lt k) (n+i) = n(n+1) cdots (n+k-1)`. 下降幂与二项式系数只相差系数 `k!`: `n^(ul k)/(k!) = (n;k)`. 下面考虑负指数的下降幂. 设 `k` 为正整数, 定义 `n^(ul(-k)) = prod_(1 le i le k) (n+i)^-1` `= 1/((n+1)(n+2) cdots (n+k))`. 特别地 `n^(ul(-1)) = 1/(n+1)`.

下降幂的差分公式 `k` 为非负整数时, 由 Pascal 公式 `(n+1; k) = (n; k) + (n; k-1)` 知 `Delta (n;k) = (n; k-1)`. 写成下降幂形式就是 `Delta n^(ul k) = k n^(ul(k-1))`. 可以验证上式对于负指数下降幂仍成立. 因此, 下降幂的差分仍是一个下降幂, 这对应于微分学的公式 `"d"/dx x^n = n x^(n-1)`. 然而 `n^k` 没有简单的差分公式, 所以在离散数学的世界里, 用下降幂代替普通的幂, 往往能简化计算.

求和

1. 如同积分与微分一样, 求和与差分具有线性性: `sum (a x_n + b y_n) = a sum x_n + b sum y_n`, `quad Delta (a x_n + b y_n) = a Delta x_n + b Delta y_n`.
2. 如同积分与微分互逆一样, 求和与差分互为逆运算: `sum_(0 le k lt n) Delta x_k = [x_k]_0^n = x_n - x_0`, `quad Delta sum_(0 le k lt n) x_k = x_n`. 特别当 `x_0 = 0` 时, 两式结果相等.
3. 分部求和公式 (Abel 变换) 因为 `Delta (x_k y_k) = x_(k+1) y_(k+1) - x_k y_(k+1) + x_k y_(k+1) - x_k y_k` `= y_(k+1) Delta x_k + x_k Delta y_k`, 所以 `sum_(0 le k lt n) x_k Delta y_k` `= [x_k y_k]_0^n - sum_(0 le k lt n) y_(k+1) Delta x_k`. 分部求和的几何意义见下图: 左上阶梯形 (蓝) = 大矩形 - 左下小矩形 (黄) - 右下阶梯形 (红).

类似可证 `sum_(k=0)^(n-m) x_k (y_(k+m) - y_k)` `= [x_k y_k]_0^n - sum_(k=0)^(n-m) y_(k+m) (x_(k+m)-x_k)`.

求 `sum k 2^k`, 其中求和范围是 `0 le k lt n`.

`sum k 2^k = sum k Delta 2^k = [k 2^k]_0^n - sum 2^(k+1) Delta k` `= n 2^n - 2 sum Delta 2^k = n 2^n - 2 [2^k]_0^n` `= n 2^n - 2(2^n - 1)`.

求 `sum H_k`, 求和范围是 `0 le k lt n`.

`sum H_k = sum H_k Delta k = [k H_k]_0^n - sum (k+1) Delta H_k` `= n(H_n - 1)`.

[来自 我是的的是我] `sum_(k=1)^n (H(k))/k = 1/2(H_n^2 + H_n^((2)))`, 其中 `H_n^((r)) = sum_(k=1)^n 1/k^r` 是广义调和数.

`n = 0` 时等式显然成立; 右边差分得 `1/2[H_n^2 - H_(n-1)^2 + H_n^((2)) - H_(n-1)^((2))]` `= 1/2[1/n(2H_(n-1) + 1/n) + 1/n^2]` `= (H_n)/n`, 等于左边的差分; 等式得证.

记左边为 `S`, 于是 `S = [H_k H_(k-1)]_1^(n+1) - sum_(k=1)^n H_k/(k+1)`
`= H_(n+1) H_n - sum_(k=1)^n H_k/(k+1)`
`= H_n^2 - sum_(k=1)^(n-1) H_k/(k+1)`
`= H_n^2 + sum_(k=1)^(n-1) 1/(k+1)^2 - sum_(k=1)^(n-1) H_(k+1)/(k+1)`
`= H_n^2 + H_n^((2)) - S`.

一般地, 对任意正整数 `a, b`, 运用分部求和得 `sum_(k=1)^n H_k^((a))/k^b = H_n^((a)) H_n^((b)) + H_n^((a+b)) - sum_(k=1)^n H_k^((b))/k^a`. 特别当 `a = b` 时, `sum_(k=1)^n H_k^((a))/k^a = 1/2((H_n^((a)))^2 + H_n^((2a)))`.

应用离散微积分的方法导出朱世杰恒等式: `sum_(k=0)^n (k;p)` `= sum_(k=0)^n Delta (k; p+1)` `= (n+1; p+1) - (0; p+1)` `= (n+1; p+1)`.

使用下降幂的差分求 `sum_(k=1)^n k^2`, `sum_(k=1)^n k^3`.

1. `sum_(k=1)^(n-1) k^2` `= sum k + sum k (k-1)` `= (n(n-1))/2 + (n(n-1)(n-2))/3` `= (n(n-1)(2n-1))/6`.
2. `sum_(k=1)^(n-1) k^3` `= sum k + 3sum k^(ul 2) + sum k^(ul 3)` `= n^(ul 2)/2 + n^(ul 3) + n^(ul 4)/4` `= ((n(n-1))/2)^2`.

几个公式

本节的最后导出几个有用的公式.

组合反演公式 关于数列 `{x_n}, {y_n}`, 若 `y_n = sum_(k=0)^n (n;k) x_k`, 则 `x_n = sum_(k=0)^n (n;k) (-1)^(n-k) y_k`. 如果记 `x_n = x^n`, `y_n = y^n`, 以上结论形式地记为 `y^n = (x+1)^n` `rArr x^n = (y-1)^n`. 或对称地, `y^n = (1-x)^n` `rArr x^n = (1-y)^n`. 用矩阵记为 `(a_(i j))^-1 = (-1)^(i-j) (a_(i j))`, `quad a_(i j) = (i;j)`.

代入验证, `sum_(k=0)^n (n;k) (-1)^(n-k) sum_(j=0)^k (k;j) x_j` `= sum_(j=0)^n x_j sum_(k=j)^n (n;k)(k;j) (-1)^(n-k)` `= sum_(j=0)^n x_j A_j`. 其中 `A_j = sum_(k=0)^(n-j) (n;k)(n-k;j) (-1)^k` `= sum_(k=0)^(n-j) (n;j)(n-j;k) (-1)^k` `= (n;j) (1-1)^(n-j)` `= {1, if j = n; 0, if j lt n:}` 因此 `sum_(j=0)^n x_j A_j = x_n`.

`k` 阶差分公式 `Delta^k x_n = sum_(j=0)^k (k;j) (-1)^(k-j) x_(n+j)`. 如果记 `x_n = x^n`, 上式形式地写为 `Delta^k x^n = (x-1)^k x^n`, 亦即 `Delta = x - 1`. 该表达式与指数函数的 `k` 阶差分一致.
反演得到离散 Maclaurin 公式: `x_(n+k) = sum_(j=0)^k (k;j) Delta^j x_n` `= sum_(j=0)^k (Delta^j x_n)/(j!) k^(ul j)`. 形式记为 `x^(n+k) = (Delta+1)^k x^n`, 亦即 `x = Delta + 1`. 向后差分、中心差分也有类似的表达式.

`k = 0` 时, `Delta^0 x_n = x_n` 成立; 假设等式对 `k-1` 成立, 则 `Delta^k x_n = Delta^(k-1) x_(n+1) - Delta^(k-1) x_n`
`= sum_(j=0)^(k-1) (k-1; j) (-1)^(k-1-j) x_(n+1+j)` `- sum_(j=0)^(k-1) (k-1; j) (-1)^(k-1-j) x_(n+j)`
`= x_(n+k) - (-1)^(k-1) x_n` `+ sum_(j=1)^(k-1) (-1)^(k-j) x_(n+j) ((k-1; j-1) + (k-1; j))`
`= sum_(j=0)^k (k;j) (-1)^(k-j) x_(n+j)`.

离散 Maclaurin 公式 设 `x_n` 为数列通项, `m` 为非负整数, 且对任意 `k gt m` 有 `Delta^k x_0 = 0`, 则 `x_n = sum_(k=0)^m (n;k) Delta^k x_0` `= sum_(k=0)^m (Delta^k x_0)/(k!) n^(ul k)`.

首先由大于 `m` 阶的差分都等于零知, `x_n` 是关于 `n` 的多项式.
因为 `(n; k)` 是关于 `n` 的 `k` 次多项式, `k = 0, 1, cdots, m`, 所以 `(n;k)` 线性无关, 它们组成多项式空间 `bbb P{::}_m[x]` 的基.
下面计算 `(n;k)` 的差分. 因为 `Delta_n (n;k) = (n;k-1)`, 所以 `{: Delta_n^j (n;k) |_(n=0) = {: (n; k-j) |_(n=0)` `= delta_(j k) = { 1, if j = k; 0, if j != k :}`. 现在将多项式 `x_n` 用基底线性表示为 `x_n = sum_(k=0)^m c_k (n; k)`, 则由差分的线性性知, `Delta_n^j x_0 = sum_(k=0)^m c_k {:Delta_n^j (n;k)|_(n=0)` `= sum_(k=0)^m c_k delta_(j k)` `= c_j`. 所以公式成立.

带余项的离散 Maclaurin 公式 `x_n = sum_(k=0)^m (n;k) Delta^k x_0 + R_m(n)`, 其中 `R_m(n) = sum_(0 le k lt n) (n-k-1; m) Delta^(m+1) x_k`.

设 `m gt 0`, 对 `R_m(n)` 分部求和, `{: R_m(n) ,= sum_(0 le k lt n-1) (n-k-1; m) Delta^(m+1) x_k; ,= [(n-k-1; m) Delta^m x_k]_(k=0)^(n-1) + sum_(0 le k lt n-1) (n-k-2; m-1) Delta^m x_(k+1); ,= -(n-1; m) Delta^m x_0 + sum_(1 le k lt n) (n-k-1; m-1) Delta^m x_k; ,= -(n-1; m) Delta^m x_0 - (n-1; m-1) Delta^m x_0 + sum_(0 le k lt n) (n-k-1; m-1) Delta^m x_k; ,= -(n; m) Delta^m x_0 + R_(m-1) (n).:}` 于是 `R_0(n) = sum_(1 le k le m) (R_(k-1)(n) - R_k(n)) + R_m(n)`, 即 `x_n - x_0 = sum_(1 le k le m) (n;k) Delta^k x_0 + R_m(n)`.

Stirling 数

Stirling 数 从线性代数的角度, `{x^k}_(k=0)^n` 与 `{x^(ul k)}_(k=0)^n` 是次数不超过 `n` 的多项式组成的线性空间 `bbb P{::}_n[x]` 的两个基底. 为了方便离散微积分运算, 需要将多项式的基由普通幂换为下降幂. 下面定义的两类 Stirling 数分别给出了两组基底之间互化的系数 (即过渡矩阵): `x^(ul n) = sum_(k=0)^n (-1)^(n-k) s(n, k) x^k`, `quad x^n = sum_(k=0)^n S(n, k) x^(ul k)`. 其中 `s(n, k)` 称为第一类 Stirling 数, 也记为 `[n;k]`, `S(n, k)` 称为第二类 Stirling 数, 也记为 `{n;k}`.

两类 Stirling 数的反演 由两类 Stirling 数构成的系数矩阵互逆知, `y_n = sum_(k=0)^n S(n, k) x_n` `iff` `x_n = sum_(k=0)^n (-1)^(n-k) s(n, k) y_n`.

我们先讨论第二类 Stirling 数.

第二类 Stirling 数

          1
        0   1
      0   1   1
    0   1   3   1
  0   1   7   6   1
0   1  15  25  10   1

递推公式 `S(n, k) = { 1, if n = k; 0, if n gt 0 and k = 0; S(n-1, k-1) + k S(n-1, k), otherwise; :}`

从定义式 `x^n = sum_(k=0)^n S(n, k) x^(ul k)` 出发. 考虑两边的 `n` 次项系数和常数项分别可得 `S(n, n) = 1` (`n ge 0`) 和 `S(m, 0) = 0` (`m gt 0`). 现在设 `n gt 0`, 两边对 `x` 差分得到 `(x+1)^n - x^n = sum_(k=1)^n S(n, k) k x^(ul (k-1))`. 另一方面 `(x+1)^n - x^n` `= 1/(x+1) (x+1)^(n+1) - x^n`
`= 1/(x+1) sum_(k=1)^(n+1) S(n+1, k) (x+1)^(ul k)` `- sum_(k=1)^n S(n, k) x^(ul k)`
`= sum_(k=1)^(n+1) S(n+1, k) x^(ul(k-1))` `- sum_(k=2)^(n+1) S(n, k-1) x^(ul(k-1))`
`= sum_(k=1)^n S(n+1, k) x^(ul(k-1))` `- sum_(k=1)^n S(n, k-1) x^(ul(k-1))`.
比较得 `k S(n, k) = S(n+1, k) - S(n, k-1)`.

通项公式 `S(n,k) = 1/k! sum_(j=0)^k (k;j) (-1)^j (k-j)^n`. 令 `S^k = k^n`, 上式形式写为 `S(n, k) = 1/k! (S-1)^k`.

只需验证这个通项公式适合边界条件和递推公式, 过程略.

在定义式 `x^n = sum_(k=0)^n S(n, k) x^(ul k)` 中, 考虑数列 `1^n, 2^n, 3^n, cdots, n^n`, 由离散 Maclaurin 公式和 `k` 阶差分公式 `S(n,k) = 1/k! {:Delta_x^k x^n|_(x=0)` `= 1/k! sum_(i=0)^k (k;j) (-1)^j (k-j)^n`.

组合意义 `n` 元集到 `k` 元集的满射数等于 `k! S(n, k)`. 换言之, 将 `n` 元集划分为 `k` 份, 每份不空的方法数为 `S(n, k)`.

应用容斥原理, `n` 元集到 `k` 元集的满射数目等于全部映射 `k^n`, 减去有一个元素无原像的情形 `(k;1) (k-1)^n`, 再加上有两个元素无原像的情形 `(k;2) (k-2)^n`, 再减去有三个元素无原像的情形... 最终得到 `sum_(j=0)^k (k;j) (-1)^j (k-j)^n = k! S(n, k)`.

假设将 `n` 元集划分为 `k` 份, 每份不空的方法数为 `T(n, k)`, 下证 `T(n, k)` 和 `S(n, k)` 满足相同的递推公式, 从而 `T(n, k) = S(n, k)`. 事实上, 显然有 `T(n, n) = 1`, `T(n, 1) = 1`. 接下来考虑 `n` 元集中的固定元素 `x_0`,
1. 若 `x_0` 在划分中独占一份 (住单人间, 没有室友), 只需再将剩下的 `n-1` 个元素分为 `k-1` 份, 方法数为 `T(n-1, k-1)`;
2. 若 `x_0` 在划分中不是独占一份, 只需将剩下的 `n-1` 个元素分为 `k` 份, 然后选择一份将 `x_0` 加入其中 (选择一个房间入住), 方法数为 `k T(n-1, k)`.
综上 `T(n, k) = T(n-1, k-1) + k T(n-1, k)`, 恰好满足 `S(n, k)` 的递推公式.

`n` 元集到 `k` 元集的映射, 相当于 `n` 元集到 `j` 元集的满射总和, `j = 0, 1, cdots, k`: `k^n = sum_(j=0)^k (k;j) S'(n, j)`, 反演得 `S'(n, k) = sum_(j=0)^k (k;j) (-1)^(k-j) j^n`, 因此 `S'(n, k) = k! S(n, k)`.

等幂求和 记 `(:n;k:) = S(n,k) k!`, `S_n(x) = sum_(i=1)^x i^n`. 则 `S_n(x) = sum_(k=0)^n (:n;k:) (x+1; k+1)`, `quad n ge 1`. 关于第一类和第二类 Bernoulli 数, 分别有 `B_n^(-) = sum_(k=0)^n (:n;k:) (-1)^k/(k+1)`, `quad n ge 0`,
`B_n^(+) = sum_(k=0)^n (:n;k:) (-1)^(k-1)/(k(k+1))`, `quad n ge 1`.

由第二类 Stirling 数定义 `i^n = sum_(k=0)^n (:n;k:) (i;k)`, 两边令 `i` 从 `0` 到 `x` 求和就得到第一式. 再和 Bernoulli 数的公式 `S_n(x) = 1/(n+1) sum_(k=0)^n (n+1;k) B_k^(+) x^(n+1-k)` 比较一次项系数, 就得到 `B_n^(+)` 的表达式. 同理比较 `S_n(x-1) = sum_(k=0)^n (:n;k:) (x;k+1)`, `quad n ge 1`,
`S_n(x-1) = 1/(n+1) sum_(k=0)^n (n+1;k) B_k^(-) x^(n+1-k)`, `quad n ge 1` 的一次项系数得 `B_n^(-) = sum_(k=0)^n (:n;k:) (-1)^k/(k+1)`, `quad n ge 1`. 容易验证上式对 `n=0` 也成立.

第一类 Stirling 数

          1
        0   1
      0   1   1
    0   2   3   1
  0   6   11  6   1
0   24  50  35  10   1

验算: `[ 1; 0, 1; 0, 1, 1; 0, 1, 3, 1; 0, 1, 7, 6, 1; 0, 1, 15, 25, 10, 1; ] [ 1; 0, 1; 0, -1, 1; 0, 2, -3, 1; 0, -6, 11, -6, 1; 0, 24, -50, 35, -10, 1; ] = bm I`

递推公式 `s(n, k) = { 1, if n = k; 0, if n gt 0 and k = 0; s(n-1, k-1) + (n-1) s(n-1, k), otherwise; :}` 第一类 Stirling 数没有实用的通项公式.

Euler 数

Euler 数 对正整数 `n` 和 正整数 `k le n`, `A_(n, k) = { 1, if n = 1 or k = n; (n-k+1) A_(n-1,k-1) + k A_(n-1, k), otherwise; :}` 表示 `1~n` 的排列中, 恰有 `k-1` 个相邻逆序对的排列数. 比如 `1, 3, 5, 4, 2` 中, `(5, 4)`, `(4, 2)` 是两个相邻逆序对. 由定义 `sum_(k=1)^n A_(n, k) = n!`.

          1
        1   1
      1   4   1
    1  11  11   1
  1  26  66  26   1
1  57 302 302  57   1

[来自 诗许] 多项式的基底表示 `x^n = sum_(k=1)^n (x+k-1; n) A_(n, k)`. 此公式两边对 `n` 求和即得到又一个等幂求和公式.

归纳证明. `n = 1` 时 `x = (x+1-1; 1) A_(1, k)` 成立. 假设命题对 `n` 成立, 考虑 `n+1` 的情形: `sum_(k=1)^(n+1) (x+k-1; n+1) A_(n+1, k)`
`= (x; n+1) + (x+n; n+1)` `+ sum_(k=2)^n (x+k-1; n+1)[(n-k+2) A_(n,k-1) + k A_(n, k)]`
`= (x; n+1) + (x+n; n+1)` `+ sum_(k=1)^(n-1) (x+k; n+1) (n-k+1) A_(n,k) + sum_(k=2)^n (x+k-1; n+1) k A_(n,k)`
`= (x; n+1) + (x+n; n+1) + (x+1; n+1) n + (x+n-1; n+1) n` `+ sum_(k=2)^(n-1) A_(n,k) (x+k-1; n)((x+k)(n-k+1) + k(x+k-n-1))/(n+1)`
`= sum_(k=1)^n A_(n, k) (x+k-1; n) x`
`= x^(n+1)`.

Euler 数的显式表达式为 `A_(n, k) = sum_(j=0)^k (n+1; j) (-1)^j (k-j)^n`.

杂例

哥隆尺 (Golomb ruler) [OEIS A004137] [科学网 一把数学的尺子 哥隆尺]
尺子上面有 `m` 个整数刻度 `x_1, cdots, x_m in ZZ`, 现在给出两个条件:
1. 这尺子可以量出 `0` 到 `n` 之间的所有整数长度, 即 `[0..n] sube {|x_i - x_j|: i, j in [1..m]}`;
2. 每个整数长度没有多于一种测量方法, 即 `AA i, j, k, l in [1..m]`, `x_i - x_j = x_k - x_l iff i = k and j = l`;
如果尺子满足条件 2, 则称它为哥隆尺, 如果同时满足两个条件, 则称为完美哥隆尺. 尺子长度大于 6 时不存在完美哥隆尺.
问题: 若尺子满足条件 1, 求尺子长度 `n` 的最大值.

from itertools import combinations

def check(t, n):
    s = set(abs(t[i] - t[j]) for i in range(len(t)) for j in range(i, len(t)))
    s = s.union(t)
    s = s.union(n-v for v in t)
    return len(s) == n

# 输入尺子长度 n, 求最小的 k 使得尺子能量出 [0..n] 的所有长度
def golomb(n, mink=1):
    if n < 3:
        return n+1
    for k in range(mink, n-1):
        for t in combinations(range(1,n), k):
            if check(t, n):
                print(0, *t, n)
                return len(t)+2
    return -1

print(golomb(6))