java - Il modo più veloce per determinare se la radice quadrata di un intero è un numero intero



math optimization perfect-square (25)

Per quanto riguarda il metodo Carmac, sembra che sarebbe abbastanza semplice ripetere ancora una volta, il che dovrebbe raddoppiare il numero di cifre di precisione. È, dopo tutto, un metodo iterativo estremamente troncato - Newton, con un'ottima prima ipotesi.

Per quanto riguarda il tuo attuale migliore, vedo due micro-ottimizzazioni:

  • sposta il segno di spunta contro 0 dopo il controllo usando mod255
  • riorganizzare i poteri di divisione di quattro per saltare tutti i controlli per il solito caso (75%).

Vale a dire:

// Divide out powers of 4 using binary search

if((n & 0x3L) == 0) {
  n >>=2;

  if((n & 0xffffffffL) == 0)
    n >>= 32;
  if((n & 0xffffL) == 0)
      n >>= 16;
  if((n & 0xffL) == 0)
      n >>= 8;
  if((n & 0xfL) == 0)
      n >>= 4;
  if((n & 0x3L) == 0)
      n >>= 2;
}

Ancora meglio potrebbe essere un semplice

while ((n & 0x03L) == 0) n >>= 2;

Ovviamente, sarebbe interessante sapere quanti numeri vengono abbattuti ad ogni checkpoint - dubito piuttosto che i controlli siano veramente indipendenti, il che rende le cose complicate.

Sto cercando il modo più veloce per determinare se un valore long è un quadrato perfetto (cioè la sua radice quadrata è un altro intero):

  1. L'ho fatto in modo semplice, usando la funzione Math.sqrt () integrata, ma mi chiedo se c'è un modo per farlo più velocemente limitando te stesso al dominio di solo intero.
  2. Il mantenimento di una tabella di ricerca è impraticabile (poiché vi sono circa 2 31.5 numeri interi il cui quadrato è inferiore a 2 63 ).

Ecco il modo molto semplice e diretto che sto facendo ora:

public final static boolean isPerfectSquare(long n)
{
  if (n < 0)
    return false;

  long tst = (long)(Math.sqrt(n) + 0.5);
  return tst*tst == n;
}

Note: sto usando questa funzione in molti problemi di Project Euler . Quindi nessun altro dovrà mai mantenere questo codice. E questo tipo di micro-ottimizzazione potrebbe davvero fare la differenza, dato che parte della sfida consiste nel fare ogni algoritmo in meno di un minuto, e questa funzione dovrà essere chiamata milioni di volte in alcuni problemi.

Una nuova soluzione pubblicata da A. Rex ha dimostrato di essere ancora più veloce. In una corsa sui primi 1 miliardo di interi, la soluzione richiedeva solo il 34% del tempo di utilizzo della soluzione originale. Mentre l'hack di John Carmack è un po 'migliore per piccoli valori di n , il vantaggio rispetto a questa soluzione è piuttosto ridotto.

Ecco la soluzione di A. Rex, convertita in Java:

private final static boolean isPerfectSquare(long n)
{
  // Quickfail
  if( n < 0 || ((n&2) != 0) || ((n & 7) == 5) || ((n & 11) == 8) )
    return false;
  if( n == 0 )
    return true;

  // Check mod 255 = 3 * 5 * 17, for fun
  long y = n;
  y = (y & 0xffffffffL) + (y >> 32);
  y = (y & 0xffffL) + (y >> 16);
  y = (y & 0xffL) + ((y >> 8) & 0xffL) + (y >> 16);
  if( bad255[(int)y] )
      return false;

  // Divide out powers of 4 using binary search
  if((n & 0xffffffffL) == 0)
      n >>= 32;
  if((n & 0xffffL) == 0)
      n >>= 16;
  if((n & 0xffL) == 0)
      n >>= 8;
  if((n & 0xfL) == 0)
      n >>= 4;
  if((n & 0x3L) == 0)
      n >>= 2;

  if((n & 0x7L) != 1)
      return false;

  // Compute sqrt using something like Hensel's lemma
  long r, t, z;
  r = start[(int)((n >> 3) & 0x3ffL)];
  do {
    z = n - r * r;
    if( z == 0 )
      return true;
    if( z < 0 )
      return false;
    t = z & (-z);
    r += (z & t) >> 1;
    if( r > (t  >> 1) )
    r = t - r;
  } while( t <= (1L << 33) );
  return false;
}

private static boolean[] bad255 =
{
   false,false,true ,true ,false,true ,true ,true ,true ,false,true ,true ,true ,
   true ,true ,false,false,true ,true ,false,true ,false,true ,true ,true ,false,
   true ,true ,true ,true ,false,true ,true ,true ,false,true ,false,true ,true ,
   true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,false,true ,false,
   true ,true ,true ,false,true ,true ,true ,true ,false,true ,true ,true ,false,
   true ,false,true ,true ,false,false,true ,true ,true ,true ,true ,false,true ,
   true ,true ,true ,false,true ,true ,false,false,true ,true ,true ,true ,true ,
   true ,true ,true ,false,true ,true ,true ,true ,true ,false,true ,true ,true ,
   true ,true ,false,true ,true ,true ,true ,false,true ,true ,true ,false,true ,
   true ,true ,true ,false,false,true ,true ,true ,true ,true ,true ,true ,true ,
   true ,true ,true ,true ,true ,false,false,true ,true ,true ,true ,true ,true ,
   true ,false,false,true ,true ,true ,true ,true ,false,true ,true ,false,true ,
   true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,false,true ,true ,
   false,true ,false,true ,true ,false,true ,true ,true ,true ,true ,true ,true ,
   true ,true ,true ,true ,false,true ,true ,false,true ,true ,true ,true ,true ,
   false,false,true ,true ,true ,true ,true ,true ,true ,false,false,true ,true ,
   true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,false,false,
   true ,true ,true ,true ,false,true ,true ,true ,false,true ,true ,true ,true ,
   false,true ,true ,true ,true ,true ,false,true ,true ,true ,true ,true ,false,
   true ,true ,true ,true ,true ,true ,true ,true ,false,false,true ,true ,false,
   true ,true ,true ,true ,false,true ,true ,true ,true ,true ,false,false,true ,
   true ,false,true ,false,true ,true ,true ,false,true ,true ,true ,true ,false,
   true ,true ,true ,false,true ,false,true ,true ,true ,true ,true ,true ,true ,
   true ,true ,true ,true ,true ,false,true ,false,true ,true ,true ,false,true ,
   true ,true ,true ,false,true ,true ,true ,false,true ,false,true ,true ,false,
   false,true ,true ,true ,true ,true ,false,true ,true ,true ,true ,false,true ,
   true ,false,false,true ,true ,true ,true ,true ,true ,true ,true ,false,true ,
   true ,true ,true ,true ,false,true ,true ,true ,true ,true ,false,true ,true ,
   true ,true ,false,true ,true ,true ,false,true ,true ,true ,true ,false,false,
   true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,
   false,false,true ,true ,true ,true ,true ,true ,true ,false,false,true ,true ,
   true ,true ,true ,false,true ,true ,false,true ,true ,true ,true ,true ,true ,
   true ,true ,true ,true ,true ,false,true ,true ,false,true ,false,true ,true ,
   false,true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,true ,false,
   true ,true ,false,true ,true ,true ,true ,true ,false,false,true ,true ,true ,
   true ,true ,true ,true ,false,false,true ,true ,true ,true ,true ,true ,true ,
   true ,true ,true ,true ,true ,true ,false,false,true ,true ,true ,true ,false,
   true ,true ,true ,false,true ,true ,true ,true ,false,true ,true ,true ,true ,
   true ,false,true ,true ,true ,true ,true ,false,true ,true ,true ,true ,true ,
   true ,true ,true ,false,false
};

private static int[] start =
{
  1,3,1769,5,1937,1741,7,1451,479,157,9,91,945,659,1817,11,
  1983,707,1321,1211,1071,13,1479,405,415,1501,1609,741,15,339,1703,203,
  129,1411,873,1669,17,1715,1145,1835,351,1251,887,1573,975,19,1127,395,
  1855,1981,425,453,1105,653,327,21,287,93,713,1691,1935,301,551,587,
  257,1277,23,763,1903,1075,1799,1877,223,1437,1783,859,1201,621,25,779,
  1727,573,471,1979,815,1293,825,363,159,1315,183,27,241,941,601,971,
  385,131,919,901,273,435,647,1493,95,29,1417,805,719,1261,1177,1163,
  1599,835,1367,315,1361,1933,1977,747,31,1373,1079,1637,1679,1581,1753,1355,
  513,1539,1815,1531,1647,205,505,1109,33,1379,521,1627,1457,1901,1767,1547,
  1471,1853,1833,1349,559,1523,967,1131,97,35,1975,795,497,1875,1191,1739,
  641,1149,1385,133,529,845,1657,725,161,1309,375,37,463,1555,615,1931,
  1343,445,937,1083,1617,883,185,1515,225,1443,1225,869,1423,1235,39,1973,
  769,259,489,1797,1391,1485,1287,341,289,99,1271,1701,1713,915,537,1781,
  1215,963,41,581,303,243,1337,1899,353,1245,329,1563,753,595,1113,1589,
  897,1667,407,635,785,1971,135,43,417,1507,1929,731,207,275,1689,1397,
  1087,1725,855,1851,1873,397,1607,1813,481,163,567,101,1167,45,1831,1205,
  1025,1021,1303,1029,1135,1331,1017,427,545,1181,1033,933,1969,365,1255,1013,
  959,317,1751,187,47,1037,455,1429,609,1571,1463,1765,1009,685,679,821,
  1153,387,1897,1403,1041,691,1927,811,673,227,137,1499,49,1005,103,629,
  831,1091,1449,1477,1967,1677,697,1045,737,1117,1737,667,911,1325,473,437,
  1281,1795,1001,261,879,51,775,1195,801,1635,759,165,1871,1645,1049,245,
  703,1597,553,955,209,1779,1849,661,865,291,841,997,1265,1965,1625,53,
  1409,893,105,1925,1297,589,377,1579,929,1053,1655,1829,305,1811,1895,139,
  575,189,343,709,1711,1139,1095,277,993,1699,55,1435,655,1491,1319,331,
  1537,515,791,507,623,1229,1529,1963,1057,355,1545,603,1615,1171,743,523,
  447,1219,1239,1723,465,499,57,107,1121,989,951,229,1521,851,167,715,
  1665,1923,1687,1157,1553,1869,1415,1749,1185,1763,649,1061,561,531,409,907,
  319,1469,1961,59,1455,141,1209,491,1249,419,1847,1893,399,211,985,1099,
  1793,765,1513,1275,367,1587,263,1365,1313,925,247,1371,1359,109,1561,1291,
  191,61,1065,1605,721,781,1735,875,1377,1827,1353,539,1777,429,1959,1483,
  1921,643,617,389,1809,947,889,981,1441,483,1143,293,817,749,1383,1675,
  63,1347,169,827,1199,1421,583,1259,1505,861,457,1125,143,1069,807,1867,
  2047,2045,279,2043,111,307,2041,597,1569,1891,2039,1957,1103,1389,231,2037,
  65,1341,727,837,977,2035,569,1643,1633,547,439,1307,2033,1709,345,1845,
  1919,637,1175,379,2031,333,903,213,1697,797,1161,475,1073,2029,921,1653,
  193,67,1623,1595,943,1395,1721,2027,1761,1955,1335,357,113,1747,1497,1461,
  1791,771,2025,1285,145,973,249,171,1825,611,265,1189,847,1427,2023,1269,
  321,1475,1577,69,1233,755,1223,1685,1889,733,1865,2021,1807,1107,1447,1077,
  1663,1917,1129,1147,1775,1613,1401,555,1953,2019,631,1243,1329,787,871,885,
  449,1213,681,1733,687,115,71,1301,2017,675,969,411,369,467,295,693,
  1535,509,233,517,401,1843,1543,939,2015,669,1527,421,591,147,281,501,
  577,195,215,699,1489,525,1081,917,1951,2013,73,1253,1551,173,857,309,
  1407,899,663,1915,1519,1203,391,1323,1887,739,1673,2011,1585,493,1433,117,
  705,1603,1111,965,431,1165,1863,533,1823,605,823,1179,625,813,2009,75,
  1279,1789,1559,251,657,563,761,1707,1759,1949,777,347,335,1133,1511,267,
  833,1085,2007,1467,1745,1805,711,149,1695,803,1719,485,1295,1453,935,459,
  1151,381,1641,1413,1263,77,1913,2005,1631,541,119,1317,1841,1773,359,651,
  961,323,1193,197,175,1651,441,235,1567,1885,1481,1947,881,2003,217,843,
  1023,1027,745,1019,913,717,1031,1621,1503,867,1015,1115,79,1683,793,1035,
  1089,1731,297,1861,2001,1011,1593,619,1439,477,585,283,1039,1363,1369,1227,
  895,1661,151,645,1007,1357,121,1237,1375,1821,1911,549,1999,1043,1945,1419,
  1217,957,599,571,81,371,1351,1003,1311,931,311,1381,1137,723,1575,1611,
  767,253,1047,1787,1169,1997,1273,853,1247,413,1289,1883,177,403,999,1803,
  1345,451,1495,1093,1839,269,199,1387,1183,1757,1207,1051,783,83,423,1995,
  639,1155,1943,123,751,1459,1671,469,1119,995,393,219,1743,237,153,1909,
  1473,1859,1705,1339,337,909,953,1771,1055,349,1993,613,1393,557,729,1717,
  511,1533,1257,1541,1425,819,519,85,991,1693,503,1445,433,877,1305,1525,
  1601,829,809,325,1583,1549,1991,1941,927,1059,1097,1819,527,1197,1881,1333,
  383,125,361,891,495,179,633,299,863,285,1399,987,1487,1517,1639,1141,
  1729,579,87,1989,593,1907,839,1557,799,1629,201,155,1649,1837,1063,949,
  255,1283,535,773,1681,461,1785,683,735,1123,1801,677,689,1939,487,757,
  1857,1987,983,443,1327,1267,313,1173,671,221,695,1509,271,1619,89,565,
  127,1405,1431,1659,239,1101,1159,1067,607,1565,905,1755,1231,1299,665,373,
  1985,701,1879,1221,849,627,1465,789,543,1187,1591,923,1905,979,1241,181
};

Ho provato le diverse soluzioni presentate di seguito.

  • Dopo approfonditi test, ho trovato che aggiungere 0.5 al risultato di Math.sqrt () non è necessario, almeno non sulla mia macchina.
  • L' hack di John Carmack è stato più veloce, ma ha dato risultati errati a partire da n = 410881. Tuttavia, come suggerito da BobbyShaftoe , possiamo usare l'hack Carmack per n <410881.
  • Il metodo di Newton era leggermente più lento di Math.sqrt() . Ciò è probabilmente dovuto al fatto che Math.sqrt() utilizza qualcosa di simile al metodo di Newton, ma implementato nell'hardware in modo che sia molto più veloce rispetto a Java. Inoltre, il metodo di Newton richiedeva ancora l'uso del doppio.
  • Un metodo di Newton modificato, che utilizzava alcuni trucchi per coinvolgere solo la matematica intera, richiedeva alcuni hack per evitare l'overflow (voglio che questa funzione funzioni con tutti gli interi con segno positivo a 64 bit) ed era ancora più lento di Math.sqrt() .
  • Il chop binario era ancora più lento. Questo ha senso perché il chop binario richiede in media 16 passaggi per trovare la radice quadrata di un numero a 64 bit.

L'unico suggerimento che ha mostrato miglioramenti è stato fatto da John D. Cook . È possibile osservare che l'ultima cifra esadecimale (ovvero gli ultimi 4 bit) di un quadrato perfetto deve essere 0, 1, 4 o 9. Ciò significa che il 75% dei numeri può essere immediatamente eliminato come possibile quadrato. L'implementazione di questa soluzione ha comportato una riduzione del 50% in runtime.

Sulla base del suggerimento di John, ho studiato le proprietà degli ultimi n bit di un quadrato perfetto. Analizzando gli ultimi 6 bit, ho trovato che solo 12 valori su 64 sono possibili per gli ultimi 6 bit. Ciò significa che l'81% dei valori può essere eliminato senza l'uso della matematica. L'implementazione di questa soluzione ha comportato una riduzione dell'8% in runtime (rispetto all'algoritmo originale). Analizzando più di 6 bit si ottiene un elenco di possibili bit di fine che è troppo grande per essere pratico.

Ecco il codice che ho usato, che viene eseguito nel 42% del tempo richiesto dall'algoritmo originale (basato su una corsa sui primi 100 milioni di interi). Per valori di n inferiori a 410881, viene eseguito solo nel 29% del tempo richiesto dall'algoritmo originale.

private final static boolean isPerfectSquare(long n)
{
  if (n < 0)
    return false;

  switch((int)(n & 0x3F))
  {
  case 0x00: case 0x01: case 0x04: case 0x09: case 0x10: case 0x11:
  case 0x19: case 0x21: case 0x24: case 0x29: case 0x31: case 0x39:
    long sqrt;
    if(n < 410881L)
    {
      //John Carmack hack, converted to Java.
      // See: http://www.codemaestro.com/reviews/9
      int i;
      float x2, y;

      x2 = n * 0.5F;
      y  = n;
      i  = Float.floatToRawIntBits(y);
      i  = 0x5f3759df - ( i >> 1 );
      y  = Float.intBitsToFloat(i);
      y  = y * ( 1.5F - ( x2 * y * y ) );

      sqrt = (long)(1.0F/y);
    }
    else
    {
      //Carmack hack gives incorrect answer for n >= 410881.
      sqrt = (long)Math.sqrt(n);
    }
    return sqrt*sqrt == n;

  default:
    return false;
  }
}

Note :

  • Secondo i test di John, l'uso or istruzioni sono più veloci in C ++ rispetto all'utilizzo di uno switch , ma in Java e C # non sembra esserci differenza tra or e switch .
  • Ho anche provato a creare una tabella di ricerca (come array statico privato di 64 valori booleani). Quindi al posto di uno switch o di una istruzione, vorrei solo dire if(lookup[(int)(n&0x3F)]) { test } else return false; . Con mia sorpresa, questo è stato (solo leggermente) più lento. Non sono sicuro del perché. Questo perché i limiti dell'array sono controllati in Java .

"Sto cercando il modo più veloce per determinare se un valore lungo è un quadrato perfetto (cioè la sua radice quadrata è un altro intero)."

Le risposte sono impressionanti, ma non sono riuscito a vedere un semplice controllo:

controlla se il primo numero a destra della lunga è un membro del set (0,1,4,5,6,9). Se non lo è, allora non può essere un "quadrato perfetto".

per esempio.

4567 - non può essere un quadrato perfetto.


Metodo di Newton con aritmetica intera

Se si desidera evitare operazioni non interi, è possibile utilizzare il metodo seguente. In pratica utilizza il metodo di Newton modificato per l'aritmetica dei numeri interi.

/**
 * Test if the given number is a perfect square.
 * @param n Must be greater than 0 and less
 *    than Long.MAX_VALUE.
 * @return <code>true</code> if n is a perfect
 *    square, or <code>false</code> otherwise.
 */
public static boolean isSquare(long n)
{
    long x1 = n;
    long x2 = 1L;

    while (x1 > x2)
    {
        x1 = (x1 + x2) / 2L;
        x2 = n / x1;
    }

    return x1 == x2 && n % x1 == 0L;
}

Questa implementazione non può competere con le soluzioni che usano Math.sqrt. Tuttavia, le sue prestazioni possono essere migliorate utilizzando i meccanismi di filtraggio descritti in alcuni degli altri post.


Project Euler è menzionato nei tag e molti dei problemi in esso richiedono il controllo dei numeri >> 2 ^ 64. La maggior parte delle ottimizzazioni di cui sopra non funzionano facilmente quando si lavora con un buffer da 80 byte.

Ho usato java BigInteger e una versione leggermente modificata del metodo di Newton, uno che funziona meglio con gli interi. Il problema era che i quadrati esatti n ^ 2 convergevano in (n-1) invece di n perché n ^ 2-1 = (n-1) (n + 1) e l'errore finale era solo un gradino sotto il divisore finale e il algoritmo terminato. È stato facile risolvere aggiungendo uno all'argomento originale prima di calcolare l'errore. (Aggiungi due per le radici dei cubi, ecc.)

Un bell'attributo di questo algoritmo è che puoi immediatamente capire se il numero è un quadrato perfetto - l'errore finale (non la correzione) nel metodo di Newton sarà zero. Una semplice modifica consente inoltre di calcolare rapidamente floor (sqrt (x)) anziché il numero intero più vicino. Questo è utile con diversi problemi di Eulero.


Stavo pensando ai tempi orribili che ho trascorso nel corso di Analisi Numerica.

E poi ricordo, c'era questa funzione che girava attorno alla rete dal codice sorgente di Quake:

float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;  // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 ); // wtf?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
  #endif
  #endif
  return y;
}

Che calcola fondamentalmente una radice quadrata, usando la funzione di approssimazione di Newton (non ricordo il nome esatto).

Dovrebbe essere utilizzabile e potrebbe anche essere più veloce, proviene da uno dei fenomenali giochi del software di identificazione!

È scritto in C ++ ma non dovrebbe essere troppo difficile riutilizzare la stessa tecnica in Java una volta ottenuta l'idea:

L'ho trovato originariamente su: http://www.codemaestro.com/reviews/9

Il metodo di Newton spiegato a wikipedia: http://en.wikipedia.org/wiki/Newton%27s_method

Puoi seguire il link per ulteriori spiegazioni su come funziona, ma se non ti interessa molto, allora questo è più o meno quello che ricordo leggendo il blog e prendendo il corso di Analisi Numerica:

  • in * (long*) &ypratica si tratta di una funzione di conversione rapida che consente di applicare le operazioni sui numeri interi ai byte non elaborati.
  • la 0x5f3759df - (i >> 1);linea è un valore di seme precalcolato per la funzione di approssimazione.
  • la * (float*) &iconverte il valore a virgola mobile.
  • la y = y * ( threehalfs - ( x2 * y * y ) )linea esegue di nuovo il valore sulla funzione in modo basale.

La funzione di approssimazione fornisce valori più precisi, più si itera la funzione sul risultato. Nel caso di Quake, un'iterazione è "abbastanza buona", ma se non fosse per te ... allora potresti aggiungere tutta l'iterazione di cui hai bisogno.

Questo dovrebbe essere più veloce perché riduce il numero di operazioni di divisione eseguite in un ingenuo quadrato di rooting fino a una semplice divisione per 2 (in realtà un'operazione di * 0.5Fmoltiplicazione) e sostituirlo con un numero limitato di operazioni di moltiplicazione.


Considerando per la lunghezza di bit generale (anche se qui ho usato un tipo specifico), ho cercato di progettare algo semplicistico come di seguito. Inizialmente è richiesto un controllo semplice e ovvio per 0,1,2 o <0. Di seguito è semplice, nel senso che non tenta di utilizzare le funzioni matematiche esistenti. La maggior parte dell'operatore può essere sostituita con operatori bit-saggio. Non ho ancora provato con i dati benchmark. Non sono né esperto di matematica o di progettazione di algoritmi informatici in particolare, mi piacerebbe vederti a sottolineare il problema. So che ci sono molte possibilità di miglioramento lì.

int main()
{
    unsigned int c1=0 ,c2 = 0;  
    unsigned int x = 0;  
    unsigned int p = 0;  
    int k1 = 0;  
    scanf("%d",&p);  
    if(p % 2 == 0) {  
        x = p/2; 
    }  
    else {  
        x = (p/2) +1;  
    }  
    while(x) 
    {
        if((x*x) > p) {  
            c1 = x;  
            x = x/2; 
        }else {  
            c2 = x;  
            break;  
        }  
    }  
    if((p%2) != 0)  
        c2++;

    while(c2 < c1) 
    {  
        if((c2 * c2 ) == p) {  
            k1 = 1;  
            break;  
        }  
        c2++; 
    }  
    if(k1)  
        printf("\n Perfect square for %d", c2);  
    else  
        printf("\n Not perfect but nearest to :%d :", c2);  
    return 0;  
}  

Non so se questo è stato menzionato prima. Ma ho trovato una soluzione here :

int result = (int)(floor(sqrt(b)) - ceil(sqrt(a)) + 1);

Un problema intero merita una soluzione intera. così

Fai una ricerca binaria sugli interi (non negativi) per trovare il più grande intero t tale che t**2 <= n. Quindi verifica se r**2 = nesattamente. Questo richiede tempo O (log n).

Se non sai come eseguire la ricerca binaria negli interi positivi perché il set è illimitato, è facile. Stai iniziando calcolando la tua crescente funzione f (sopra f(t) = t**2 - n) su potenze di due. Quando lo vedi diventa positivo, hai trovato un limite superiore. Quindi puoi eseguire una ricerca binaria standard.


Ho controllato tutti i possibili risultati quando si osservano gli ultimi n bit di un quadrato. Esaminando successivamente più bit, è possibile eliminare fino al 5/6 degli ingressi. In realtà ho progettato questo per implementare l'algoritmo Factorization di Fermat, ed è molto veloce lì.

public static boolean isSquare(final long val) {
   if ((val & 2) == 2 || (val & 7) == 5) {
     return false;
   }
   if ((val & 11) == 8 || (val & 31) == 20) {
     return false;
   }

   if ((val & 47) == 32 || (val & 127) == 80) {
     return false;
   }

   if ((val & 191) == 128 || (val & 511) == 320) {
     return false;
   }

   // if((val & a == b) || (val & c == d){
   //   return false;
   // }

   if (!modSq[(int) (val % modSq.length)]) {
        return false;
   }

   final long root = (long) Math.sqrt(val);
   return root * root == val;
}

L'ultimo bit di pseudocodice può essere utilizzato per estendere i test per eliminare più valori. I test sopra riportati sono per k = 0, 1, 2, 3

  • a è della forma (3 << 2k) - 1
  • b è del modulo (2 << 2k)
  • c è nella forma (2 << 2k + 2) - 1
  • d è del formato (2 << 2k - 1) * 10

    In primo luogo verifica se ha un residuo quadrato con moduli di potenza di due, quindi verifica in base a un modulo finale, quindi utilizza Math.sqrt per eseguire un test finale. Mi è venuta l'idea dal primo posto e ho tentato di estenderla. Apprezzo qualsiasi commento o suggerimento.

    Aggiornamento: utilizzando il test con un modulo, (modSq) e una base del modulo di 44352, il mio test viene eseguito nel 96% del tempo di quello nell'aggiornamento dell'OP per numeri fino a 1.000.000.000.


  • Dovrai fare un po 'di benchmarking. Il miglior algoritmo dipenderà dalla distribuzione dei tuoi input.

    Il tuo algoritmo potrebbe essere quasi ottimale, ma potresti voler fare un rapido controllo per escludere alcune possibilità prima di chiamare la tua routine radice quadrata. Ad esempio, guarda l'ultima cifra del tuo numero in esadecimale facendo un po '"e". I quadrati perfetti possono terminare solo a 0, 1, 4 o 9 nella base 16, quindi per il 75% dei tuoi ingressi (supponendo che siano distribuiti in modo uniforme) puoi evitare una chiamata alla radice quadrata in cambio di un po 'di twittling molto veloce.

    Kip ha confrontato il seguente codice implementando il trucco esadecimale. Durante il test dei numeri da 1 a 100.000.000, questo codice ha funzionato due volte più velocemente dell'originale.

    public final static boolean isPerfectSquare(long n)
    {
        if (n < 0)
            return false;
    
        switch((int)(n & 0xF))
        {
        case 0: case 1: case 4: case 9:
            long tst = (long)Math.sqrt(n);
            return tst*tst == n;
    
        default:
            return false;
        }
    }
    

    Quando ho testato il codice analogo in C ++, in realtà è più lento dell'originale. Tuttavia, quando ho eliminato l'istruzione switch, il trucco hex ancora una volta rende il codice due volte più veloce.

    int isPerfectSquare(int n)
    {
        int h = n & 0xF;  // h is the last hex "digit"
        if (h > 9)
            return 0;
        // Use lazy evaluation to jump out of the if statement as soon as possible
        if (h != 2 && h != 3 && h != 5 && h != 6 && h != 7 && h != 8)
        {
            int t = (int) floor( sqrt((double) n) + 0.5 );
            return t*t == n;
        }
        return 0;
    }
    

    L'eliminazione dell'istruzione switch ha avuto scarso effetto sul codice C #.


    Questa è una rielaborazione da decimale a binario del vecchio algoritmo della calcolatrice di Marchant (mi dispiace, non ho un riferimento), in Ruby, adattato specificamente per questa domanda:

    def isexactsqrt(v)
        value = v.abs
        residue = value
        root = 0
        onebit = 1
        onebit <<= 8 while (onebit < residue)
        onebit >>= 2 while (onebit > residue)
        while (onebit > 0)
            x = root + onebit
            if (residue >= x) then
                residue -= x
                root = x + onebit
            end
            root >>= 1
            onebit >>= 2
        end
        return (residue == 0)
    end
    

    Ecco un workup di qualcosa di simile (per favore non votarmi per lo stile di codifica / odori o clunky O / O - è l'algoritmo che conta, e C ++ non è la mia lingua madre). In questo caso, stiamo cercando residui == 0:

    #include <iostream>  
    
    using namespace std;  
    typedef unsigned long long int llint;
    
    class ISqrt {           // Integer Square Root
        llint value;        // Integer whose square root is required
        llint root;         // Result: floor(sqrt(value))
        llint residue;      // Result: value-root*root
        llint onebit, x;    // Working bit, working value
    
    public:
    
        ISqrt(llint v = 2) {    // Constructor
            Root(v);            // Take the root 
        };
    
        llint Root(llint r) {   // Resets and calculates new square root
            value = r;          // Store input
            residue = value;    // Initialise for subtracting down
            root = 0;           // Clear root accumulator
    
            onebit = 1;                 // Calculate start value of counter
            onebit <<= (8*sizeof(llint)-2);         // Set up counter bit as greatest odd power of 2 
            while (onebit > residue) {onebit >>= 2; };  // Shift down until just < value
    
            while (onebit > 0) {
                x = root ^ onebit;          // Will check root+1bit (root bit corresponding to onebit is always zero)
                if (residue >= x) {         // Room to subtract?
                    residue -= x;           // Yes - deduct from residue
                    root = x + onebit;      // and step root
                };
                root >>= 1;
                onebit >>= 2;
            };
            return root;                    
        };
        llint Residue() {           // Returns residue from last calculation
            return residue;                 
        };
    };
    
    int main() {
        llint big, i, q, r, v, delta;
        big = 0; big = (big-1);         // Kludge for "big number"
        ISqrt b;                            // Make q sqrt generator
        for ( i = big; i > 0 ; i /= 7 ) {   // for several numbers
            q = b.Root(i);                  // Get the square root
            r = b.Residue();                // Get the residue
            v = q*q+r;                      // Recalc original value
            delta = v-i;                    // And diff, hopefully 0
            cout << i << ": " << q << " ++ " << r << " V: " << v << " Delta: " << delta << "\n";
        };
        return 0;
    };
    

    Sono in ritardo per la festa, ma spero di fornire una risposta migliore; più corto e (supponendo che il mio benchmark sia corretto) anche molto faster .

    long goodMask; // 0xC840C04048404040 computed below
    {
        for (int i=0; i<64; ++i) goodMask |= Long.MIN_VALUE >>> (i*i);
    }
    
    public boolean isSquare(long x) {
        // This tests if the 6 least significant bits are right.
        // Moving the to be tested bit to the highest position saves us masking.
        if (goodMask << x >= 0) return false;
        final int numberOfTrailingZeros = Long.numberOfTrailingZeros(x);
        // Each square ends with an even number of zeros.
        if ((numberOfTrailingZeros & 1) != 0) return false;
        x >>= numberOfTrailingZeros;
        // Now x is either 0 or odd.
        // In binary each odd square ends with 001.
        // Postpone the sign test until now; handle zero in the branch.
        if ((x&7) != 1 | x <= 0) return x == 0;
        // Do it in the classical way.
        // The correctness is not trivial as the conversion from long to double is lossy!
        final long tst = (long) Math.sqrt(x);
        return tst * tst == x;
    }
    

    Il primo test cattura la maggior parte dei non quadrati rapidamente. Utilizza una tabella di 64 elementi compressa in un lungo, quindi non c'è alcun costo di accesso alla matrice (controlli di indirezione e limiti). Per un long uniformemente casuale, c'è una probabilità dell'81.25% di finire qui.

    Il secondo test cattura tutti i numeri con un numero dispari di due nella loro fattorizzazione. Il metodo Long.numberOfTrailingZeros è molto veloce in quanto ottiene JIT-ed in una singola istruzione i86.

    Dopo aver lasciato cadere gli zeri finali, il terzo test gestisce i numeri che terminano con 011, 101 o 111 in binario, che non sono quadrati perfetti. Si preoccupa anche dei numeri negativi e gestisce anche 0.

    Il test finale ritorna alla double aritmetica. Poiché il double ha solo mantissa a 53 bit, la conversione da long a double include l'arrotondamento per valori grandi. Tuttavia, il test è corretto (a meno che la proof non sia corretta).

    Cercare di incorporare l'idea mod255 non ha avuto successo.


    La chiamata sqrt non è perfettamente precisa, come è stato menzionato, ma è interessante e istruttivo che non elimini le altre risposte in termini di velocità. Dopo tutto, la sequenza delle istruzioni di linguaggio assembly per un sqrt è minuscola. Intel ha un'istruzione hardware, che non è utilizzata da Java, credo perché non è conforme all'IEEE.

    Allora perché è lento? Perché Java sta effettivamente chiamando una routine C attraverso JNI, ed è in realtà più lento a farlo che chiamare una subroutine Java, che a sua volta è più lenta di farlo in linea. Questo è molto fastidioso, e Java avrebbe dovuto trovare una soluzione migliore, ad esempio la creazione di chiamate in libreria a virgola mobile, se necessario. Oh bene.

    In C ++, sospetto che tutte le alternative complesse vadano perse in velocità, ma non le ho controllate tutte. Quello che ho fatto e ciò che le persone Java troveranno utile, è un semplice trucco, un'estensione del test dei casi speciali suggerito da A. Rex. Usa un singolo valore lungo come array di bit, che non è controllato. In questo modo, hai una ricerca booleana a 64 bit.

    typedef unsigned long long UVLONG
    UVLONG pp1,pp2;
    
    void init2() {
      for (int i = 0; i < 64; i++) {
        for (int j = 0; j < 64; j++)
          if (isPerfectSquare(i * 64 + j)) {
        pp1 |= (1 << j);
        pp2 |= (1 << i);
        break;
          }
       }
       cout << "pp1=" << pp1 << "," << pp2 << "\n";  
    }
    
    
    inline bool isPerfectSquare5(UVLONG x) {
      return pp1 & (1 << (x & 0x3F)) ? isPerfectSquare(x) : false;
    }
    

    La routine isPerfectSquare5 viene eseguita in circa 1/3 del tempo sulla mia macchina core2 duo. Sospetto che ulteriori ritocchi lungo la stessa linea potrebbero ridurre il tempo in media, ma ogni volta che controlli, stai negoziando altri test per eliminarli di più, quindi non puoi andare più lontano su quella strada.

    Certamente, piuttosto che avere un test separato per il negativo, è possibile controllare i 6 bit alti allo stesso modo.

    Nota che tutto quello che sto facendo è eliminare i possibili quadrati, ma quando ho un potenziale caso devo chiamare l'originale, isPerfectSquare.

    La routine init2 viene chiamata una volta per inizializzare i valori statici di pp1 e pp2. Nota che nella mia implementazione in C ++, sto usando unsigned long long, quindi visto che sei firmato, dovresti usare l'operatore >>>.

    Non c'è alcuna necessità intrinseca che i limiti controllino l'array, ma l'ottimizzatore di Java deve capire questa cosa abbastanza rapidamente, quindi non li biasimo per questo.


    Non sono sicuro se sarebbe più veloce, o addirittura preciso, ma potresti usare http://www.codemaestro.com/reviews/9 , l'algoritmo per risolvere la radice quadrata più velocemente. Probabilmente potresti testarlo facilmente per tutti i possibili numeri interi a 32 bit e convalidare che hai effettivamente ottenuto risultati corretti, in quanto è solo un appunto. Tuttavia, ora che ci penso, anche l'uso del doppio è approssimativo, quindi non sono sicuro di come entrerebbe in gioco.


    Per le prestazioni, molto spesso devi fare alcune compromissioni. Altri hanno espresso vari metodi, tuttavia, hai notato che l'attacco di Carmack era più veloce fino a certi valori di N. Quindi, dovresti controllare la "n" e se è inferiore a quel numero N, usa l'hack di Carmack, altrimenti usa qualche altro metodo descritto nelle risposte qui.


    Mi piace l'idea di utilizzare un metodo quasi corretto su alcuni degli input. Ecco una versione con un "offset" più alto. Il codice sembra funzionare e passa il mio semplice caso di test.

    Basta sostituire il tuo:

    if(n < 410881L){...}
    

    codice con questo:

    if (n < 11043908100L) {
        //John Carmack hack, converted to Java.
        // See: http://www.codemaestro.com/reviews/9
        int i;
        float x2, y;
    
        x2 = n * 0.5F;
        y = n;
        i = Float.floatToRawIntBits(y);
        //using the magic number from 
        //http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
        //since it more accurate
        i = 0x5f375a86 - (i >> 1);
        y = Float.intBitsToFloat(i);
        y = y * (1.5F - (x2 * y * y));
        y = y * (1.5F - (x2 * y * y)); //Newton iteration, more accurate
    
        sqrt = Math.round(1.0F / y);
    } else {
        //Carmack hack gives incorrect answer for n >= 11043908100.
        sqrt = (long) Math.sqrt(n);
    }
    

    Voglio che questa funzione funzioni con tutti gli interi con segno positivo a 64 bit

    Math.sqrt()funziona con i doppi come parametri di input, quindi non otterrai risultati accurati per numeri interi maggiori di 2 ^ 53 .


    Ho scoperto un metodo che funziona ~ 35% più veloce del tuo 6bits + codice Carmack + sqrt, almeno con la mia CPU (x86) e il linguaggio di programmazione (C / C ++). I risultati possono variare, soprattutto perché non so come si svilupperà il fattore Java.

    Il mio approccio è triplice:

    1. In primo luogo, filtrare risposte ovvie. Questo include numeri negativi e guardando gli ultimi 4 bit. (Ho trovato che guardare gli ultimi sei non mi ha aiutato.) int64 x anche sì per 0. (Leggendo il codice qui sotto, nota che il mio input è int64 x .)
      if( x < 0 || (x&2) || ((x & 7) == 5) || ((x & 11) == 8) )
          return false;
      if( x == 0 )
          return true;
    2. Successivamente, controlla se è un modulo quadrato 255 = 3 * 5 * 17. Poiché questo è un prodotto di tre numeri primi distinti, solo circa 1/8 dei residui mod 255 sono quadrati. Tuttavia, secondo la mia esperienza, chiamare l'operatore modulo (%) costa più del beneficio che si ottiene, quindi uso i trucchi bit che coinvolgono 255 = 2 ^ 8-1 per calcolare il residuo. (Nel bene o nel male, non sto usando il trucco di leggere singoli byte di una parola, solo bit per bit e spostamenti.)
      int64 y = x;
      y = (y & 4294967295LL) + (y >> 32); 
      y = (y & 65535) + (y >> 16);
      y = (y & 255) + ((y >> 8) & 255) + (y >> 16);
      // At this point, y is between 0 and 511.  More code can reduce it farther.
      
      Per verificare effettivamente se il residuo è un quadrato, cerco la risposta in una tabella precalcolata.
      if( bad255[y] )
          return false;
      // However, I just use a table of size 512
      
    3. Infine, prova a calcolare la radice quadrata usando un metodo simile al lemma di Hensel . (Non penso sia applicabile direttamente, ma funziona con alcune modifiche.) Prima di farlo, divido tutti i poteri di 2 con una ricerca binaria:
      if((x & 4294967295LL) == 0)
          x >>= 32;
      if((x & 65535) == 0)
          x >>= 16;
      if((x & 255) == 0)
          x >>= 8;
      if((x & 15) == 0)
          x >>= 4;
      if((x & 3) == 0)
          x >>= 2;
      A questo punto, affinché il nostro numero sia quadrato, deve essere 1 mod 8.
      if((x & 7) != 1)
          return false;
      La struttura di base del lemma di Hensel è la seguente. (Nota: codice non testato, se non funziona, prova t = 2 o 8.)
      int64 t = 4, r = 1;
      t <<= 1; r += ((x - r * r) & t) >> 1;
      t <<= 1; r += ((x - r * r) & t) >> 1;
      t <<= 1; r += ((x - r * r) & t) >> 1;
      // Repeat until t is 2^33 or so.  Use a loop if you want.
      L'idea è che ad ogni iterazione, si aggiunge un bit su r, la radice quadrata "corrente" di x; ogni radice quadrata è accurata modulo una potenza sempre più grande di 2, vale a dire t / 2. Alla fine, r e t / 2-r saranno le radici quadrate di x modulo t / 2. (Si noti che se r è una radice quadrata di x, allora anche -r è vero anche i numeri di modulo, ma attenzione, modulo alcuni numeri, le cose possono avere anche più di 2 radici quadrate, in particolare, questo include poteri di 2. ) Perché la nostra radice quadrata effettiva è inferiore a 2 ^ 32, a quel punto possiamo effettivamente controllare se r o t / 2-r sono radici quadrate reali. Nel mio codice attuale, utilizzo il seguente ciclo modificato:
      int64 r, t, z;
      r = start[(x >> 3) & 1023];
      do {
          z = x - r * r;
          if( z == 0 )
              return true;
          if( z < 0 )
              return false;
          t = z & (-z);
          r += (z & t) >> 1;
          if( r > (t >> 1) )
              r = t - r;
      } while( t <= (1LL << 33) );
      L'accelerazione qui viene ottenuta in tre modi: valore iniziale precalcolato (equivalente a ~ 10 iterazioni del ciclo), uscita precedente del ciclo e salto di alcuni valori t. Per l'ultima parte, guardo z = r - x * x , e impostiamo t come la più grande potenza di 2 divisioni z con un po 'di trucco. Questo mi consente di saltare i valori che non avrebbero comunque influenzato il valore di r. Il valore iniziale precalcolato nel mio caso seleziona la radice quadrata "minimo più piccolo" modulo 8192.

    Anche se questo codice non funziona più velocemente per te, spero che ti piacciano alcune delle idee che contiene. Segue un codice completo e testato, comprese le tabelle precalcolate.

    typedef signed long long int int64;
    
    int start[1024] =
    {1,3,1769,5,1937,1741,7,1451,479,157,9,91,945,659,1817,11,
    1983,707,1321,1211,1071,13,1479,405,415,1501,1609,741,15,339,1703,203,
    129,1411,873,1669,17,1715,1145,1835,351,1251,887,1573,975,19,1127,395,
    1855,1981,425,453,1105,653,327,21,287,93,713,1691,1935,301,551,587,
    257,1277,23,763,1903,1075,1799,1877,223,1437,1783,859,1201,621,25,779,
    1727,573,471,1979,815,1293,825,363,159,1315,183,27,241,941,601,971,
    385,131,919,901,273,435,647,1493,95,29,1417,805,719,1261,1177,1163,
    1599,835,1367,315,1361,1933,1977,747,31,1373,1079,1637,1679,1581,1753,1355,
    513,1539,1815,1531,1647,205,505,1109,33,1379,521,1627,1457,1901,1767,1547,
    1471,1853,1833,1349,559,1523,967,1131,97,35,1975,795,497,1875,1191,1739,
    641,1149,1385,133,529,845,1657,725,161,1309,375,37,463,1555,615,1931,
    1343,445,937,1083,1617,883,185,1515,225,1443,1225,869,1423,1235,39,1973,
    769,259,489,1797,1391,1485,1287,341,289,99,1271,1701,1713,915,537,1781,
    1215,963,41,581,303,243,1337,1899,353,1245,329,1563,753,595,1113,1589,
    897,1667,407,635,785,1971,135,43,417,1507,1929,731,207,275,1689,1397,
    1087,1725,855,1851,1873,397,1607,1813,481,163,567,101,1167,45,1831,1205,
    1025,1021,1303,1029,1135,1331,1017,427,545,1181,1033,933,1969,365,1255,1013,
    959,317,1751,187,47,1037,455,1429,609,1571,1463,1765,1009,685,679,821,
    1153,387,1897,1403,1041,691,1927,811,673,227,137,1499,49,1005,103,629,
    831,1091,1449,1477,1967,1677,697,1045,737,1117,1737,667,911,1325,473,437,
    1281,1795,1001,261,879,51,775,1195,801,1635,759,165,1871,1645,1049,245,
    703,1597,553,955,209,1779,1849,661,865,291,841,997,1265,1965,1625,53,
    1409,893,105,1925,1297,589,377,1579,929,1053,1655,1829,305,1811,1895,139,
    575,189,343,709,1711,1139,1095,277,993,1699,55,1435,655,1491,1319,331,
    1537,515,791,507,623,1229,1529,1963,1057,355,1545,603,1615,1171,743,523,
    447,1219,1239,1723,465,499,57,107,1121,989,951,229,1521,851,167,715,
    1665,1923,1687,1157,1553,1869,1415,1749,1185,1763,649,1061,561,531,409,907,
    319,1469,1961,59,1455,141,1209,491,1249,419,1847,1893,399,211,985,1099,
    1793,765,1513,1275,367,1587,263,1365,1313,925,247,1371,1359,109,1561,1291,
    191,61,1065,1605,721,781,1735,875,1377,1827,1353,539,1777,429,1959,1483,
    1921,643,617,389,1809,947,889,981,1441,483,1143,293,817,749,1383,1675,
    63,1347,169,827,1199,1421,583,1259,1505,861,457,1125,143,1069,807,1867,
    2047,2045,279,2043,111,307,2041,597,1569,1891,2039,1957,1103,1389,231,2037,
    65,1341,727,837,977,2035,569,1643,1633,547,439,1307,2033,1709,345,1845,
    1919,637,1175,379,2031,333,903,213,1697,797,1161,475,1073,2029,921,1653,
    193,67,1623,1595,943,1395,1721,2027,1761,1955,1335,357,113,1747,1497,1461,
    1791,771,2025,1285,145,973,249,171,1825,611,265,1189,847,1427,2023,1269,
    321,1475,1577,69,1233,755,1223,1685,1889,733,1865,2021,1807,1107,1447,1077,
    1663,1917,1129,1147,1775,1613,1401,555,1953,2019,631,1243,1329,787,871,885,
    449,1213,681,1733,687,115,71,1301,2017,675,969,411,369,467,295,693,
    1535,509,233,517,401,1843,1543,939,2015,669,1527,421,591,147,281,501,
    577,195,215,699,1489,525,1081,917,1951,2013,73,1253,1551,173,857,309,
    1407,899,663,1915,1519,1203,391,1323,1887,739,1673,2011,1585,493,1433,117,
    705,1603,1111,965,431,1165,1863,533,1823,605,823,1179,625,813,2009,75,
    1279,1789,1559,251,657,563,761,1707,1759,1949,777,347,335,1133,1511,267,
    833,1085,2007,1467,1745,1805,711,149,1695,803,1719,485,1295,1453,935,459,
    1151,381,1641,1413,1263,77,1913,2005,1631,541,119,1317,1841,1773,359,651,
    961,323,1193,197,175,1651,441,235,1567,1885,1481,1947,881,2003,217,843,
    1023,1027,745,1019,913,717,1031,1621,1503,867,1015,1115,79,1683,793,1035,
    1089,1731,297,1861,2001,1011,1593,619,1439,477,585,283,1039,1363,1369,1227,
    895,1661,151,645,1007,1357,121,1237,1375,1821,1911,549,1999,1043,1945,1419,
    1217,957,599,571,81,371,1351,1003,1311,931,311,1381,1137,723,1575,1611,
    767,253,1047,1787,1169,1997,1273,853,1247,413,1289,1883,177,403,999,1803,
    1345,451,1495,1093,1839,269,199,1387,1183,1757,1207,1051,783,83,423,1995,
    639,1155,1943,123,751,1459,1671,469,1119,995,393,219,1743,237,153,1909,
    1473,1859,1705,1339,337,909,953,1771,1055,349,1993,613,1393,557,729,1717,
    511,1533,1257,1541,1425,819,519,85,991,1693,503,1445,433,877,1305,1525,
    1601,829,809,325,1583,1549,1991,1941,927,1059,1097,1819,527,1197,1881,1333,
    383,125,361,891,495,179,633,299,863,285,1399,987,1487,1517,1639,1141,
    1729,579,87,1989,593,1907,839,1557,799,1629,201,155,1649,1837,1063,949,
    255,1283,535,773,1681,461,1785,683,735,1123,1801,677,689,1939,487,757,
    1857,1987,983,443,1327,1267,313,1173,671,221,695,1509,271,1619,89,565,
    127,1405,1431,1659,239,1101,1159,1067,607,1565,905,1755,1231,1299,665,373,
    1985,701,1879,1221,849,627,1465,789,543,1187,1591,923,1905,979,1241,181};
    
    bool bad255[512] =
    {0,0,1,1,0,1,1,1,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,1,0,1,1,1,1,0,1,
     1,1,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,1,1,1,0,1,1,1,
     0,1,0,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,1,1,0,0,1,1,1,1,1,1,1,1,0,1,
     1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,0,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,1,
     1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,0,1,1,1,1,1,
     1,1,1,1,1,1,0,1,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,
     1,1,1,0,0,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,
     1,0,1,1,1,0,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,
     0,0,1,1,0,1,1,1,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,1,0,1,1,1,1,0,1,
     1,1,0,1,0,1,1,1,1,1,1,1,1,1,1,1,1,0,1,0,1,1,1,0,1,1,1,1,0,1,1,1,
     0,1,0,1,1,0,0,1,1,1,1,1,0,1,1,1,1,0,1,1,0,0,1,1,1,1,1,1,1,1,0,1,
     1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,0,1,1,1,0,1,1,1,1,0,0,1,1,1,1,1,1,
     1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,0,0,1,1,1,1,1,0,1,1,0,1,1,1,1,1,
     1,1,1,1,1,1,0,1,1,0,1,0,1,1,0,1,1,1,1,1,1,1,1,1,1,1,0,1,1,0,1,1,
     1,1,1,0,0,1,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,1,1,1,1,0,0,1,1,1,
     1,0,1,1,1,0,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,0,1,1,1,1,1,1,1,1,
     0,0};
    
    inline bool square( int64 x ) {
        // Quickfail
        if( x < 0 || (x&2) || ((x & 7) == 5) || ((x & 11) == 8) )
            return false;
        if( x == 0 )
            return true;
    
        // Check mod 255 = 3 * 5 * 17, for fun
        int64 y = x;
        y = (y & 4294967295LL) + (y >> 32);
        y = (y & 65535) + (y >> 16);
        y = (y & 255) + ((y >> 8) & 255) + (y >> 16);
        if( bad255[y] )
            return false;
    
        // Divide out powers of 4 using binary search
        if((x & 4294967295LL) == 0)
            x >>= 32;
        if((x & 65535) == 0)
            x >>= 16;
        if((x & 255) == 0)
            x >>= 8;
        if((x & 15) == 0)
            x >>= 4;
        if((x & 3) == 0)
            x >>= 2;
    
        if((x & 7) != 1)
            return false;
    
        // Compute sqrt using something like Hensel's lemma
        int64 r, t, z;
        r = start[(x >> 3) & 1023];
        do {
            z = x - r * r;
            if( z == 0 )
                return true;
            if( z < 0 )
                return false;
            t = z & (-z);
            r += (z & t) >> 1;
            if( r > (t  >> 1) )
                r = t - r;
        } while( t <= (1LL << 33) );
    
        return false;
    }

    Ho eseguito la mia analisi di alcuni degli algoritmi in questo thread e ho trovato alcuni nuovi risultati. Puoi vedere quei vecchi risultati nella cronologia delle modifiche di questa risposta, ma non sono accurati, come ho fatto un errore, e ho perso tempo ad analizzare diversi algoritmi che non sono vicini. Tuttavia, tirando le lezioni da diverse risposte diverse, ora ho due algoritmi che schiacciano il "vincitore" di questo thread. Ecco la cosa principale che faccio in modo diverso rispetto a tutti gli altri:

    // This is faster because a number is divisible by 2^4 or more only 6% of the time
    // and more than that a vanishingly small percentage.
    while((x & 0x3) == 0) x >>= 2;
    // This is effectively the same as the switch-case statement used in the original
    // answer. 
    if((x & 0x7) != 1) return false;
    

    Tuttavia, questa semplice linea, che per la maggior parte delle volte aggiunge una o due istruzioni molto rapide, semplifica enormemente la switch-casedichiarazione in una sola istruzione. Tuttavia, può aggiungere al runtime se molti dei numeri testati hanno significativi fattori di potenza di due.

    Gli algoritmi seguenti sono i seguenti:

    • Internet : risposta di Kip
    • Durron - La mia risposta modificata utilizzando la risposta one-pass come base
    • DurronTwo - La mia risposta modificata usando la risposta in due passaggi (di @JohnnyHeggheim), con alcune altre leggere modifiche.

    Ecco un runtime di esempio se i numeri sono generati usando Math.abs(java.util.Random.nextLong())

     0% Scenario{vm=java, trial=0, benchmark=Internet} 39673.40 ns; ?=378.78 ns @ 3 trials
    33% Scenario{vm=java, trial=0, benchmark=Durron} 37785.75 ns; ?=478.86 ns @ 10 trials
    67% Scenario{vm=java, trial=0, benchmark=DurronTwo} 35978.10 ns; ?=734.10 ns @ 10 trials
    
    benchmark   us linear runtime
     Internet 39.7 ==============================
       Durron 37.8 ============================
    DurronTwo 36.0 ===========================
    
    vm: java
    trial: 0
    

    Ed ecco un runtime di esempio se viene eseguito solo sul primo milione di long:

     0% Scenario{vm=java, trial=0, benchmark=Internet} 2933380.84 ns; ?=56939.84 ns @ 10 trials
    33% Scenario{vm=java, trial=0, benchmark=Durron} 2243266.81 ns; ?=50537.62 ns @ 10 trials
    67% Scenario{vm=java, trial=0, benchmark=DurronTwo} 3159227.68 ns; ?=10766.22 ns @ 3 trials
    
    benchmark   ms linear runtime
     Internet 2.93 ===========================
       Durron 2.24 =====================
    DurronTwo 3.16 ==============================
    
    vm: java
    trial: 0
    

    Come puoi vedere, DurronTwofa meglio per gli input di grandi dimensioni, perché usa molto spesso il trucco magico, ma viene scombinato rispetto al primo algoritmo e Math.sqrtperché i numeri sono molto più piccoli. Nel frattempo, il più semplice Durronè un grande vincitore perché non deve mai dividere per 4 molte volte nel primo milione di numeri.

    Ecco Durron:

    public final static boolean isPerfectSquareDurron(long n) {
        if(n < 0) return false;
        if(n == 0) return true;
    
        long x = n;
        // This is faster because a number is divisible by 16 only 6% of the time
        // and more than that a vanishingly small percentage.
        while((x & 0x3) == 0) x >>= 2;
        // This is effectively the same as the switch-case statement used in the original
        // answer. 
        if((x & 0x7) == 1) {
    
            long sqrt;
            if(x < 410881L)
            {
                int i;
                float x2, y;
    
                x2 = x * 0.5F;
                y  = x;
                i  = Float.floatToRawIntBits(y);
                i  = 0x5f3759df - ( i >> 1 );
                y  = Float.intBitsToFloat(i);
                y  = y * ( 1.5F - ( x2 * y * y ) );
    
                sqrt = (long)(1.0F/y);
            } else {
                sqrt = (long) Math.sqrt(x);
            }
            return sqrt*sqrt == x;
        }
        return false;
    }
    

    E DurronTwo

    public final static boolean isPerfectSquareDurronTwo(long n) {
        if(n < 0) return false;
        // Needed to prevent infinite loop
        if(n == 0) return true;
    
        long x = n;
        while((x & 0x3) == 0) x >>= 2;
        if((x & 0x7) == 1) {
            long sqrt;
            if (x < 41529141369L) {
                int i;
                float x2, y;
    
                x2 = x * 0.5F;
                y = x;
                i = Float.floatToRawIntBits(y);
                //using the magic number from 
                //http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
                //since it more accurate
                i = 0x5f375a86 - (i >> 1);
                y = Float.intBitsToFloat(i);
                y = y * (1.5F - (x2 * y * y));
                y = y * (1.5F - (x2 * y * y)); //Newton iteration, more accurate
                sqrt = (long) ((1.0F/y) + 0.2);
            } else {
                //Carmack hack gives incorrect answer for n >= 41529141369.
                sqrt = (long) Math.sqrt(x);
            }
            return sqrt*sqrt == x;
        }
        return false;
    }
    

    E il mio punto di riferimento: (Richiede Google caliper 0.1-rc5)

    public class SquareRootBenchmark {
        public static class Benchmark1 extends SimpleBenchmark {
            private static final int ARRAY_SIZE = 10000;
            long[] trials = new long[ARRAY_SIZE];
    
            @Override
            protected void setUp() throws Exception {
                Random r = new Random();
                for (int i = 0; i < ARRAY_SIZE; i++) {
                    trials[i] = Math.abs(r.nextLong());
                }
            }
    
    
            public int timeInternet(int reps) {
                int trues = 0;
                for(int i = 0; i < reps; i++) {
                    for(int j = 0; j < ARRAY_SIZE; j++) {
                        if(SquareRootAlgs.isPerfectSquareInternet(trials[j])) trues++;
                    }
                }
    
                return trues;   
            }
    
            public int timeDurron(int reps) {
                int trues = 0;
                for(int i = 0; i < reps; i++) {
                    for(int j = 0; j < ARRAY_SIZE; j++) {
                        if(SquareRootAlgs.isPerfectSquareDurron(trials[j])) trues++;
                    }
                }
    
                return trues;   
            }
    
            public int timeDurronTwo(int reps) {
                int trues = 0;
                for(int i = 0; i < reps; i++) {
                    for(int j = 0; j < ARRAY_SIZE; j++) {
                        if(SquareRootAlgs.isPerfectSquareDurronTwo(trials[j])) trues++;
                    }
                }
    
                return trues;   
            }
        }
    
        public static void main(String... args) {
            Runner.main(Benchmark1.class, args);
        }
    }
    

    AGGIORNAMENTO: ho creato un nuovo algoritmo più veloce in alcuni scenari, più lento in altri, ho ottenuto benchmark diversi basati su diversi input. Se calcoliamo il modulo 0xFFFFFF = 3 x 3 x 5 x 7 x 13 x 17 x 241, possiamo eliminare il 97,82% dei numeri che non possono essere quadrati. Questo può essere (fatto) in una riga, con 5 operazioni bit a bit:

    if (!goodLookupSquares[(int) ((n & 0xFFFFFFl) + ((n >> 24) & 0xFFFFFFl) + (n >> 48))]) return false;
    

    L'indice risultante è 1) il residuo, 2) il residuo + 0xFFFFFFo 3) il residuo + 0x1FFFFFE. Ovviamente, abbiamo bisogno di una tabella di ricerca per i residui modulo 0xFFFFFF, che si tratta di un file 3mb (in questo caso memorizzato come numeri decimali del testo ASCII, non ottimale ma chiaramente migliorabile con un ByteBuffere così via. Ma dal momento che è precalcolo non lo fa t molto importa. Puoi trovare il file qui (o generarlo da solo):

    public final static boolean isPerfectSquareDurronThree(long n) {
        if(n < 0) return false;
        if(n == 0) return true;
    
        long x = n;
        while((x & 0x3) == 0) x >>= 2;
        if((x & 0x7) == 1) {
            if (!goodLookupSquares[(int) ((n & 0xFFFFFFl) + ((n >> 24) & 0xFFFFFFl) + (n >> 48))]) return false;
            long sqrt;
            if(x < 410881L)
            {
                int i;
                float x2, y;
    
                x2 = x * 0.5F;
                y  = x;
                i  = Float.floatToRawIntBits(y);
                i  = 0x5f3759df - ( i >> 1 );
                y  = Float.intBitsToFloat(i);
                y  = y * ( 1.5F - ( x2 * y * y ) );
    
                sqrt = (long)(1.0F/y);
            } else {
                sqrt = (long) Math.sqrt(x);
            }
            return sqrt*sqrt == x;
        }
        return false;
    }
    

    Lo carico in un booleanarray come questo:

    private static boolean[] goodLookupSquares = null;
    
    public static void initGoodLookupSquares() throws Exception {
        Scanner s = new Scanner(new File("24residues_squares.txt"));
    
        goodLookupSquares = new boolean[0x1FFFFFE];
    
        while(s.hasNextLine()) {
            int residue = Integer.valueOf(s.nextLine());
            goodLookupSquares[residue] = true;
            goodLookupSquares[residue + 0xFFFFFF] = true;
            goodLookupSquares[residue + 0x1FFFFFE] = true;
        }
    
        s.close();
    }
    

    Esempio di runtime. Ha battuto Durron(versione uno) in ogni prova che ho corso.

     0% Scenario{vm=java, trial=0, benchmark=Internet} 40665.77 ns; ?=566.71 ns @ 10 trials
    33% Scenario{vm=java, trial=0, benchmark=Durron} 38397.60 ns; ?=784.30 ns @ 10 trials
    67% Scenario{vm=java, trial=0, benchmark=DurronThree} 36171.46 ns; ?=693.02 ns @ 10 trials
    
      benchmark   us linear runtime
       Internet 40.7 ==============================
         Durron 38.4 ============================
    DurronThree 36.2 ==========================
    
    vm: java
    trial: 0
    

    La seguente semplificazione della soluzione di maaartinus sembra radere alcuni punti percentuali dal runtime, ma non sono abbastanza bravo nel benchmarking per produrre un benchmark di cui mi possa fidare:

    long goodMask; // 0xC840C04048404040 computed below
    {
        for (int i=0; i<64; ++i) goodMask |= Long.MIN_VALUE >>> (i*i);
    }
    
    public boolean isSquare(long x) {
        // This tests if the 6 least significant bits are right.
        // Moving the to be tested bit to the highest position saves us masking.
        if (goodMask << x >= 0) return false;
        // Remove an even number of trailing zeros, leaving at most one.
        x >>= (Long.numberOfTrailingZeros(x) & (-2);
        // Repeat the test on the 6 least significant remaining bits.
        if (goodMask << x >= 0 | x <= 0) return x == 0;
        // Do it in the classical way.
        // The correctness is not trivial as the conversion from long to double is lossy!
        final long tst = (long) Math.sqrt(x);
        return tst * tst == x;
    }
    

    Vale la pena controllare come omettere il primo test,

    if (goodMask << x >= 0) return false;
    

    influenzerebbe le prestazioni.


    Dovresti liberarti della parte 2-power di N fin dall'inizio.

    2nd Edit Dovrebbe essere l'espressione magica per m sotto

    m = N - (N & (N-1));
    

    e non come scritto

    Fine della seconda modifica

    m = N & (N-1); // the lawest bit of N
    N /= m;
    byte = N & 0x0F;
    if ((m % 2) || (byte !=1 && byte !=9))
      return false;
    

    1a modifica:

    Miglioramento secondario:

    m = N & (N-1); // the lawest bit of N
    N /= m;
    if ((m % 2) || (N & 0x07 != 1))
      return false;
    

    Fine della prima modifica

    Ora continua come al solito. In questo modo, quando arrivi alla parte in virgola mobile, ti sei già liberato di tutti i numeri la cui parte a 2 potenze è dispari (circa la metà), e quindi consideri solo 1/8 di quello che è rimasto. Cioè si esegue la parte a virgola mobile sul 6% dei numeri.


    Se fai un binario per cercare di trovare la radice quadrata "giusta", puoi facilmente riconoscere se il valore che hai ottenuto è abbastanza vicino da dire:

    (n+1)^2 = n^2 + 2n + 1
    (n-1)^2 = n^2 - 2n + 1
    

    Quindi, dopo aver calcolato n^2, le opzioni sono:

    • n^2 = target : fatto, restituisce true
    • n^2 + 2n + 1 > target > n^2 : sei vicino, ma non è perfetto: restituisci falso
    • n^2 - 2n + 1 < target < n^2 : idem
    • target < n^2 - 2n + 1 : binary chop su un lower n
    • target > n^2 + 2n + 1 : binary chop su un livello più alto n

    (Scusa, questo usa ncome tua ipotesi corrente, e targetper il parametro. Scusa per la confusione!)

    Non so se sarà più veloce o meno, ma vale la pena provarci.

    EDIT: Il binary chop non deve prendere l'intero range di interi, (2^x)^2 = 2^(2x)quindi, una volta che hai trovato il bit del set più alto nel tuo target (che può essere fatto con un trucchetto, ho dimenticato esattamente come) puoi ottenere rapidamente una serie di potenziali risposte. Intendiamoci, un ingenuo binario sta ancora andando a prendere 31 o 32 iterazioni.


    Ecco il modo più semplice e conciso, anche se non so come si confronta in termini di cicli della CPU. Funziona alla grande se vuoi solo sapere se la radice è un numero intero. Se ti interessa davvero se è un numero intero, puoi anche capirlo. Ecco una funzione semplice (e pura):

    public static boolean isRootWhole(double number) {
        return Math.sqrt(number) % 1 == 0;
    }
    

    Se non hai bisogno di micro-ottimizzazione, questa risposta è migliore in termini di semplicità e manutenibilità. Se si ottengono numeri negativi, forse si vorrà utilizzare Math.abs () sull'argomento numero come argomento Math.sqrt ().

    Sulla mia CPU Intel i7-4790 da 3,6 Ghz, una corsa di questo algoritmo su 0 - 10.000.000 ha richiesto una media di 35 - 37 nanosecondi per calcolo. Ho eseguito 10 sequenze sequenziali, stampando il tempo medio trascorso su ciascuno dei dieci milioni di calcoli sqrt. Per completare la corsa totale, sono necessari poco più di 600 ms.

    Se si sta eseguendo un numero minore di calcoli, i calcoli precedenti richiedono più tempo.


    Dovrebbe essere molto più veloce usare http://en.wikipedia.org/wiki/Newton%27s_method per calcolare l' Integer Square Root , quindi quadrare questo numero e controllare, come nella tua attuale soluzione. Il metodo di Newton è la base per la soluzione di Carmack menzionata in altre risposte. Dovresti essere in grado di ottenere una risposta più rapida poiché ti interessa solo la parte intera della radice, permettendoti di interrompere l'algoritmo di approssimazione prima.

    Un'altra ottimizzazione che puoi provare: se la radice digitale di un numero non finisce in 1, 4, 7 o 9, il numero non è un quadrato perfetto. Questo può essere usato come un modo rapido per eliminare il 60% dei tuoi input prima di applicare l'algoritmo a radice quadrata più lenta.


    Questo è uguale al tuo, ma con meno operazioni. (Ho contato 9 op nel ciclo del tuo codice, incluso test e incremento i nel ciclo for e 3 assegnazioni, ma forse alcuni di quelli scompaiono quando sono codificati in ASM? Ci sono 6 opzioni nel codice qui sotto, se contate g*g>n come due (nessun incarico)).

    int isqrt(int n) {
      int g = 0x8000;
      int c = 0x8000;
      for (;;) {
        if (g*g > n) {
          g ^= c;
        }
        c >>= 1;
        if (c == 0) {
          return g;
        }
        g |= c;
      }
    }
    

    L'ho preso here È possibile eliminare un confronto se si srotola il loop e si passa al punto appropriato in base al bit non zero più alto nell'input.

    Aggiornare

    Ho pensato di più sull'uso del metodo di Newton. In teoria, il numero di bit di accuratezza dovrebbe raddoppiare per ciascuna iterazione. Ciò significa che il metodo di Newton è molto peggiore di altri suggerimenti quando ci sono pochi bit corretti nella risposta; tuttavia, la situazione cambia dove ci sono molti bit corretti nella risposta. Considerando che molti suggerimenti sembrano prendere 4 cicli per bit, ciò significa che un'iterazione del metodo di Newton (16 cicli per divisione + 1 per aggiunta + 1 per spostamento = 18 cicli) non vale la pena a meno che non fornisca più di 4 bit.

    Quindi, il mio suggerimento è di costruire 8 bit della risposta con uno dei metodi suggeriti (8 * 4 = 32 cicli), quindi eseguire un'iterazione del metodo di Newton (18 cicli) per raddoppiare il numero di bit a 16. Questo è un totale di 50 cicli (più forse un extra di 4 cicli per ottenere 9 bit prima di applicare il metodo di Newton, più forse 2 cicli per superare l'overshoot occasionalmente sperimentato dal metodo di Newton). Questo è un massimo di 56 cicli che, per quanto posso vedere, rivaleggia con gli altri suggerimenti.

    Secondo aggiornamento

    Ho codificato l'idea dell'algoritmo ibrido. Il metodo stesso di Newton non ha costi generali; basta applicare e raddoppiare il numero di cifre significative. Il problema è di avere un numero prevedibile di cifre significative prima di applicare il metodo di Newton. Per questo, abbiamo bisogno di capire dove verrà visualizzato il bit più significativo della risposta. Usando una modifica del metodo di sequenza DeBruijn veloce dato da un altro poster, possiamo eseguire tale calcolo in circa 12 cicli nella mia stima. D'altra parte, conoscere la posizione del msb della risposta accelera tutti i metodi (media, non il caso peggiore), quindi sembra che valga la pena non importa quale.

    Dopo aver calcolato il msb della risposta, eseguo un numero di round dell'algoritmo suggerito sopra, quindi lo concludo con uno o due round del metodo di Newton. Decidiamo quando eseguire il metodo di Newton con il seguente calcolo: un bit della risposta impiega circa 8 cicli in base al calcolo nei commenti; un round del metodo di Newton richiede circa 18 cicli (divisione, addizione e maiuscole e forse assegnazione), quindi dovremmo eseguire il metodo di Newton solo se ne ricaveremo almeno tre bit. Quindi per le risposte a 6 bit, possiamo eseguire il metodo lineare 3 volte per ottenere 3 bit significativi, quindi eseguire il metodo di Newton 1 volta per ottenere un altro 3. Per le risposte a 15 bit, eseguiamo il metodo lineare 4 volte per ottenere 4 bit, quindi Newton metodo due volte per ottenere un altro 4 poi un altro 7. E così via.

    Questi calcoli dipendono dal sapere esattamente quanti cicli sono necessari per ottenere un po 'dal metodo lineare rispetto a quanti sono richiesti dal metodo di Newton. Se la "economia" cambia, ad esempio, scoprendo un modo più rapido per costruire i bit in modo lineare, la decisione su quando invocare il metodo di Newton cambierà.

    Ho srotolato i loop e ho implementato le decisioni come interruttori, che spero si traducano in ricerche di tabelle veloci in assembly. Non sono assolutamente sicuro di avere il numero minimo di cicli in ogni caso, quindi forse è possibile effettuare ulteriori regolazioni. Ad esempio, per s = 10, puoi provare a ottenere 5 bit, quindi applicare il metodo di Newton una volta anziché due.

    Ho testato l'algoritmo accuratamente per correttezza. Alcuni ulteriori piccoli aumenti sono possibili se in alcuni casi sei disposto ad accettare risposte leggermente errate. Almeno due cicli vengono utilizzati dopo aver applicato il metodo di Newton per correggere un errore di off-by-one che si verifica con i numeri del modulo m^2-1 . E un ciclo viene utilizzato per verificare l'input 0 all'inizio, poiché l'algoritmo non può gestire quell'input. Se sai che non prenderai mai la radice quadrata di zero, puoi eliminare quel test. Infine, se nella risposta sono necessari solo 8 bit significativi, è possibile eliminare uno dei calcoli del metodo di Newton.

    #include <inttypes.h>
    #include <stdint.h>
    #include <stdbool.h>
    #include <stdio.h>
    
    uint32_t isqrt1(uint32_t n);
    
    int main() {
      uint32_t n;
      bool it_works = true;
      for (n = 0; n < UINT32_MAX; ++n) {
        uint32_t sr = isqrt1(n);
        if ( sr*sr > n || ( sr < 65535 && (sr+1)*(sr+1) <= n )) {
          it_works = false;
          printf("isqrt(%" PRIu32 ") = %" PRIu32 "\n", n, sr);
        }
      }
      if (it_works) {
        printf("it works\n");
      }
      return 0;
    }
    
    /* table modified to return shift s to move 1 to msb of square root of x */
    /*
    static const uint8_t debruijn32[32] = {
        0, 31, 9, 30, 3,  8, 13, 29,  2,  5,  7, 21, 12, 24, 28, 19,
        1, 10, 4, 14, 6, 22, 25, 20, 11, 15, 23, 26, 16, 27, 17, 18
    };
    */
    
    static const uint8_t debruijn32[32] = {
      15,  0, 11, 0, 14, 11, 9, 1, 14, 13, 12, 5, 9, 3, 1, 6,
      15, 10, 13, 8, 12,  4, 3, 5, 10,  8,  4, 2, 7, 2, 7, 6
    };
    
    /* based on CLZ emulation for non-zero arguments, from
     * http://.com/questions/23856596/counting-leading-zeros-in-a-32-bit-unsigned-integer-with-best-algorithm-in-c-pro
     */
    uint8_t shift_for_msb_of_sqrt(uint32_t x) {
      x |= x >>  1;
      x |= x >>  2;
      x |= x >>  4;
      x |= x >>  8;
      x |= x >> 16;
      x++;
      return debruijn32 [x * 0x076be629 >> 27];
    }
    
    uint32_t isqrt1(uint32_t n) {
      if (n==0) return 0;
    
      uint32_t s = shift_for_msb_of_sqrt(n);
      uint32_t c = 1 << s;
      uint32_t g = c;
    
      switch (s) {
      case 9:
      case 5:
        if (g*g > n) {
          g ^= c;
        }
        c >>= 1;
        g |= c;
      case 15:
      case 14:
      case 13:
      case 8:
      case 7:
      case 4:
        if (g*g > n) {
          g ^= c;
        }
        c >>= 1;
        g |= c;
      case 12:
      case 11:
      case 10:
      case 6:
      case 3:
        if (g*g > n) {
          g ^= c;
        }
        c >>= 1;
        g |= c;
      case 2:
        if (g*g > n) {
          g ^= c;
        }
        c >>= 1;
        g |= c;
      case 1:
        if (g*g > n) {
          g ^= c;
        }
        c >>= 1;
        g |= c;
      case 0:
        if (g*g > n) {
          g ^= c;
        }
      }
    
      /* now apply one or two rounds of Newton's method */
      switch (s) {
      case 15:
      case 14:
      case 13:
      case 12:
      case 11:
      case 10:
        g = (g + n/g) >> 1;
      case 9:
      case 8:
      case 7:
      case 6:
        g = (g + n/g) >> 1;
      }
    
      /* correct potential error at m^2-1 for Newton's method */
      return (g==65536 || g*g>n) ? g-1 : g;
    }
    

    Nella prova leggera sulla mia macchina (che certamente non è come la tua), la nuova routine isqrt1 funziona in media circa il 40% più velocemente rispetto alla precedente routine isqrt ho dato.





    java math optimization perfect-square