Radice quadrata intera in python

Question

Esiste una radice quadrata intera da qualche parte in python o nelle librerie standard? Voglio che sia esatto (cioè restituisca un intero) e abbaia se non c'è soluzione.

Al momento ho arrotolato la mia ingenua uno:

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

Ma è brutto e non mi fido per i grandi numeri interi. Potrei scorrere i quadrati e rinunciare se ho superato il valore, ma presumo che sarebbe un po 'lento fare qualcosa del genere. Anche Credo che probabilmente sarei reinventare la ruota, qualcosa come questo deve sicuramente esistere in pitone già ...

Answer 1

il metodo di Newton funziona perfettamente su interi:

def isqrt(n): 
    x = n 
    y = (x + 1) // 2 
    while y < x: 
     x = y 
     y = (x + n // x) // 2 
    return x 

Questo restituisce il più grande intero x per i quali x * x non supera n.Se vuoi verificare se il risultato è esattamente la radice quadrata, esegui semplicemente la moltiplicazione per verificare se n è un quadrato perfetto.

Discuto di questo algoritmo e di altri tre algoritmi per il calcolo delle radici quadrate, a my blog.

Answer 2

Prova questa condizione (nessun calcolo aggiuntivo):

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

Se avete bisogno di tornare un int (non uno float con zero finale) quindi assegnare una seconda variabile o calcolare int(i) due volte.

Answer 3

Sembra che si potrebbe verificare in questo modo:

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

Aggiornamento:

Come lei ha sottolineato quanto sopra non riesce per grandi valori di n. Per quelli che seguono sembra promettente, che è un adattamento del codice C di esempio, di Martin Guy @ UKC, giugno 1985, per il metodo di calcolo numerico cifra per cifra relativamente semplice dell'aspetto menzionato nell'articolo di Wikipedia Methods of computing square roots:

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

uscita:

    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 4

la funzione non riesce per i grandi ingressi:

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

ValueError: input was not a perfect square 

C'è una recipe on the ActiveState site che sho si spera che sia più affidabile dato che usa solo matematica intera. Si basa su una precedente domanda StackOverflow: Writing your own square root function

Answer 5

Una possibilità sarebbe quella di utilizzare il modulo decimal, e farlo in carri sufficientemente preciso:

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 

che credo dovrebbe funzionare:

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

I galleggianti non possono essere rappresentati con precisione sui computer. Puoi testare l'impostazione di prossimità desiderata epsilon su un valore piccolo all'interno della precisione dei float di 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 7

Ci scusiamo per la risposta molto tarda; Sono appena incappato in questa pagina. Nel caso in cui qualcuno visiti questa pagina in futuro, il modulo python gmpy2 è progettato per funzionare con input molto grandi e include tra le altre cose una funzione radice quadrata intera.

Esempio:

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

Certo, tutto avrà il tag "MPZ", ma mpz di sono compatibili con int di:

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

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

Vedi my other answer per una discussione di esecuzione da parte di questo metodo rispetto ad un certo altre risposte a questa domanda.

download: https://code.google.com/p/gmpy/ algoritmo radice quadrata

Answer 8

lunga mano

Si scopre che esiste un algoritmo per calcolare le radici quadrate che si può calcolare a mano, qualcosa di simile a lungo divisione. Ogni iterazione dell'algoritmo produce esattamente una cifra della radice quadrata risultante mentre consuma due cifre del numero di cui si cerca la radice quadrata. Mentre la versione "long hand" dell'algoritmo è specificata in decimale, funziona in qualsiasi base, con il binario che è il più semplice da implementare e forse il più veloce da eseguire (a seconda della sottostante rappresentazione bignum).

Poiché questo algoritmo opera su numeri cifra per cifra, produce risultati esatti per quadrati perfettamente grandi e arbitrariamente grandi, e per quadratini non perfetti, può produrre tante cifre di precisione (a destra della cifra decimale) come desiderato.

Ci sono due belle writeups sul sito "Dr. Math" che spiegano l'algoritmo:

• Square Roots in Binary
• Longhand Square Roots

Ed ecco un'implementazione in 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) 

Si potrebbe facilmente modificare questa funzione t o conduci ulteriori iterazioni per calcolare la parte frazionaria della radice quadrata. Mi interessava soprattutto il calcolo delle radici di grandi quadrati perfetti.

Non sono sicuro di come questo si colleghi all'algoritmo "intero di metodo di Newton". Sospetto che il metodo di Newton sia più veloce, dal momento che in teoria può generare più bit della soluzione in un'unica iterazione, mentre l'algoritmo della "mano lunga" genera esattamente un bit della soluzione per iterazione.

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

Answer 9

Ecco un'implementazione molto semplice:

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 

Questa è solo una ricerca binaria. Inizializza il valore m in modo che sia la più grande potenza di 2 che non superi la radice quadrata, quindi controlla se è possibile impostare ciascun bit più piccolo mantenendo il risultato non più grande della radice quadrata. (. Controllare i bit uno alla volta, in ordine discendente)

Per ragionevolmente grandi valori di n (per esempio, circa 10**6000, o nelle vicinanze 20000 bit), questo sembra essere:

• veloce rispetto alla Implementazione del metodo di Newton described by user448810.
• Molto, molto più lento del metodo integrato gmpy2 in my other answer.
• Paragonabile a, ma leggermente più lento, della radice quadrata lunga described by nibot.

Tutti questi approcci riescono su input di queste dimensioni, ma sulla mia macchina, questa funzione richiede circa 1,5 secondi, mentre @ Nibot di dura circa 0,9 secondi, @ user448810 di dura circa 19 secondi, e il gmpy2 built-in il metodo richiede meno di un millisecondo (!). Esempio:

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

Questa funzione può essere generalizzato facilmente, anche se non è così bello perché non ho abbastanza precisa di un'ipotesi iniziale per 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 

Tuttavia, si noti che gmpy2 ha anche un metodo i_root.

In effetti questo metodo potrebbe essere adattato e applicato a qualsiasi funzione (non negativa, crescente) f per determinare un "intero inverso di f". Tuttavia, per scegliere un valore iniziale efficiente di m, vorresti comunque sapere qualcosa su f.

Modifica: Grazie a @Greggo per aver sottolineato che la funzione i_sqrt può essere riscritta per evitare l'uso di moltiplicazioni. Questo dà un notevole incremento delle prestazioni!

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 

noti che per costruzione, il k esimo bit m << 1 non è impostato, in modo bit-o può essere utilizzato per implementare l'aggiunta di (m<<1) + (1<<k). In definitiva ho (2*m*(2**k) + 2**(2*k)) scritto come (((m<<1) | (1<<k)) << k), quindi è tre turni e un bit per bit o (seguito da una sottrazione per ottenere new_diff). Forse c'è ancora un modo più efficiente per ottenere questo? Indipendentemente da ciò, è molto meglio che moltiplicare m*m! Confronta con sopra:

>>> 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 10

ho trovato questa discussione un paio di giorni fa e ri-scritto nibot's solution, e tagliando il numero di iterazioni a metà e facendo alcune altre modifiche di prestazioni minori, sono stato in grado di migliorare le prestazioni da un fattore di ~ 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) 

Ecco i risultati da timeit:

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

Poi, per il confronto, ho implementato l'algoritmo radice più comunemente usato intero quadrato: Newton's method. Questa definizione è molto più compatta.

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

Risulta che anche una versione ottimizzata di longhand radici quadrate è più lento rispetto al metodo di Newton, impiegando circa 1,5 volte più a lungo.

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

Quindi, in conclusione, se avete bisogno di una funzione veloce puro Python intero radice quadrata, non guardare oltre quello fornito in precedenza.

Modifica: ho corretto il bug nel metodo di Newton sopra. Sulla mia macchina, funziona ~ 10% più veloce di user448810's solution.

Answer 11

ho confrontato i diversi metodi qui riportati con un ciclo:

for i in range (1000000): # 700 msec 
    r=int(123456781234567**0.5+0.5) 
    if r**2==123456781234567:rr=r 
    else:rr=-1 

scoprendo che questo è più veloce e non hanno bisogno di matematica-import. Molto lungo potrebbe non riuscire, ma guarda questo

15241576832799734552675677489**0.5 = 123456781234567.0