精确比较非常大整数(十亿级)的两个指数的算法

Question

我们要比较 a^b 和 c^d，并判断第一个是更小、更大还是相等(其中 ^ 表示求幂)。

显然，对于非常大的数字，我们无法明确计算这些值。

在这种情况下，最常见的做法是在两边都应用 log 并将 b * log(a) 与 d * log(c) 进行比较。这里的问题是日志是浮点运算，因此我们不能 100% 相信我们的答案(可能有一些值非常接近，并且由于浮点错误，我们得到了错误的答案)。

有解决这个问题的算法吗？我一直在为此搜索内部网，但我只能找到仅适用于特定情况的解决方案(例如，其中一个指数是另一个指数的倍数)，或者以某种方式使用 float (对数、除法)等。

Answer 1

这是两个合二为一的问题:

1. 它们相等吗？

2. 如果不是，哪个更大？

正如 Peter O. 所观察到的，使用提供任意精度小数类型的语言进行构建是最简单的。我将使用 Python 3。

我们假设 a ≤ c 不失一般性(必要时交换)和 b与 d 互质(两者除以最大公约数)。

为了触及问题的核心，我假设 a, c > 0和 b, d ≥ 0 .去除这个假设很乏味但并不困难。

等式检验

在一些简单的情况下 a = 1或 b = 0或 c = 1或 d = 0 .

另外，a^b = c^d 的必要条件是

我。 b ≥ d , 因为否则 b < d , 连同 a ≤ c暗示 a^b < c^d ;

二。 a是 c 的除数，因为我们从 (i) 知道 a^b = c^d是 c^b = c^(b−d) c^d 的除数.

当这些条件成立时，我们可以除以 a^d将问题减少到测试是否a^(b−d) = (c/a)^d .

在 Python 3 中:

def equal_powers(a, b, c, d):
    while True:
        lhs_is_one = a == 1 or b == 0
        rhs_is_one = c == 1 or d == 0
        if lhs_is_one or rhs_is_one:
            return lhs_is_one and rhs_is_one
        if a > c:
            a, b, c, d = c, d, a, b
        if b < d:
            return False
        q, r = divmod(c, a)
        if r != 0:
            return False
        b -= d
        c = q


def test_equal_powers():
    for a in range(1, 25):
        for b in range(25):
            for c in range(1, 25):
                for d in range(25):
                    assert equal_powers(a, b, c, d) == (a ** b == c ** d)


test_equal_powers()

不等式检验

一旦我们确定这两个量不相等，就该找出哪个更大了。 (如果没有相等性测试，这里的代码可能永远运行下去。)

如果您真的要这样做，您应该查阅有关计算初等函数的实际引用资料。我将尝试做最简单的事情。

复习微积分的时间到了。我们有泰勒级数

−log x = (1−x) + (1−x)^2/2 + (1−x)^3/3 + (1−x)^4/4 + ...

要获得下限，请截断序列。为了获得上限，我们可以截断但替换最后一项 (1−x)^n/n与 (1−x)^n/n (1/x) , 因为

(1−x)^n/n (1/x)
= (1−x)^n/n (1 + (1−x) + (1−x)^2 + ...)
= (1−x)^n/n + (1−x)^(n+1)/n + (1−x)^(n+2)/n + ...
> (1−x)^n/n + (1−x)^(n+1)/(n+1) + (1−x)^(n+2)/(n+2) + ...

为了获得良好的收敛速度，我们需要 0.5 ≤ x < 1 ，我们可以通过除以 x 来实现2 的幂。

在 Python 中，我们将实数表示为包含真值的收缩区间的无限生成器。一旦间隔为b log a和 d log c是不相交的，我们可以确定它们的比较方式。

import fractions


def minus(x, y):
    while True:
        x_lo, x_hi = next(x)
        y_lo, y_hi = next(y)
        yield x_lo - y_hi, x_hi - y_lo


def times(b, x):
    for lo, hi in x:
        yield b * lo, b * hi


def restricted_log(a):
    series = 0
    n = 0
    numerator = 1
    while True:
        n += 1
        numerator *= 1 - a
        series += fractions.Fraction(numerator, n)
        yield -(series + fractions.Fraction(numerator * (1 - a), (n + 1) * a)), -series


def log(a):
    n = 0
    while a >= 1:
        a = fractions.Fraction(a, 2)
        n += 1
    return minus(restricted_log(a), times(n, restricted_log(fractions.Fraction(1, 2))))


def less_powers(a, b, c, d):
    lhs = times(b, log(a))
    rhs = times(d, log(c))
    while True:
        lhs_lo, lhs_hi = next(lhs)
        rhs_lo, rhs_hi = next(rhs)
        if lhs_hi < rhs_lo:
            return True
        if rhs_hi < lhs_lo:
            return False


def test_less_powers():
    for a in range(1, 10):
        for b in range(10):
            for c in range(1, 10):
                for d in range(10):
                    if a ** b != c ** d:
                        assert less_powers(a, b, c, d) == (a ** b < c ** d)


test_less_powers()