跳到主内容

「译」整数的快速开方

原文

http://www.nuprl.org/MathLibrary/integer_sqrt/

示例代码由 Python 重写。

简介

对于一个自然数 \(x(x \in {0, 1, 2, 3, ...})\),它的自然数平方根 \(r\) 定义为满足 \(r^2 \le x < (r + 1)^2\) 的自然数。下面这个表显示了在一个较小范围内整数平方根的分布:

/images/isqrt-1024.svg

对于给定的自然数 \(x\),我们发现其平方根与 \(x-1\) 的平方根有关:有时是相同的,有时是其前驱。因此可以通过测试 \(x-1\) 的平方根来编写递归的算法:

def isqrt(x: int) -> int:
    if x == 0:
        return 0
    else:
        r2 = isqrt(x - 1)
        r3 = r2 + 1
        if x < r3 * r3:
            return r2
        else:
            return r3

根据自然数归纳法的基本法则,要证明一个任意自然数 \(x\) 的陈述 \(P(x)\),则需要证明基例 \(P(0)\) 并证明 \(P(x)\) 可以通过 \(P(x-1)\) 来派生。上面的递归函数实际上通过归纳法证明了下面这个定理:

\begin{equation*} \forall x: \N. (\exists r: \{\N | (r^2 \le x) \land (x < (r+1)^2)\}) \tag{ht1:\text{自然数标准规范}} \end{equation*}

当我们在 Nuprl 中证明该定理时,我们进行了构造性的证明,这意味着为了证明整数平方根的存在,我们必须证明如何构造它。 然后我们可以从该构造性证明中提取程序,并且这样做实际上证明了该程序是正确的 —— 我们知道结果将满足 th1 的规范。我们说该程序在构造上是正确的。

但是,尽管该算法是正确的,但它并不是最有效的方法。 为了以这种方式求解 \(x\) 的整数平方根,必须首先求解 \(x-1\) 的根。 同样,为了求解 \(x-1\) 的根,必须首先找到 \(x-2\) 的根。 递减 1 的过程一直持续到 0 为止。所以现在,我们需要另一个可以证明 th1 的,更快的算法。

快速归纳

为了找到一个快速的算法,我们转向自然数的完全归纳法。完全归纳法告诉我们,要对任意自然数 \(x\) 证明一个陈述 \(P(x)\),需要证明对于任意 \(y < x\),都能从 \(P(y)\) 派生出 \(P(x)\)。这是个更强的假设,形式上可以表示为以下定理( \(\N{}x\) 表示比 \(x\) 小的自然数):

\begin{align*} \forall [P:\N \rightarrow \mathbb{P}]. ((\forall x: \N. ((\forall y: \N{}x. P[y]) \Rightarrow P[x]) \Rightarrow (\forall x: \N. P[x]))) \\ \tag{th2:\text{完全自然归纳}} \end{align*}

标准归纳法是完全归纳法在 \(y=(x-1)\) 时的特例。另外尽管基本陈述在 th2 中未指出,它仍然是证明的一步。例如,当对 \(x-1\) 使用完全归纳法时,由于 \(x=0\) 时的 \((x-1)\) 并不是自然数,需要单独考虑 \(x=0\) 时的情况。

这个强归纳原理打开了用比 \(1\) 更大的步长递减 \(x\) 的大门。事实上,我们可以证明对于任意 \(b \in \N \land b > 1\)\(x \div b\) 也是一个可接受的步长。我们可以将其想象成完全归纳法的又一个特例。对于此原理,我们需要 \((x\div b) < x\)。因此我们说明 \(b>1\) 而且基例 \(P(0)\) 被单独证明。

\begin{align*} \forall b:\{b:\mathbb{Z}| 1 < b\}. \forall [P:\mathbb{N} \rightarrow \mathbb{P}]. (P[0] \Rightarrow (\forall x:\mathbb{N}^{+}. (P[x \div b] \Rightarrow P[x])) \Rightarrow (\forall x:\mathbb{N}. P[x])) \\ \tag{th3:\text{除法自然归纳}} \end{align*}

所以为了改进整数平方根算法的速度,我们可能需要使用整数除法来代替减一来步进到 \(0\) 。但是我们应该除以哪个数字?使用 th3 描述的归纳原理,我们可以假设 \(x\div b\) 存在一个根 \(r_0\),然后用 \(r_0\) 来构造 \(x\) 的一个根。更特殊地,归纳假设存在一个自然数 \(r_0\),满足 \(r_0^2 \le x \div b < (r_0 + 1)^2\),并且我们需要构造一个新的 \(r\),满足 \(r^2 \le x < (r+1)^2\)。在这个时候忽略余数并且只使用一些代数,我们可以发现 \(\sqrt{b} \cdot r_0\) 可能会表现为 \(x\) 的平方根。因此,我们需要尝试 \(b\) 的完全平方,这很正常,因为我们在求平方根。如下章节所属,如果我们在归纳步骤中除以 \(4\),那么 \(x\) 的平方根将为 \(2r_0, (2 r_0 + 1)\) 之一,取决于 \(x\) 是否小于 \((2r_0 + 1)^2\) 。当我们证明了这个结论,我们将其提取为一个高效的算法:

def isqrt(x: int) -> int:
    if x == 0:
        return 0
    else:
        z = x // 4
        r2 = 2 * isqrt(z)
        r3 = r2 + 1
        if x < r3 * r3:
            return r2
        else:
            return r3

该结果可以推广到更普遍的情况,即找到任意正整数 \(n\) 的自然数的第 \(n\) 个整数根,例如找到自然数 \(r\) 使得 \(r^n \le x <(r+1)^n\)。 可以使用 th3 将除数设为 \(2^n\) 来归纳证明更普遍的根的存在。可以从证明中提取以下算法:

def nroot(n: int):
    b = 2**n

    def nroot_(x: int):
        if x == 0:
            return x
        else:
            z = x // b
            r2 = 2* nroot_(z)
            r3 = r2 + 1
            if r3**n < x+1:
                return r3
            else:
                return r2

    return nroot_

使用快速归纳证明整数平方根(略)

原文第三章节使用 NuPRL 语言来证明完全归纳法的细节,由于非专业,忽略不译。

查找整数平方根的其他方法

如果我们要走到大街上随便找一个人并要求他或她算出一个数字的整数平方根(并解释我们的整数平方根的意思),我们认为那个人将如何解决?

对于小于 100 的数字,人们很可能会使用简单的猜测和检查方法来找到解决方案。例如,如果提示找到 56 的根,人们可能会认为 \(7^2 = 49\) 以及 \(8^2 = 64\),所以答案必须是 7 。

这种猜测检查方法是我们可以正式编写的搜索算法的缩写。 要找到数字 \(x\) 的根 \(r\) ,我们可以从 \(r = 0\) 开始,然后将 \(r\) 加 1,直到找到 \((r + 1)^2>x\) 。 或者,我们可以编写一种二分搜索算法:首先以 \((x \div 2)\) 作为根进行测试,然后仅根据结果在上方或下方进行搜索。 然后,每次都测试搜索区域的中点以继续。 通常,搜索算法涉及测试输出 \(r\) 作为潜在解决方案。 相反,前面描述的所有方法都涉及对输入变量 \(x\) 的递归。

现在,超出我们记忆或愿意测试的更大数量的数字呢? 好吧,我们可能只是拿出一个计算器来找到十进制答案,然后四舍五入到最接近的整数。 这个想法促使我们探索使用实数,十进制平方根查找整数平方根的方法。 此方法将涉及两个主要步骤:

  1. 找到实平方根;

  2. 将结果截断为整数。 假设我们找到一种计算实平方根的算法,那么第二步可能会遇到问题。 直观地,4.000000001和3.999999999可能是等效的实数,但是截断显然会导致两个不同的整数。

在研究这种情况时,我们意识到,以相反的方向处理整数平方根与实平方根之间的连接可能会更有趣。 也就是说,我们意识到使用整数平方根对计算实平方根很重要,而不是使用实平方根来查找整数平方根。 这个方向为您可能正在问自己的问题提供了见解,整数平方根起什么作用?

我们发现,在计算实平方根时,可以使用整数平方根来提高效率。 在下一节中,我们将介绍实数的一些概念,并讨论实数平方根。(后略)

总结

函数

复杂度

测试 \(x-1\) 的 isqrt 函数

\(O(x)\)

测试 \(x \div 4\) 的 isqrt 函数

\(O(\log{x})\)

二分查找

\(O(\log{x})\)

测试 \(x \div 2^n\) 的 nroot 函数

\(O(\log{x})\)

这里有一个集成的脚本,可以让你运行并显示文中两种整数平方根和一种整数任意根算法的效率:

tr-fast-isqrt.py (源文件)

def count(title: str):
    class Counter:
        def __init__(self, fn):
            self._fn = fn
            self._c = 0
            self.r = None

        def __call__(self, *args, **kwargs):
            self._c += 1
            self.r = self._fn(*args, **kwargs)
            return self.r

        def show(self):
            print(f"{title}[{self._c}]: {self.r}")

    return Counter


@count("isqrt1")
def isqrt1(x: int) -> int:
    if x == 0:
        return 0
    else:
        r2 = isqrt1(x - 1)
        r3 = r2 + 1
        if x < r3 * r3:
            return r2
        else:
            return r3


@count("isqrt2")
def isqrt2(x: int) -> int:
    if x == 0:
        return 0
    else:
        z = x // 4
        r2 = 2 * isqrt2(z)
        r3 = r2 + 1
        if x < r3 * r3:
            return r2
        else:
            return r3

def nroot(n: int):
    "Curring 化,以便携带 show 方法"
    b = 2**n

    @count("nroot")
    def nroot_(x: int):
        if x == 0:
            return x
        else:
            z = x // b
            r2 = 2* nroot_(z)
            r3 = r2 + 1
            if r3**n < x+1:
                return r3
            else:
                return r2

    return nroot_

if __name__ == "__main__":
    n = 156
    isqrt1(n)
    isqrt2(n)

    x, b = 128, 3
    nroot_ = nroot(b)
    nroot_(x)

    isqrt1.show()
    isqrt2.show()
    nroot_.show()