Ускорение вычисления квадратных чисел [дубликат]

Я нашел эту ветку несколько дней назад и повторно написал решение nibot , и, сократив количество итераций пополам и выполняя некоторые другие незначительные улучшения производительности, я смог улучшить производительность с помощью фактор ~ 2.4:

def isqrt(n):
    a = 0 # a is the current answer.
    r = 0 # r is the current remainder.
    for s in reversed(range(0, n.bit_length(), 2)): # Shift n by s bits.
        t = n >> s & 3 # t is the two next most significant bits of n.
        r = r << 2 | t # Increase the remainder as if no new bit is set.
        c = a << 2 | 1 # c is an intermediate value used for comparison.
        b = r >= c     # b is the next bit in the remainder.
        if b: 
            r -= c     # b has been set, so reduce the remainder.
        a = a << 1 | b # Update the answer to include b.
    return (a, r)

Вот результаты из timeit:

>>> timeit('isqrt(12345678901234567890)', setup='from __main__ import isqrt')
8.862877120962366

Затем для сравнения я применил наиболее часто используемый алгоритм с квадратным корнем: Метод Ньютона . Это определение гораздо более компактно.

def isqrt(n):
    x, y = n, n >> 1
    while x > y:
        x, y = y, (y + n//y) >> 1
    return (x, n - x*x)

Оказывается, что даже оптимизированная версия длинномерных квадратных корней медленнее метода Ньютона, занимая примерно в 1,5 раза.

>>> timeit('isqrt(12345678901234567890)', setup='from __main__ import isqrt')
5.74083631898975

Итак, в заключение, если вам нужна быстрая чистая функция квадратного корня Python, посмотрите не дальше, чем приведенная выше.

Изменить: я исправил ошибку в методе Ньютона выше. На моей машине он работает на ~ 10% быстрее, чем user448810. .

Один из вариантов - использовать модуль 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

Попробуйте это условие (без дополнительных вычислений):

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

Похоже, вы можете проверить это следующим образом:

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

Обновление:

Как вы указали выше, сбой при больших значениях n. Для тех, которые выглядят многообещающе, что является адаптацией примера кода C Мартином Гаем @ UKC, июнь 1985 года, для относительно просто выглядящего двоичного численного метода расчета по цифре, упомянутого в статье Википедии . Методы вычисление квадратных корней :

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)

Длинный алгоритм квадратного корня

Оказывается, существует алгоритм вычисления квадратных корней, который вы можете вычислить вручную, что-то вроде длинного деления. Каждая итерация алгоритма производит ровно одну цифру полученного квадратного корня, потребляя две цифры числа, квадратный корень которого вы ищете. В то время как «длинная рука» версия алгоритма указана в десятичной системе, она работает в любой базе, причем бинарность является самой простой для реализации и, возможно, самой быстрой для выполнения (в зависимости от базового представления bignum).

этот алгоритм работает с цифрами по цифрам, он дает точные результаты для произвольно больших совершенных квадратов, а для не-идеальных квадратов может потребоваться столько цифр точности (справа от десятичной точки), сколько желательно.

На сайте «Dr. Math» есть два замечательных отчета, которые объясняют алгоритм:

• Квадратные корни в двоичном файле
• 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)

Вы можете легко изменить эту функцию, чтобы провести дополнительные итерации до вычислить дробную часть квадратного корня. Меня больше всего интересовали вычисления корней больших совершенных квадратов.

Я не уверен, как это сравнивается с алгоритмом «целочисленный метод Ньютона». Я подозреваю, что метод Ньютона быстрее, поскольку он может в принципе генерировать несколько бит решения на одной итерации, тогда как алгоритм «длинной руки» генерирует ровно один бит решения на итерацию.

Source repo: https://gist.github.com/tobin/11233492

Я сравнил различные методы, приведенные здесь, с циклом:

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