python中的整数平方根-Java 学习之路

在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 ......

11 回答

15
牛顿方法在整数上运行得非常好：
```
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是否为完美平方 .

我在my blog讨论了这个算法以及另外三种计算平方根的算法 .
回复于 2024-05-09T01:44:50+08:00
6
很抱歉很晚才回复;我只是偶然发现了这个页面 . 如果将来有人访问此页面，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
```
有关此方法相对于此问题的其他答案的性能的讨论，请参阅my other answer .

下载：https://code.google.com/p/gmpy/
回复于 2024-05-09T01:44:50+08:00
2
Long-hand square root algorithm

事实证明，有一种计算平方根的算法，您可以手动计算，例如长除法 . 算法的每次迭代都会生成所得到的平方根的正好一位数，同时消耗您寻找的平方根的数字的两位数 . 虽然算法的“长手”版本以十进制形式指定，但它适用于任何基础，二进制最简单的实现，也许执行速度最快（取决于底层的bignum表示） .

因为该算法逐个数字地对数字进行操作，所以它为任意大的正方形产生精确的结果，对于非完美的正方形，可以根据需要产生尽可能多的精度数字（在小数点的右侧） .

“数学博士”网站上有两篇很好的文章解释了这个算法：
- Square Roots in Binary
- Longhand Square Roots
这是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)
```
您可以轻松修改此函数以执行其他迭代来计算平方根的小数部分 . 我最感兴趣的是计算大型完美正方形的根 .

我不确定这与“整数牛顿方法”算法相比如何 . 我怀疑Newton的方法更快，因为它原则上可以在一次迭代中生成解决方案的多个位，而“长手”算法每次迭代生成解决方案的一个位 .

来源回购：https://gist.github.com/tobin/11233492
回复于 2024-05-09T01:44:50+08:00

这是一个非常简单的实现：

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 位），这似乎是：

比牛顿的方法实现更快described by user448810 .
比my other answer中的 gmpy2 内置方法慢得多 .
与Longhand Square Root described by nibot相比，但有点慢 .

所有这些方法都在这个尺寸的输入上成功，但在我的机器上，这个功能需要大约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

这个函数可以很容易地推广，尽管's not quite as nice because I don'对于 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 .

Edit: Thanks to @Greggo for pointing out that the i_sqrt function can be rewritten to avoid using any multiplications. 这会带来令人印象深刻的性能提升！

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

注意，通过构造， m << 1 的 k 未设置，所以按位 - 或者可以用来实现 (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

回复于 2024-05-09T01:44:50+08:00

-2

一种选择是使用 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

回复于 2024-05-09T01:44:50+08:00

-4
我在这里对小（0 ... 222）和大（250001）输入的每个（正确）函数进行了基准测试 . 两个案例中明显的获胜者首先是gmpy2.isqrt suggested by mathmandan，其次是ActiveState recipe linked by NPE . ActiveState配方有一堆可以由shift替换的分区，这使得它更快（但仍落后于 gmpy2.isqrt ）：
```
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 ms大
- int(gmpy2.isqrt()) *：0.3μs小，0.07 ms大
- ActiveState (optimized as above) ：0.6μs小，17.0 ms大
- ActiveState (npe)：1.0μs小，17.3 ms大
- castlebravo long-hand：4μs小，80 ms大
- mathmandan improved：2.7μs小，120 ms大
- martineau（带this correction）：2.3μs小，140 ms大
- nibot：8μs小，1000 ms大
- mathmandan：1.8μs小，2200 ms大
- castlebravo Newton’s method：1.5μs小，19000 ms大
- user448810：1.4μs小，20000 ms大
（*由于 gmpy2.isqrt 返回 gmpy2.mpz 对象，其行为大多数但不完全像 int ，您可能需要将其转换回 int 用于某些用途 . ）
回复于 2024-05-09T01:44:50+08:00

好像你可以像这样检查：

if int(math.sqrt(n))**2 == n:
    print n, 'is a perfect square'

更新：

正如您所指出的那样，上面的 n 值很大 . 对于那些看起来很有希望的人，这是1985年6月Martin Guy @ UKC的示例C代码的改编，用于维基百科文章_2474957中提到的相对简单的二进制数字逐位计算方法：

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)

回复于 2024-05-09T01:44:50+08:00

-2
您的功能因大输入而失败：
```
In [26]: isqrt((10**100+1)**2)

ValueError: input was not a perfect square
```
有一个recipe on the ActiveState site应该更可靠，因为它只使用整数数学 . 它基于早期的StackOverflow问题：Writing your own square root function
回复于 2024-05-09T01:44:50+08:00
70
浮动无法在计算机上精确显示 . 您可以在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')
```
回复于 2024-05-09T01:44:50+08:00
6
我将这里给出的不同方法与循环进行了比较：
```
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
```
回复于 2024-05-09T01:44:50+08:00
5
尝试这种情况（无需额外计算）：
```
def isqrt(n):
  i = math.sqrt(n)
  if i != int(i):
    raise ValueError('input was not a perfect square')  
  return i
```
如果你需要它返回一个 int （不是带有尾随零的 float ），那么要么分配第二个变量，要么计算 int(i) 两次 .
回复于 2024-05-09T01:44:50+08:00

python中的整数平方根

11 回答

相关问题