Racine carrée entière en python
y a-t-il une racine carrée entière quelque part dans python, ou dans les bibliothèques standard? Je veux qu'il soit exact (c.-à-d. retourner un entier), et l'écorce s'il n'y a pas de solution.
au moment où j'ai roulé mon propre naïf:
def isqrt(n):
i = int(math.sqrt(n) + 0.5)
if i**2 == n:
return i
raise ValueError('input was not a perfect square')
mais c'est laid et je n'y crois pas vraiment pour les grands entiers. J'ai pu parcourir les places et si j'ai dépassé la valeur, mais je présume qu'il est un peu lent à faire quelque chose comme ça. Aussi, Je je suppose que je serais probablement en train de réinventer la roue, quelque chose comme ça doit sûrement déjà exister en python...
11 réponses
la méthode de Newton fonctionne parfaitement sur les entiers:
def isqrt(n):
x = n
y = (x + 1) // 2
while y < x:
x = y
y = (x + n // x) // 2
return x
renvoie le plus grand entier x pour qui x * x ne pas dépasser n . Si vous voulez vérifier si le résultat est exactement la racine carrée, simplement effectuer la multiplication pour vérifier si n est un carré parfait.
je discute de cet algorithme, et trois autres algorithmes pour calculer des racines carrées, à mon blog .
désolé pour la réponse très tardive, je viens de tomber sur cette page. Dans le cas où quelqu'un visite cette page à l'avenir, le module python gmpy2 est conçu pour fonctionner avec de très grandes entrées, et inclut entre autres une fonction racine carrée entière.
exemple:
>>> import gmpy2
>>> gmpy2.isqrt((10**100+1)**2)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001L)
>>> gmpy2.isqrt((10**100+1)**2 - 1)
mpz(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000L)
accordé, tout aura l'étiquette "mpz", mais mpz sont compatibles avec int's:
>>> gmpy2.mpz(3)*4
mpz(12)
>>> int(gmpy2.mpz(12))
12
voir mon autre réponse pour une analyse du rendement de cette méthode par rapport à d'autres réponses à cette question.
télécharger: https://code.google.com/p/gmpy/
la main à la racine carrée de l'algorithme de
il s'avère qu'il existe un algorithme pour calculer les racines carrées que vous pouvez calculer à la main, quelque chose comme la division longue. Chaque itération de l'algorithme produit un chiffre exactement de la racine carrée tout en consommant deux chiffres du nombre dont la racine carrée que vous cherchez. Alors que la version "longue main" de l'algorithme est spécifié en décimal, il fonctionne dans n'importe quelle base, avec être binaire la plus simple à mettre en œuvre et peut-être la plus rapide à exécuter (selon la représentation de bignum sous-jacente).
parce que cet algorithme fonctionne sur les nombres chiffre par chiffre, il produit des résultats exacts pour les carrés parfaits arbitrairement grands, et pour les non-carrés parfaits, peut produire autant de chiffres de précision (à la droite de la place décimale) que désiré.
il y a deux belles Écritures sur le site" Dr. Math "qui expliquent l'algorithme:
et voici une implémentation en 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)
, Vous pouvez facilement modifier cette fonction pour effectuer des itérations pour calculer la partie fractionnaire de la racine carrée. J'étais le plus intéressé par des racines de calcul de carrés parfaits grands.
Je ne suis pas sûr de comparer cela à l'algorithme de "méthode de Newton entier". Je soupçonne que la méthode de Newton est plus rapide, car elle peut en principe générer plusieurs bits de la solution en une itération, tandis que l'algorithme "main longue" génère exactement un bit de la solution par itération.
source repo: https://gist.github.com/tobin/11233492
Voici une mise en œuvre très simple:
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
ce n'est qu'une recherche binaire. Initialiser la valeur m pour être la plus grande puissance de 2 qui ne dépasse pas la racine carrée, puis vérifier si chaque plus petit bit peut être réglé tout en gardant le résultat pas plus grand que la racine carrée. (Vérifier les bits un à un, dans l'ordre décroissant.)
pour des valeurs raisonnablement élevées de n (disons, autour de 10**6000 , ou autour de 20000 bits), ce qui semble être:
- plus rapide que L'implémentation de la méthode de Newton décrit par l'utilisateur 448810 .
- beaucoup, Beaucoup plus lent que le
gmpy2méthode intégrée dans mon autre réponse . - Comparable à, mais un peu plus lent que, la racine carrée de Longhand décrit par nibot .
toutes ces approches réussissent sur des entrées de cette taille, mais sur ma machine, Cette fonction prend environ 1,5 seconde, alors que @Nibot prend environ 0,9 seconde, @user448810 prend environ 19 secondes, et la méthode intégrée gmpy2 prend moins d'une milliseconde(!). Exemple:
>>> 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
cette fonction peut être généralisée facilement, bien qu'elle ne soit pas tout à fait aussi agréable parce que je n'ai pas tout à fait aussi précis d'une supposition initiale pour 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
cependant, notez que gmpy2 a aussi une méthode i_root .
en fait, cette méthode pourrait être adaptée et appliquée à n'importe quelle fonction (non négative, croissante) f pour déterminer un "entier inverse de f ". Cependant, pour choisir une valeur initiale efficace de m , vous voulez tout de même savoir quelque chose sur f .
La fonction i_sqrt peut être réécrite pour éviter toute multiplication. cela donne un coup de pouce de performance impressionnant!
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
noter que par construction, le k TH bit de m << 1 n'est pas réglé, donc bitwise-ou peut être utilisé pour mettre en œuvre l'ajout de (m<<1) + (1<<k) . Finalement j'ai (2*m*(2**k) + 2**(2*k)) écrit comme (((m<<1) | (1<<k)) << k) , donc c'est trois quarts et un bitwise-or (suivi d'une soustraction pour obtenir new_diff ). Peut-être qu'il est encore un moyen plus efficace de faire cela? Quoi qu'il en soit, c'est bien mieux que de multiplier m*m ! Comparer avec ci-dessus:
>>> 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
une option serait d'utiliser le module decimal et de le faire sur des flotteurs suffisamment précis:
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
qui je pense devrait fonctionner:
>>> 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
semble que vous pourriez vérifier comme ceci:
if int(math.sqrt(n))**2 == n:
print n, 'is a perfect square'
mise à jour:
comme vous l'avez fait remarquer ci-dessus échoue pour les grandes valeurs de n . Pour ceux-ci, ce qui suit semble prometteur, ce qui est une adaptation du code exemple C, par Martin Guy @ UKC, juin 1985, pour la méthode de calcul numérique binaire à l'apparence relativement simple mentionnée dans L'article de Wikipedia méthodes de calcul des racines carrées :
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 '')
sortie:
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)
votre fonction échoue pour les entrées importantes:
In [26]: isqrt((10**100+1)**2)
ValueError: input was not a perfect square
il y a une recette sur le site actif qui devrait être, espérons-le, plus fiable puisqu'elle utilise des maths entières seulement. Il est basé sur une question antérieure de StackOverflow: écrire votre propre fonction de racine carrée
j'ai trouvé ce fil il y a quelques jours et j'ai réécrit la solution de nibot , et en coupant le nombre d'itérations en deux et en faisant quelques autres petits ajustements de performance, j'ai été capable d'améliorer la performance d'un facteur de ~ 2,4: "15198090920"
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)
Voici les résultats de timeit :
>>> timeit('isqrt(12345678901234567890)', setup='from __main__ import isqrt')
8.862877120962366
ensuite, pour comparaison, j'ai mis en œuvre l'algorithme de racine carrée entier le plus couramment utilisé: Newton la méthode . Cette définition est beaucoup plus compacte.
def isqrt(n):
x, y = n, n >> 1
while x > y:
x, y = y, (y + n//y) >> 1
return (x, n - x*x)
il s'avère que même une version optimisée de racines carrées longhand est plus lente que la méthode de Newton, prenant environ 1,5 fois plus de temps.
>>> timeit('isqrt(12345678901234567890)', setup='from __main__ import isqrt')
5.74083631898975
donc, en conclusion, Si vous avez besoin d'une fonction de racine carrée Python pure rapide, ne cherchez pas plus loin que celui fourni ci-dessus.
Edit: j'ai corrigé le bug dans la méthode de Newton ci-dessus. Sur mon machine, il fonctionne ~10% plus rapide que la solution user448810 .
ne peuvent pas être représentés avec précision sur les ordinateurs. Vous pouvez tester pour un réglage de proximité epsilon à une petite valeur dans la précision des flotteurs de 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')
j'ai comparé les différentes méthodes données ici avec une boucle:
for i in range (1000000): # 700 msec
r=int(123456781234567**0.5+0.5)
if r**2==123456781234567:rr=r
else:rr=-1
constatant que celui-ci est le plus rapide et n'a pas besoin d'importation mathématique. Très long pourrait échouer, mais regardez ce
15241576832799734552675677489**0.5 = 123456781234567.0
essayer cette condition (pas de calcul supplémentaire):
def isqrt(n):
i = math.sqrt(n)
if i != int(i):
raise ValueError('input was not a perfect square')
return i
si vous en avez besoin pour retourner un int (pas un float avec un zéro de queue), alors Assignez une 2e variable ou calculez int(i) deux fois.