python 或标准库中是否有整数平方根？我希望它是准确的（即返回一个整数），如果没有解决方案就吠叫。
此刻我推出了我自己的天真之一：
def isqrt(n):
i = int(math.sqrt(n) + 0.5)
if i**2 == n:
return i
raise ValueError('input was not a perfect square')
但它很丑陋，我真的不相信它对于大整数。如果我超过了该值，我可以遍历方块并放弃，但我认为做这样的事情会有点慢。另外我想我可能会重新发明轮子，像这样的东西肯定已经存在于 python 中......
牛顿法对整数非常有效：
def isqrt(n):
x = n
y = (x + 1) // 2
while y < x:
x = y
y = (x + n // x) // 2
return x
这将返回x * x不超过n的最大整数x。如果要检查结果是否正好是平方根，只需执行乘法以检查n是否为完美平方。
我在我的博客上讨论了这个算法，以及其他三种计算平方根的算法。
更新： Python 3.8math.isqrt在标准库中有一个函数！
我在小 (0…2 22 ) 和大 (2 50001 ) 输入上对每个（正确的）函数进行了基准测试。在这两种情况下，最明显的赢家是gmpy2.isqrtmathmandan 排名第一，其次是 Python 3.8 math.isqrt，其次是由 NPE 链接的 ActiveState recipe排名第三。ActiveState 配方有一堆可以用班次代替的部门，这使它更快一点（但仍然落后于本机函数）：
def isqrt(n):
if n > 0:
x = 1 << (n.bit_length() + 1 >> 1)
while True:
y = (x + n // x) >> 1
if y >= x:
return x
x = y
elif n == 0:
return 0
else:
raise ValueError("square root not defined for negative numbers")
基准测试结果：
gmpy2.isqrt()(mathmandan)：小 0.08 µs，大 0.07 msint(gmpy2.isqrt())*：小 0.3 µs，大 0.07 msmath.isqrt：小 0.13 µs，大 0.9 ms（* 由于gmpy2.isqrt返回一个gmpy2.mpz对象，它的行为主要但不完全像 an int，您可能需要将其转换回 anint以用于某些用途。）
很抱歉回复太晚了；我只是偶然发现了这个页面。如果将来有人访问此页面，python 模块 gmpy2 旨在处理非常大的输入，其中包括整数平方根函数。
例子：
>>> import gmpy2
>>> gmpy2.isqrt((10**100+1)**2)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001L)
>>> gmpy2.isqrt((10**100+1)**2 - 1)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000L)
当然，一切都会有“mpz”标签，但 mpz 与 int 兼容：
>>> gmpy2.mpz(3)*4
mpz(12)
>>> int(gmpy2.mpz(12))
12
有关此方法相对于该问题的其他一些答案的性能的讨论，请参阅我的其他答案。
长手平方根算法
事实证明，有一种计算平方根的算法，你可以手动计算，比如长除法。算法的每次迭代都会产生结果平方根的一个数字，同时消耗您寻找平方根的数字的两个数字。虽然算法的“长手”版本以十进制指定，但它适用于任何基础，二进制最容易实现，也许执行速度最快（取决于底层的 bignum 表示）。
因为该算法对数字逐位进行运算，所以它可以为任意大的完全平方产生精确的结果，而对于非完全平方，它可以产生任意多的精度（小数点右侧）。
“Dr. Math”网站上有两篇很好的文章解释了该算法：
这是 Python 中的一个实现：
def exact_sqrt(x):
"""Calculate the square root of an arbitrarily large integer.
The result of exact_sqrt(x) is a tuple (a, r) such that a**2 + r = x, where
a is the largest integer such that a**2 <= x, and r is the "remainder". If
x is a perfect square, then r will be zero.
The algorithm used is the "long-hand square root" algorithm, as described at
http://mathforum.org/library/drmath/view/52656.html
Tobin Fricke 2014-04-23
Max Planck Institute for Gravitational Physics
Hannover, Germany
"""
N = 0 # Problem so far
a = 0 # Solution so far
# We'll process the number two bits at a time, starting at the MSB
L = x.bit_length()
L += (L % 2) # Round up to the next even number
for i in xrange(L, -1, -1):
# Get the next group of two bits
n = (x >> (2*i)) & 0b11
# Check whether we can reduce the remainder
if ((N - a*a) << 2) + n >= (a<<2) + 1:
b = 1
else:
b = 0
a = (a << 1) | b # Concatenate the next bit of the solution
N = (N << 2) | n # Concatenate the next bit of the problem
return (a, N-a*a)
您可以轻松修改此函数以进行额外的迭代来计算平方根的小数部分。我最感兴趣的是计算大完全平方的根。
我不确定这与“整数牛顿法”算法相比如何。我怀疑牛顿的方法更快，因为它原则上可以在一次迭代中生成多位解，而“长手”算法每次迭代只生成一位解。
这是一个非常简单的实现：
def i_sqrt(n):
i = n.bit_length() >> 1 # i = floor( (1 + floor(log_2(n))) / 2 )
m = 1 << i # m = 2^i
#
# Fact: (2^(i + 1))^2 > n, so m has at least as many bits
# as the floor of the square root of n.
#
# Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2)
# >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED.
#
while m*m > n:
m >>= 1
i -= 1
for k in xrange(i-1, -1, -1):
x = m | (1 << k)
if x*x <= n:
m = x
return m
这只是一个二分搜索。将值初始化为m不超过平方根的 2 的最大幂，然后检查是否可以设置每个较小的位，同时保持结果不大于平方根。（按降序逐个检查位。）
对于相当大的值n（例如，大约10**6000或大约20000位），这似乎是：
gmpy2中的内置方法慢得多。所有这些方法都在这种大小的输入上成功，但在我的机器上，这个函数大约需要 1.5 秒，而 @Nibot 大约需要 0.9 秒，@user448810 大约需要 19 秒，而 gmpy2 内置方法需要不到一毫秒（！）。例子：
>>> import random
>>> import timeit
>>> import gmpy2
>>> r = random.getrandbits
>>> t = timeit.timeit
>>> t('i_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. # This function
1.5102493192883117
>>> t('exact_sqrt(r(20000))', 'from __main__ import *', number = 5)/5. # Nibot
0.8952787937686366
>>> t('isqrt(r(20000))', 'from __main__ import *', number = 5)/5. # user448810
19.326695976676184
>>> t('gmpy2.isqrt(r(20000))', 'from __main__ import *', number = 5)/5. # gmpy2
0.0003599147067689046
>>> all(i_sqrt(n)==isqrt(n)==exact_sqrt(n)[0]==int(gmpy2.isqrt(n)) for n in (r(1500) for i in xrange(1500)))
True
这个函数可以很容易地推广，虽然它不是很好，因为我没有对 的初始猜测非常精确m：
def i_root(num, root, report_exactness = True):
i = num.bit_length() / root
m = 1 << i
while m ** root < num:
m <<= 1
i += 1
while m ** root > num:
m >>= 1
i -= 1
for k in xrange(i-1, -1, -1):
x = m | (1 << k)
if x ** root <= num:
m = x
if report_exactness:
return m, m ** root == num
return m
但是，请注意，gmpy2也有一个i_root方法。
事实上，这种方法可以适用于任何（非负的、递增的）函数f来确定“整数倒数f”。但是，要选择一个有效的初始值，m您仍然需要了解有关f.
编辑：感谢@Greggo 指出i_sqrt可以重写该函数以避免使用任何乘法。这产生了令人印象深刻的性能提升！
def improved_i_sqrt(n):
assert n >= 0
if n == 0:
return 0
i = n.bit_length() >> 1 # i = floor( (1 + floor(log_2(n))) / 2 )
m = 1 << i # m = 2^i
#
# Fact: (2^(i + 1))^2 > n, so m has at least as many bits
# as the floor of the square root of n.
#
# Proof: (2^(i+1))^2 = 2^(2i + 2) >= 2^(floor(log_2(n)) + 2)
# >= 2^(ceil(log_2(n) + 1) >= 2^(log_2(n) + 1) > 2^(log_2(n)) = n. QED.
#
while (m << i) > n: # (m<<i) = m*(2^i) = m*m
m >>= 1
i -= 1
d = n - (m << i) # d = n-m^2
for k in xrange(i-1, -1, -1):
j = 1 << k
new_diff = d - (((m<<1) | j) << k) # n-(m+2^k)^2 = n-m^2-2*m*2^k-2^(2k)
if new_diff >= 0:
d = new_diff
m |= j
return m
请注意，通过构造，未设置 的k第 位m << 1，因此可以使用按位或来实现 的加法(m<<1) + (1<<k)。最终我(2*m*(2**k) + 2**(2*k))写成(((m<<1) | (1<<k)) << k)，所以它是三个班次和一个按位或（然后是减法得到new_diff）。也许还有更有效的方法来获得这个？无论如何，这比乘法要好得多m*m！与上面的比较：
>>> t('improved_i_sqrt(r(20000))', 'from __main__ import *', number = 5)/5.
0.10908999762373242
>>> all(improved_i_sqrt(n) == i_sqrt(n) for n in xrange(10**6))
True
一种选择是使用该decimal模块，并以足够精确的浮点数进行：
import decimal
def isqrt(n):
nd = decimal.Decimal(n)
with decimal.localcontext() as ctx:
ctx.prec = n.bit_length()
i = int(nd.sqrt())
if i**2 != n:
raise ValueError('input was not a perfect square')
return i
我认为应该工作：
>>> isqrt(1)
1
>>> isqrt(7**14) == 7**7
True
>>> isqrt(11**1000) == 11**500
True
>>> isqrt(11**1000+1)
Traceback (most recent call last):
File "<ipython-input-121-e80953fb4d8e>", line 1, in <module>
isqrt(11**1000+1)
File "<ipython-input-100-dd91f704e2bd>", line 10, in isqrt
raise ValueError('input was not a perfect square')
ValueError: input was not a perfect square
好像你可以这样检查：
if int(math.sqrt(n))**2 == n:
print n, 'is a perfect square'
更新：
正如您所指出的，对于较大的n. 对于那些看起来很有希望的人，这是对示例 C 代码的改编，Martin Guy @ UKC，1985 年 6 月，用于维基百科文章计算平方根的方法中提到的相对简单的二进制数字逐位计算方法：
from math import ceil, log
def isqrt(n):
res = 0
bit = 4**int(ceil(log(n, 4))) if n else 0 # smallest power of 4 >= the argument
while bit:
if n >= res + bit:
n -= res + bit
res = (res >> 1) + bit
else:
res >>= 1
bit >>= 2
return res
if __name__ == '__main__':
from math import sqrt # for comparison purposes
for i in range(17)+[2**53, (10**100+1)**2]:
is_perfect_sq = isqrt(i)**2 == i
print '{:21,d}: math.sqrt={:12,.7G}, isqrt={:10,d} {}'.format(
i, sqrt(i), isqrt(i), '(perfect square)' if is_perfect_sq else '')
输出：
0: math.sqrt= 0, isqrt= 0 (perfect square)
1: math.sqrt= 1, isqrt= 1 (perfect square)
2: math.sqrt= 1.414214, isqrt= 1
3: math.sqrt= 1.732051, isqrt= 1
4: math.sqrt= 2, isqrt= 2 (perfect square)
5: math.sqrt= 2.236068, isqrt= 2
6: math.sqrt= 2.44949, isqrt= 2
7: math.sqrt= 2.645751, isqrt= 2
8: math.sqrt= 2.828427, isqrt= 2
9: math.sqrt= 3, isqrt= 3 (perfect square)
10: math.sqrt= 3.162278, isqrt= 3
11: math.sqrt= 3.316625, isqrt= 3
12: math.sqrt= 3.464102, isqrt= 3
13: math.sqrt= 3.605551, isqrt= 3
14: math.sqrt= 3.741657, isqrt= 3
15: math.sqrt= 3.872983, isqrt= 3
16: math.sqrt= 4, isqrt= 4 (perfect square)
9,007,199,254,740,992: math.sqrt=9.490627E+07, isqrt=94,906,265
100,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,020,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001: math.sqrt= 1E+100, isqrt=10,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,000,001 (perfect square)
Python 的默认math库有一个整数平方根函数：
math.isqrt(n)返回非负整数n的整数平方根。这是n的精确平方根的下限，或者等效地是满足a² ≤ n的最大整数 a 。
下面的脚本提取整数平方根。它不使用除法，只使用位移，因此速度非常快。它在反平方根上使用牛顿法，该技术由Quake III Arena出名，在维基百科文章Fast inverse square root中提到。
算法的计算策略s = sqrt(Y)如下。
B = 2**(2*k)对于某个整数 k。我们想要找到 X，其中 x = X/B，并且 x = 1 / sqrt(y)。s = X*Y/(2**(3*k))。我们实际上不创建分数或执行任何除法。所有的算术都是用整数完成的，我们使用位移来除以 B 的各种幂。
范围缩减让我们找到了一个很好的初始近似值来提供给牛顿方法。这是区间 [1/4, 1) 中平方根反比的 2 次极小极大多项式逼近的一个版本：
（对不起，我在这里颠倒了 x & y 的含义，以符合通常的约定）。该近似值的最大误差约为 0.0355 ~= 1/28。这是显示错误的图表：
使用这个 poly，我们的初始 x 以至少 4 或 5 位的精度开始。牛顿方法的每一轮都将精度翻倍，因此如果我们想要的话，不需要很多轮就可以获得数千位。
""" Integer square root
Uses no divisions, only shifts
"Quake" style algorithm,
i.e., Newton's method for 1 / sqrt(y)
Uses a quadratic minimax polynomial for the first approximation
Written by PM 2Ring 2022.01.23
"""
def int_sqrt(y):
if y < 0:
raise ValueError("int_sqrt arg must be >= 0, not %s" % y)
if y < 2:
return y
# print("\n*", y, "*")
# Range reduction.
# Find k such that 1/4 <= y/b < 1, where b = 2 ** (k*2)
j = y.bit_length()
# Round k*2 up to the next even number
k2 = j + (j & 1)
# k and some useful multiples
k = k2 >> 1
k3 = k2 + k
k6 = k3 << 1
kd = k6 + 1
# b cubed
b3 = 1 << k6
# Minimax approximation: x/b ~= 1 / sqrt(y/b)
x = (((463 * y * y) >> k2) - (896 * y) + (698 << k2)) >> 8
# print(" ", x, h)
# Newton's method for 1 / sqrt(y/b)
epsilon = 1 << k
for i in range(1, 99):
dx = x * (b3 - y * x * x) >> kd
x += dx
# print(f" {i}: {x} {dx}")
if abs(dx) <= epsilon:
break
# s == sqrt(y)
s = x * y >> k3
# Adjust if too low
ss = s + 1
return ss if ss * ss <= y else s
def test(lo, hi, step=1):
for y in range(lo, hi, step):
s = int_sqrt(y)
ss = s + 1
s2, ss2 = s * s, ss * ss
assert s2 <= y < ss2, (y, s2, ss2)
print("ok")
test(0, 100000, 1)
这段代码肯定比and慢。它的目的只是为了说明算法。看看如果用 C 实现它会有多快会很有趣......math.isqrtdecimal.Decimal.sqrt
这是一个实时版本，在 SageMathCell 服务器上运行。设置hi<= 0 以计算并显示 中设置的单个值的结果lo。您可以在输入框中输入表达式，例如设置为hi0 和lo获取。2 * 10**100sqrt(2) * 10**50
您的功能因大量输入而失败：
In [26]: isqrt((10**100+1)**2)
ValueError: input was not a perfect square
ActiveState 网站上有一个配方，希望它更可靠，因为它只使用整数数学。它基于早期的 StackOverflow 问题：编写您自己的平方根函数
浮点数不能在计算机上精确表示。您可以在 python 的浮点数精度范围内将所需的接近度设置 epsilon 测试为一个小值。
def isqrt(n):
epsilon = .00000000001
i = int(n**.5 + 0.5)
if abs(i**2 - n) < epsilon:
return i
raise ValueError('input was not a perfect square')
我用循环比较了这里给出的不同方法：
for i in range (1000000): # 700 msec
r=int(123456781234567**0.5+0.5)
if r**2==123456781234567:rr=r
else:rr=-1
发现这个是最快的，不需要数学导入。很长可能会失败，但看看这个
15241576832799734552675677489**0.5 = 123456781234567.0
试试这个条件（没有额外的计算）：
def isqrt(n):
i = math.sqrt(n)
if i != int(i):
raise ValueError('input was not a perfect square')
return i
如果您需要它返回一个int（不是float带有尾随零的 a），则分配第二个变量或计算int(i)两次。