Najszybszy sposób ustalenia, czy pierwiastek kwadratowy z liczby całkowitej jest liczbą całkowitą

1453

Szukam najszybszego sposobu ustalenia, czy longwartość jest idealnym kwadratem (tzn. Jej pierwiastek kwadratowy jest inną liczbą całkowitą):

  1. Zrobiłem to w prosty sposób, korzystając z wbudowanej Math.sqrt() funkcji, ale zastanawiam się, czy istnieje sposób, aby to zrobić szybciej, ograniczając się do domeny zawierającej tylko liczby całkowite.
  2. Utrzymywanie tabeli odnośników jest niepraktyczne (ponieważ istnieje około 2 31,5 liczb całkowitych, których kwadrat jest mniejszy niż 2 63 ).

Oto bardzo prosty i bezpośredni sposób, w jaki to teraz robię:

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

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

Uwaga: Używam tej funkcji w wielu problemach z Project Euler . Dlatego nikt inny nie będzie musiał utrzymywać tego kodu. Tego rodzaju mikrooptymalizacja może w rzeczywistości coś zmienić, ponieważ częścią wyzwania jest wykonanie każdego algorytmu w mniej niż minutę, a funkcja ta będzie musiała zostać wywołana miliony razy w przypadku niektórych problemów.

Próbowałem różnych rozwiązań tego problemu:

  • Po wyczerpujących testach odkryłem, że dodawanie 0.5do wyniku Math.sqrt () nie jest konieczne, przynajmniej nie na moim komputerze.
  • Szybko odwrotność pierwiastka kwadratowego był szybszy, ale to dało nieprawidłowych wyników dla n> = 410881. Jednak, jak sugeruje BobbyShaftoe , możemy użyć hack FISR dla n <410881.
  • Metoda Newtona była o wiele wolniejsza niż Math.sqrt(). Jest tak prawdopodobnie dlatego, że Math.sqrt()używa czegoś podobnego do metody Newtona, ale zaimplementowanej w sprzęcie, dzięki czemu jest znacznie szybsza niż w Javie. Ponadto Metoda Newtona nadal wymagała użycia podwójnych.
  • Zmodyfikowana metoda Newtona, która wykorzystywała kilka sztuczek, aby zaangażować tylko matematykę całkowitą, wymagała kilku hacków, aby uniknąć przepełnienia (chcę, aby ta funkcja działała ze wszystkimi dodatnimi liczbami całkowitymi ze znakiem 64-bitowym), i nadal była wolniejsza niż Math.sqrt().
  • Binarny kotlet był jeszcze wolniejszy. Ma to sens, ponieważ przecięcie binarne wymaga średnio 16 przejść, aby znaleźć pierwiastek kwadratowy liczby 64-bitowej.
  • Według testów Johna używanie orinstrukcji w C ++ jest szybsze niż używanie a switch, ale w Javie i C # wydaje się, że nie ma różnicy między ori switch.
  • Próbowałem także utworzyć tabelę odnośników (jako prywatną statyczną tablicę 64 wartości boolowskich). Następnie zamiast przełącznika lub orinstrukcji, powiedziałbym tylko if(lookup[(int)(n&0x3F)]) { test } else return false;. Ku mojemu zaskoczeniu było to (tylko nieco) wolniejsze. Wynika to z faktu, że granice tablic są sprawdzane w Javie .
java math optimization perfect-square Kip
źródło
21
To jest kod Java, gdzie int == 32 bity i długi == 64 bity, i oba są podpisane.
Kip
14
@Shreevasta: Przeprowadziłem testy na dużych wartościach (większych niż 2 ^ 53), a twoja metoda daje fałszywie pozytywne wyniki. Pierwszy napotkany jest dla n = 9007199326062755, który nie jest kwadratem idealnym, ale jest zwracany jako jeden.
Kip
37
Nie nazywaj tego „hackem Johna Carmacka”. Nie wymyślił tego.
user9282,
84
@mama - Być może, ale przypisuje mu się to. Henry Ford nie wynalazł samochodu, Wright Bros. nie wynalazł samolotu, a Galleleo nie był pierwszym, który odkrył, że Ziemia obraca się wokół Słońca ... świat składa się ze skradzionych wynalazków (i miłość).
Robert Fraser
4
Możesz uzyskać niewielki wzrost prędkości w „szybkim błędzie”, używając czegoś podobnego ((1<<(n&15))|65004) != 0, zamiast trzech osobnych kontroli.
Nabb

Odpowiedzi:

735

Opracowałem metodę, która działa ~ 35% szybciej niż twój 6-bitowy kod + Carmack + kod sqrt, przynajmniej z moim procesorem (x86) i językiem programowania (C / C ++). Twoje wyniki mogą się różnić, szczególnie dlatego, że nie wiem, jak będzie się grał czynnik Java.

Moje podejście jest trojakie:

  1. Najpierw odfiltruj oczywiste odpowiedzi. Obejmuje to liczby ujemne i patrząc na ostatnie 4 bity. (Odkryłem, że spojrzenie na ostatnie sześć nie pomogło.) Odpowiadam również na „tak” (czytając poniższy kod, zauważ, że moje dane wejściowe są int64 x). 
    if( x < 0 || (x&2) || ((x & 7) == 5) || ((x & 11) == 8) )
    return false;
if( x == 0 )
    return true;
  2. Następnie sprawdź, czy jest to moduł kwadratowy 255 = 3 * 5 * 17. Ponieważ jest to iloczyn trzech różnych liczb pierwszych, tylko około 1/8 reszt mod 255 to kwadraty. Jednak z mojego doświadczenia wynika, że ​​wywołanie operatora modulo (%) kosztuje więcej niż korzyść, którą się otrzymuje, więc używam sztuczek bitowych obejmujących 255 = 2 ^ 8-1 do obliczenia pozostałości. (Na lepsze lub gorsze, nie używam sztuczki polegającej na odczytywaniu poszczególnych bajtów ze słowa, tylko bitowe i zmiany). 
    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.
    Aby faktycznie sprawdzić, czy pozostałość jest kwadratem, szukam odpowiedzi w obliczonej wcześniej tabeli. 
    if( bad255[y] )
    return false;
// However, I just use a table of size 512
  3. Na koniec spróbuj obliczyć pierwiastek kwadratowy przy użyciu metody podobnej do lematu Hensela . (Nie sądzę, aby można go było zastosować bezpośrednio, ale działa z pewnymi modyfikacjami.) Przedtem dzielę wszystkie moce 2 za pomocą wyszukiwania binarnego: 
    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;
    W tym momencie, aby nasza liczba była kwadratem, musi to być 1 mod 8. 
    if((x & 7) != 1)
    return false;
    Podstawowa struktura lematu Hensela jest następująca. (Uwaga: nieprzetestowany kod; jeśli nie działa, spróbuj t = 2 lub 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.
    Chodzi o to, że przy każdej iteracji dodajesz jeden bit do r, „bieżący” pierwiastek kwadratowy z x; każdy pierwiastek kwadratowy jest dokładnym modułem o coraz większej potędze 2, mianowicie t / 2. Na końcu r i t / 2-r będą pierwiastkami kwadratowymi x modulo t / 2. (Zauważ, że jeśli r jest pierwiastkiem kwadratowym z x, to tak samo jest z -r. To prawda, nawet liczby modulo, ale uwaga, modulo niektóre liczby, rzeczy mogą mieć nawet więcej niż 2 pierwiastki kwadratowe; w szczególności obejmuje to potęgi 2. ) Ponieważ nasz rzeczywisty pierwiastek kwadratowy jest mniejszy niż 2 ^ 32, w tym momencie możemy faktycznie sprawdzić, czy r lub t / 2-r są rzeczywistymi pierwiastkami kwadratowymi. W moim aktualnym kodzie używam następującej zmodyfikowanej pętli: 
    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) );
    Przyspieszenie tutaj uzyskuje się na trzy sposoby: wstępnie obliczoną wartość początkową (równoważną ~ 10 iteracjom pętli), wcześniejsze wyjście z pętli i pominięcie niektórych wartości t. Na ostatnią część patrzę z = r - x * xi ustawiam t na największą potęgę 2 dzielących z odrobiną sztuczki. To pozwala mi pominąć wartości t, które i tak nie wpłynęłyby na wartość r. Wstępnie obliczona wartość początkowa w moim przypadku wybiera moduł „najmniejszego dodatniego” pierwiastka kwadratowego 8192.

Nawet jeśli ten kod nie działa szybciej dla Ciebie, mam nadzieję, że podoba Ci się niektóre zawarte w nim pomysły. Następuje pełny, przetestowany kod, w tym wstępnie obliczone tabele.

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;
}
A. Rex
źródło
5
Łał! Spróbuję przekonwertować to na Javę i przeprowadzić porównanie, a także sprawdzić dokładność wyników. Dam ci znać, co znajdę.
Kip
79
Wow, to jest piękne. Widziałem wcześniej Hensela podnoszącego (obliczanie pierwiastków wielomianów modulo a prime), ale nawet nie zdawałem sobie sprawy, że lemat można ostrożnie obniżyć do samego końca, obliczając pierwiastki kwadratowe liczb; to jest ... podnoszący na duchu :)
ShreevatsaR
3
@nightcracker Nie ma. 9 < 0 => false, 9&2 => 0, 9&7 == 5 => false, 9&11 == 8 => false.
primo
53
Maartinus opublikował 2x szybsze rozwiązanie (i znacznie krótsze) poniżej, nieco później, które nie wydaje się zbytnio kochane.
Jason C
3
Wygląda na to, że znaczną przewagę prędkości w różnych rozwiązaniach uzyskuje się przez odfiltrowanie oczywistych kwadratów. Czy ktoś ocenił sytuację, w której odfiltrowano za pomocą rozwiązania Maartinus, a następnie po prostu użył funkcji sqrt, ponieważ jest to funkcja wbudowana?
user1914292,
377

Jestem spóźniony na przyjęcie, ale mam nadzieję, że udzielę lepszej odpowiedzi; krótszy i (zakładając, że mój test jest poprawny) również znacznie szybszy .

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

Pierwszy test szybko wychwytuje większość elementów innych niż kwadraty. Używa 64-elementowej tabeli zapakowanej w długi, więc nie ma żadnych kosztów dostępu do tablicy (pośrednie i sprawdzanie granic). Dla równomiernie losowego longprawdopodobieństwo, że skończy się tutaj, wynosi 81,25%.

Drugi test wyłapuje wszystkie liczby o nieparzystej liczbie dwójkowej w rozkładzie na czynniki. Metoda Long.numberOfTrailingZerosjest bardzo szybka, ponieważ przekształca JIT w pojedynczą instrukcję i86.

Po usunięciu końcowych zer trzeci test obsługuje liczby kończące się na 011, 101 lub 111 w systemie binarnym, które nie są idealnymi kwadratami. Dba również o liczby ujemne, a także obsługuje 0.

Ostatni test powraca do doublearytmetyki. Podobnie jak doublema tylko 53 bity mantysy, konwersja z longna doubleobejmuje zaokrąglanie dużych wartości. Niemniej jednak test jest poprawny (chyba że dowód jest błędny).

Próba wprowadzenia pomysłu mod255 nie powiodła się.

maaartinus
źródło
3
To ukryte maskowanie wartości przesunięcia jest trochę ... złe. Czy masz pojęcie, dlaczego tak jest w specyfikacji Java?
dfeuer
6
@dfeuer Wydaje mi się, że istnieją dwa powody: 1. Przesunięcie o więcej nie ma sensu. 2. To tak, jakby HW działa i każdy, kto używa operacji bitowych, jest zainteresowany wydajnością, więc robienie czegokolwiek innego byłoby niewłaściwe. -goodMask badanie to robi, ale robi to przed właściwym przesunięciem. Musisz więc to powtórzyć, ale w ten sposób jest to prostsze i AFAIK trochę szybsze i równie dobre.
maaartinus,
3
@dfeuer Dla testu porównawczego ważne jest jak najszybsze udzielenie odpowiedzi, a sama końcowa liczba zerowa nie daje odpowiedzi; to tylko krok przygotowawczy. i86 / amd64 to zrobić. Nie mam pojęcia o małych procesorach w telefonach komórkowych, ale w najgorszym przypadku Java musi wygenerować dla nich instrukcję AND, która z pewnością jest prostsza niż na odwrót.
maaartinus
2
@Sebastian chyba lepiej Test: if ((x & (7 | Integer.MIN_VALUE)) != 1) return x == 0;.
maaartinus
4
„Ponieważ podwójny ma tylko 56-bitową mantysę” -> Powiedziałbym, że bardziej prawdopodobne jest 53-bitowe . Również
chux
132

Będziesz musiał przeprowadzić testy porównawcze. Najlepszy algorytm będzie zależeć od rozkładu twoich danych wejściowych.

Twój algorytm może być prawie optymalny, ale możesz zrobić szybkie sprawdzenie, aby wykluczyć pewne możliwości przed wywołaniem procedury pierwiastka kwadratowego. Na przykład spójrz na ostatnią cyfrę swojego numeru szesnastkowego, wykonując nieco „i”. Idealne kwadraty mogą kończyć się na 0, 1, 4 lub 9 na podstawie 16, więc dla 75% twoich danych wejściowych (zakładając, że są one równomiernie rozmieszczone) możesz uniknąć wezwania do pierwiastka kwadratowego w zamian za bardzo szybkie kręcenie bitów.

Kip przeprowadził analizę porównawczą następującego kodu implementującego sztuczkę szesnastkową. Podczas testowania liczb od 1 do 100 000 000 ten kod działał dwa razy szybciej niż oryginał.

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

Kiedy testowałem analogiczny kod w C ++, faktycznie działał wolniej niż oryginał. Kiedy jednak wyeliminowałem instrukcję switch, sztuczka szesnastkowa ponownie sprawia, że ​​kod jest dwa razy szybszy.

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

Wyeliminowanie instrukcji switch miało niewielki wpływ na kod C #.

John D. Cook
źródło
to całkiem sprytne ... nie pomyślałbym o tym
warren
Dobra uwaga na temat bitów końcowych. Spróbowałbym połączyć ten test z kilkoma innymi uwagami tutaj.
PeterAllenWebb
3
Doskonałe rozwiązanie. Zastanawiasz się, jak to wymyśliłeś? Czy jest to ustalona zasada, czy po prostu coś wymyśliłeś? : D
Jeel Shah,
3
@ LarsH Nie ma potrzeby dodawania 0,5, zobacz moje rozwiązanie, aby uzyskać link do dowodu.
maaartinus
2
@JerryGoyal To zależy od kompilatora i wartości spraw. W idealnym kompilatorze przełącznik jest zawsze co najmniej tak szybki, jak w innym przypadku. Ale kompilatory nie są idealne, więc najlepiej jest wypróbować to, co zrobił John.
fishinear
52

Myślałem o okropnych czasach, które spędziłem na kursie analizy numerycznej.

A potem pamiętam, że ta funkcja krążyła po sieci z kodu źródłowego 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;
}

Który zasadniczo oblicza pierwiastek kwadratowy, używając funkcji aproksymacji Newtona (nie pamiętam dokładnej nazwy).

Powinien być użyteczny, a może nawet szybszy, pochodzi z jednej z fenomenalnych gier oprogramowania id!

Jest napisany w C ++, ale nie powinno być zbyt trudne ponowne użycie tej samej techniki w Javie, gdy tylko wpadniesz na pomysł:

Pierwotnie znalazłem na: http://www.codemaestro.com/reviews/9

Metoda Newtona wyjaśniona na wikipedii: http://en.wikipedia.org/wiki/Newton%27s_method

Możesz skorzystać z linku, aby uzyskać więcej informacji o tym, jak to działa, ale jeśli nie przejmujesz się tym, to mniej więcej to pamiętam z czytania bloga i z kursu analizy numerycznej:

  • * (long*) &yjest w zasadzie fast funkcja konwersji do operacji tak długo całkowita może być stosowany na surowych bajtów.
  • 0x5f3759df - (i >> 1);jest to wartość nasion uprzednio obliczone przez funkcję aproksymacji.
  • * (float*) &ikonwertuje wartość z powrotem do zmiennoprzecinkowych.
  • y = y * ( threehalfs - ( x2 * y * y ) )linia bascially iteracje wartości ponad funkcję ponownie.

Funkcja aproksymacji podaje bardziej precyzyjne wartości, im bardziej iterujesz funkcję nad wynikiem. W przypadku Quake'a jedna iteracja jest „wystarczająco dobra”, ale jeśli to nie było dla ciebie ... to możesz dodać tyle iteracji, ile potrzebujesz.

Powinno to być szybsze, ponieważ zmniejsza liczbę operacji dzielenia wykonanych naiwnym kwadratowym rootowaniu do zwykłego dzielenia przez 2 (w rzeczywistości * 0.5Foperację mnożenia) i zastępuje ją kilkoma ustalonymi liczbami operacji mnożenia.

czakryt
źródło
9
Należy zauważyć, że zwraca 1 / sqrt (liczba), a nie sqrt (liczba). Przeprowadziłem pewne testy, ale nie udaje się to przy n = 410881: magiczna formuła Johna Carmacka zwraca 642.00104, gdy faktyczny pierwiastek kwadratowy to 641.
Kip
11
Możesz spojrzeć na artykuł Chrisa Lomontsa o szybkich odwrotnych pierwiastkach kwadratowych: lomont.org/Math/Papers/2003/InvSqrt.pdf Używa tej samej techniki co tutaj, ale z inną liczbą magiczną. Artykuł wyjaśnia, dlaczego wybrano magiczną liczbę.
4
Ponadto, beyond3d.com/content/articles/8 i beyond3d.com/content/articles/15 rzuciły nieco światła na pochodzenie tej metody. Często przypisuje się to Johnowi Carmackowi, ale wygląda na to, że oryginalny kod został (być może) napisany przez Gary'ego Tarolli, Grega Walsha i prawdopodobnie innych.
3
W Javie nie można także pisać na komputerze typu float i ints.
Antymon
10
@Antimony, kto mówi? FloatToIntBits i IntToFloatBits istnieją już od wersji Java 1.0.2.
corsiKa
38

Nie jestem pewien, czy byłoby to szybsze, czy nawet dokładne, ale możesz użyć algorytmu Magical Square Root Johna Carmacka , aby szybciej rozwiązać pierwiastek kwadratowy. Prawdopodobnie mógłbyś łatwo przetestować to dla wszystkich możliwych 32-bitowych liczb całkowitych i sprawdzić, czy faktycznie masz poprawne wyniki, ponieważ jest to tylko przybliżenie. Jednak teraz, gdy o tym myślę, użycie podwójnych jest również przybliżone, więc nie jestem pewien, jak to by się stało.

Kibbee
źródło
10
Uważam, że sztuczka Carmacka jest obecnie bezcelowa. Wbudowana instrukcja sqrt jest o wiele szybsza niż kiedyś, więc lepiej jest po prostu wykonać regularny pierwiastek kwadratowy i sprawdzić, czy wynikiem jest liczba całkowita. Jak zawsze, sprawdź to.
lipiec
4
Ta przerwa zaczyna się od n = 410881, magiczna formuła Johna Carmacka zwraca 642.00104, gdy faktyczny pierwiastek kwadratowy wynosi 641.
Kip
11
Niedawno użyłem sztuczki Carmacka w grze Java i była ona bardzo skuteczna, dając przyspieszenie około 40%, więc nadal jest przydatna, przynajmniej w Javie.
finnw
3
@Robert Fraser Tak + 40% w ogólnej liczbie klatek na sekundę. Gra posiadała układ fizyki cząstek, który zajmował prawie wszystkie dostępne cykle procesora, zdominowany przez funkcję pierwiastka kwadratowego i funkcję zaokrąglania do najbliższej liczby całkowitej (którą zoptymalizowałem również przy użyciu podobnego nieco hackowania).
finnw
5
Link jest zepsuty.
Pixar
36

Jeśli wykonasz dwójkę binarną, aby znaleźć „właściwy” pierwiastek kwadratowy, możesz dość łatwo wykryć, czy wartość, którą masz, jest wystarczająco bliska, aby powiedzieć:

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

Po obliczeniu n^2opcje są następujące:

  • n^2 = target: gotowe, zwróć wartość true
  • n^2 + 2n + 1 > target > n^2 : jesteś blisko, ale to nie jest idealne: zwróć false
  • n^2 - 2n + 1 < target < n^2 : to samo
  • target < n^2 - 2n + 1 : binarny kotlet na niższym poziomie n
  • target > n^2 + 2n + 1 : dwójkowy na wyższym n

(Przepraszamy, ten parametr jest używany njako bieżące przypuszczenie i targetza parametr. Przepraszamy za zamieszanie!)

Nie wiem, czy to będzie szybsze, czy nie, ale warto spróbować.

EDYCJA: Binarny rąbek nie musi również przyjmować całego zakresu liczb całkowitych, (2^x)^2 = 2^(2x)więc gdy już znajdziesz bit najwyższego zestawu w swoim celu (można to zrobić za pomocą sztuczki polegającej na kręceniu się; zapominam dokładnie jak) możesz szybko uzyskać szereg potencjalnych odpowiedzi. Pamiętaj, że naiwny binarny kotlet nadal będzie wymagał tylko 31 lub 32 iteracji.

Jon Skeet
źródło
Moje pieniądze są na tego rodzaju podejściu. Unikaj wywoływania funkcji sqrt (), ponieważ oblicza ona pierwiastek kwadratowy, a potrzebujesz tylko kilku pierwszych cyfr.
PeterAllenWebb
3
Z drugiej strony, jeśli zmiennoprzecinkowa jest wykonywana w dedykowanej jednostce FP, może ona używać wszelkiego rodzaju zabawnych sztuczek. Nie chciałbym stawiać na to bez testu porównawczego :) (mogę spróbować tego wieczoru jednak w C #, tylko po to, żeby zobaczyć ...)
Jon Skeet,
8
Sprzętowe sqrts są obecnie dość szybkie.
Adam Rosenfield
24

Przeprowadziłem własną analizę kilku algorytmów w tym wątku i opracowałem kilka nowych wyników. Możesz zobaczyć te stare wyniki w historii edycji tej odpowiedzi, ale nie są one dokładne, ponieważ popełniłem błąd i straciłem czas na analizę kilku algorytmów, które nie są blisko. Jednak wyciągając wnioski z kilku różnych odpowiedzi, mam teraz dwa algorytmy, które miażdżą „zwycięzcę” tego wątku. Oto podstawowa rzecz, którą robię inaczej niż wszyscy inni:

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

Jednak ten prosty wiersz, który przez większość czasu dodaje jedną lub dwie bardzo szybkie instrukcje, znacznie upraszcza switch-caseinstrukcję w jedną instrukcję if. Może jednak zwiększyć czas działania, jeśli wiele z testowanych liczb ma znaczącą moc dwóch czynników.

Poniższe algorytmy są następujące:

  • Internet - opublikowana odpowiedź Kipa
  • Durron - Moja zmodyfikowana odpowiedź przy użyciu odpowiedzi jednoprzebiegowej jako podstawy
  • DurronTwo - Moja zmodyfikowana odpowiedź przy użyciu odpowiedzi dwuprzebiegowej (autor: @JohnnyHeggheim), z kilkoma drobnymi modyfikacjami.

Oto przykładowe środowisko wykonawcze, jeśli liczby są generowane przy użyciu 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

Oto przykładowe środowisko uruchomieniowe, jeśli działa tylko na pierwszych milionach długości:

 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

Jak widać, DurronTwolepiej sprawdza się w przypadku dużych nakładów, ponieważ bardzo często korzysta z magicznej sztuczki, ale staje się nieczytelny w porównaniu z pierwszym algorytmem i Math.sqrtponieważ liczby są znacznie mniejsze. Tymczasem prostsze Durronjest ogromnym zwycięzcą, ponieważ nigdy nie musi dzielić 4 razy wiele razy w pierwszym milionie liczb.

Oto 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;
}

I 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;
}

I moja uprząż porównawcza: (wymaga suwmiarki Google 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);
    }
}

AKTUALIZACJA: Stworzyłem nowy algorytm, który jest szybszy w niektórych scenariuszach, wolniejszy w innych, otrzymałem różne testy porównawcze na podstawie różnych danych wejściowych. Obliczając modulo 0xFFFFFF = 3 x 3 x 5 x 7 x 13 x 17 x 241, możemy wyeliminować 97,82% liczb, które nie mogą być kwadratami. Można to zrobić (w pewnym sensie) w jednym wierszu za pomocą 5 operacji bitowych:

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

Otrzymany wskaźnik to albo 1) reszta, 2) reszta + 0xFFFFFF, albo 3) reszta + 0x1FFFFFE. Oczywiście musimy mieć tabelę przeglądową dla reszt modulo 0xFFFFFF, która ma rozmiar około 3 MB (w tym przypadku jest przechowywana jako liczby dziesiętne tekstu ascii, nie jest optymalna, ale można ją łatwo poprawić za pomocą ByteBuffera itd. Ale ponieważ jest to wstępne obliczenie , nie robi tego) to ważne. Możesz znaleźć plik tutaj (lub sam go wygenerować):

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

Ładuję go do booleantablicy takiej jak ta:

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();
}

Przykładowe środowisko wykonawcze. Pokonał Durron(wersja pierwsza) w każdym badaniu, które prowadziłem.

 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
durron597
źródło
3
Olbrzymi stół przeglądowy nie wydaje się dobrym pomysłem. Brak pamięci podręcznej jest wolniejszy (~ 100 do 150 cykli) niż sprzętowa instrukcja sqrt x86 (~ 20 cykli). Jeśli chodzi o przepustowość, możesz wytrzymać wiele wyjątkowych braków pamięci podręcznej, ale nadal eksmitujesz inne przydatne dane. Ogromna tablica przeglądowa byłaby tego warta tylko wtedy, gdyby była DUŻO szybsza niż jakakolwiek inna opcja, a ta funkcja była głównym czynnikiem wpływającym na wydajność całego programu.
Peter Cordes,
1
@ SwissFrank: czy sprawdzanie kwadratów jest jedyną rzeczą, którą robi Twój program? Tabela odnośników może wyglądać dobrze na mikrobenchmycie, który wywołuje ją wielokrotnie w ciasnej pętli, ale w prawdziwym programie, który ma inne dane w zestawie roboczym, nie jest dobrze.
Peter Cordes
1
Bitmapę bitów 0x1FFFFFE trwa 4 mega bajtów , jeśli jest przechowywany w postaci upakowanej bitmapy. Pamięć podręczna L3 trafiona na nowoczesny pulpit Intela ma> 40 cykli opóźnienia i gorzej na dużym Xeonie; dłuższe niż sprzętowe sqrt + opóźnienie wielokrotne. Jeśli jest przechowywany jako mapa bajtowa z 1 bajtem na wartość, ma około 32 MB; większy niż pamięć podręczna L3 czegokolwiek poza wielordzeniowym Xeonem, gdzie wszystkie rdzenie dzielą jedną ogromną pamięć podręczną. Więc jeśli twoje dane wejściowe mają jednolity losowy rozkład w wystarczająco dużym zakresie danych wejściowych, otrzymasz wiele braków pamięci podręcznej L2, nawet w ciasnej pętli. (prywatny procesor L2 na procesor Intel ma tylko 256 KB przy opóźnieniu
Peter Cordes
1
@SwissFrank: Och, jeśli wszystko, co robisz, to sprawdzanie roota, to jest to możliwe dzięki bitmapie, aby uzyskać trafienia w L3. Patrzyłem na opóźnienia, ale wiele chybień może być w locie jednocześnie, więc przepustowość jest potencjalnie dobra. Przepustowość OTOH, SIMD, sqrtpsa nawet sqrtpd(podwójna precyzja) nie jest tak zła na Skylake, ale niewiele lepsza niż opóźnienie na starych procesorach. W każdym razie 7-cpu.com/cpu/Haswell.html ma fajne liczby eksperymentalne i strony dla innych procesorów. Pdf przewodnik mikroprocesora Agner Fog zawiera pewne opóźnienia w pamięci podręcznej dla uarches Intela i AMD: agner.org/optimize
Peter Cordes
1
Korzystanie z SIMD x86 z Java jest problemem, a do czasu dodania kosztu konwersji int-> fp i fp-> int jest prawdopodobne, że bitmapa mogłaby być lepsza. Potrzebujesz doubleprecyzji, aby uniknąć zaokrąglania liczby całkowitej poza zakresem + -2 ^ 24 (więc liczba całkowita 32-bitowa może znajdować się poza tym), i sqrtpdjest wolniejsza niż, sqrtpsa także przetwarza tylko połowę liczby elementów na instrukcję (na wektor SIMD) .
Peter Cordes
18

Zastosowanie metody Newtona powinno być znacznie szybsze , aby obliczyć pierwiastek kwadratowy z liczby całkowitej , a następnie wyprostować tę liczbę i sprawdzić, tak jak ma to miejsce w obecnym rozwiązaniu. Metoda Newtona jest podstawą rozwiązania Carmack wspomnianego w kilku innych odpowiedziach. Powinieneś być w stanie uzyskać szybszą odpowiedź, ponieważ interesuje Cię tylko całkowita liczba części root, co pozwala wcześniej zatrzymać algorytm aproksymacji.

Kolejna optymalizacja, którą możesz wypróbować: Jeśli pierwiastek cyfrowy liczby nie kończy się na 1, 4, 7 lub 9, liczba nie jest idealnym kwadratem. Można to wykorzystać jako szybki sposób na wyeliminowanie 60% danych wejściowych przed zastosowaniem algorytmu wolniejszego pierwiastka kwadratowego.

Bill jaszczurka
źródło
1
Cyfrowy pierwiastek jest ściśle obliczeniowo równoważny modulo, dlatego należy go rozważyć wraz z innymi metodami modulo, takimi jak mod 16 i mod 255.
Christian Oudard
1
Czy jesteś pewien, że cyfrowy root jest równoważny modulo? Wydaje się, że jest to coś zupełnie innego, jak wyjaśnia link. Zauważ, że lista zawiera 1,4,7,9, a nie 1,4,5,9.
Fractaly,
1
Cyfrowy pierwiastek w systemie dziesiętnym jest równoważny z użyciem modulo 9 (studnia dr (n) = 1 + ((n-1) mod 9); więc również niewielkie przesunięcie). Liczby 0,1,4,5,9 odnoszą się do modułu 16, a 0, 1, 4, 7 odnoszą się do modułu 9 - co odpowiada 1, 4, 7, 9 cyfrowemu pierwiastkowi.
Hans Olsson,
16

Chcę, aby ta funkcja działała ze wszystkimi dodatnimi liczbami całkowitymi ze znakiem 64-bitowym

Math.sqrt()działa z podwójnymi jako parametrami wejściowymi, więc nie uzyskasz dokładnych wyników dla liczb całkowitych większych niż 2 ^ 53 .

mrzl
źródło
5
Właściwie przetestowałem odpowiedź na wszystkich idealnych kwadratach większych niż 2 ^ 53, a także na wszystkich liczbach od 5 poniżej każdego idealnego kwadratu do 5 powyżej każdego idealnego kwadratu i otrzymuję prawidłowy wynik. (błąd zaokrąglania jest korygowany, gdy zaokrąglam odpowiedź sqrt do długiej, a następnie kwadrat tej wartości i porównuję)
Kip
2
@Kip: Chyba udowodniłem, że to działa .
maaartinus
Wyniki nie są idealnie dokładne, ale dokładniejsze niż mogłoby się wydawać. Jeśli założymy co najmniej 15 dokładnych cyfr po konwersji na podwójną i po pierwiastku kwadratowym, to wystarczy, ponieważ nie potrzebujemy więcej niż 11: 10 cyfr dla 32-bitowego pierwiastka kwadratowego i mniej niż 1 dla miejsca dziesiętnego, ponieważ +0,5 zaokrągla do najbliższego.
mwfearnley
3
Math.sqrt () nie jest całkowicie dokładne, ale nie musi. W pierwszym poście tst jest liczbą całkowitą zbliżoną do sqrt (N). Jeśli N nie jest kwadratem, to tst * tst! = N, bez względu na wartość tst. Jeśli N jest kwadratem idealnym, to sqrt (N) <2 ^ 32, i dopóki sqrt (N) jest obliczany z błędem <0,5, mamy się dobrze.
gnasher729
13

Dla przypomnienia, innym podejściem jest wykorzystanie pierwotnego rozkładu. Jeśli każdy czynnik rozkładu jest parzysty, liczba jest idealnym kwadratem. Chcecie więc sprawdzić, czy liczbę można rozłożyć na iloczyn kwadratów liczb pierwszych. Oczywiście nie trzeba uzyskiwać takiego rozkładu, aby sprawdzić, czy on istnieje.

Najpierw zbuduj tabelę kwadratów liczb pierwszych, które są mniejsze niż 2 ^ 32. Jest to znacznie mniej niż tabela wszystkich liczb całkowitych do tego limitu.

Rozwiązanie byłoby wtedy takie:

boolean isPerfectSquare(long number)
{
    if (number < 0) return false;
    if (number < 2) return true;

    for (int i = 0; ; i++)
    {
        long square = squareTable[i];
        if (square > number) return false;
        while (number % square == 0)
        {
            number /= square;
        }
        if (number == 1) return true;
    }
}

Myślę, że to trochę tajemnicze. Na każdym kroku sprawdza, czy kwadrat liczby pierwszej dzieli liczbę wejściową. Jeśli tak, dzieli liczbę przez kwadrat tak długo, jak to możliwe, aby usunąć ten kwadrat z głównego rozkładu. Jeśli w wyniku tego procesu dojdziemy do 1, to liczbą wejściową był rozkład kwadratu liczb pierwszych. Jeśli kwadrat staje się większy niż sama liczba, to nie ma możliwości, aby ten kwadrat lub jakiekolwiek większe kwadraty mogły go podzielić, więc liczba nie może być rozkładem kwadratów liczb pierwszych.

Biorąc pod uwagę dzisiejszy sqrt wykonywany sprzętowo i potrzebę obliczenia liczb pierwszych tutaj, myślę, że to rozwiązanie jest znacznie wolniejsze. Ale powinno to dać lepsze wyniki niż rozwiązanie z sqrt, które nie będzie działać powyżej 2 ^ 54, jak mówi mrzl w swojej odpowiedzi.

Cyrille Ka
źródło
1
podział liczb całkowitych jest wolniejszy niż FP sqrt na obecnym sprzęcie. Ten pomysł nie ma szans. >. <Nawet w 2008 r. sqrtsdPrzepustowość Core2 wynosi 1 na 6-58c. Jest idivto jeden na 12-36 motocykli. (opóźnienia podobne do przepływności: żadna jednostka nie jest potokowana).
Peter Cordes,
sqrt nie musi być idealnie dokładny. Dlatego sprawdzasz, podnosząc do kwadratu wynik i porównując liczbę całkowitą, aby zdecydować, czy wejściowa liczba całkowita ma dokładną liczbę całkowitą sqrt.
Peter Cordes,
11

Wskazano, że ostatnie dcyfry idealnego kwadratu mogą przyjmować tylko określone wartości. Ostatnie dcyfry (w bazie b) liczby nsą takie same jak reszta, gdy njest podzielona przez bd, tj. w notacji C n % pow(b, d).

Można to uogólnić na dowolny moduł m, tj. n % mmożna użyć, aby wykluczyć pewien procent liczb z idealnych kwadratów. Moduł, którego obecnie używasz, wynosi 64, co pozwala na 12, tj. 19% pozostałych, jak to możliwe kwadratów. Przy odrobinie kodowania znalazłem moduł 110880, który pozwala tylko 2016, tj. 1,8% pozostałych jako możliwe kwadraty. Tak więc w zależności od kosztu operacji modułu (tj. Podziału) i wyszukiwania tabeli względem pierwiastka kwadratowego na komputerze, użycie tego modułu może być szybsze.

Nawiasem mówiąc, jeśli Java ma sposób na przechowywanie spakowanej tablicy bitów dla tabeli odnośników, nie używaj jej. 110880 32-bitowe słowa w dzisiejszych czasach to niewiele pamięci RAM, a pobranie słowa maszynowego będzie szybsze niż pobranie jednego bitu.

Hugh Allen
źródło
Miły. Czy pracowałeś nad tym algebraicznie, czy metodą prób i błędów? Rozumiem, dlaczego jest tak skuteczny - wiele kolizji między idealnymi kwadratami, np. 333 ^ 2% 110880 == 3 ^ 2, 334 ^ 2% 110880 == 26 ^ 2, 338 ^ 2% 110880 == 58 ^ 2 .. ,
finnw
IIRC była to brutalna siła, ale zauważ, że 110880 = 2 ^ 5 * 3 ^ 2 * 5 * 7 * 11, co daje 6 * 3 * 2 * 2 * 2 - 1 = 143 właściwych dzielników.
Hugh Allen,
Odkryłem, że z powodu ograniczeń wyszukiwania, 44352 działa lepiej, z przepustowością 2,6%. Przynajmniej w mojej realizacji.
Fractaly
1
Dzielenie liczb całkowitych ( idiv) jest równe lub gorsze kosztem FP sqrt ( sqrtsd) na obecnym sprzęcie x86. Ponadto całkowicie nie zgadzaj się z unikaniem pól bitowych. Współczynnik trafień w pamięci podręcznej będzie o wiele lepszy w przypadku pola bitowego, a testowanie nieco w polu bitowym to tylko jedna lub dwie prostsze instrukcje niż testowanie całego bajtu. (W przypadku małych tabel, które mieszczą się w pamięci podręcznej nawet jako pola niebędące bitami, najlepsza byłaby tablica bajtów, a nie 32-bitowe int. X86 ma dostęp do jednego bajtu z prędkością równą 32-bitowemu dwordowi.)
Peter Cordes
11

Problem liczb całkowitych zasługuje na rozwiązanie liczb całkowitych. A zatem

Wykonaj wyszukiwanie binarne na (nieujemnych) liczbach całkowitych, aby znaleźć największą liczbę całkowitą taką, że t**2 <= n. Następnie sprawdź, czy r**2 = ndokładnie. To zajmuje czas O (log n).

Jeśli nie wiesz, jak przeszukiwać binarnie liczby całkowite dodatnie, ponieważ zestaw jest nieograniczony, jest to łatwe. Zaczynasz od obliczenia rosnącej funkcji f (powyżej f(t) = t**2 - n) na potęgach dwóch. Kiedy zobaczysz, że zmienia się na dodatnią, znalazłeś górną granicę. Następnie możesz wykonać standardowe wyszukiwanie binarne.

Pułkownik Panic
źródło
W rzeczywistości czas byłby co najmniej O((log n)^2)dlatego, że mnożenie nie jest czasem stałym, ale w rzeczywistości ma dolną granicę O(log n), co staje się widoczne podczas pracy z dużymi liczbami o wielu dokładnościach. Ale zasięg tej wiki wydaje się być 64-bitowy, więc może jest to nbd.
10

Następujące uproszczenie rozwiązania maaartinus wydaje się zmniejszać o kilka punktów procentowych czas działania, ale nie jestem wystarczająco dobry w testowaniu porównawczym, aby stworzyć benchmark, któremu mogę zaufać:

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

Warto sprawdzić, jak pominąć pierwszy test,

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

wpłynie na wydajność.

dfeuer
źródło
2
Wyniki są tutaj . Usunięcie pierwszego testu jest złe, ponieważ rozwiązuje większość przypadków dość tanio. Źródło jest w mojej odpowiedzi (zaktualizowane).
maaartinus
9

Aby uzyskać wydajność, bardzo często trzeba iść na kompromisy. Inni wyrażali różne metody, jednak zauważyłeś, że hack Carmacka był szybszy do pewnych wartości N. Następnie powinieneś sprawdzić „n”, a jeśli jest on mniejszy niż liczba N, użyj hacka Carmacka, w przeciwnym razie użyj innej opisanej metody w odpowiedziach tutaj.

BobbyShaftoe
źródło
Włączyłem również twoją sugestię do rozwiązania. Również ładny uchwyt. :)
Kip
8

Jest to najszybsza implementacja Java, jaką mogłem wymyślić, wykorzystując kombinację technik sugerowanych przez innych w tym wątku.

  • Test Mod-256
  • Niedokładny test mod-3465 (pozwala uniknąć dzielenia liczb całkowitych kosztem niektórych fałszywych trafień)
  • Pierwiastek kwadratowy zmiennoprzecinkowy, zaokrąglić i porównać z wartością wejściową

Eksperymentowałem również z tymi modyfikacjami, ale nie pomogły one w wydajności:

  • Dodatkowy test mod-255
  • Dzielenie wartości wejściowej przez potęgi 4
  • Szybki odwrotny pierwiastek kwadratowy (aby pracować dla wysokich wartości N, potrzebuje 3 iteracji, co wystarcza, aby był wolniejszy niż sprzętowa funkcja pierwiastka kwadratowego).

public class SquareTester {

    public static boolean isPerfectSquare(long n) {
        if (n < 0) {
            return false;
        } else {
            switch ((byte) n) {
            case -128: case -127: case -124: case -119: case -112:
            case -111: case -103: case  -95: case  -92: case  -87:
            case  -79: case  -71: case  -64: case  -63: case  -60:
            case  -55: case  -47: case  -39: case  -31: case  -28:
            case  -23: case  -15: case   -7: case    0: case    1:
            case    4: case    9: case   16: case   17: case   25:
            case   33: case   36: case   41: case   49: case   57:
            case   64: case   65: case   68: case   73: case   81:
            case   89: case   97: case  100: case  105: case  113:
            case  121:
                long i = (n * INV3465) >>> 52;
                if (! good3465[(int) i]) {
                    return false;
                } else {
                    long r = round(Math.sqrt(n));
                    return r*r == n; 
                }
            default:
                return false;
            }
        }
    }

    private static int round(double x) {
        return (int) Double.doubleToRawLongBits(x + (double) (1L << 52));
    }

    /** 3465<sup>-1</sup> modulo 2<sup>64</sup> */
    private static final long INV3465 = 0x8ffed161732e78b9L;

    private static final boolean[] good3465 =
        new boolean[0x1000];

    static {
        for (int r = 0; r < 3465; ++ r) {
            int i = (int) ((r * r * INV3465) >>> 52);
            good3465[i] = good3465[i+1] = true;
        }
    }

}
finnw
źródło
7

Powinieneś pozbyć się 2-częściowej mocy N od samego początku.

2. edycja Magiczne wyrażenie dla m poniżej powinno być

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

a nie jak napisano

Koniec 2. edycji

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

1. edycja:

Niewielka poprawa:

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

Koniec 1. edycji

Teraz kontynuuj jak zwykle. W ten sposób, zanim dotrzesz do części zmiennoprzecinkowej, pozbyłeś się już wszystkich liczb, których część 2-mocowa jest nieparzysta (około połowa), a następnie rozważasz tylko 1/8 pozostałej części. Czyli uruchamiasz część zmiennoprzecinkową na 6% liczb.

David Lehavi
źródło
7

Project Euler jest wymieniony w tagach i wiele problemów w nim wymaga sprawdzania liczb >> 2^64 . Większość wyżej opisanych optymalizacji nie działa łatwo, gdy pracujesz z 80-bajtowym buforem.

Użyłem java BigInteger i nieco zmodyfikowanej wersji metody Newtona, która działa lepiej z liczbami całkowitymi. Problem polegał na tym, że dokładne kwadraty były n^2zbieżne (n-1)zamiast nponieważn^2-1 = (n-1)(n+1) a błąd końcowy był tylko o krok poniżej dzielnika końcowego i algorytm zakończył się. Łatwo było to naprawić, dodając jeden do oryginalnego argumentu przed obliczeniem błędu. (Dodaj dwa dla pierwiastków kostki itp.)

Jedną z fajnych cech tego algorytmu jest to, że można natychmiast stwierdzić, czy liczba jest idealnym kwadratem - błąd końcowy (nie korekta) w metodzie Newtona wyniesie zero. Prosta modyfikacja pozwala również szybko obliczyć floor(sqrt(x))zamiast najbliższej liczby całkowitej. Jest to przydatne w przypadku kilku problemów Eulera.

bgiles
źródło
1
Myślałem o tym samym, że te algorytmy nie tłumaczą się dobrze na bufory wieloprecyzyjne. Więc pomyślałem, że przykleję to tutaj ... Znalazłem probabilistyczny test kwadratury o lepszej asymptotycznej złożoności dla wielkich liczb ..... gdzie zastosowania teorii liczb często się . Choć nie zna Project Euler ... wygląda interesująco.
6

To przeróbka z dziesiętnej na dwójkową starego algorytmu kalkulatora Marchanta (przepraszam, nie mam referencji) w Rubim, dostosowanego specjalnie do tego pytania:

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

Oto kilka podobnych rzeczy (proszę nie głosować na mnie za kodowanie stylu / zapachów lub niezdarnego O / O - liczy się algorytm, a C ++ nie jest moim językiem ojczystym). W tym przypadku szukamy pozostałości == 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;
};
Brent.Longborough
źródło
Liczba iteracji wygląda na O (ln n), gdzie n jest długością bitu v, więc wątpię, że zaoszczędzi to dużo dla większych v. Zmienna zmiennoprzecinkowa sqrt jest powolna, może 100-200 cykli, ale matematyka na liczbach całkowitych nie jest albo za darmo. Kilkanaście iteracji po 15 cykli każdy i byłoby to pranie. Mimo to +1 za bycie interesującym.
Tadmas
Właściwie wierzę, że dodawanie i odejmowanie może być wykonane przez XOR.
Brent.Longborough
To był głupi komentarz - tylko XOR może dokonać dodania; odejmowanie jest arytmetyczne.
Brent.Longborough
1
Czy tak naprawdę jest jakakolwiek istotna różnica między czasem działania XOR a dodawaniem?
Tadmas
1
@Tadmas: prawdopodobnie nie wystarczy, aby złamać zasadę „optymalizuj później”. (:-)
Brent.Longborough
6

Jak już wspomniano, wywołanie sqrt nie jest idealnie dokładne, ale jest interesujące i pouczające, że nie rozwala innych odpowiedzi pod względem szybkości. W końcu sekwencja instrukcji języka asemblera dla sqrt jest niewielka. Intel ma instrukcję sprzętową, która, jak sądzę, nie jest używana przez Javę, ponieważ nie jest zgodna z IEEE.

Dlaczego więc jest wolny? Ponieważ Java w rzeczywistości wywołuje procedurę C za pośrednictwem JNI, i jest to wolniejsze niż wywoływanie podprogramu Java, co samo w sobie jest wolniejsze niż wykonywanie go bezpośrednio. Jest to bardzo denerwujące i Java powinna była wymyślić lepsze rozwiązanie, tj. Wbudować w zmiennoprzecinkowe wywołania biblioteki, jeśli to konieczne. No cóż.

Podejrzewam, że w C ++ wszystkie złożone alternatywy straciłyby na szybkości, ale nie sprawdziłem ich wszystkich. To, co zrobiłem i co ludzie Javy będą przydatni, to prosty hack, rozszerzenie specjalnych testów przypadków sugerowanych przez A. Rexa. Użyj pojedynczej długiej wartości jako tablicy bitowej, która nie jest zaznaczona granicami. W ten sposób masz 64-bitowe wyszukiwanie boolowskie.

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

Procedura toPerfectSquare5 działa w około 1/3 czasu na mojej maszynie Core2 Duo. Podejrzewam, że dalsze poprawki w tej samej linii mogą średnio skrócić czas, ale za każdym razem, gdy sprawdzasz, wymieniasz więcej testów na więcej eliminacji, więc nie możesz iść zbyt daleko na tej drodze.

Z pewnością zamiast osobnego testu na wynik ujemny, możesz sprawdzić wysokie 6 bitów w ten sam sposób.

Zauważ, że wszystko, co robię, to eliminowanie możliwych kwadratów, ale gdy mam potencjalny przypadek, muszę zadzwonić do oryginału, wstawiony jestPerfectSquare.

Procedura init2 jest wywoływana raz, aby zainicjować wartości statyczne pp1 i pp2. Zauważ, że w mojej implementacji w C ++ używam niepodpisanego długiego, więc ponieważ jesteś zalogowany, będziesz musiał użyć operatora >>>.

Nie ma wewnętrznej potrzeby sprawdzania tablicy, ale optymalizator Javy musi szybko to rozgryźć, więc nie winię ich za to.

hydrodog
źródło
3
Założę się, że dwa razy się mylisz. 1. Intel sqrt jest zgodny z IEEE. Jedynymi niezgodnymi instrukcjami są instrukcje goniometryczne dla argumentów lange. 2. Java używa funkcji wewnętrznych dla Math.sqrt, bez JNI .
maaartinus
1
Nie zapomniałeś użyć pp2? Rozumiem, że pp1służy do testowania sześciu najmniej znaczących bitów, ale nie sądzę, aby testowanie kolejnych sześciu bitów miało jakikolwiek sens.
maaartinus
6

Podoba mi się pomysł użycia prawie poprawnej metody na niektórych danych wejściowych. Oto wersja z wyższym „przesunięciem”. Kod wydaje się działać i przekazuje mój prosty przypadek testowy.

Wystarczy wymienić:

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

kod z tym:

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);
}
Jonny Heggheim
źródło
6

Biorąc pod uwagę ogólną długość bitów (chociaż użyłem tutaj określonego typu), próbowałem zaprojektować uproszczone algo, jak poniżej. Początkowo wymagana jest prosta i oczywista kontrola dla wartości 0,1,2 lub <0. Poniższe jest proste w tym sensie, że nie próbuje używać żadnych istniejących funkcji matematycznych. Większość operatorów można zastąpić operatorami bitowymi. Nie testowałem jednak żadnych danych porównawczych. Nie jestem ekspertem w matematyce ani w szczególności w projektowaniu algorytmów komputerowych. Chciałbym zobaczyć, jak zwracasz uwagę na problem. Wiem, że istnieje wiele szans na poprawę.

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;  
}
nabam serbang
źródło
@Kip: Problem z przeglądarką.
nabam serbang
1
Potrzebujesz wcięcia.
Steve Kuo,
5

Sprawdziłem wszystkie możliwe wyniki, gdy zaobserwowano ostatnie n bitów kwadratu. Poprzez sukcesywne badanie większej liczby bitów można wyeliminować do 5/6 wejść. Właściwie zaprojektowałem to, aby zaimplementować algorytm faktoryzacji Fermata i jest tam bardzo szybki.

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

Ostatni bit pseudokodu można wykorzystać do rozszerzenia testów w celu wyeliminowania większej liczby wartości. Powyższe testy dotyczą k = 0, 1, 2, 3

  • a ma postać (3 << 2k) - 1
  • b ma postać (2 << 2k)
  • c ma postać (2 << 2k + 2) - 1
  • d ma postać (2 << 2k - 1) * 10

    Najpierw sprawdza, czy ma kwadratową resztę z modułami potęgi dwóch, następnie testuje na podstawie końcowego modułu, a następnie używa Math.sqrt, aby wykonać test końcowy. Wpadłem na pomysł z pierwszego postu i próbowałem go rozwinąć. Doceniam wszelkie komentarze lub sugestie.

    Aktualizacja: Używając testu według modułu (modSq) i podstawy modułu 44352, mój test działa w 96% czasu w porównaniu z aktualizacją OP dla liczb do 1 000 000 000.

    • Fractaly
    źródło
    2

    Oto rozwiązanie dziel i zwyciężaj.

    Jeśli pierwiastek kwadratowy z liczby naturalnej ( number) jest liczbą naturalną ( solution), możesz łatwo określić zakres na solutionpodstawie liczby cyfr number:

    • numberma 1 cyfrę: solutionw zakresie = 1 - 4
    • numberma 2 cyfry: solutionw zakresie = 3 - 10
    • numberma 3 cyfry: solutionw zakresie = 10 - 40
    • numberma 4 cyfry: solutionw zakresie = 30 - 100
    • numberma 5 cyfr: solutionw zakresie = 100 - 400

    Zwróć uwagę na powtórzenie?

    Możesz użyć tego zakresu w podejściu do wyszukiwania binarnego, aby sprawdzić, czy istnieje solution:

    number == solution * solution

    Oto kod

    Oto moja klasa SquareRootChecker

    public class SquareRootChecker {

    private long number;
    private long initialLow;
    private long initialHigh;

    public SquareRootChecker(long number) {
        this.number = number;

        initialLow = 1;
        initialHigh = 4;
        if (Long.toString(number).length() % 2 == 0) {
            initialLow = 3;
            initialHigh = 10;
        }
        for (long i = 0; i < Long.toString(number).length() / 2; i++) {
            initialLow *= 10;
            initialHigh *= 10;
        }
        if (Long.toString(number).length() % 2 == 0) {
            initialLow /= 10;
            initialHigh /=10;
        }
    }

    public boolean checkSquareRoot() {
        return findSquareRoot(initialLow, initialHigh, number);
    }

    private boolean findSquareRoot(long low, long high, long number) {
        long check = low + (high - low) / 2;
        if (high >= low) {
            if (number == check * check) {
                return true;
            }
            else if (number < check * check) {
                high = check - 1;
                return findSquareRoot(low, high, number);
            }
            else  {
                low = check + 1;
                return findSquareRoot(low, high, number);
            }
        }
        return false;
    }

}

    A oto przykład, jak go używać.

    long number =  1234567;
long square = number * number;
SquareRootChecker squareRootChecker = new SquareRootChecker(square);
System.out.println(square + ": " + squareRootChecker.checkSquareRoot()); //Prints "1524155677489: true"

long notSquare = square + 1;
squareRootChecker = new SquareRootChecker(notSquare);
System.out.println(notSquare + ": " + squareRootChecker.checkSquareRoot()); //Prints "1524155677490: false"
    MWB
    źródło
    2
    Uwielbiam tę koncepcję, ale chciałbym uprzejmie zwrócić uwagę na poważną wadę: liczby są w bazie dwójkowej. Konwersja bazy 2 na bazę 10 za pomocą toStringjest niezwykle kosztowną operacją w porównaniu do operatorów bitowych. Tak więc, aby spełnić cel pytania - wydajność - musisz użyć operatorów bitowych zamiast 10 ciągów podstawowych. Znowu bardzo podoba mi się twoja koncepcja. Niezależnie od tego, twoja implementacja (w obecnej formie) jest zdecydowanie najwolniejsza ze wszystkich możliwych rozwiązań opublikowanych dla tego pytania.
    Jack Giffin
    1

    Jeśli chodzi o szybkość, to dlaczego nie podzielić na partycje najczęściej używanego zestawu danych wejściowych i ich wartości, a następnie wykonać zoptymalizowany algorytm magiczny, który wymyśliłeś w wyjątkowych przypadkach?

    Eliasz
    źródło
    Problem polega na tym, że nie ma „powszechnie używanego zestawu danych wejściowych” - zwykle iteruję listę, więc nie użyję tych samych danych wejściowych dwa razy.
    Kip
    1

    Powinno być możliwe spakowanie „nie może być idealnym kwadratem, jeśli ostatnie X cyfr to N” znacznie wydajniej! Użyję 32-bitowych liczb całkowitych java i wygeneruję wystarczającą ilość danych, aby sprawdzić ostatnie 16 bitów liczby - to 2048 liczb szesnastkowych int.

    ...

    Ok. Albo natknąłem się na teorię liczb, która jest trochę poza mną, albo w moim kodzie jest błąd. W każdym razie oto kod:

    public static void main(String[] args) {
    final int BITS = 16;

    BitSet foo = new BitSet();

    for(int i = 0; i< (1<<BITS); i++) {
        int sq = (i*i);
        sq = sq & ((1<<BITS)-1);
        foo.set(sq);
    }

    System.out.println("int[] mayBeASquare = {");

    for(int i = 0; i< 1<<(BITS-5); i++) {
        int kk = 0;
        for(int j = 0; j<32; j++) {
            if(foo.get((i << 5) | j)) {
                kk |= 1<<j;
            }
        }
        System.out.print("0x" + Integer.toHexString(kk) + ", ");
        if(i%8 == 7) System.out.println();
    }
    System.out.println("};");
}

    a oto wyniki:

    (ed: elided za niską wydajność w prettify.js; zobacz historię wersji, aby zobaczyć.)

    paulmurray
    źródło
    1

    Metoda Newtona z arytmetyką liczb całkowitych

    Jeśli chcesz uniknąć operacji na liczbach całkowitych, możesz skorzystać z poniższej metody. Zasadniczo wykorzystuje Metodę Newtona zmodyfikowaną dla arytmetyki liczb całkowitych.

    /**
 * 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;
}

    Ta implementacja nie może konkurować z używanymi rozwiązaniami Math.sqrt. Jednak jego wydajność można poprawić za pomocą mechanizmów filtrujących opisanych w niektórych innych postach.

    awenturyn
    źródło
    1

    Obliczanie pierwiastków kwadratowych metodą Newtona jest niesamowicie szybkie ... pod warunkiem, że wartość początkowa jest rozsądna. Jednak nie ma rozsądnej wartości początkowej, aw praktyce kończymy zachowaniem dzielenia i logowania (2 ^ 64).
    Aby być naprawdę szybkim, potrzebujemy szybkiego sposobu na uzyskanie rozsądnej wartości początkowej, a to oznacza, że ​​musimy zejść do języka maszynowego. Jeśli procesor dostarcza instrukcję taką jak POPCNT w Pentium, która zlicza zera wiodące, możemy użyć tej wartości, aby uzyskać wartość początkową z połową znaczących bitów. Ostrożnie możemy znaleźć stałą liczbę kroków Newtona, które zawsze będą wystarczające. (W ten sposób zrezygnowano z konieczności wykonywania pętli i wykonywania bardzo szybko).

    Drugim rozwiązaniem jest funkcja zmiennoprzecinkowa, która może mieć szybkie obliczenia sqrt (takie jak koprocesor i87). Nawet wycieczka przez exp () i log () może być szybsza niż Newton zdegenerowana do wyszukiwania binarnego. Jest to trudny aspekt, zależna od procesora analiza tego, co i jeśli konieczne jest późniejsze udoskonalenie.

    Trzecie rozwiązanie rozwiązuje nieco inny problem, ale warto o nim wspomnieć, ponieważ sytuacja jest opisana w pytaniu. Jeśli chcesz obliczyć bardzo wiele pierwiastków kwadratowych dla liczb, które różnią się nieznacznie, możesz użyć iteracji Newtona, jeśli nigdy nie ponownie zainicjujesz wartości początkowej, ale po prostu zostaw ją tam, gdzie poprzednio zostało obliczone. Użyłem tego z powodzeniem w co najmniej jednym problemie Eulera.

    Albert van der Horst
    źródło
    Uzyskanie dobrego oszacowania nie jest zbyt trudne. Możesz użyć liczby cyfr liczby, aby oszacować dolną i górną granicę rozwiązania. Zobacz także moją odpowiedź, w której proponuję rozwiązanie dziel i zwyciężaj.
    MWB
    Jaka jest różnica między POPCNT a liczeniem cyfr? Tyle że możesz wykonać POPCNT w ciągu jednej nanosekundy.
    Albert van der Horst
    1

    Pierwiastek kwadratowy liczby, biorąc pod uwagę, że liczba ta jest kwadratem idealnym.

    Złożoność to log (n)

    /**
 * Calculate square root if the given number is a perfect square.
 * 
 * Approach: Sum of n odd numbers is equals to the square root of n*n, given 
 * that n is a perfect square.
 *
 * @param number
 * @return squareRoot
 */

public static int calculateSquareRoot(int number) {

    int sum=1;
    int count =1;
    int squareRoot=1;
    while(sum<number) {
        count+=2;
        sum+=count;
        squareRoot++;
    }
    return squareRoot;
}
    Sajjad Ali Vayani
    źródło
    0

    Jeśli chcesz prędkości, biorąc pod uwagę, że liczby całkowite mają skończony rozmiar, podejrzewam, że najszybszym sposobem byłoby (a) podzielenie parametrów według rozmiaru (np. Na kategorie według największego zestawu bitów), a następnie sprawdzenie wartości względem tablicy idealnych kwadratów w tym zakresie.

    Niebiański Łasica
    źródło
    2
    W odległości długiej znajdują się 2 ^ 32 idealne kwadraty. Ten stół byłby ogromny. Ponadto przewaga obliczania wartości nad dostępem do pamięci może być ogromna.
    PeterAllenWebb
    O nie, nie ma, są 2 ^ 16. 2 ^ 32 to 2 ^ 16 do kwadratu. Jest 2 ^ 16.
    Celestial M Weasel
    3
    tak, ale zasięg długiego wynosi 64 bity, a nie 32 bity. sqrt (2 ^ 64) = 2 ^ 32. (Ignoruję bit znaku, aby matematyka była trochę łatwiejsza ... w rzeczywistości są (długie) (2 ^ 31,5) = 3037000499 idealne kwadraty)
    Kip
    0

    Jeśli chodzi o metodę Carmaca, wydaje się, że dość łatwo byłoby powtórzyć jeszcze raz, co powinno podwoić liczbę cyfr dokładności. Jest to w końcu niezwykle skrócona metoda iteracyjna - metoda Newtona, z bardzo dobrym pierwszym odgadnięciem.

    Jeśli chodzi o twoje obecne najlepsze, widzę dwie mikrooptymalizacje:

    • przesuń czek vs 0 po czeku za pomocą mod255
    • zmienić kolejność dzielących mocy czterech, aby pominąć wszystkie kontrole w zwykłej (75%) sprawie.

    To znaczy:

    // 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;
}

    Jeszcze lepiej może być proste

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

    Oczywiście byłoby interesujące wiedzieć, ile liczb zostaje wyzerowanych w każdym punkcie kontrolnym - raczej wątpię, czy kontrole są naprawdę niezależne, co sprawia, że ​​wszystko jest trudne.

    Ben
    źródło