Chengyuan Ma's Blog

幽居默处而观万物之变,尽其自然之理,而断之于中

0%

自从上一次2月6日写关于SAM的笔记之后似乎这个博客就停更了好久~

(**为什么在归档里有许多本篇之前SAM之后的?因为出于以下的原因某人把本地的笔记整理传了上来,通过改变date属性穿越时空把一些笔记送回了时间线的前端(即本地笔记的创建日期),如果看修改日期还是2019/12/8,因为是今天传的=_=**)

在这大半年里还是发生了很多事情。

比如说省选粗心失掉100分导致没进队。

比如说CSP再心态爆炸导致分数不是最理想。

比如自己的探究课题结题了

比如自己暑假里参加实习去GIX之类的

比如自己SAT也非常翻车的样子

比如说自己仍然在纠结~

回忆起来记得自己当初创设这个博客,很大程度上是觉得别人大佬的博客都很酷!

然后传了一些文章。

之后立flag发奋要搞OI,刚开始肯定是非常兴奋的,于是遇到很多题都会写文章,处于一个兴奋期吧。

当时也把去年在狗社搞得一些机器人的成果传了上去。

之后频率慢慢就往下降了。

一方面自己去年NOIP没考好,所以后半学期就偏向于划水了在咸鱼学校不脱产的高中生是这个样子的。

另一方面有时候也觉得维护博客有点烦,因为是公网上的总想着自己写得要尽善尽美一些,久而久之写一篇文章的时间甚至超过了做题目本身的时间。“啊刷题要紧反正做过的题目我大致都有印象博客这种浪费时间的东西还是停一停吧”有很长一段时间都这么想过对自己记忆的谜之自信然而最后也没有刷特别多的题,真是一个悲伤的故事,我们学校的事情实在是太多了。

自己倒是在本地维护了零零散散的一些做题心得。但是因为写得比较散,并没有上传。

导致一路咕到现在的这种状况。

最近才意识到博客不就是分享自己心得的吗……其实不必要过于正式。

感觉自己朋友圈里面的很多东西都可以发到博客里来的样子。(反正朋友圈里看懂的人似乎也不多233)

发现的很多好玩的也可以发的样子。

脑洞也可以发的样子。

这样想来维护博客似乎也没那么累,对于就是在本地写好Markdown,然后git同步一下,可能也就图片传图床可能会比较烦的样子。

果然还是自己比较懒啊

这个样子怎么行呢!?

果然这个博客还是要有点生气啊!

恢复更新!(Flag预定)

先整理一下本地的一些零零碎碎的做题笔记,然后之后再看看有啥好玩的以后也发在这上面吧!

然而深知自己尿性,但愿这不要又是一次短暂的心血来潮吧!

本地笔记整理~原稿写于11/10,2019/12/8上传

题意

若干年后。

主持魔法学院重建工作的 \omega坐在围墙上,与一旁的少女交谈。

「老师老师,你知道咱们学院是怎么传递信息的吗?」

「当然知道啊!」\omega 笑着摸了摸少女的头说道,「魔法学院有 n处地点,它们依次编号为 1n 。我们搭设了n-1条线路,第i条线路连接 u_i号地点和 v_i号地点。当然,任意两处地点间都存在一条由线路构成的路径。」

少女调皮地歪了下脑袋。「但是老师,为什么我在1号地点,经常收不到 n号地点的信息呢?」

「因为怕狡猾的敌人窃听我们啊。」\omega回答道,「学院内部只使用 m个频道,每处地点分别会选择性地屏蔽一些频道。两处地点能够直接通信,当且仅当它们被一条线路所连接,且存在两处地点都没有屏蔽的频道。为了测评通讯质量,我们定义『通讯值』为最多能把n处地点分成多少非空组,满足不同组内的地点无法直接通信。」

「这样一成不变的话,行动不是迟早会暴露在敌人面前吗?」

「是啊。所以我们打算轮流执行q次修改,第i次修改会改变x_i号地点的频道屏蔽状态。现在我来考考你,在修改之前以及每次修改后,你能分别求出通讯值是多少吗?」

少女微微一笑,仿佛心中已有了答案。

这题目的题面,有点尬啊。

解法

考虑u,v如果S_u \cap S_v \neq \emptyset那么u,v必在同一组内,因此若定义该条件为连边条件的话问题就变成了动态计算树子图的连通分量个数。

因为每一次都是单点修改,所以每次只要计算少连了多少边就可以得知答案增加了多少。

因此,问题被转化为如何快速计算一个节点所有相邻节点中与其集合交集非空的节点数,或者交集为空的节点数。

对于一般的树,暴力大约是足够了,时间复杂度就是\mathcal O(q\deg(u)),但是对于菊花图\deg(u) = \mathcal{O}(n),于是乎就没了。

考虑转化为有根树,每个节点只有一个父亲,对于这个父亲我们现场算,我们接着维护一个数组f[u][S]表示u的儿子当中集合为S的节点个数,然后在询问时变为\mathcal O\left(2^m\right)查询。在修改时\mathcal O(1)修改父亲的数组值。

这当然是可以的,但是m\le 20。还是不够。

注意到如果A \cap B = \emptyset,那么B \subseteq \overline A ,因此,我们还可以维护数组f[u][S]记录u的儿子中集合为S的子集的节点数,询问便是f[u][\overline{S_u}],变成了\mathcal O(1)的复杂度,倒是修改变为了\mathcal O\left(2^m\right)。依然不行。

考虑合并两种解法,对于集合的前\left\lfloor \frac{m}{2} \right\rfloor种元素,我们用第一种方法,对于集合的后\left\lceil \frac{m}{2} \right\rceil种元素,我们用第二种方法,那么单次的查询修改就都到了\mathcal O\left(2^{m\over 2}\right)。时间复杂度上是能通过该题了。

但是这张表要\mathcal O\left(n2^m\right)的空间不就炸了吗?

为了解决炸空间的问题,我们只选择度数大于d使用这种方法,对于小度数的节点暴力就完事了。

总体复杂度就是\mathcal O \left(nd + \frac{n}{d}2^{m\over 2}\right)

代码

#include <cstdio>
#include <algorithm>
using namespace std;
const int N = 100010;
const int MS = 1 << 10;
const int D = 200;
const int T = 500;
int n, m, fa[N], w[N], all, tp, ans = 0;
int fst[N], nxt[2 * N], to[2 * N], tot = 1;
int f[D][MS][MS], id[N], cnt = 0 , deg[N];
void link(int u, int v) {
    deg[u]++; deg[v]++;
    nxt[++tot] = fst[u];
    fst[u] = tot; to[tot] = v;
}
void dfs(int u, int p) {
    fa[u] = p;
    for (int e = fst[u]; e; e = nxt[e])
        if (to[e] != p) dfs(to[e], u);
}
inline int readbin() {
    char s[22]; scanf("%s", s);
    int ans = 0;
    for (int i = 0; i < m; i++)
        ans = (ans << 1) | (s[i] - '0');
    return ans;
}
void modify(int u, int x, int d) {
    int k = id[u];
    int mask = (1 << tp) - 1;
    int fr = x >> tp, bk = x & mask, c = mask ^ bk;
    for (int s = c; s; s = c & (s - 1))
        f[k][fr][mask ^ s] += d;
    f[k][fr][mask] += d;
}
int query(int u, int x) {
    int k = id[u];
    x = all ^ x;
    int mask = (1 << tp) - 1;
    int fr = x >> tp, bk = x & mask, ans = 0;
    for (int s = fr; s; s = fr & (s - 1))
        ans += f[k][s][bk];
    return ans + f[k][0][bk];
}
void solve() {
    all = (1 << m) - 1;
    tp = (m + 1) / 2;
    for (int i = 1; i <= n; i++) 
        if (deg[i] > T) id[i] = ++cnt;
    for (int i = 1; i <= n; i++)
        if (id[fa[i]]) modify(fa[i], w[i], 1);
    int Q; scanf("%d", &Q);
    while (Q--) {
        int x, c; scanf("%d ", &x);
        c = readbin();
        int last = 0, now = 0;
        if (id[x]) {
            last = query(x, w[x]);
            now = query(x, c);
            if (fa[x]) {
                last += (w[x] & w[fa[x]]) == 0;
                now += (c & w[fa[x]]) == 0;
            }
        } else {
            for (int e = fst[x]; e; e = nxt[e]) {
                last += (w[x] & w[to[e]]) == 0;
                now += (c & w[to[e]]) == 0;
            }
        }
        printf("%d\n", (ans += now - last));
        if (id[fa[x]]) {
            modify(fa[x], w[x], -1);
            modify(fa[x], c, 1);
        }
        w[x] = c; 
    }
}
int main() {
    scanf("%d%d", &n, &m);
    for (int i = 1; i < n; i++) {
        int u, v; scanf("%d%d", &u, &v);
        link(u, v); link(v, u);
    }
    for (int i = 1; i <= n; i++) w[i] = readbin(); 
    dfs(1, 0); ans = n;
    for (int i = 2; i <= n; i++) if (w[i] & w[fa[i]]) ans--;
    printf("%d\n", ans);
    solve();
    return 0;
}

本地笔记整理重发,原文写于2019/10/28,2019/12/8上传

题面

题意:求使得关于\{x_i\}的方程

\sum_{i=1}^n a_ix_i = b
存在非负整数解的的b \in [l, r]的个数。

数据范围:n \le 12, a_i \le 5\times 10^5, l, r \in 10^{12}

神题。不妨设a_1\{a_i\}当中的最小值。则若对于任意的b \equiv x \pmod{a_1}当中的b如果符合题意有解,那么b+a_1也一定有解。令满足这个性质的最小的bb_x

考虑在此条件下如何计算解的数目,设m以下满足题意的bf(m)个,则可以得出计算式:

f(m) = \sum_{x=0}^{a_1-1} [b_x \le m]\left(\left\lfloor\frac{m-b_x}{a_1}\right\rfloor+1\right)
原理显然。

答案就是f(r) - f(l - 1)

接下来考虑如何求取b_x。考虑构建具有a_1个点,编号分别为0,1,2,\cdots,a_1 - 1的图,对于编号u的点,向(u + a_i) \bmod a_1, \ (i=2,3,\cdots, n)连边。表示在模a_1意义下,如果u可以有解,那么u+a_i也一定有解。

由于0是一定有平凡解的,我们从0开始跑最短路,就可以求出b_x。这背后的思想无疑是非常巧妙的。

复杂度就是\mathcal{O}(na_1(\log n + \log a_1)),即为Dijkstra的复杂度。

代码

#include <cstdio>
#include <algorithm>
#include <queue>

using namespace std;
typedef long long ll;
typedef pair<ll, int> pli;
const int N = 13;
const int V = 500010;
const int E = N * V;
const ll INF = 1ll << 60;
ll n, Bl, Br, a[N], mn, len[E];
int fst[V], nxt[E], to[E], tot = 1;
ll dis[V];
void link(int u, int v, ll w) {
    nxt[++tot] = fst[u];
    fst[u] = tot; to[tot] = v; len[tot] = w;
}
void dijkstra() {
    static priority_queue<pli, vector<pli>, greater<pli> > q;
    fill(dis, dis + mn, INF);
    dis[0] = 0; q.push(make_pair(0, 0));
    while (!q.empty()) {
        pli p = q.top(); q.pop();
        int u = p.second;
        if (dis[u] != p.first) continue;
        for (int e = fst[u]; e; e = nxt[e]) {
            int v = to[e];
            if (dis[v] > dis[u] + len[e]) {
                dis[v] = dis[u] + len[e];
                q.push(make_pair(dis[v], v));
            }
        }
    }
}
ll calc(ll x) {
    ll ans = 0;
    for (int i = 0; i < mn; i++)
        if (x >= dis[i]) ans += (x - dis[i]) / mn + 1;
    return ans;
}
int main() {
    scanf("%lld%lld%lld", &n, &Bl, &Br);
    for (int i = 1; i <= n; i++) scanf("%lld", &a[i]);
    sort(a + 1, a + 1 + n); mn = a[1];
    for (int i = 0; i < mn; i++)
        for (int j = 2; j <= n; j++)
            link(i, (i + a[j]) % mn, a[j]);
    dijkstra();
    printf("%lld\n", calc(Br) - calc(Bl - 1));  
    return 0;
}

本地笔记整理第一发,比较散,2019/12/8上传

分析

题目要求求取带修改的树上最大独立集,朴素\BigO{n^2}DP如下:

\begin{aligned}
    \dpa[u][0] &= \sum_{v \in \mathrm{ch}(u)} \max\left\{\dpa[v][1], \dpa[v][0]\right\} \\
    \dpa[u][1] &= \mathrm{val}[u] + \sum_{v \in \mathrm{ch}(u)} \dpa[v][0]
\end{aligned}

但是这么一次转移是\BigO{n}的,考虑优化,定义\mathrm{lch}(u)u的轻儿子集合。\mathrm{hch}(u)u的重儿子,稍微修改一下,定义

\begin{aligned}
\ldp[u][0] &= \sum_{v \in \lch{u}} \max\left\{\dpa[v][1], \dpa[v][0]\right\} \\
\ldp[u][1] &= \mathrm{val}[u] + \sum_{v \in \lch{u}} \dpa[v][0] \\
\end{aligned}

则有

\begin{aligned}
    \dpa[u][0] &= \max\left\{\dpa[\hch{u}][0],\dpa[\hch{u}][1]\right\} + \ldp[u][0]\\
    \dpa[u][1] &= \dpa[\hch{u}][0] + \ldp[u][1]
\end{aligned}

使用Max-Plus Algebra的方式重写:

\begin{pmatrix}
    \dpa[u][0] \\
    \dpa[u][1]
\end{pmatrix} = 
\begin{pmatrix}
    \ldp[u][0] & \ldp[u][0] \\
    \ldp[u][1] & -\infty
\end{pmatrix}
\begin{pmatrix}
    \dpa[\hch{u}][0] \\
    \dpa[\hch{u}][1]
\end{pmatrix}
接下来只需要树链剖分维护矩阵转移就行了,单次修改的时间复杂度为\BigO{\log^2n}

代码

#include <cstdio>
#include <algorithm>
#include <cstring>

using namespace std;
const int N = 100010;
const int T = 4 * N;
const int INF = 0x3f3f3f3f;
int n;
int fst[N], nxt[2 * N], to[2 * N], tot = 1;
struct vec {
    int a[2];
    vec() { memset(a, 0, sizeof a); }
    int &operator [](int x) { return a[x]; }
    const int &operator [](int x) const { return a[x]; }
} eye;
struct mat {
    int a[2][2];
    mat() { memset(a, 0, sizeof a); }   
    int *operator [](int x) { return a[x]; }
    const int *operator [](int x) const { return a[x]; }
    mat operator *(const mat &b) const {
        mat c;
        for (int i = 0; i < 2; i++) {
            for (int j = 0; j < 2; j++) {
                c[i][j] = -INF;
                for (int k = 0; k < 2; k++)
                    c[i][j] = max(c[i][j], a[i][k] + b[k][j]);
            }
        }
        return c;
    }
    vec operator *(const vec &b) const {
        vec c;
        for (int i = 0; i < 2; i++) {
            c[i] = -INF;
            for (int j = 0; j < 2; j++)
                c[i] = max(c[i], a[i][j] + b[j]);
        }
        return c;
    }
};
void link(int u, int v) {
    nxt[++tot] = fst[u];
    fst[u] = tot; to[tot] = v;
}
int cnt = 0, id[N], top[N], dep[N], bot[N], sz[N], son[N], fa[N];
int a[N], ldp[N][2], dp[N][2];
mat val[T], init[N];
void dfs1(int u, int p) {
    dep[u] = dep[p] + 1; fa[u] = p;
    sz[u] = 1; int mx = -1;
    for (int e = fst[u]; e; e = nxt[e]) {
        int v = to[e];
        if (v == p) continue;
        dfs1(v, u); sz[u] += sz[v];
        if (sz[v] > mx) mx = sz[v], son[u] = v;
    }
}
void dfs2(int u, int tp) {
    id[u] = ++cnt, top[u] = tp;
    bot[top[u]] = u;
    if (son[u]) dfs2(son[u], tp);
    for (int e = fst[u]; e; e = nxt[e]) {
        int v = to[e];
        if (v == fa[u] || v == son[u]) continue;
        dfs2(v, v); 
    }
}
#define lch (rt << 1)
#define rch (lch | 1)
#define larg lch, l, mid
#define rarg rch, mid + 1, r
void pushup(int rt) {
    val[rt] = val[lch] * val[rch];
}
void build(int rt, int l, int r) {
    if (l == r) return;
    int mid = (l + r) >> 1;
    build(larg); build(rarg);
    pushup(rt);
}
mat query(int rt, int l, int r, int ql, int qr) {
    if (ql <= l && r <= qr) return val[rt];
    int mid = (l + r) >> 1;
    if (qr <= mid) return query(larg, ql, qr);
    if (ql > mid) return query(rarg, ql, qr);
    return query(larg, ql, qr) * query(rarg, ql, qr);
}
void update(int rt, int l, int r, int p, const mat &v) {
    if (l == r) { val[rt] = v; return; }
    int mid = (l + r) >> 1;
    if (p <= mid) update(larg, p, v);
    else update(rarg, p, v);
    pushup(rt);
}
void updateu(int u) {
    mat m; m[0][0] = m[0][1] = ldp[u][0];
    m[1][0] = ldp[u][1]; m[1][1] = -INF;
    update(1, 1, n, id[u], m);
}
void dfs3(int u) {
    ldp[u][0] = 0, ldp[u][1] = a[u];
    for (int e = fst[u]; e; e = nxt[e]) {
        int v = to[e];
        if (v == fa[u]) continue;
        dfs3(v);
        if (v == son[u]) continue;
        ldp[u][0] += max(dp[v][0], dp[v][1]);
        ldp[u][1] += dp[v][0];
    }
    updateu(u);
    vec res = query(1, 1, n, id[u], id[bot[top[u]]]) * eye;
    dp[u][0] = res[0]; dp[u][1] = res[1];
}
void modify(int u, int p) {
    ldp[u][1] += p - a[u]; a[u] = p;
    updateu(u);
    while (true) {
        u = top[u];
        vec res = query(1, 1, n, id[u], id[bot[u]]) * eye;
        if (u != 1) {
            ldp[fa[u]][0] += max(res[0], res[1]) - max(dp[u][0], dp[u][1]);
            ldp[fa[u]][1] += res[0] - dp[u][0];
        }
        dp[u][0] = res[0]; dp[u][1] = res[1];
        if (u == 1) break;
        updateu(u = fa[u]);
    }
}
int main() {
    int q; scanf("%d%d", &n, &q);
    for (int i = 1; i <= n; i++) scanf("%d", &a[i]);
    for (int i = 1; i < n; i++) {
        int x, y; scanf("%d%d", &x, &y);
        link(x, y); link(y, x);
    }
    dfs1(1, 0); dfs2(1, 1); build(1, 1, n); dfs3(1);
    while (q--) {
        int x, y; scanf("%d%d", &x, &y);
        modify(x, y);
        printf("%d\n", max(dp[1][0], dp[1][1]));
    }
    return 0;
}

本地稿整理,原稿供班级公众号,写于2019-09-12,2019-12-08上传

最近数学课上在讲数列。众所周知,我们主要探寻数列的两个属性:

  1. 通项公式
  2. n项和

关于第二点,我们已经在课上学习了等差数列和等比数列前n项和的计算方法,但是更为复杂的数列呢?例如数列\{a_n\},a_n=n^2,n\in\mathbb{N}^{\ast}的前n项和:

S_n = \sum_{x=1}^n x^2
或是更复杂一点,a_n=2^nn^2,n\in\mathbb{N}^{\ast}
S_n = \sum_{x=1}^{n}2^xx^2
或者更复杂呢?

诚然,你可能知道平方和的公式,对于第二个式子也许也有巧办法解决。但是很多方法仍存在不少局限性。有没有更通用的办法呢?

如果你自学过一点基础的微积分,你会不会想到,既然

\int_{1}^n2^xx^2\mathrm{d}x
可以通过分部积分比较套路地演算,变成数列之后有没有类似的方法呢?

这篇文章想分享的,就是在通过类比微积分构造出的这样一种系统化的求和方法。

下面的文章如果自学过微积分的话可能会比较好读,没学过也没关系,反正本身不用到任何微积分概念。

这个方法是我在《具体数学》里面看到的,也似乎只有这一本书里讲到这个方法,网上都没有。但是真的很有效,也很成系统,很震撼。这里对记号和名词做了一点小修改并避免了一些概念。

以下统一把数列看做函数,然后因为下面有非常多的函数,一定会牵涉到定义域的问题。但是好在这些函数的定义域都是平凡且显然的,所以就不写了。

基础记号与理念

定义f(x)对于x差分为且记作:

\frac{\mathrm{\delta} f(x)}{\mathrm{\delta} x} = f(x+1)-f(x)
易证差分对加减法有分配率:
\frac{\mathrm{\delta} \left[f(x) \pm g(x)\right]}{\mathrm{\delta} x} = \frac{\mathrm{\delta} f(x)}{\mathrm{\delta} x} \pm \frac{\mathrm{\delta} g(x)}{\mathrm{\delta} x}
且易证其线性:
\frac{\mathrm{\delta} kf(x)}{\mathrm{\delta} x} = k\frac{\mathrm{\delta} f(x)}{\mathrm{\delta} x}
又易证常数的差分为0(显然地)。

F(x)的差分为f(x),我们观察到对于a,b\in\mathbb{Z}, b > a,有:

F(b)-F(a) = \sum_{x=a}^{b-1} f(x)
那么,一个求和问题就可以被规约为求同一函数两个函数值差的问题。

因此,在这个框架下,如果要计算\sum_{x=a}^b f(x),我们就要:

  1. 找到合适的F(x),其差分为f(x)
  2. 计算F(b+1)-F(a)

(前几天课上讲的“累减法”是以上步骤的反向简化版)

当中最核心,也是最有意思的一个步骤,自然就是寻找F(x)。寻找F(x)的过程本质上是差分的逆运算。不妨把这个逆运算称作”逆差分“。因为使用F(x)易产生指代不清的问题,我们为f(x)的逆差分定义如下记号:

F(x) = \sum f(x)\mathrm{\delta} x
则由定义,可得:
\frac{\mathrm{\delta} \left(\sum f(x) \mathrm{\delta} x\right)}{\mathrm{\delta} x} = f(x), \sum \frac{\mathrm{\delta} f(x)}{\mathrm{\delta} x} \mathrm{\delta} x = f(x)
举些例子:
\begin{aligned}
\frac{\mathrm{\delta} x}{\mathrm{\delta} x} = 1 &\Leftrightarrow \sum1\mathrm{\delta} x = x + C\\
\frac{\mathrm{\delta} x^2}{\mathrm{\delta} x} = 2x+1 &\Leftrightarrow \sum (2x+1)\mathrm{\delta} x = x^2 + C
\end{aligned}
为什么要+C呢?因为f(x)的逆差分不唯一:对于任意已知的f(x)逆差分加上任意的常数依然是f(x)的逆差分(差分的时候常数就抵消掉了)。这里的C便是这个意思。

逆差分由于是差分的逆运算,其对于加法的分配率以及线性亦显然。

接下来的这篇文章的中心便是:如何计算\sum f(x) \mathrm{\delta} x

以下默认x,n \in \mathbb Z

“幂”法则

对微积分有些许了解的同学一定知道积分当中的幂法则:

\int x^m \mathrm{d}x=\frac{x^{m+1}}{m+1}+C
它来源于微分的幂法则:
\frac{\mathrm{d}x^m}{\mathrm{d}x} = mx^{m-1}
类似的幂法则在我们刚刚建立的框架当中的类比是什么呢?

直接把微分换成差分?

\frac{\mathrm{\delta} x^m}{\mathrm{\delta} x} = (x+1)^m-x^m = ?
我们发现好像没有什么优雅的结果诶!

会不会是"幂"出了问题呢?

定义下降阶乘幂x^{\underline{m}}

x^{\underline m} = x(x-1)(x-2)\cdots(x-m+1) = \frac{x!}{(x - m)!} , \quad x>0,m\in \mathbb{Z}
(当然,最右边那个有阶乘的式子要有意义,应该x\ge m,如果x<m,那么不妨令整个阶乘幂为0

为什么要定义那么奇怪的幂呢?考查下降阶乘幂的差分:

\begin{aligned}
\frac{\mathrm{\delta} x^{\underline m}}{\mathrm{\delta} x} &= (x + 1)^{\underline m} - x^{\underline m} \\
&= \frac{(x + 1)!}{(x + 1 - m)!} - \frac{x!}{(x - m)!} \\
&= \frac{(x + 1)! - (x + 1 - m)x!}{(x + 1 - m)!} \\
&= \frac{x![(x + 1) - (x + 1 - m)]}{(x + 1 - m)!} \\
&= \frac{mx!}{[x - (m - 1)]!} \\
&= mx^{\underline{m - 1}}
\end{aligned}
我们发现在使用下降阶乘幂的情况下,我们就可以把微分当中的幂法则类比到差分当中了。

同时,我们不妨探索一下常规幂的性质在下降阶乘幂中的类比:

x^{\underline{m + n}} = x^{\underline m}(x - m)^{\underline n}, \quad m,n\in \mathbb{Z}
以及当m<0时:
x^{\underline{-m}} = \frac{1}{(x + 1)(x + 2)\cdots(x + m)}
这些类比得来的性质有助于我们加深对于下降阶乘幂的认识。

既然我们得到了差分的幂法则,逆差分中的幂法则也就呼之欲出了:

\sum x^{\underline m} \mathrm{\delta} x = \frac{x^{\underline{m+1}}}{m+1}+C,\quad m\in\mathbb{Z}, m\neq -1
这些法则有什么用呢?我们先试试看用我们的方法推导等差数列求和公式吧!

欲求式为

S_n = \sum_{x=1}^n [a_1+(x-1)d]
那么我们计算:
\begin{aligned}
    \sum [a_1+(x-1)d]\mathrm{\delta} x &= \sum[a_1-d+dx]\mathrm{\delta} x \\
    &= \sum(a_1 - d)\mathrm{\delta} x+\sum dx\mathrm{\delta} x \\
    &= (a_1-d)x + d\sum x^{\underline 1} \delta x \\
    &= (a_1-d)x + d\cdot\frac{x^{\underline 2}}{2} + C
\end{aligned}
然后根据之前提到的步骤分别代入x=1x=n+1作差相减,就是答案:
\begin{aligned}
\sum_{x=1}^n [a_1+(x-1)d] &= \left[(a_1-d)(n+1)+d\cdot\frac{(n+1)^{\underline 2}}{2}\right] - \left[(a_1 - d)\cdot 1 + d\cdot\frac{1^{\underline 2}}{2} \right] \\
&= (a_1-d)n + \frac{n(n+1)}{2}d \\
&= na_1+\frac{n(n-1)}{2}d
\end{aligned}
是不是很简单?再比如说我们要求平方和
\sum_{x=1}^n x^2
然后我们发现
x^2 = x^{\underline 2} + x^{\underline 1}
那我们可以借此求取x^2的逆差分:
\begin{aligned}
\sum x^2 \mathrm{\delta} x &= \sum x^{\underline 2}\mathrm{\delta} x + \sum x^{\underline 1}\mathrm{\delta} x \\
&= \frac{x^{\underline 3}}{3} + \frac{x^{\underline 2}}{2} + C \\
&= \frac{2x(x - 1)(x - 2) + 3x(x - 1)}{6} + C \\
&= \frac{x(x - 1)(2x - 1)}{6} + C
\end{aligned}

代入上下界:

\begin{aligned}
\sum_{x=1}^n x^2 &= \frac{(n+1)[(n+1) - 1][2(n+1) - 1]}{6} - \frac{1(1- 1)(2\cdot1 - 1)}{6} \\
&= \frac{n(n+1)(2n+1)}{6}
\end{aligned}
我们一不小心就推出了平方和公式?

再来一个:

\sum_{x=1}^n \frac{1}{x(x+1)}
我们先把它转变成:
\sum_{x=0}^{n-1}\frac{1}{(x+1)(x+2)}
因为这样凑到了形式:
\frac{1}{(x+1)(x+2)} = x^{\underline{-2}}
然后用幂法则逆差分:
\sum x^{\underline{-2}} \mathrm{\delta} x = \frac{x^{\underline{-1}}}{-1} + C= -\frac{1}{x+1} + C
然后代入x=0x=n相减:
\begin{aligned}
\sum_{x=0}^{n-1}\frac{1}{(x+1)(x+2)} &= \left(-\frac{1}{n+1}\right) - \left(-\frac{1}{0+1}\right) \\
&= 1-\frac{1}{n+1}
\end{aligned}
这不就是裂项吗?

其实,通过归纳可以证明,任何整式多项式都可以被表示为若干个下降阶乘幂的线性组合(且方法是显然的)。且裂项公式的另类证明也说明幂法则在指数为负数的时候也可以使用。(对于复杂的分式,我们也可以待定系数拆着做)。

指数法则

自然而然地,我们接下来探究指数函数(包括底数为负数)的差分与逆差分:

\frac{\mathrm{\delta} a^{x}}{\mathrm{\delta} x} = a^{x+1}-a^x=(a-1)a^x
因此
\sum a^x\mathrm{\delta} x=\frac{a^x}{a-1} + C ,\quad a\neq 1
(要验证?再差分一下)

因此

\sum_{x=0}^{n-1}a^x = \frac{a^n}{a-1} - \frac{a^0}{a-1} = \frac{a^n-1}{a-1} ,\quad a\neq 1
啊,是非平凡等比数列的求和公式呢。

刚刚的结论也可以稍作拓展:

\begin{aligned}
\frac{\mathrm{\delta} a^{x+k}}{\mathrm{\delta} x} &= (a-1)a^{x+k}\\ \sum a^{x+k}\mathrm{\delta} x&=\frac{a^{x+k}}{a-1} + C
\end{aligned}
同时可以注意下:
\frac{\mathrm{\delta} 2^x}{\mathrm{\delta} x} = 2^x, \sum2^x\mathrm{\delta} x=2^x+C
这是一个有意思的现象:2^x(或更为广义来讲,2^{x+k})的差分与逆差分等于其自身。这与微积分中e^x的地位是相似的。

分部求和

如果仅有指数法则和幂法则,那么这套方法其实并没有那么优秀(计算的式子我们常规方法也能算啊!)

真正让这套方法威力大幅提升的,是分部求和。

你也许知道分部积分是从微分的积法则变化而来的:

\begin {aligned}
&\frac{\mathrm{d}f(x)g(x)}{\mathrm{d}x} = f(x)\frac{\mathrm{d}g(x)}{\mathrm{d}x} + g(x)\frac{\mathrm{d}f(x)}{\mathrm{d}x} \\
\Longleftrightarrow& f(x)\frac{\mathrm{d}g(x)}{\mathrm{d}x} = \frac{\mathrm{d}f(x)g(x)}{\mathrm{d}x}- g(x)\frac{\mathrm{d}f(x)}{\mathrm{d}x}\\
\Longleftrightarrow  &\int f(x)\mathrm{d}g(x) = f(x)g(x) - \int g(x)\mathrm{d}f(x) \\
\end{aligned}
(出于统一风格的需要没有用u,v,写得繁复了一点)

我们类比着探寻差分下的积法则:

\begin{aligned}
    \frac{\mathrm{\delta} f(x)g(x)}{\mathrm{\delta} x} &= f(x+1)g(x+1) - f(x)g(x) \\
    &= f(x+1)g(x+1)-f(x)g(x+1)+f(x)g(x+1)-f(x)g(x) \\
    &= f(x)[g(x+1) - g(x)] + g(x+1)[f(x+1)-f(x)] \\
    &= f(x)\frac{\mathrm{\delta} g(x)}{\mathrm{\delta} x} + g(x+1)\frac{\mathrm{\delta} f(x)}{\mathrm{\delta} x}
\end{aligned}
因此
\sum f(x) \mathrm{\delta} g(x) = f(x)g(x)-\sum g(x+1)\mathrm{\delta} f(x)
(\sum f(x) \mathrm{\delta} g(x) = \sum f(x)\frac{\mathrm{\delta} g(x)}{\cancel{\mathrm{\delta} x}} \cancel{\mathrm{\delta} x},虽然不是直接“约分”,但是是一种简便的记号)

这便是逆差分当中的分部求和方法。

怎么用?

比如说我们要计算:

\sum x^22^x \mathrm{\delta} x
那么令f(x)=x^2,\frac{\mathrm{\delta} g(x)}{\mathrm{\delta} x} = 2^x \Rightarrow g(x)=2^x:
\begin{aligned}
\sum x^22^x \mathrm{\delta} x = \sum x^2 \mathrm{\delta}2^x&=x^22^x- \sum 2^{x+1}\mathrm{\delta} x^2 \\
&= x^22^x- \sum 2^{x+1}(2x+1)\mathrm{\delta} x \\
&= x^22^x - \sum 2^{x+1}\mathrm{\delta} x - 2\sum x2^{x+1}\mathrm{\delta} x \\
&= x^22^x - 2^{x+1} - 4\sum x2^x\mathrm{\delta} x
\end{aligned}
再来,令f(x) = x,\frac{\mathrm{\delta} g(x)}{\mathrm{\delta} x} = 2^x \Rightarrow g(x)=2^x
\begin{aligned}
\sum x2^x\mathrm{\delta} x &= x2^x- \sum 2^{x+1} \mathrm{\delta} x \\
&= x2^x - 2^{x+1} + C \\
\end{aligned}
代回去:
\begin{aligned}
\sum x^22^x \mathrm{\delta} x &= x^22^x - 2^{x+1} - 4\left(x2^x - 2^{x+1}\right) + C \\
&=(x^2-4x+6)2^x+C
\end{aligned}
然后上下界就按需求随意加好了~

g(x)=2^x+C的时候推导还成立吗?要不要考虑呢?)

换元求和

你也许会想

既然我们都把分部积分类比到了分部求和,那么再把换元积分类比到换元求和岂不美哉?

很遗憾,这在通常意义上是不可以的。

为什么?

因为换元积分逆用的是微分的链式法则。

差分有链式法则吗?

\begin{aligned}
    \frac{\mathrm{\delta} f[g(x)]}{\mathrm{\delta} x} &= f[g(x+1)] - f[g(x)] \\
    &\stackrel{?}{=} \left\{f[g(x)+1]-f[g(x)]\right\}[g(x+1)-g(x)] \\
    &= \frac{\mathrm{\delta} f[g(x)]}{\mathrm{\delta} g(x)} \cdot \frac{\mathrm{\delta} g(x)}{\mathrm{\delta} x}
\end{aligned}
打问号的那一步是推不过来的,除非g(x+1)=g(x)+1 \Longleftrightarrow \frac{\mathrm{\delta} g(x)}{\mathrm{\delta} x} = 1。在这种特例条件下,g(x)=x+C,也就是说我们只能在这样形式的g(x)上应用链式法则/换元求和——这个形式不就是平移吗?

而在其他形式的g(x)当中,正如我们所看到的一样,推不过去。

其实,我们没有必要将所有微积分里面的要素全部类比到我们的体系里面。比如说,微分的链式法则与商法则之所以必要是因为微分按照定义计算十分繁琐,而差分的定义本身就可以用来进行计算,因而不需要链式法则与商法则。类似地,微积分里面一些关于三角的技巧也不用,也不能类比到我们的框架当中,因为在离散的情况很难对待三角。

小结

这篇文章给大家分享了一整套可以用来机械化一部分求和计算的方法,并在此方法框架内推导了我们所学过的许多求和公式。这个方法精彩且难能可贵的是,这个方法完全是借助于类比的思想借鉴微积分体系建立起来的。因而如果同学们有关于微积分的了解,那么会发现知识的迁移可以很快完成。这种类比的方法是我们数学学习过程中应该提倡的。

这套方法有名字吗?

《具体数学》当中称该方法为“有限微积分”,我个人认为“离散微积分”更为恰当,但是既然是离散的,那么“微”和“积”自然无从说起,我个人一直想回避给这个方法命名(就是一个微积分的类比又不是什么开拓性创新有啥好命名的)。

考试/作业可以用吗?

我觉得肯定是不行的。首先我确信考试或者作业大概率不会出这种求和。其次我确定很多人并没有听说过这个方法。再其次这套方法和其他技巧还不大一样,是成体系的,所以如果你要用也要像这篇文章一样从头开始定义差分写一大堆。然而如果真的碰到了这样的题目,你当然……可以……在草稿纸上……假装是灵光一闪得到答案……然后归纳证明对吧……我觉得老师大概也不会拿你怎么样?

调和级数与自然对数

这里只是对于幂法则的一个简要补充,可以跳过。

你也许注意到,正如积分的幂法则一样,我们逆差分的幂法则规定指数不能为-1。原因显然。

那指数-1怎么办呢?

你也许知道

\int_1^n x^{-1}\mathrm{d}x=\ln n
那类比着不妨看看
\sum_{x=0}^{n-1} x^{\underline{-1}} = \frac{1}{1}+\frac{1}{2}+\cdots+\frac{1}{n}=?
这不是调和级数吗?这和\ln n 看起来一点都不像啊!然而,可以证明,
\lim_{n\to\infty} \left(- \ln n + \sum_{x=1}^n \frac{1}{x}\right) = \gamma \approx 0.57722
\gamma被称作欧拉-马歇罗尼常数。当然这和这篇文章的主题已经一点关系都没有了,权当是对于指数为-1情况的探究与补全吧。

Bonus

暑假里刷计算机竞赛的时候小马碰到一道神奇的数列题。当时所有人都是打表找规律做的(计算机竞赛日常,你写出程序就行),但是我一直很好奇有没有严谨的推导。数列是这个样子的:

\begin{aligned}
    a_1 &= 1 \\
    a_n &= \left(\sum_{i=1}^{n-1} i\cdot a_i\right) \bmod n ,\quad n\ge2, n\in \mathbb{N^{\ast}}
\end{aligned}
a_n的通项公式。我趁这篇文章的机会问问各位大佬有没有想法,有的话可以和我分享一下。

回顾

微积分是数学分析中一类极为有力的工具,其中微分的核心是微分算子D,定义为:

Df(x) = \lim_{h \to 0} \frac{f(x + h) - f(x)}{h}
DF(x) = f(x)我们称F(x)f(x)的原函数,他们两个一起构成了积分的核心——牛顿-莱布尼茨公式:
\int_a^b f(x)dx = F(b) - F(a)
在微积分当中我们有一整套系统的规则来方便地计算各种形式的微分和积分,而依稀记得一句话:积分是黎曼和的极限,即积分是求和的连续情况,求和是积分的离散情况,既然我们可以系统地计算积分,那么我们可不可以系统地计算和式呢?答案是肯定的,《具体数学》中为我们引入的一类工具便是——有限微积分。(个人觉得叫离散微积分更好?)

定义和记号

仿照微分算子D的定义我们首先定义差分算子D

\Delta f(x) = f(x + 1) - f(x)
可以发现,差分算子本质是微分算子定义中h = 1的特例,只不过在微分中h无限接近于0,而在差分当中最小只能为1。

类似地,若\Delta F(x) = f(x)我们也可以称F(x)f(x)的“原函数”,显然我们有:

\sum_a^b f(x) \delta x  = F(b) - F(a)
我们在这里仿照定积分引入了一个新的记号\sum_a^b f(x)\delta x,不难发现:
\sum_a^b f(x) \delta x = \sum_{x = a}^{b - 1}f(x)
类似地,仿照不定积分的定义,我们还可以定义一个更为抽象的概念”不定和式“,即
\sum f(x) \delta x = F(x) + C
C也起类似积分常数的作用。

注入灵魂

照葫芦画瓢到此为止,然而如果只有这么些记号是没有卵用的,就像你不能只靠微分的极限定义和莱布尼茨公式计算微积分一样,我们也不能指望靠这些记号就可以快速计算和式了,我们还需要类似于微积分当中的幂法则等一系列规则的支持。

首先不证自明的一点是:在有限微积分中,差分算子对于加法和减法的分配率成立

幂法则

我们先考察微积分当中的幂法则在这里有没有等价的替代品,众所周知,幂法则是:

Dx^m = mx^{m-1}
如果简单地把微分算子换成差分我们会得到:
\Delta x^m = (x + 1)^m - x^m
接下来的展开势必就要用到二项式定理了,事到如今只要是正常人都知道最后是得不到$ mx^{m-1} $的,我们发现直接把幂法则照搬过来是行不通的,如果问题不是出在差分算子上那就是出在“幂”上了,那么现在问题变成了:我们能不能找到一种“幂”的定义,使它在差分算子作用下表现出类似幂法则的性质呢?

答案是肯定的,这种“幂”被称作下降阶乘幂,记作x^{\underline m},《具体数学》中对它的定义是:

x^{\underline m} = x(x - 1)(x - 2)\cdots(x - m + 1)
显然无论是从形式还是从名字来看它都有一个更简洁的阶乘定义:
x^{\underline m} = \frac{x!}{(x - m)!}
在《具体数学》中作者花了不少笔墨来把它最上面的定义推广到m < 0的情形,但是一切在阶乘定义中就很显然了,对于m > 0,我们有:
x^{\underline{-m}} = \frac{1}{(x + 1)(x + 2)\cdots(x + m)}
我们也有类似x^{m + n} = x^mx^n的规则:
x^{\underline{m + n}} = x^{\underline m}(x - m)^{\underline n}
好的,让我们回到幂法则,那么如果我们结合下降阶乘幂和差分算子,我们会得到什么呢?
\begin{aligned}
\Delta x^{\underline m} &= (x + 1)^{\underline m} - x^{\underline m} \\
&= \frac{(x + 1)!}{(x + 1 - m)!} - \frac{x!}{(x - m)!} \\
&= \frac{(x + 1)! - (x + 1 - m)x!}{(x + 1 - m)!} \\
&= \frac{x![(x + 1) - (x + 1 - m)]}{(x + 1 - m)!} \\
&= \frac{mx!}{[x - (m - 1)]!} \\
&= mx^{\underline{m - 1}}
\end{aligned}
我们成功地推出了“幂法则”!

顺带一提,除了下降阶乘幂之外还有上升阶乘幂,定义为x^{\overline m} = \frac{(x + m - 1)!}{(x - 1)!},也有类似幂法则的性质,但是用起来比较麻烦,故我们使用下降阶乘幂。

既然我们已经成功导出了有限微积分中差分算子的幂法则,那么对于不定和式,其幂法则也呼之欲出:

\sum x^{\underline m}\delta x = \frac{x^{\underline{m + 1}}}{m + 1} + C
正如积分中对于幂法则在m = -1的时候需要特殊处理,即\int x^{-1} dx = \ln x + C,在求和当中,我们也有类似的:
\sum x^{\underline{-1}}\delta x = H_x + C
其中H_x表示调和级数的第x项部分和,我们知道当x \to \infty时有H_x = \ln x + \gamma,其中\gamma是欧拉-马歇罗尼常数。从此我们也可以看出上述规则确实是积分中对应公式的一个离散模拟。

那么这些法则有什么用呢?我们可以再试试看求解如下和式:

\sum_{x = 1}^n x^2
我在先前成套方法的介绍中曾利用成套方法对它进行过求解,过程相对来说比较繁琐,我们看看怎么用有限微积分秒杀这道题。

首先,注意到:

x^2 = x^{\underline 2} + x^{\underline 1}
所以:
\begin{aligned}
\sum x^2 \delta x &= \sum x^{\underline 2}\delta x + \sum x^{\underline 1}\delta x \\
&= \frac{x^{\underline 3}}{3} + \frac{x^{\underline 2}}{2} + C \\
&= \frac{2x(x - 1)(x - 2) + 3x(x - 1)}{6} + C \\
&= \frac{x(x - 1)(2x - 1)}{6} + C
\end{aligned}
接下来就是展现真正技术的时刻:
\begin{aligned}
\sum_{x = 1}^n x^2 &= \sum_1^{n + 1} x^2 \delta x \\
&= \frac{n(n + 1)(2n + 1)}{6} + C - \frac{1 \cdot 0 \cdot 1}{6} - C \\
&= \frac{n(n + 1)(2n + 1)}{6}
\end{aligned}
就非常舒服。通过归纳法可以证明,对于任何多项式,我们总是可以把它转换为若干个下降阶乘幂的线性组合,因此,对于一切多项式的和式,我们都能用有限微积分系统地求取!

运用指数为-2的幂法则,你就会发现你得到了所谓的裂项公式。

e^x的对等

在微积分中,e^x的具有非常好的性质,即导数和积分都等于自身,那么在有限微积分当中有没有类似的函数呢,幸运的是,找到这样的函数并不困难:

\begin{aligned}
\Delta f(x) &= f(x) \\ 
f(x + 1) - f(x) &= f(x) \\
f(x + 1) &= 2f(x)
\end{aligned}
这就很明显了,f(x) = 2^x,或者对其平移的结果都满足差分等于自身的性质。

因此,无惑乎我们有:

\begin{aligned}
\sum_{x = 1}^n 2^x &= \sum_1^{n + 1} 2^x \delta x \\
&= 2^{n + 1} - 2^0 \\
&= 2^{n + 1} - 1
\end{aligned}

指数函数

我们知道在微积分当中:

\int a^xdx = \frac{a^x}{\ln a} + C
我们接下来探究在有限微积分当中是否有类似结论的出现,我们先猜想:
\sum a^x\delta x = ka^x + C
两边取差分:
\begin{aligned}
a^x &= ka^{x + 1} - ka^x \\
a^x &= ka^x(a - 1) \\
k &= \frac{1}{a - 1}
\end{aligned}
因此:
\sum a^x\delta x = \frac{a^x}{a - 1} + C
使用这一规则可以直接推出等比数列求和公式

其他法则?

我们要不要把链式法则,积法则,商法则之类的全部搬到有限微积分里面了呢?其实是不用的,因为这些法则都是辅助我们计算微分的,其存在的原因便是使用微分的极限定义计算导数实在是太困难了,而差分算子则不然,计算一个函数的差分只要代进去就行了。

但是,为了后文的方便,在此我们着重来讨论积法则在有限微积分当中的对等,在微积分当中,积法则是:

D(uv) = uDv + vDu
那么在有限微积分当中:
\begin{aligned}
\Delta(u(x)v(x)) &= u(x + 1)v(x + 1) - u(x)v(x) \\
&= u(x + 1)v(x + 1) - u(x)v(x + 1) + u(x)v(x + 1) - u(x)v(x) \\
&= u(x)\Delta v(x) + v(x + 1)\Delta u(x)
\end{aligned}
记移进算子Ef(x) = f(x + 1),则有限微积分的积法则可以被表示为:
\begin{aligned}
\Delta(uv) &= u\Delta v + Ev\Delta u \\
&= Eu\Delta v + v\Delta u
\end{aligned}

分部求和

对于积分,我们有分部积分公式:

\int udv = uv - \int vdu
其本质是对微分积法则移项并两边积分得到的成果,我们能不能也用上面推导出的有限微积分的“积法则”推导出分部求和呢?
\begin{aligned}
\Delta(uv) &= u\Delta v + Ev\Delta u \\
\Rightarrow u\Delta v &= \Delta(uv) - Ev\Delta u \\
\Rightarrow \sum u \Delta v &= uv - \sum Ev\Delta u
\end{aligned}
最后的结果就是分部求和的公式,分部积分的经典例题是\int xe^xdx,对应到有限微积分当中就是\sum x2^x\delta x,运用我们刚得到的分部求和公式,令u = x, \Delta u = 1, v = 2^x, \Delta v = 2^x我们就得到了:
\begin{aligned}
\sum x2^x\delta x &= x2^x - \sum 2^{x + 1}\delta x + C\\
&=x2^x - 2^{x + 1} + C
\end{aligned}
因此
\begin{aligned}
\sum_{i = 0}^nx2^x &= \sum_0^{n + 1} x2^x\delta x \\
&= (n + 1)2^{n + 1} - 2^{n + 2} - 0\cdot 1 + 2 \\
&= (n - 1)2^{n + 1} + 2
\end{aligned}
计算这个和式甚至不需要思考!

我们来试试上次成套方法博客中我说需要一定技巧的\sum_{x = 1}^n (-1)^xx^2吧!

u = x^2, \Delta v = (-1)^x可得\Delta u = 2x + 1, v = -\frac{ (-1)^x}{2}

\begin{aligned}
\sum (-1)^xx^2\delta x &= -\frac{x^2(-1)^x}{2} - \sum-\frac{ (-1)^{x + 1}}{2}(2x + 1)\delta x + C\\
&= -\frac{x^2(-1)^x}{2} - \sum\frac{(-1)^x(2x + 1)}{2}\delta x + C\\
&= -\frac{x^2(-1)^x}{2} - \sum\frac{(2x + 1)}{2}(-1)^x\delta x + C
\end{aligned}
观察这个式子的后项,我们故技重施,令u = \frac{(2x + 1)}{2}, \Delta u = 1, \Delta v = (-1)^x
\begin{aligned}
\sum\frac{(2x + 1)}{2}(-1)^x\delta x &= -\frac{(2x + 1)(-1)^x}{4} - \sum-\frac{ (-1)^{x + 1}}{2}\delta x + C\\
&= -\frac{(2x + 1)(-1)^x}{4} - \frac{1}{2}\sum(-1)^x\delta x + C\\
&= -\frac{(2x + 1)(-1)^x}{4} + \frac{1}{4}(-1)^x + C\\
&= -\frac{(-1)^xx}{2} + C
\end{aligned}
回到上式:
\begin{aligned}
\sum (-1)^xx^2\delta x &= -\frac{x^2(-1)^x}{2} + \frac{(-1)^xx}{2} + C\\
&= -\frac{(-1)^xx(x - 1)}{2} + C
\end{aligned}
给这个和式加上下界:
\begin{aligned}
\sum_{x = 1}^n (-1)^xx^2 &= -\frac{(-1)^{n + 1}n(n + 1)}{2} + -\frac{(-1)^1\cdot 0 \cdot 1}{2} \\
&= \frac{(-1)^nn(n + 1)}{2}
\end{aligned}
以上的过程如行云流水一般地系统化和顺畅,只要熟悉我们拆分出指数项,然后对多项式项不断使用差分降次——这其实是分部积分的技巧,但是得益于有限微积分和微积分的联系我们得以直接拿来用了。

总结

综上所述,基于整数建立起来的有限微积分虽然不一定有传统的微积分那么好的性质,但是依然可以系统地高效地处理很多的求和问题,我们在有限微积分的公式推导中还顺便导出了平方和公式,裂项公式和等比数列公式等,而这一切公式的导出在有限微积分的背景下是那么地自然,以至于让人无法相信第一个属于小学奥数,第二个属于初中数学,而最后一个要在高中才能教到。虽然这也和我们独到的证明方法有关(一般来说这三个公式分别需要用到归纳法,叠缩和扰动法证明),但是可以预见的是,有限微积分的确在和式的处理上具有提纲掣领的作用,不容小觑。

本地笔记整理重发,原文写于2019/8/28,2019/12/8上传

分析

首先使用哈希计算字符串两两之间的overlap,这个的时间复杂度不是特别显然:

\sum_{i=1}^n\sum_{j=1}^n\min\{a_i,a_j\}
其中a_i, a_j分别是第i个和第j个字符串的长度。假设a_i从小到大排了序,那么就变为:
\sum_{i=1}^n \left[\sum_{j=1}^ia_j+(n-i)a_i\right]
这样每一个a_i都会被计算n-i+1+n-i=2n-2i+1次,则时间复杂度为:
\begin{aligned}
\sum_{i=1}^n (2n-2i+1)a_i
\end{aligned}
显然这个时候越前面的a_i越大越好,但是a_i是递增的,所以极限情况是各a_i相等取最大值,时间复杂度就是O\left(n\sum_{i=1}^na_i\right)

然后Floyd 矩阵乘法即可。

双哈希大法好

卡常Nice

#include <cstdio>
#include <algorithm>
#include <vector>
#include <cstring>

using namespace std;
const int MOD1 = 998244353, MOD2 = 1e9 + 7;
const int EPS = 26;
const int N = 210;
const int S = 100010;
typedef long long ll;
typedef pair<int, int> hsh;
int n, m;
vector<hsh> h[N];
hsh p[S];
inline hsh operator +(const hsh &a, const hsh &b) {
    hsh ret = make_pair(a.first + b.first, a.second + b.second);
    if (ret.first >= MOD1) ret.first -= MOD1;
    if (ret.second >= MOD2) ret.second -= MOD2;
    return ret;
}
inline hsh operator -(const hsh &a, const hsh &b) {
    hsh ret = make_pair(a.first - b.first, a.second - b.second);
    if (ret.first < 0) ret.first += MOD1;
    if (ret.second < 0) ret.second += MOD2;
    return ret;
}
inline hsh operator +(const hsh &a, char b) {
    int c = b - 'a';
    return make_pair(((ll)a.first * EPS + c) % MOD1, 
            ((ll)a.second * EPS + c) % MOD2);
}
inline hsh operator *(const hsh &a, const hsh &b) {
    return make_pair((ll)a.first * b.first % MOD1, 
            (ll)a.second * b.second % MOD2);
}
inline hsh cut(int i, int l, int r) {
    return h[i][r] - h[i][l] * p[r - l];    
}
inline int intersect(int a, int b) {
    int la = h[a].size() - 1, lb = h[b].size() - 1;
    int lim = a == b ? la - 1 : min(la, lb);
    int ret = 0;
    for (int i = 1; i <= lim; i++)
        if (cut(a, la - i, la) == cut(b, 0, i)) ret = i;
    return ret;
}
struct floyd {
    ll d[N][N];
    floyd() { memset(d, 0x3f, sizeof d); }
    ll *operator [](int x) { return d[x]; }
    const ll *operator [](int x) const { return d[x]; }
    floyd operator *(const floyd &b) const {
        floyd ret;
        for (int k = 1; k <= n; k++)
            for (int i = 1; i <= n; i++)
                for (int j = 1; j <= n; j++) 
                    ret[i][j] = min(ret[i][j], d[i][k] + b[k][j]);
        return ret;
    }
} f;
int main() {
    scanf("%d%d", &n, &m);
    int maxlen = 0;
    for (int i = 1; i <= n; i++) {
        static char tmp[S]; scanf("%s", tmp + 1);
        int len = strlen(tmp + 1);
        maxlen = max(maxlen, len);
        h[i].push_back(make_pair(0, 0));
        for (int j = 1; j <= len; j++)
            h[i].push_back(h[i].back() + tmp[j]);
    }
    p[0] = make_pair(1, 1);
    for (int i = 1; i <= maxlen; i++) p[i] = p[i - 1] * make_pair(EPS, EPS);
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) 
            f[i][j] = h[j].size() - 1 - intersect(i, j);
    floyd ans; for (int i = 1; i <= n; i++) ans[i][i] = h[i].size() - 1;
    int y = m - 1;
    for (; y; y >>= 1) {
        if (y & 1) ans = ans * f;
        f = f * f;
    }   
    ll mn = 1ll << 60;
    for (int i = 1; i <= n; i++)
        for (int j = 1; j <= n; j++) mn = min(mn, ans[i][j]);
    printf("%lld\n", mn);
    return 0;
}

本地笔记整理,原文写于2019/7/31,2019/12/8上传

这个枚举子集啊从\mathcal O(4^n)\mathcal O(3^n)只优化了一半,最后的目标一定是\mathcal O(2^n)! ——MR

考虑加速在序列\{a\}, \{b\}上的集合并卷积:

c_k = \sum_{i \cup j = k} a_ib_j
利用和FFT类似的思想,考虑构造变换\operatorname{FMT}\{a\},使得:
\operatorname{FMT}\{a\}_i\times\operatorname{FMT}\{b\}_i = \operatorname{FMT}\{c\}_i
观察到可以构造:
\operatorname{FMT}\{a\}_i = \sum_{j \subseteq i} a_j
即前缀子集和。该变换成立的原因显然。

然而朴素构造子集和的时间复杂度是O(3^n)的,有没有办法优化呢?

显然是有的,我们采用增量的方法,逐个考虑集合的每个元素x,如果一个集合i应该包含x,则其子集可以包含也可以不包含x,因此进行操作:a_i = a_i + a_{i\setminus x}

代码是很短的:

for (int i = 1; i < (1 << n); i <<= 1) // 枚举元素
    for (int s = i << 1, l = 0; l < (1 << n); l += s) // [l, l + s) 是一段区间
        for (int k = 0; k < i; k++) // 其中前半段没有i,后半段有i
            a[l + i + k] += a[l + k]

逆向变换的代码也很简单:

for (int i = 1; i < (1 << n); i <<= 1)
    for (int s = i << 1, l = 0; l < (1 << n); l += s)
        for (int k = 0; k < i; k++) 
            a[l + i + k] -= a[l + k]

时间复杂度为

\sum_{i = 0}^{n - 1} 2^{n - i - 1}\cdot 2^i = O(n2^n)

本地笔记整理,原文写于2019/7/8,2019/12/8上传

Double Dabble算法是二进制转BCD的常用算法。

整数Double Dabble算法

(下文当中,"BCD位"指一个十进制位对应的四位二进制)

首先,如果我们要计算一个二进制数的十进制值,除了位权展开以外,可以从左往右如此迭代:

x_{n + 1} = 2x_n + \text{该二进制位}
这是非常显然的。

在Double Dabble的图示当中,左边一列存储的是BCD的结果:

Image result for double dabble algorithm

Double Dabble算法本身其实就是在这BCD的结果列进行如上的迭代操作(乘2加位)。

为了这样做,直觉上来讲,我们可以按照二进制的方法将BCD列中的结果乘2(左移),然后进行一定的矫正(之前的左移可能破坏了BCD的性质),然后在加上原数的一个二进制位(最右边移进)。

左移是平凡的,加上一位也是平凡的(因为左移之后最右边一位永远是2,不会产生进位),问题在乘2后为了维护BCD性质所进行的矫正上。

首先,我们应该注意到,忽略进位,如果一个BCD位大于10,我们需要减10来获得正确的位,而减10在四位二进制意义下就等同于加6,即110_2

Double Double算法的首先高在它甚至在乘2之前就进行了矫正。我们考虑哪些BCD位乘2之后会出现问题,答案非常平凡:满5的。而因为是在左移之前,加上的110_2也就变成11_2,即加3

这就是“满五加三”的由来。

不仅如此,Double Dabble算法高在因为这是对于左移之前的四位进行的操作,因而加3之后会进位到四位当中的最高位,其左移之后就到左边的BCD位上去了...相当于巧妙handle进位的问题。

小数Double Dabble算法

小数Double Dabble似乎没多少人提到……但是本质确实是和整数算法一致的。

我们进行小数二进制-十进制转化可以利用如下的从右往左的迭代过程:

x_{n + 1} = \frac{x_n + \text{该二进制位}}{2}
我们同样要对结果列进行这样的操作。

首先是x_n + \text{该二进制位}的部分,由于结果列现在表示小数部分,加上去的二进制列就在整个结果列的左侧(整数部分),这无疑是平凡的。

接下来是除2的部分。基于二进制的思想,我们通过右移来完成除2的操作。同时我们需要进行一定的矫正来维护BCD的性质。

我们观察到,出现问题的BCD位,一定是在右移之前是奇数的位,而奇数最低位上的1,在右移之后,就到了右边的最高位上,也就是右边的哪一位右移之后会大于等于8

而回到十进制,如果一个奇数位除以二,结果应该给右边的位加上5,而我们之前算法中右移过来的1让它加了8,所以要减去3

这就是“满八减三”的由来。

后缀自动机(Suffix Automaton)本质上就是接受字符串集合等于一个给定字符串S后缀集合的最小化DFA。SAM可以在O(|S|)的时间内构造出来。

等价关系

在下文中统一设母串为S

定义1:\operatorname{endpos}(s)表示子串sS中的结束位置集合。如S = abcababc, s = ab,则\operatorname{endpos}(s) = \{2, 5, 7\}。特别地,\operatorname{endpos}(\emptyset) = \{1,2,3, \cdots, |S|\}

引理1:两个子串u, w\operatorname{endpos}等价的,且假设|u| < |w|,则uw的一个后缀。

证明:脑补一下,不证自明。

引理2:考虑两个子串u, w,假设|u| < |w|\operatorname{endpos}(u)\operatorname{endpos}(v)有且仅有以下两种关系:

  1. \operatorname{endpos}(w) \subseteq \operatorname{endpos}(u):此时u一定是w的后缀。
  2. \operatorname{endpos}(w) \cap \operatorname{endpos}(u) = \emptyset:其他。

证明:如果\operatorname{endpos}(w) \cap \operatorname{endpos}(u) \neq \emptyset,那么它们至少都会在原串的某一点结束,结合|u| < |w|u一定就是w的后缀,如此一来w结束的地方u也一定结束,因此\operatorname{endpos}(w) \subseteq \operatorname{endpos}(u)

引理3:考虑一个\operatorname{endpos}等价类,如果将其中的串按照长度从小到大排序,则后一个串的长度一定比前一个串的长度严格大1,或者说,这个等价类的字符串取遍了某一区间内的所有长度。

证明:记这个等价类当中最短的串为u,最长的串为w。首先由引理1,一个等价类必定由w和它的一系列后缀组成。考虑长度在u, w之间的w的任意后缀x,因为xw的后缀,由引理2,\operatorname{endpos}(w) \subseteq \operatorname{endpos}(x),同时又因为ux的后缀,\operatorname{endpos}(x) \subseteq \operatorname{endpos}(u),结合\operatorname{endpos}(w) = \operatorname{endpos}(u),则xu, w\operatorname{endpos}等价。

后缀链接

\operatorname{endpos}集合描述了一个串在母串当中的出现情况。我们希望构建出来的DFA要最小化,因此我们自然希望一个节点表示一整个\operatorname{endpos}等价类即SAM当中每个节点的\operatorname{endpos}互不相同。

定义2:对于一个节点v,令\operatorname{endpos}(v)表示其表示的\operatorname{endpos}等价类。

我们又知道一个节点表示的字符串是S的某一个前缀的前几个后缀(引理3)(假设把后缀长度从大到小排序)。

定义3:既然某个节点v​表示某前缀的前几个后缀,令\operatorname{link}(v)​表示第一个不在等价类v​中的后缀所在的等价类节点,称为v​后缀链接

定义4:\min(v)表示节点v表示的最短子串,\max(v)表示节点v表示的最长子串。

观察1:结合定义3,显然\left|\min(v)\right| = \left|\max(\operatorname{link}(v))\right| + 1

定义4:设SAM初始节点为v_0,令\operatorname{link}(v_0) = \emptyset, \min(v_0) = \max(v_0) = \emptyset

引理4:节点与后缀链接构成一棵以初始节点v_0​为根的树。

证明:后缀链接\operatorname{link}(v)的表示的子串长度是严格小于v。顺着后缀链接总可以退回到表示子串长度为0v_0

观察2:由引理2,定义3,\operatorname{endpos}(v) \subsetneq \operatorname{endpos}(\operatorname{link}(v))

观察3:因此,后缀链接组成的树也完全对应\operatorname{endpos}​等价类的子集关系。

构建算法

我们采用增量式的方法构建后缀自动机,即考虑在S的SAM的基础上构建S|c的后缀自动机:

在算法中,我们对于一个节点v只存储以下信息:

  1. \operatorname{link}(v):后缀链接
  2. \operatorname{len}(v) = |\max(v)|​
  3. 转移边

算法流程如下:

  1. v_{\text{last}}表示原来表示整个S串的节点。
  2. 创建新节点v_{\text{cur}}用于表示串S|c,令\operatorname{len}(v_{\text{cur}}) = \operatorname{len}(v_{\text{last}}) + 1
  3. v_{\text{last}}​开始顺着后缀链接一路向上,对于没有c​转移边的节点添加到v_{\text{cur}}​的后缀转移。如果一路到了v_0​,那么令\operatorname{link}(v_{\text{cur}}) = v_0​,否则在第一个存在c​转移边的节点p​出停住。记p​c​转移到q​
  4. 判断:
    1. 如果\operatorname{len}(q) = \operatorname{len}(p) + 1​,则令\operatorname{link}(v_{\text{cur}}) = q​
    2. 否则,将q复制一份得到q',令\operatorname{len}(q') = \operatorname{len}(p) + 1。令\operatorname{link}(q) = q'\operatorname{link}(v_{\text{cur}}) = q'。同时再顺着p的后缀链接向上,把所有指向qc转移边重定向到q',直到到达初始节点或者遇到一个非q的转移为止。

正确性证明

首先,算法从创建一个新状态v_{\text{cur}}​开始,这个节点表示的字符串是S|c​,同时也开辟一个新的等价类。

接下来对于之前的后缀我们都要做相应的更新(加上c才算新的后缀),最简单的方法就是给这些节点加上往v_{\text{cur}}的字符c转移边,但是这些转移边必须不与之前的转移冲突才行,因此我们到节点p就停止了,再加就冲突了。(注意我们在这里每给一个节点v(设表示串为S_v)加一条c转移边,我们就把S_v|c加到了v_{\text{cur}}表示的串集合当中。这决定了v_{\text{cur}}的后缀链接应该链接到哪里。)

首先如果我们找不到节点p,也就是我们一路到了初始节点,那说明这个c就根本没有S中出现过,也说明v_{\text{cur}}一下子表示了S|c的所有非空后缀了,因此由后缀链接的定义\operatorname{link}(v_{\text{cur}}) = v_0

如果我们确实在p节点停下了,令p代表的串集为S_p,那q的串集S_q至少包含S_p|c,而此时S_p|c\operatorname{endpos}集合是被拓展的(因为现在S_p|c也会在S|c的末尾结束),因此:

  1. 如果S_q只包含S_p|c,也就是\operatorname{len}(q) = \operatorname{len}(p)+1那么我们可以直接拓展q的等价类。此时\operatorname{link}(v_{\text{cur}}) = q
  2. 如果S_q不只包含S_p|c,也就是\operatorname{len}(q) > \operatorname{len}(p) + 1的话,那么S_q中只有一部分串的\operatorname{endpos}集合被拓展了,此时我们有必要重新开一个节点来表示拓展得到的新等价类。那些比S_p|c长的,没有拓展的串仍然保留在qS_p|c被单独移到了q',由定义那么\operatorname{link}(q) = q'。最后我们需要做一些收尾工作,也就是把一些节点从转移到q变成转移到q',我们发现这些节点代表的串一定是S_p的后缀,因此只要一路p的后缀链接向上就行了。

关于时空复杂度,首先这个算法本身的保证了状态数的上界是2|S| - 1 = O(|S|)​的。

其次这个算法的转移数也是线性的,这两个性质一起确保了整个算法的均摊复杂度是线性的。具体证明可以看这里

代码实现

个人认为SAM的代码实现只要理解了SAM的原理其实是比SA简单的:

namespace sam {
    int cnt = 1, last = 1;
    const int EPSILON = 26;
    int len[N], link[N], nxt[N][EPSILON];
    
    void extend(int c) {
        int cur = ++cnt, p = last;
        len[cur] = len[last] + 1;
        for (; p && !nxt[p][c]; p = link[p]) nxt[p][c] = cur;
        if (!p) link[cur] = 1;
        else {
            int q = nxt[p][c];
            if (len[q] == len[p] + 1) link[cur] = q;
            else {
                int nq = ++cnt;
                len[nq] = len[p] + 1;
                link[nq] = link[q];
                copy(nxt[q], nxt[q] + EPSILON, nxt[nq]);
                for (; nxt[p][c] == q; p = link[p]) nxt[p][c] = nq;
                link[q] = link[cur] = nq;
            }
        }
        last = cur;
    }
}

SAM绝对是代码中“微言精义”的典范233。

其次在DFS遍历SAM的后缀链接树的时候其实有一个小技巧,就是把节点按照\operatorname{len}从大到小的顺序进行访问(这一步可以桶排),这样可以确保儿子节点一定在父亲节点之前访问到。