整数平方根在python整数平方根在python(Integer square root in py

2019-05-06 11:36发布

站内文章 / 移动开发

33 0

男人必须洒脱

女 | 书童

私信

是否有一个整数平方根蟒蛇的地方，或在标准库？我想这是准确的（即返回一个整数），和树皮如果没有解决方案。

在我推出我自己的一个幼稚的时刻：

def isqrt(n):
    i = int(math.sqrt(n) + 0.5)
    if i**2 == n:
        return i
    raise ValueError('input was not a perfect square')

但它的丑陋，我真的不相信它的大整数。我可以通过迭代广场和放弃，如果我已经超过了价值，但我相信它会有点慢，做这样的事情。此外，我想我可能会被重新发明轮子，这样的事情必须肯定蟒蛇已经存在...

Answer 1:

牛顿的方法很完善的整数：

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是一个完美的正方形。

我讨论这个算法，和其他三种算法来计算平方根，在我的博客。

Answer 2:

很抱歉的响应到很晚; 我只是偶然到这个页面。如果有人访问这个页面在未来，Python模块gmpy2设计有非常大的投入工作，除其他事项外包括一个整数平方根函数。

例：

>>> import gmpy2
>>> gmpy2.isqrt((10**100+1)**2)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001L)
>>> gmpy2.isqrt((10**100+1)**2 - 1)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000L)

当然，一切都将有“MPZ”的标签，但MPZ的是与诠释的兼容：

>>> gmpy2.mpz(3)*4
mpz(12)

>>> int(gmpy2.mpz(12))
12

见我的其他答案的表现相对这种方法的的一些其他的答案对这个问题的讨论。

下载： https://code.google.com/p/gmpy/

Answer 3:

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

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

比牛顿方法实现更快的通过user448810描述。
多，比慢得多gmpy2在内置的方法我其他的答案。
相媲美，但比略慢，在手写平方根由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

该功能可以很容易地一概而论，但不是因为我还没有经历那样精确的初始猜测的并不像不错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

Answer 4:

长手平方根算法

事实证明，没有计算平方根，您可以通过手工计算，像长划分的算法。该算法的每一次迭代产生所产生的平方根恰好有一个数字，而消耗它的平方根你追求数量的两位数。虽然在小数指定的算法的“长手”的版本，它可以在任何基地，拥有二进制感最简单的实现，或许以最快的速度执行（具体取决于基础BIGNUM表示）。

因为此算法上的数字操作的数字逐位，根据需要它产生用于任意大的完全平方准确的结果，而对于非理想平方，能产生尽可能多的位精度（至小数点的右边）。

有关于“数学博士”网站两个漂亮的writeups解释算法：

在二进制平方根
手写平方根

下面是用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)

你可以很容易地修改此功能进行更多的迭代计算平方根的小数部分。我最感兴趣的是计算的大完全平方根。

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

来源回购： https://gist.github.com/tobin/11233492

Answer 5:

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

Answer 6:

我在这里基准上都小（0 ... 2 ^月22日）和大（2 ^50001）每投入（正确的）功能。在这两种情况下明显的赢家是gmpy2.isqrt通过mathmandan建议排在首位，其次是由NPE相连的ActiveState食谱排在第二。该ActiveState的配方有一堆师可以通过轮班替换，这使得它快一点（但仍落后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微秒小，0.07毫秒大
int(gmpy2.isqrt()) *：0.3微秒小，0.07毫秒大
ActiveState的（如以上优化的）：0.6微秒小，17.0毫秒大
的ActiveState（NPE）：1.0微秒小，17.3毫秒大
castlebravo长手：4微秒小，80毫秒的大
mathmandan改善：2.7微秒小，120毫秒大
蒂诺（具有这种校正）：2.3微秒小，140毫秒大
nibot ：8微秒小，1000毫秒大
mathmandan ：1.8微秒小，2200毫秒大
castlebravo牛顿法：1.5微秒小，19000毫秒大
user448810 ：1.4微秒小，20000毫秒大

（*由于gmpy2.isqrt返回gmpy2.mpz对象，其行为大多但不完全像一个int ，你可能需要将其转换回int对于某些用途。）

Answer 7:

好像你可以检查是这样的：

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

更新：

正如你指出了上述故障为大型值n 。对于那些以下看起来很有希望，这是例如C代码的适配，由Martin盖伊@ 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)

Answer 8:

对大投入的功能失败：

In [26]: isqrt((10**100+1)**2)

ValueError: input was not a perfect square

有一个在ActiveState的网站食谱应该希望更可靠，因为它仅使用整数数学。它是基于一个较早的StackOverflow问题：编写自己的平方根函数

Answer 9:

花车不能准确代表的计算机上。您可以测试所需的接近设定小量Python的浮点数的精度范围内的较小值。

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')

Answer 10:

我比较了循环这里给出的不同的方法：

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

Answer 11:

试试这个条件（没有额外计算）：

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)的两倍。

文章来源: Integer square root in python

标签： python math integer sqrt

男人必须洒脱

女 | 书童

私信

收藏的人(0)

Ta的文章更多文章

0条评论

还没有人评论过~