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が完全な平方根かどうかを確認します。
このアルゴリズムと、平方根を計算するための他の 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
これは単なる二分探索です。平方根を超えない最大の 2 の累乗に値mを初期化し、結果が平方根を超えないようにしながら、小さいビットごとに設定できるかどうかを確認します。(一度に 1 ビットずつ、降順にチェックしてください。)
かなり大きな値n(たとえば、約10**6000、または約20000ビット) の場合、次のようになります。
gmpy2の組み込みメソッドよりもはるかに遅いです。これらのアプローチはすべてこのサイズの入力で成功しますが、私のマシンでは、この関数には約 1.5 秒かかりますが、@ Nibot の場合は約 0.9 秒、@ user448810 の場合は約 19 秒かかり、gmpy2 組み込みメソッドは 1 ミリ秒未満かかります。 (!)。例：
>>> 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。
i_sqrt編集：乗算を使用しないように関数を書き直すことができることを指摘してくれた@Greggoに感謝します。これにより、パフォーマンスが大幅に向上します。
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
構造上、 のkth ビットm << 1は設定されていないため、ビットごとの OR を使用して の加算を実装できることに注意してください(m<<1) + (1<<k)。最終的に私は(2*m*(2**k) + 2**(2*k))として書いた(((m<<1) | (1<<k)) << k)ので、それは 3 つのシフトと 1 つのビット単位の OR (その後に get への減算が続きます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
1つのオプションは、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。それらについては、ウィキペディアの記事「平方根の計算方法」で言及されている比較的単純に見える 2 進数の桁ごとの計算方法について、1985 年 6 月に Martin Guy @ UKC によって作成されたサンプル C コードの適応である次のものが有望に見えます。:
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)
あなたの関数は大きな入力に対して失敗します:
In [26]: isqrt((10**100+1)**2)
ValueError: input was not a perfect square
ActiveState サイトには、整数演算のみを使用するため、より信頼性の高いレシピが掲載されています。これは、以前の StackOverflow の質問:自分の平方根関数を書くことに基づいています。
ここに示すさまざまな方法をループと比較しました。
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
float はコンピューター上で正確に表現できません。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')
この条件を試してください（追加の計算はありません）：
def isqrt(n):
i = math.sqrt(n)
if i != int(i):
raise ValueError('input was not a perfect square')
return i
int（末尾にゼロが付いていない）を返す必要がある場合はfloat、2番目の変数を割り当てるか、2int(i)回計算します。