经典算法|分治策略与求解 | SLie's Blog|琴弦之轮

基本思想

在计算机科学中，分治法是一种很重要的算法。

字面上的解释是“分而治之”，就是把一个复杂的问题分成两个或更多的相同或相似的子问题，再把子问题分成更小的子问题……直到最后子问题可以简单的直接求解，原问题的解即子问题的解的合并。

不难总结出分治策略的三大步骤：

分解|Divide 将原始问题划分或者归结为规模较小的子问题
解决|Conquer 递归或迭代求解每个子问题
合并|Combine 将子问题的解综合得到原问题的解

注意：

子问题应与原始问题性质完全一样
子问题之间可彼此独立地求解
递归停止时子问题可直接求解

对于分治法求解问题的分析，通常我们可以用递推式，又叫 递归式|Recurrence 来进行刻画。并且一般地，其递归式形如：

$T(n)=\sum_i^ka_iT(n-i)+f(n)\tag{1}$

或：

$T(n)=cT(n/b)+g(n)\tag{2}$

即规模为 $n$ 的问题，其时间可以用若干个小规模的子问题时间的线性组合，再加上非递归代价 $f(n)$ ，这里可以代表分治算法中用于分解和合并的代价，组成。

对于递归式的分析与计算，已在下面的本站文章中给出，此处从略。

算法分析基础|Algorithm-Overview

https://qslie.top/posts/4ea9ab3c

改进措施

减少子问题个数
增加预处理

经典实例一览

排序问题类运用

二分查找问题
二分排序问题
归并排序问题
快速排序问题

移步👉🏻 排序算法全TIPS收集企划

快速选择类运用

最大最小选择问题
第 $k$ 大选择问题
线性时间选择问题

移步 👉🏻 选择算法也要收录

数学运用

最大子序和问题
最近点对问题
快速幂算法
大整数乘法问题
Stassen算法
FFT快速傅里叶变换

其他生活运用

汉诺塔问题
芯片测试问题
棋盘覆盖问题
循环赛日程表

最大子序和| Maximum Subarray

给定一个整数序列 $L=(a_1,a_2,...,a_n)$ ，要求找出一个具有最大和的连续子序列（子序列最少包含一个元素），返回其最大和。

不难写出，最大和的表示：

$\max_{i,j}\left\{0,\max_{1\leq i\leq j\leq n}\sum_{k=i}^ja_k\right\}$

输入：整数序列 $L=(a_1,a_2,...,a_n)$
输出：最大和
实例：最大子数组和|LeetCode

算法思想

假设给定序列 $L$ 在计算机中以数组 nums 进行存储，考虑如下的分治策略.

将数组 $\boldsymbol{nums[}left..right\boldsymbol ]$ 从中心点 $mid$ 处划分为左右两个子数组：
$\boldsymbol{nums[}left..mid\boldsymbol ],\quad \boldsymbol{nums[}mid+1..right\boldsymbol ]$

则原始问题的解只有3种情况：

由左子数组内的某个最大连续子数组组成；
由右子数组内的某个最大连续子数组组成；
横跨中心点由某个最大连续子数组 $\boldsymbol{nums[}i..mid]+\boldsymbol{nums[}mid..j\boldsymbol ]$ 组成。

对于1，2两点，相当于是规模更小的同类子问题，可以直接递归实现。
对于第3点，可以考虑从中心点出发往左右两边扩展，直到子数组的和最大。

算法伪码如下：

$\begin{aligned} &\text{Algorithm: }\;\text{Find-Max-in-Mid}(nums,left,mid,right)\\\\ 1.&\;left\text{-}sum\leftarrow-\infty\\ 2.&\;sum\leftarrow0\\ 3.&\;\mathbf{for}\;i=mid\;\mathbf{downto}\;left\;\mathbf{do}\\ 4.&\;\qquad sum\leftarrow sum+nums[i]\\ 5.&\;\qquad \mathbf{if}\;sum\gt left\text{-}sum\;\mathbf{then}\\ 6.&\;\qquad\qquad left\text{-}sum\leftarrow sum\\ 7.&\;\qquad\qquad max\text{-}left\leftarrow i\\ 8.&\;right\text{-}sum\leftarrow-\infty\\ 9.&\;sum\leftarrow0\\ 10.&\;\mathbf{for}\;j=mid+1\;\mathbf{to}\;right\;\mathbf{do}\\ 11.&\;\qquad sum\leftarrow sum+nums[j]\\ 12.&\;\qquad \mathbf{if}\;sum\gt right\text{-}sum\;\mathbf{then}\\ 13.&\;\qquad\qquad right\text{-}sum\leftarrow sum\\ 14.&\;\qquad\qquad max\text{-}right\leftarrow j\\ 15.&\;\mathbf{return}\;(max\text{-}left,max\text{-}right,left\text{-}sum+right\text{-}sum) \end{aligned}$

$\begin{aligned} &\text{Algorithm: }\;\text{Maximum-Subarray}(nums,left,right)\\\\ 1.&\;\mathbf{if}\;left=right\;\mathbf{then}\\ 2.&\;\qquad\mathbf{return}\;nums[left]\\ 3.&\;\mathbf{else}\\ 4.&\;\qquad mid\leftarrow \lfloor(left+right)/2\rfloor\\ 5.&\;\qquad left\text{-}max\leftarrow \text{Maximum-Subarray}(nums,left,mid)\\ 6.&\;\qquad right\text{-}max\leftarrow \text{Maximum-Subarray}(nums,mid+1,left)\\ 7.&\;\qquad mid\text{-}max\leftarrow \text{Find-Max-in-Mid}(nums,left,mid,right)\\ 8.&\;\qquad\mathbf{return}\;\max{(left\text{-}max,right\text{-}max,mid\text{-}max)} \end{aligned}$

时间复杂度
容易得出，该算法的递归式为 $T(n)=2T(n/2)+\Theta(n)$ 。利用主定理可以得出其时间复杂度为 $\Theta(n\log n)$ .

编程实现

int maxInMid(vector<int>& nums,int left, int mid, int right){
	int left_sum = INT_MIN;
	int right_sum = INT_MIN;
	int sum = 0, max_left = 0, max_right = 0;
	for(int i = mid; i >= left; i--){
		sum += nums[i];
		if(sum > left_sum){
			left_sum = sum;
			max_left = i;
		}
	}
	sum = 0;
	for(int j = mid+1; j <= right; j++){
		sum += nums[j];
		if(sum > right_sum){
			right_sum = sum;
			max_right = j;
		}
	}
	return left_sum+right_sum;
}
int maxSubArrayReal(vector<int>& nums, int left, int right) {
	if(left == right)
		return nums[left];
	int mid = (left + right)/2;
	int left_max = maxSubArrayReal(nums, left, mid);
	int right_max = maxSubArrayReal(nums, mid+1, right);
	int cross_max = maxInMid(nums, left, mid, right);
	return max(max(left_max,right_max),cross_max);
}
int maxSubArray(vector<int>& nums){
	return maxSubArrayReal(nums, 0, nums.size()-1);
}

快速幂算法| Quick Power

移步 👉🏻 斐波那契与快速幂

最近点对问题| Closest pair of points

给定 $d$ 维空间中的 $n$ 个点（ $n\gt1$ ），输出两个点，要求这两个点的欧式距离是此空间中所有点对距离中最小的。

可参考👉🏻Closest pair of points|UCSB大学课件

采用暴力算法需要对 $C_n^2$ 对点进行依次比较，显然其时间复杂度在 $O(dn^2)$ 。而如果采用分治思想，可以将复杂度降到 $O(n\log n)$ 。

对于一维空间，我们可以用目前最好是排序算法先对数据进行排序，然后进行一次相邻点距离搜索即可得到结果，复杂度也在 $O(n\log n)$ 级别。当然，我们也可以使用分治思想，我们将在二维平面里详细阐述，思想是一样的。而对于 $d$ 维空间也是一样，不同的细节就在分治之后的合并部分，详情见课件。

这里我们只以二维空间为例，探讨如何利用分支方法处理二维平面的最近点对问题。

算法思想

将二维平面 $P$ 划分为点个数大致相等的两个平面 $P_L,P_R$ 然后分别求取两个子问题，再考虑合并子问题。

Divide：作中垂线 $l$ ，使得 $P_L,P_R$ 分别划分出均有 $n/2$ 个点的点集 $S_1,S_2$ ；
Conquer：分别计算两个平面里面最近的点对 $\delta_1,\delta_2$ ；
Combine：计算两个平面之间的最近点对 $(p,q)\quad p\in S_1,q\in S_2$ ，得到 $\delta_{cross}$ ；
取三个值中的最小值返回对应点对.

事实上，取 $\delta=\min(\delta_1,\delta_2)$ ，然后以中垂线为中心往左右两边分别扩张 $\delta$ 的范围后，处于这个范围内的点对才有可能使得它们的距离 $\delta_{cross}\lt\delta$ 。

如图，以其中一个 $p\in S_1$ 为例，所有与之距离在 $\delta$ 以内的 $q\in S_2$ 只会坐落在这样的尺寸为 $\delta\times2\delta$ 的矩形 $R$ 中：

最近点对的合并部分

值得注意的是，如果将这个矩形按 $\delta/2\times2\delta/3$ 划分为6个小矩形后（如下图所示），每个小矩形的对角线长 $d=5\delta/6\lt\delta$ 。那么就不会有两个点都落在这个小矩形里（如果有，那么右半平面的最近距离就不可能是现在的这个 $\delta$ ，而是更小的了）。

化简为6个小矩形

也就是说， $p$ 最多只需要和6个 $q\in S_2$ 进行比较！那么接下来问题就是如何快速地选到对应的 $p$ 和 $q$ 了。

事实上，如果我们对 $(l-\delta,l+\delta)$ 区域内的点进行排序：以 $x$ 轴为第一关键字、 $y$ 轴为第二关键字。那么我们其实可以直接遍历该区域内所有的点（不再考虑是属于 $P_L$ 还是 $P_R$ 了），然后直接计算该点与此后7个点的距离。如果有小于 $\delta$ 的点对，那么这两个点一定是跨中垂线的（因为纯左或者纯右的点对的最短距离一定不超过 $\delta$ ）。

为什么是“此后”，不考虑“此前”？因为是按顺序遍历的，所以当前点之前的点已经把当前点与“此前”的距离计算一遍了。
为什么是7个点？与前面我们推出 $S_2$ 最多有6个点与 $p$ 点计算距离的思想类似，这里是考虑一个点下方矩形，这个矩形里最多还有7个点使得距离在 $\delta$ 以内。

另外，其实我们可以不必每一次都对点进行排序，只需在分治之前一次性先排序一次即可。这样一来我们便可以求出其时间复杂度： $W(n)=T(n)+\Theta(n\log n)$ 。其中 $T(n)$ 为递归过程，有：

$T(n)=2T(n/2)+O(n)\Rightarrow\;T(n)=O(n\log n)$

编程实现

#include <iostream>
#include <string.h>
#include <vector>
#include <cstdlib>
#include <climits>
#include <algorithm>

using namespace std;

long long min_d = INT64_MAX;

struct Point
{
    int id;
    int x,y;
    Point(){}
    Point(int _x, int _y){
        x = _x;
        y = _y;
    }
};

long long distence(Point a, Point b){
    if(a.id == b.id) return INT_MAX; //处理边界问题
    return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}

int cmp_with_x(Point a, Point b){
    if(a.x != b.x) return (a.x)<(b.x);
    else return (a.y)<(b.y);
}

//处理合并部分
Point CrossPart(Point *S, int mid, Point P_L, Point P_R, int l, int r){
    long long Lmin = distence(S[P_L.x],S[P_L.y]);
    long long Rmin = distence(S[P_R.x],S[P_R.y]);
    long long delta = min(Lmin, Rmin);
    //初始化结果
    Point P;
    if(Lmin > Rmin) P = P_R;
    if(Lmin < Rmin) P = P_L;
    if(Lmin == Rmin){
        if(S[P_L.x].id < S[P_R.x].id || (S[P_L.x].id == S[P_R.x].id && S[P_L.y].id < S[P_R.y].id)){
            P = P_L;
        }else{
            P = P_R;
        }
    }
  
    //将处于(mid-delta,mid+delta)这个带状区域的下标收集到 BarSet中
    vector<int> BarSet;
    for(int i = l; i <= r; i++){
        if((S[i].x-S[mid].x)*(S[i].x-S[mid].x) < delta){
            BarSet.push_back(i);
        }
    }

    min_d = delta;
    long long temp_d;
    int a, b;
    // 遍历带状区域的所有点，依次往下计算与7个点的距离
    for(int i = 0; i < BarSet.size(); i++){
        for(int j = i+1; j < BarSet.size() && S[BarSet[j]].y <= S[BarSet[i]].y+delta; j++){
            temp_d = distence(S[BarSet[i]],S[BarSet[j]]);
            //更新点
            if(temp_d > min_d) continue;
            a = BarSet[i]; b = BarSet[j];
            if(S[a].id > S[b].id) swap(a, b);
            //满足同距离时序号最小
            if(temp_d == min_d){
                if(S[a].id < S[P.x].id || (S[a].id == S[P.x].id && S[b].id < S[P.y].id)){
                    P.x = a; P.y = b;
                }
            }else{
                min_d = temp_d;
                P.x = a; P.y = b;
            }
        }
    }
    return P;
}

Point PartCacula(Point *S, int l, int r){
    // 只有两个点时，直接返回
    if(r - l + 1 < 3){
        if(S[l].id > S[r].id) return Point(r, l);
        return Point(l, r); //注：此处 Point的 x,y只是表示经过排序后得到的下标，而不是编号
    }
    int mid = l+(r-l+1)/2;
    Point P_L = PartCacula(S, l, mid);
    Point P_R = PartCacula(S, mid+1, r);
    return CrossPart(S, mid, P_L, P_R, l, r);
}


Point ClosestPair(Point* S, int n){
    //预排序
    sort(S, S+n, cmp_with_x);
    //计算结果
    return PartCacula(S, 0, n-1);
}

#define CPC
int main(void){
#ifdef CPC
    char bufin[40012],bufout[40012];
    for(int k = 2; k <= 2; k++){
        sprintf(bufin,"%02d.in", k);
        sprintf(bufout,"%02d.out2",k);
        freopen(bufin,"rb", stdin);
        freopen(bufout, "wb", stdout);
#endif
    int n, m;
    int i, j;
    while(scanf("%d", &n) != EOF){
        Point* S = new Point[n];
        for(int i = 0; i < n; i++){
            S[i].id = i;
            scanf("%d%d", &S[i].x, &S[i].y);
        }

        Point res = ClosestPair(S, n);
        int id1 = S[res.x].id;
        int id2 = S[res.y].id;
        cout << id1 << " " << id2 << endl;
        cout << "=> distance: " << distence(S[res.x],S[res.y]) << endl;
    }

#ifdef CPC
    }
#endif
    return 0;
}

大整数乘法

#include<stdio.h>
#include<string.h>
#include<stdlib.h>
#include<iostream>
#include<algorithm>
#include<vector>
using namespace std;

// 整数位乘
struct BigBinary
{
    vector<int> x;     // 由低位到高位保存二进制位
    bool neg;               // 标记数的正负
    void Init(){x.clear(); neg = false;}
    BigBinary(){Init();}
    int is_zero(){ //消除高位的0
        for(int i = (int)x.size() - 1; i >= 0; i --) {
            if(x[i] == 0) x.pop_back();
            else break;
        }
        return x.empty();
    }
    void Print(){ //规格化输出
        if(neg && !x.empty()) printf("-");
        for(int i = x.size() - 1; i >= 0; i--)
            printf("%d", x[i]);
        if(x.empty()) printf("0");
        printf("\n");
    }
};
BigBinary Add(BigBinary &a, BigBinary &b)
{
    // TODO: 返回 a + b 的结果，小心两数异号情况
    BigBinary c;
    int jinwei = 0, temp;
    for(int i = 0; i < std::max(a.x.size(), b.x.size()); i ++){
        if(i < a.x.size() && i < b.x.size()) temp = a.x[i] + b.x[i] + jinwei;
        else if (i < a.x.size()) temp = a.x[i] + jinwei;
        else temp = b.x[i] + jinwei;
        c.x.push_back(temp % 2);
        jinwei = temp / 2;
    }
    if(jinwei) c.x.push_back(jinwei);
    return c;
}
BigBinary Minus(BigBinary &a, BigBinary &b)
{
    // TODO: 返回 a - b 的结果，注意正负号
    BigBinary c;
    int jiewei = 0, temp;
    for(int i = 0; i < a.x.size(); i ++){
        if(i < b.x.size()) temp = a.x[i] - b.x[i] - jiewei;
        else temp = a.x[i] - jiewei;
        if(temp < 0) {c.x.push_back(temp + 2); jiewei = 1;}
        else {c.x.push_back(temp); jiewei = 0;}
    }
    return c;
}
BigBinary Mul(BigBinary &a, BigBinary &b)
{
    // TODO: 模拟竖式计算两数相乘
    BigBinary c;
    for(int i = 0; i < a.x.size() + b.x.size(); i ++) c.x.push_back(0);
    for(int i = 0; i < a.x.size(); i ++)
        for(int j = 0; j < b.x.size(); j ++)
            c.x[i + j] += a.x[i] & b.x[j];
    int temp = 0;
    for(int i = 0; i < c.x.size(); i ++){
        c.x[i] += temp;
        temp = c.x[i] / 2;
        c.x[i] = c.x[i] % 2;
    }
    c.is_zero();
    return c;
}
void MulN2(BigBinary &a, int n_2)
{
    // TODO: 为 a 添加 n_2 个二进制0，即乘以 2^(n_2)
    for(int i = 0; i < n_2; i++){
        a.x.insert(a.x.begin(),0); //补零（相当于左移）
    }
}
BigBinary FasterMul(BigBinary &a, BigBinary &b)
{
    // TODO: 分治优化位乘
    //初始：<100位时，直接Mul
    // 设 a = A * 2^(n/2) + B， b = C * 2^(n/2) + D
    // a*b = AC*2^n + [(A + B)(D + C) - AC - BD] * 2^(n/2) + BD
    // 注意：a 和 b 位数不一定相同，需要额外处理边界。
    if(a.x.size() <= 500 && b.x.size() <= 500){
        return Mul(a, b);
    }
    int i, j, k;
    BigBinary A, B, C, D;
    BigBinary AC, BD, ApB, DpC, ApB_DpC;
    //统一大小
    int n = max(a.x.size(), b.x.size());
    if(n%2!=0) n++;
    if(n-a.x.size() > 0){
        k = n-a.x.size();
        for(int i = 0; i < k; i++)
            a.x.push_back(0);
    }
    if(n-b.x.size() > 0){
        k = n-b.x.size();
        for(int i = 0; i < k; i++)
            b.x.push_back(0);
    }
    //拆分
    for(i = 0; i < n/2; i++){
        B.x.push_back(a.x[i]); 
        D.x.push_back(b.x[i]); 
    }
    for(j = n/2; j < n; j++){
        A.x.push_back(a.x[j]);
        C.x.push_back(b.x[j]);
    }

    AC = FasterMul(A, C);
    BD = FasterMul(B, D);

    //中间项
    ApB = Add(A, B); DpC = Add(D, C);
    ApB_DpC = FasterMul(ApB, DpC);
    ApB_DpC = Minus(ApB_DpC, AC);
    ApB_DpC = Minus(ApB_DpC, BD);

    // 计算
    MulN2(ApB_DpC, n/2);
    MulN2(AC, n);
    ApB_DpC = Add(ApB_DpC, BD);
    ApB_DpC = Add(ApB_DpC, AC);

    return ApB_DpC;

}

#define CPC
int main(void){
#ifdef CPC
    char bufin[40012],bufout[40012];
    for(int k = 1; k <= 1; k++){
        sprintf(bufin,"%02d.in", k);
        sprintf(bufout,"%02d.out2",k);
        freopen(bufin,"rb", stdin);
        freopen(bufout, "wb", stdout);
#endif
    int n;
    int i, j;
    char ch_a[40012], ch_b[40012];
    
    BigBinary a, b, c;
    // fgets(ch_a,40012,stdin);
    // fgets(ch_b,40012,stdin);
    cin >> ch_a >> ch_b;
    for(int i = (int)strlen(ch_a) - 1; i >= 0; i --) a.x.push_back(ch_a[i] - '0');
    for(int i = (int)strlen(ch_b) - 1; i >= 0; i --) b.x.push_back(ch_b[i] - '0');
    
    c = FasterMul(a,b);
    c.Print();
#ifdef CPC
    }
#endif
    return 0;
}

Strassen算法

目前还有一个更佳的算法：Coppersmith-Winograd 算法，该算法的时间复杂度能到达 $O(n^{2.376})$

FFT快速傅里叶变换

汉诺塔问题| Hanoi

$\begin{aligned} &\text{Algorithm: }\;\text{Hanoi}(A,C,n)\\\\ 1.&\;\mathbf{if}\;n=1\;\mathbf{then}\;\text{Move}(A,C)\;//\text{只有一个盘子时，直接 A→C}\\ 2.&\;\mathbf{else}\\ 3.&\;\qquad \text{Hanoi}(A,B,n-1)\\ 4.&\;\qquad \text{Move}(A,C)\\ 5.&\;\qquad \text{Hanoi}(B,C,n-1)\\ 6.&\;\mathbf{return} \end{aligned}$

容易分析得出，有 $T(n) = 2 T(n-1) + 1,\;T(1) = 1$ ，最终可求得该算法的时间复杂度 $T(n)=2^n-1=\Theta(2^n)$ .

芯片测试| Chip testing

移步 👉🏻 芯片测试|一道分治算法分析题

棋盘覆盖| Tromino Puzzle

这里的棋盘覆盖问题，又称 L型覆盖，Tromino谜题等。

Tromino是一个由棋盘上的三个邻接方块组成的L型骨牌。
要求，用Tromino覆盖一个缺少了一个方块（可以在棋盘上的任何位置）的 $2^n\times 2^n$ 棋盘。除了这个缺失的方块，Tromino应该覆盖棋盘上的所有方块，而且不能有重叠。
Tromino Puzzle

🎲通过该网站可以游玩一个 $8\times 8$ 的Tromino游戏👉🏻The Tromino Puzzle

算法思想
https://www3.cs.stonybrook.edu/~rezaul/Fall-2017/CSE548/CSE548-lecture-2.pdf

循环赛日程表

代码如下：

#include <iostream>
#include <string.h>
#include <vector>
#include <stack>
#include <cstdlib>
#include <climits>
#include <cmath>
#include <queue>
#include <algorithm>

using namespace std;

int main(void){
    int n, m;
    int i, j;
    while(scanf("%d", &n) != EOF){
        vector<vector<int>> a(n+1,vector<int>(n+1)); //建表

        for(i = 1; i <= n; i++){
            a[1][i] = i; //初始化第1行
        }

        // 从 m*m 的小矩阵开始, 逐渐扩张成 2m*2m 的矩阵
        for(m = 1; m <= n/2; m = m*2){
            for(i = m+1; i <= 2*m; i++){
                for(j = m+1; j <= 2*m; j++){
                    a[i][j] = a[i-m][j-m]; //右下=左上
                    a[i-m][j] = a[i-m][j-m]+m; //右上=左上+m
                    a[i][j-m] = a[i-m][j]; //左下=右上
                }
            }
        }

        //输出
        for(i = 1; i <= n; i++){
            for(j = 2; j <= n; j++){
                printf("%d ", a[i][j]);
            }
            printf("\n");
        }
    }
    return 0;
}