python中的整数平方根

Question

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

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*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与int的兼容:

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

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

请参阅我的其他答案,讨论此方法相对于此问题的其他答案的性能.

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

Answer 3

更新： Python 3.8math.isqrt在标准库中有一个功能！

我在这里用小（0…2 ²²）和大（2 ⁵⁰⁰⁰¹）输入对每个（正确）函数进行基准测试。在这两种情况下，显而易见的赢家gmpy2.isqrt由mathmandan提出，其次是Python 3.8 math.isqrt，其次是NPE链接的ActiveState食谱。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 ms
• int(gmpy2.isqrt())*：小0.3 µs，大0.07 ms
• Python 3.8 math.isqrt：小0.13 µs，大0.9 ms
• ActiveState（如上优化）：小0.6 µs，大17.0 ms
• ActiveState（NPE）：小1.0 µs，大17.3 ms
• castlebravo长手：小4 µs，大80 ms
• mathmandan改进：小2.7 µs，大120 ms
• martineau（经过此更正）：小2.3 µs，大140 ms
• nibot：小8 µs，大1000 ms
• mathmandan：小1.8 µs，大2200 ms
• castlebravo牛顿法：小1.5 µs，大19000 ms
• user448810：小1.4 µs，大20000 ms

（*由于gmpy2.isqrt返回的gmpy2.mpz对象的行为主要是但不完全类似于int，因此您可能需要将其转换回来int用于某些用途。）

Answer 4

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

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位),这似乎是:

• 比用户448810描述的牛顿方法实现更快.
• 比我的其他答案中的gmpy2内置方法慢得多.
• 与nibot 描述的Longhand Square Root相比,但有点慢.

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

长手平方根算法

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

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

"数学博士"网站上有两篇很好的文章解释了这个算法:

• 二进制中的平方根
• 手绘广场根

这是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 6

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

Python 的默认math库有一个整数平方根函数：

\n

\n

math.isqrt(n)

\n

返回非负整数n的整数平方根。这是n的精确平方根的下限，或者等效于a\xc2\xb2 \xe2\x89\xa4 n 的最大整数 a 。

\n

\n