如何算法实现加减乘除四则运算

本文记录一些基本数学运算的算法。

内容包括:

本文内容较多,需要读者静心阅读。

整数加法

不用加减乘除四则运算符,如何编程实现加法?

leetcode 题目 - 不用加减乘除做加法

回顾熟悉的十进制加法的竖式计算:

二进制加法也是类似的,满二进一:

把进位的一栏记为 c ,不考虑进位的和记为 s

如果上一位的比特值之和满二,则下一位 c1s 则只记录当前位的和的比特值。

如此一来,上面的二进制竖式计算拆解如下:

可以看到,不断迭代竖式计算,直到进位 c0 的时候,s 即最终的计算结果。

在每一栏位上,下一位的 c 和当前位的 s 的计算表格如下:

仔细观察,可知 c 的值其实是「逻辑与」的结果,s 的值其实是 「逻辑异或」的结果:

因此,c 的数值可以表达为 (a & b) << 1s 的数值可以表达为 a ^ b

每次计算出 cs 后,把 c 当做新的 as 当做新的 b , 不断循环迭代,直到 a0 ,就得出最终计算结果 b

这个方法对负整数是否成立呢?负数的最高位是符号位 1 ,下面是 8bit 负整数 -123 = 10111 的结果:

算法仍然适用,事实上,在二进制上,负数无非即大正数。

加法的位运算实现 - C 语言
// 加法
int Add(int a, int b) {
    do {
        int c = (unsigned int)(a & b) << 1;
        int s = a ^ b;
        a = c;
        b = s;
    } while (a != 0);
    return b;
}

以上算法中,每次左移一位, 步骤数不超过类型的比特位数, 因此时间复杂度是常数级别。

时间最坏情况下,当进位最多的时候。下图是 4bit 整数情况下的算法最差情况的示例图, 计算 -11 相加的过程,需要 4 步完成。

对于 32bit 有符号整数来讲,最坏情况则需要 32 步完成。

事实上,这个方法和计算机中 加法器 的实现非常类似。

一位半加器 接收两个比特的输入,经过 异或门 输出当位的和 S ,经过 与门 输出进位 C

半加器 - 图片来自 wikipedia

整数减法

减法,可以用加法来表达, a - b = a + (-b)

8bit 有符号整数为例,假如 b = 10111 ,那么 -b 可以如下得出:

又因为 -1-b 其实就是 ~b ,所以 -b = ~b+1

因此 a - b = a + (~b+1) ,即 Minus(a, b) = Add(a, Add(~b, 1))

减法的位运算实现 - C 语言
// 减法
int Minus(int a, int b) { return Add(a, Add(~b, 1)); }

整数乘法

乘法即是循环的加法,不过如此计算的话,时间复杂度是 O(a) 或者 O(b)

乘法的循环加实现 O(n) - C 语言
#define MIN(a, b) (a) > (b) ? (b) : (a)
#define MAX(a, b) (a) > (b) ? (a) : (b)

// 乘法 - 简单方法 O(min(a,b))
// 假设加法已经实现
// 以下只考虑正数实现,负数可转正数处理
int MulNaive(int a, int b) {
    int min = MIN(a, b);
    int max = MAX(a, b);
    int r = 0;

    while (min-- != 0) r += max;
    return r;
}

快一些的办法是,二分法。

举例来说,13x17 的计算过程可以拆解如下:

拆解过程中,把 a 当成了倍数,b 当成了被乘数。

如果倍数是偶数,则二分倍数,否则化倍数为偶数。

容易得出此方法的递归版本:

乘法的二分法实现 - 递归版本 - C 语言
// 乘法 - 二分法 递归版本
// 以下只考虑正数实现,负数可转正数处理
int MulRecursive(int a, int b) {
    int min = MIN(a, b);
    int max = MAX(a, b);

    if (min & 1)  // 奇数
        return MulRecursive(min - 1, max) + max;
    else {  // 偶数, 注意 min/2 即 min >> 1
        int half = MulRecursive(min >> 1, max);
        return half << 1; // 即 half x 2
    }
}

另外一种视角,是把 a x b 分解为三个部分:

  • 计算结果 result
  • 乘式因式 factor
  • 剩余倍数 remain

构成方式是 a x b = result + factor * remain 的形式。

尝试不断减少剩余倍数:

  • 如果是偶数,则分解一个 2 到左边的乘式因式 factor 中, 即 factor 增大一倍,剩余倍数减半。
  • 如果是奇数,则分解一个被乘数 factor 加到左边的 result 中,剩余倍数减一。

当剩余倍数缩减到 0 时,result 则积累到最终结果。

下面是这种算法的一种最好情况,剩余倍数不断减半:

乘法的二分法实现 - 循环版本 - C 语言
// 乘法 - 二分法循环版本
// 以下只考虑正数实现,负数可转正数处理
int MulBinary(int a, int b) {

    int min = MIN(a, b);  // 小的作倍数
    int max = MAX(a, b);  // 大的作被乘数

    int r = 0;         // 计算结果
    int factor = max;  // 乘积因子
    int remain = min;  // 剩余倍数

    while (remain != 0) {
        if (remain & 1) {
            // 剩余倍数为奇数
            r += factor;
            remain--;
        } else {
            remain >>= 1;  // 剩余倍数缩小一半
            factor <<= 1;  // 被乘数扩大两倍
        }
    }
    return r;
}

显然,时间复杂度是 O(logn) ,其中 nab 中较小的一个数。

对于 32bit 大小的整数来说,不超过 O(32)

最后一种思路,回想一下十进制数的乘法竖式:

二进制的乘法竖式,是这样的:

相比十进制乘法,它更简单,0 乘以任何数是 01 乘以任何数是自身, 用不到乘法九九算术表:

假设 ab 的比特位数是 nm , 把一个数字 b 不断左移 n 次,得到的数求和即乘法结果。

具体来说:

  1. 不断右移数字 a ,直到它为 0 时,即移动了 n 次。
  2. 每次右移数字 a 后,其当前最后一位是 a & 1,和数字 b 相乘:
    1. 0 乘以任何数是 0 , 因此不需要处理这种情况。
    2. 1 乘以任何数是自身 ,只需要处理这种情况。
    3. 将乘积左移一位,加到计算结果求和。

即得所谓的快速乘算法。

乘法的二进制竖式思路的实现 - C 语言
// 乘法 - 竖式计算方法
// 以下只考虑正数实现,负数可转正数处理
int MulFast(int a, int b) {
    int r = 0;  // 乘法结果

    while (a != 0) {
        if (a & 1) {
            // 当前位是 1, 1 乘以任何数为自身
            // 加上此时的 b
            r += b;
        }
        // 否则,当前位是 0, 0 乘以任何数是 0, 无需加上

        a >>= 1;  // 移动 n 次到 0 意味着 a 有 n 个比特位
        b <<= 1;  // 同步左移 b
    }
    return r;
}

容易知道,此算法的时间复杂度是整数的比特位数,对于 32bit 整数即 O(32)

整数除法

leetcode 题目 - 两数相除

乘法 的分析思路类似, 容易的办法是,化除法为循环的减法。

此处 O(n) 时间复杂度的算法从略。

考虑 255 / 17 的结果,和 乘法 类似,把计算式拆解为如下形式:

\[\frac {a} {b} = \frac {remain - b \times factor} {b} + {result}\]

其中:

  • 被减的因式系数 factor : 每次增大两倍,如果 remain 不足,则 factor 重置 1
  • 剩余的被除数 remain : 被除数不断被减后的剩余部分。
  • 计算结果 result :每次增长 factor 大小。

当被除数剩余部分 remain 不足 b 大小时,终止算法过程。

除法的二分法实现 - C 语言
// 以下只考虑正数实现,负数可转正数处理
int Div(int a, int b) {
    int remain = a;
    int factor = 1;
    int bfactor = b;  // b x factor.
    int result = 0;

    while (remain >= b) {
        if (remain < bfactor) {
            // 被除数不足以相减,重置 factor
            factor = 1;
            bfactor = b;
        }

        remain -= bfactor;
        result += factor;

        // factor x2
        factor = factor << 1;
        // bfactor x2
        bfactor = bfactor << 1;
    }
    return result;
}

时间复杂度是 O(logn) ,对于 32bit 整数来说,即 O(32)

另外一种思路,是源于除法竖式。

回顾十进制的除法竖式计算方法:

二进制的竖式除法也是一样的:

算法的思路则是:

  • 初始化,r 为除法结果。
  • 左对齐 ab 的有效比特起始位。
  • 如果对齐后, b 不大于 a,则更新被除数 a = a-b,同时 r 的当前位标 1
  • 否则,r 的当前位标 0
  • 到两个数右对齐为止,结束计算。

对于 32bit 整数,此算法的时间复杂度显然不超过 O(32)

除法的竖式思路实现 - C 语言
// 获取给定整数的比特位数
int GetNBits(int a) {
    int n = 0;
    while (a) {
        a >>= 1;
        n++;
    }
    return n;
}

// 除法 竖式计算
// 以下只考虑正数实现,负数可转正数处理
int DivFast(int a, int b) {
    int n = GetNBits(a);
    int m = GetNBits(b);

    int r = 0;  // 计算结果

    for (int i = n - m; i >= 0; i--) {
        // 左移 b 对齐当前的 a
        int b1 = b << i;
        if (b1 <= a) {
            // 标记结果的第 i 位为 1
            r |= (1 << i);
            a -= b1;
        }
        // 否则,r 的当前位标 0,即不用设置
    }

    return r;
}

快速幂

leetcode 题目 - pow(x, n)

快速幂的算法和前面所说的 乘法 的分析思路是类似的。

3^13 为例,展开如下:

容易得出二分法的递归版本:

幂运算的二分法递归版本 - C 语言
// 幂运算 - 二分 递归版本
// 以下只考虑正数实现,负数可转正数处理
int PowRecursive(int a, int b) {
    if (b == 0) return 1;
    // 奇数
    if (b & 1) return PowRecursive(a, b - 1) * a;
    // 偶数
    else {
        int half = PowRecursive(a, b / 2);
        return half * half;
    }
}

其时间复杂度是 O(logb) ,对于 32bit 整数来说是 O(32)

不过下面的循环迭代思路的方法更为美妙。

乘法 的拆解思路如出一辙, 把 a^b 写成 result * (factor) ^ remain 的形式:

  • 如果剩余指数 remain 是奇数,则分解一个因子 factor 给到左边的 result
  • 否则,把 remain 拆半,factor 自乘翻倍。

当剩余指数 remain 到达 0 的时候,result 就拿到了结果。

幂运算的二分法循环版本 - C 语言
// 幂运算 - 二分 循环版本
// 以下只考虑正数实现,负数可转正数处理
int PowFast(int a, int b) {
    int factor = a;
    int result = 1;
    int remain = b;
    while (remain > 0) {
        if (remain & 1) {
            // 奇数
            result = result * factor;
            remain--;
        } else {
            // 偶数
            remain = remain / 2;       // 剩余指数拆半
            factor = factor * factor;  // 因子自乘翻倍
        }
    }
    return result;
}

时间复杂度仍然是 O(b) ,对于 32bit 整数即 O(32)

平方根

如何求解 $\sqrt a$ 的值? 已知 $a >= 0$ 。

leetcode 题目 - x 的平方根

由于函数 $x^2$ 在 $x >= 0$ 区间上是单调递增的,因此可以采用二分法, 和 标准的二分查找 差不太多,介绍从略。

平方根运算 - 二分法 - C 语言
// 平方根运算 - 二分
// a 为非负数
int Sqrt(int a) {
    if (a == 0) return 0;
    if (a == 1) return 1;

    int l = 0;
    int r = a;

    while (l < r) {
        int m = (l + r) / 2;
        if (m < a / m) {  // 不采用 m * m < a 防止溢出
            l = m + 1;
        } else if (m > a / m) {
            r = m - 1;
        } else {
            return m;
        }
    }

    // 找到的 l 有可能比实际值稍大
    if (l > a / l) return l - 1;

    return l;
}

数值分析 中, 求函数近似零值的有一个常用的方法:牛顿迭代法

求解 $\sqrt a$ 的值,其实即求函数 $f(x) = x^2 - a$ 的零值。

牛顿法的思路是,先取一个值 $x = x_0$ ,然后在这一点做函数的切线,将切线和横轴的交点 $x_1$ 作为 $x$ 的新值,不断迭代,直到 $f(x)$ 和零的误差满足精度要求。

切线的斜率 $k$ 可以由函数的导数 $f’(x) = 2x$ 计算得出,即 $k = 2x_0$ 。 因此:

\[x_1 = x_0 - \frac {f(x_0)} {f'(x_0)}\]

按照这个方式迭代下去即可。

平方根运算 - 牛顿法 - C 语言
// 平方根运算 - 牛顿法
// a 为非负数
int SqrtNewton(int a) {
    if (a == 0) return 0;

    double x = a;  // 初始值起为 a 本身
    double d = 0.0;

    do {
        // delta 是 f(x) / f'(x) 即 (x^2 - a) / (2x)
        // 防溢出,写作如下方式
        d = x / 2 - a / x / 2.0;
        if (d < 1e-6) break;  // 假设误差是 1e-6
        x -= d;
    } while (x - a / x > 0);

    return x;
}

(完)

* * *
评论 首页 订阅 ↑顶部