我可能会建议四种方法来解决您的任务。首先是基于Binary Search。二是基于Newton's Method。第三是基于Shifting n-th Root Algorithm。第四种是我在picture here中描述的Chord-Tangent方法。
二分搜索已经在上面的许多答案中实现。在这里我只是介绍一下我自己对它的看法和它的实现。
作为替代方案,我还实现了优化二进制搜索方法(标记为Opt)。这个方法只是从[hi / 2, hi) 范围开始,如果我们计算第 k 个根,hi 等于 2^(num_bit_length / k)。
牛顿法在这里是新的,因为我看到它没有在其他答案中实现。它通常被认为比二分搜索更快,尽管我自己在下面的代码中的计时没有显示任何加速。因此,这里的方法仅供参考/兴趣。
移位方法比优化的二分搜索方法快 30-50%,如果用 C++ 实现应该更快,因为 C++ 具有快速的 64 位算术,该方法部分使用了该方法。
Chord-Tangent Method 是我on piece of paper 发明的,它受到启发,是对牛顿法的改进。基本上我画了一个Chord 和一个Tangent Line 并找到与水平线y = n 的交点,这两个交点形成了根解(x0, n) 位置的下限和上限近似,其中n = x0 ^ k。这种方法似乎是最快的,而所有其他方法都进行了超过 2000 次迭代,对于 8192 位数字的情况,这种方法只进行了 8 次迭代。所以这种方法比以前的(按速度)换档方法快 200-300x 倍。
作为示例,我生成了 8192 位大小的非常大的随机整数。并用这两种方法测量求三次根的时间。
在test() 函数中你可以看到我通过k = 3 作为根的幂(立方根),你可以通过任何幂而不是3。
Try it online!
def binary_search(begin, end, f, *, niter = [0]):
while begin < end:
niter[0] += 1
mid = (begin + end) >> 1
if f(mid):
begin = mid + 1
else:
end = mid
return begin
def binary_search_kth_root(n, k, *, verbose = False):
# https://en.wikipedia.org/wiki/Binary_search_algorithm
niter = [0]
res = binary_search(0, n + 1, lambda root: root ** k < n, niter = niter)
if verbose:
print('Binary Search iterations:', niter[0])
return res
def binary_search_opt_kth_root(n, k, *, verbose = False):
# https://en.wikipedia.org/wiki/Binary_search_algorithm
niter = [0]
hi = 1 << (n.bit_length() // k - 1)
while hi ** k <= n:
niter[0] += 1
hi <<= 1
res = binary_search(hi >> 1, hi, lambda root: root ** k < n, niter = niter)
if verbose:
print('Binary Search Opt iterations:', niter[0])
return res
def newton_kth_root(n, k, *, verbose = False):
# https://en.wikipedia.org/wiki/Newton%27s_method
f = lambda x: x ** k - n
df = lambda x: k * x ** (k - 1)
x, px, niter = n, 2 * n, [0]
while abs(px - x) > 1:
niter[0] += 1
px = x
x -= f(x) // df(x)
if verbose:
print('Newton Method iterations:', niter[0])
mini, minv = None, None
for i in range(-2, 3):
v = abs(f(x + i))
if minv is None or v < minv:
mini, minv = i, v
return x + mini
def shifting_kth_root(n, k, *, verbose = False):
# https://en.wikipedia.org/wiki/Shifting_nth_root_algorithm
B_bits = 64
r, y = 0, 0
B = 1 << B_bits
Bk_bits = B_bits * k
Bk_mask = (1 << Bk_bits) - 1
niter = [0]
for i in range((n.bit_length() + Bk_bits - 1) // Bk_bits - 1, -1, -1):
alpha = (n >> (i * Bk_bits)) & Bk_mask
B_y = y << B_bits
Bk_yk = (y ** k) << Bk_bits
Bk_r_alpha = (r << Bk_bits) + alpha
Bk_yk_Bk_r_alpha = Bk_yk + Bk_r_alpha
beta = binary_search(1, B, lambda beta: (B_y + beta) ** k <= Bk_yk_Bk_r_alpha, niter = niter) - 1
y, r = B_y + beta, Bk_r_alpha - ((B_y + beta) ** k - Bk_yk)
if verbose:
print('Shifting Method iterations:', niter[0])
return y
def chord_tangent_kth_root(n, k, *, verbose = False):
niter = [0]
hi = 1 << (n.bit_length() // k - 1)
while hi ** k <= n:
niter[0] += 1
hi <<= 1
f = lambda x: x ** k
df = lambda x: k * x ** (k - 1)
# https://i.stack.imgur.com/et9O0.jpg
x_begin, x_end = hi >> 1, hi
y_begin, y_end = f(x_begin), f(x_end)
for icycle in range(1 << 30):
if x_end - x_begin <= 1:
break
niter[0] += 1
if 0: # Do Binary Search step if needed
x_mid = (x_begin + x_end) >> 1
y_mid = f(x_mid)
if y_mid > n:
x_end, y_end = x_mid, y_mid
else:
x_begin, y_begin = x_mid, y_mid
# (y_end - y_begin) / (x_end - x_begin) = (n - y_begin) / (x_n - x_begin) ->
x_n = x_begin + (n - y_begin) * (x_end - x_begin) // (y_end - y_begin)
y_n = f(x_n)
tangent_x = x_n + (n - y_n) // df(x_n) + 1
chord_x = x_n + (n - y_n) * (x_end - x_n) // (y_end - y_n)
assert chord_x <= tangent_x, (chord_x, tangent_x)
x_begin, x_end = chord_x, tangent_x
y_begin, y_end = f(x_begin), f(x_end)
assert y_begin <= n, (chord_x, y_begin, n, n - y_begin)
assert y_end > n, (icycle, tangent_x - binary_search_kth_root(n, k), y_end, n, y_end - n)
if verbose:
print('Chord Tangent Method iterations:', niter[0])
return x_begin
def test():
import random, timeit
nruns = 3
bits = 8192
n = random.randrange(1 << (bits - 1), 1 << bits)
a = binary_search_kth_root(n, 3, verbose = True)
b = binary_search_opt_kth_root(n, 3, verbose = True)
c = newton_kth_root(n, 3, verbose = True)
d = shifting_kth_root(n, 3, verbose = True)
e = chord_tangent_kth_root(n, 3, verbose = True)
assert abs(a - b) <= 0 and abs(a - c) <= 1 and abs(a - d) <= 1 and abs(a - e) <= 1, (a - b, a - c, a - d, a - e)
print()
print('Binary Search timing:', round(timeit.timeit(lambda: binary_search_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
print('Binary Search Opt timing:', round(timeit.timeit(lambda: binary_search_opt_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
print('Newton Method timing:', round(timeit.timeit(lambda: newton_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
print('Shifting Method timing:', round(timeit.timeit(lambda: shifting_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
print('Chord Tangent Method timing:', round(timeit.timeit(lambda: chord_tangent_kth_root(n, 3), number = nruns) / nruns, 3), 'sec')
if __name__ == '__main__':
test()
输出:
Binary Search iterations: 8192
Binary Search Opt iterations: 2732
Newton Method iterations: 9348
Shifting Method iterations: 2752
Chord Tangent Method iterations: 8
Binary Search timing: 0.506 sec
Binary Search Opt timing: 0.05 sec
Newton Method timing: 2.09 sec
Shifting Method timing: 0.03 sec
Chord Tangent Method timing: 0.001 sec