算法初步——二分

二分查找

 　　二分查找是基于有序序列的查找方法（以下假设严格单调自增），该算法一开始令 [left,right] 为整个序列的下标区间，然后每次测试当前 [left,right] 的中间位置 mid=(left+right)/2 ，判断 A[mid]与欲查询的元素 x  的大小：

• • 如果 A[mid]  == x ，说明查找成功，退出查询
  • 如果 A[mid]  > x ，说明元素 x 在 mid 位置的左边，因此往左子区间 [left,mid-1] 继续查找
  • 如果 A[mid]  < x ，说明元素 x 在 mid 位置的右边，因此往右子区间 [mid+1,right] 继续查找

 　　因此其时间复杂度是 O(logn)。

#include <cstdio>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;

// A[] 为严格递增序列， left 为二分下界， right 为二分上界， x 为欲查询数 
int binarySearch(int A[], int left, int right, int x) {
    int mid;        // 中点
    while(left <= right) {
        mid = (left+right)/2;
        if(A[mid] == x)    return mid;
        else if(A[mid] > x) {
            right = mid-1;        // 往左区间查询 
        } else {
            left = mid+1;        // 往右区间查询 
        } 
    } 
    
    return -1;                    // 未查询到 
} 

int main() {
    const int n = 10;
    int A[n] = {1, 3, 4, 6, 7, 8, 10, 11, 12, 15};
    // 查询 3 和 5 
    printf("%d %d\n", binarySearch(A, 0, n-1, 6), binarySearch(A, 0, n-1, 5)); 

    return 0;
}

 

 　　那么，如果序列是递减的，只需把上面代码中的  A[mid] > x  改为  A[mid] < x  即可。

　　需要注意的是，如果二分上界超过 int 型数据范围的一半，那么语句  mid = (left+right)/2  有可能导致溢出，应改为  mid = left + (right-left)/2 。

 

　　

　　接下来讨论更进一步的问题：如果递增序列 A 中的元素可能重复，那么如何对给定的欲查询元素 x ，求出序列中第一个大于等于 x 的元素的位置 L 以及第一个大于 x 的元素的位置 R ，这样元素 x 在序列中的存在区间就是 [L,R) 。

　　先来考虑第一个小问：求序列中的第一个大于等于 x 的元素的位置。

• •  如果 A[mid] ≥ x ，说明第一个大于等于 x 的元素的位置一定在 mid 处或 mid 的左侧，应往左子区间 [left,mid] 继续查询，即令 right = mid 。
    
  •  如果 A[mid] < x ，说明第一个大于等于 x 的元素的位置一定在 mid 的右侧，应往右子区间 [mid+1,right] 继续查询，即令 left = mid+1 。

    // A[] 为递增序列，x 为欲查询的数，函数返回第一个大于等于 x 的元素的位置
    // left，right 为左右边界 
    int lower_bound(int A[], int left, int right, int x) {
        int mid;        // 中点
        while(left < right) {
            mid = (left+right)/2;
            if(mid >= x) {        // 中间的数大于等于 x 
                right = mid;    // 左区间 
            } else {            // 中间的数小于 x 
                left = mid+1;    // 右区间 
            }
        } 
        
        return left; 
    }

  •  考虑到欲查询元素有可能比序列中的所有元素都要大，此时应当返回 n ，因此二分上界是 n ，故二分的初始区间为 [left,right] = [0,n] 。
    

　　

　　接下来解决第二小问：求序列中第一个大于 x 的元素的位置。

• • 如果 A[mid] > x ，说明第一个大于 x 的元素的位置一定在 mid 处或 mid 的左侧，应往左子区间 [left,mid] 继续查询，即令 right = mid 。
    
  • 如果 A[mid] ≤ x ，说明第一个大于 x 的元素的位置一定在 mid 的右侧，应往右子区间 [mid+1,right] 继续查询，即令 left = mid+1 。

    // A[] 为递增序列，x 为欲查询的数，函数返回第一个大于 x 的元素的位置
    // left，right 为左右边界，传入的初值为 [0,n] 
    int upper_bound(int A[], int left, int right, int x) {
        int mid;        // 中点
        while(left < right) {
            mid = (left+right)/2;
            if(mid > x) {        // 中间的数大于 x 
                right = mid;    // 左区间 
            } else {            // 中间的数小于等于 x 
                left = mid+1;    // 右区间 
            }
        } 
        
        return left; 
    } 

     

 

　　读者会发现， lower_bound 函数 和  upper_bound 函数的代码高度相似，其实这两个函数都在解决这样一个问题：寻找有序序列中第一个满足某条件的元素的位置。此处总结了解决此类问题的固定模版：

// 解决“寻找有序序列第一个满足某条件的元素的位置”问题的固定模版
// 二分区间为左闭右闭的 [left,right]，初值必须能覆盖解的所有可能取值 
int upper_bound(int A[], int left, int right, int x) {
    int mid;        // 中点
    while(left < right) {
        mid = (left+right)/2;
        if( 条件成立 ) {        // 条件成立，第一个满足条件的元素位置 <= mid 
            right = mid;    // 左区间 
        } else {            // 条件不成立，第一个满足条件的元素位置 > mid 
            left = mid+1;    // 右区间 
        }
    } 
    
    return left;          // 返回夹出来的位置   
} 

 

 

 　　另外，如果想要寻找最有一个满足 “条件 C”的元素的位置，则可以先求第一个满足 “条件 !C”的元素的位置，然后将该位置减 1 即可

　　 而当二分区间是左开右闭区间 (left,right] 时，循环条件应当是  left+1 < right ，这样当退出循环时有  left+1 == right  成立，使得 (left,right] 才是唯一位置。且二分区间初始值应为 (-1,n] 。

 

 

二分法拓展

　　如何计算 √2 的近似值。

　　对 f(x) = x² 来说，令浮点型 left 和 right 的初值分别为 1 和 2，然后根据 left 和 right 的中点 mid 处 f(x) 的值与 2 的大小来选择子区间进行逼近：

• •  如果 f(mid) > 2，说明 mid > √2，应当在 [left,mid]的范围内继续逼近，故令 right = mid 。
  •  如果 f(mid) < 2，说明 mid < √2，应当在 [mid,right]的范围内继续逼近，故令 left= mid 。

　　上面两个步骤当 right - left < 10^-5 时结束。

const double eps = 1e-5;    // 精度为 10^-5
double f(double x) {        // 计算 f(x) 
    return x * x - 2;
} 

double calSqrt() {
    double left=1, right=2, mid;    // [left,right] = [1,2] 
    while(right - left > eps) {
        mid = (left+right)/2;
        if(f(mid) > 0) {            // mid > sqrt(2) 
            right = mid;            // 左区间 
        } else {                    // mid < sqrt(2) 
            left = mid;                // 右区间
        }
    }
    
    return mid;                        // 返回 mid 作为近似值 
}

 

 

 

　　事实上，计算 √2 的近似值的问题其实是这样一个问题的特例：给定一个定义在 [L,R] 上的单调函数 f(x) ，求方程 f(x) = 0 的根。

　　此时只需修改上述代码 f 函数部分即可，显然计算 √2 的近似值等价于求 f(x) = x² - 2 = 0 在 [1,2] 范围内的根。

 

　　装水问题：

　　在一个侧面看去是半圆的储水装置，该半圆的半径为 R ，要求往里面装入高度为 h 的水，使其从侧面看去的面积 S₁ 与半圆面积 S₂ 的比例恰好为 r 。现给定 R 和 r ，求高度 h 。

　　在这个问题中，需要寻找水面高度 h 与 面积比例 r 之间的关系。而很显然，随着水面升高，面积比例 r 一定是增大的。由此可以得到这样的思路：在 [0,R] 范围内对水面高度 h 进行二分，计算在高度下面积比例 r 的值。

#include <cstdio>
#include <string>
#include <algorithm>
#include <cmath>
using namespace std;

const double PI = acos(-1.0);     // PI
const double eps = 1e-5;        // 精度为 10^-5

double f(double R, double h) {    // 计算 r = f(h) 
    double alpha = 2 * acos((R-h)/R);
    double L = 2 * sqrt(R*R - (R-h)*(R-h));
    double S1 = alpha * R * R / 2 - L *(R-h) / 2;
    double S2 = PI * R * R /2;
    return S1 / S2; 
} 

double solve(double R, double r) {
    double left = 0, right = R, mid; 
    while(right-left > eps) {
        mid = (left+right)/2;
        if(f(R, mid) > r) {    
            right = mid;        // 左区间 [left,mid] 
        } else {
            left = mid;            // 右区间 [mid,right] 
        }
    }
    
    return mid;                    // 返回 mid 即为所求 h  
}

int main() {
    double R, r;
    scanf("%lf%lf", &R, &r);
    printf("%.4f\n", solve(R, r)); 

    return 0;
}

 

 

　　木棒切割问题：

　　给定 N 根木棒，长度均已知，现在希望通过切割它们来得到至少 K 段长度相等的木棒（长度必须为整数），问这些长度相等的木棒最长能有多长。

　　首先可以注意到一个结论：如果长度相等的木棒的长度 L 越长，那么得到的木棒段数越少。从这个角度出发便可以先到本题的算法，即二分答案，根据对当前长度 L 来说能得到的木棒段数 k 与 K 的大小关系进行二分。

 

　　最后一个问题（没有思路）：

　　给定 N 个线段的长度，试将它们头尾相接（顺序任意）地组合成一个凸多边形，使得凸多边形的外接圆的半径最大，求该最大半径。其中 N 不超过 10⁵ ，线段长度均不超过 100 ，要求算法中不涉及坐标的计算。

 

 

快速幂

　　来看一个问题：给定三个正整数 a、b、m（a<10⁹, b<10¹⁸, 1<m<10⁹），求 a^b % m 。

　　显然不能直接计算，这里要使用快速幂的做法，它基于二分的思想。快速幂基于以下事实：

• • 如果 b 是奇数，那么有 a^b = a * a^b-1 。
  • 如果 b 是偶数，那么有 a^b = a^b/2 * a^b/2 。

　　快速幂的递归写法如下：

LL binaryPow(LL a, LL b, LL m) {
    if(b == 0)    return 1;
    // b 为奇数
    if(b % 2) return a * binaryPow(a, b-1, m) % m;
    else {    // b 为偶数 
        LL mu1 = binaryPow(a, b/2, m);
        return mu1 * mu1 % m;
    }
} 

　　

　　接下来研究一下快速幂的迭代写法。

　　对 a^b 来说，如果把 b 写成二进制，那么 b 就可以写成若干二次幂之和。例如 13 的二进制是 1101 ，那么 13=2³+2²+2⁰，所以 a¹³ = a⁸ * a⁴ * a¹ 。

　　我们可以把 a^b  表示成 a^2k、... 、a⁸、 a⁴、a²、a¹ 中若干项的乘积，其中如果 b 的二进制的 i 号位为1，那么项 a²ⁱ 就被选中。

　　快速幂的迭代写法：

LL binaryPow(LL a, LL b, LL m) {
    LL ans = 1;
    while(b > 0) {
        if(b & 1) {        // 末位为 1 
            ans = ans * a % m;        // 令 ans 累乘上 a  
        }
        a = a * a % m;    // 令 a 平方
        b >>= 1;        // b 右移1，相当于除2 
    }
    return ans; 
}