Самый быстрый способ определить, является ли целочисленный квадратный корень целым числом

Я ищу самый быстрый способ определить, является ли long значение является идеальным квадратом (то есть его квадратный корень является другим целым числом):

  1. Я сделал это простым способом, используя встроенный Math.sqrt() функции, но мне интересно, есть ли способ сделать это быстрее, ограничив себя только целочисленной областью.
  2. Ведение справочной таблицы нецелесообразно (поскольку существует около 2 31,5 целых чисел, площадь которых меньше 2 63).

Вот очень простой и понятный способ сделать это сейчас:

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

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

Примечание: я использую эту функцию во многих задачах Project Euler. Так что больше никому не придется поддерживать этот код. И этот вид микрооптимизации может реально изменить ситуацию, так как одна из задач состоит в том, чтобы выполнить каждый алгоритм менее чем за минуту, и в некоторых задачах эту функцию нужно будет вызывать миллионы раз.


Я пробовал разные решения проблемы:

  • После исчерпывающего тестирования я обнаружил, что добавление 0.5 к результату Math.sqrt() нет необходимости, по крайней мере, не на моей машине.
  • Быстрый обратный квадратный корень был быстрее, но он дал неверные результаты для n >= 410881. Однако, как предполагает BobbyShaftoe, мы можем использовать хак FISR для n < 410881.
  • Метод Ньютона был немного медленнее, чем Math.sqrt(), Это наверное потому что Math.sqrt() использует что-то похожее на метод Ньютона, но реализовано в оборудовании, так что это намного быстрее, чем в Java. Кроме того, метод Ньютона все еще требовал использования двойных чисел.
  • Модифицированный метод Ньютона, который использовал несколько приемов так, чтобы была задействована только целочисленная математика, требовал некоторых хаков, чтобы избежать переполнения (я хочу, чтобы эта функция работала со всеми положительными 64-битными целыми числами со знаком), и это было все еще медленнее, чем Math.sqrt(),
  • Бинарная отбивная была еще медленнее. Это имеет смысл, потому что двоичной отбивке в среднем потребуется 16 проходов, чтобы найти квадратный корень 64-битного числа.
  • Согласно тестам Джона, используя or операторы быстрее в C++, чем с помощью switch, но в Java и C#, похоже, нет никакой разницы между or а также switch,
  • Я также попытался создать таблицу поиска (как частный статический массив из 64 логических значений). Тогда вместо того, чтобы или выключить или or Скажите, я бы просто сказал if(lookup[(int)(n&0x3F)]) { test } else return false;, К моему удивлению, это было (немного) медленнее. Это потому, что границы массивов проверяются в Java.

38 ответов

Решение

Я нашел метод, который работает на 35% быстрее, чем ваш код 6bit +Carmack+sqrt, по крайней мере, с моим процессором (x86) и языком программирования (C/C++). Ваши результаты могут отличаться, особенно потому, что я не знаю, как будет действовать фактор Java.

Мой подход тройной:

  1. Сначала отфильтруйте очевидные ответы. Это включает в себя отрицательные числа и глядя на последние 4 бита. (Я обнаружил, что просмотр последних шести не помог.) Я также отвечаю да на 0. (Читая код ниже, обратите внимание, что мой ввод int64 x.)
    if( x < 0 || (x&2) || ((x & 7) == 5) || ((x & 11) == 8) )
        return false;
    if( x == 0 )
        return true;
  2. Затем, проверьте, является ли это квадрат по модулю 255 = 3 * 5 * 17. Поскольку это произведение трех различных простых чисел, только около 1/8 из остатков по модулю 255 являются квадратами. Однако, по моему опыту, вызов оператора по модулю (%) стоит больше выгоды, которую можно получить, поэтому я использую битовые трюки с 255 = 2^8-1 для вычисления остатка. (Что бы там ни было, я не использую уловку чтения отдельных байтов из слова, только поразрядно - и сдвиги.)
    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.
    
    Чтобы на самом деле проверить, является ли остаток квадратом, я ищу ответ в предварительно вычисленной таблице.
    if( bad255[y] )
        return false;
    // However, I just use a table of size 512
    
  3. Наконец, попробуйте вычислить квадратный корень, используя метод, аналогичный лемме Хензеля. (Я не думаю, что это применимо напрямую, но работает с некоторыми изменениями.) Перед этим я делю все степени 2 с помощью двоичного поиска:
    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;
    На данный момент, чтобы наше число было квадратным, оно должно быть 1 mod 8.
    if((x & 7) != 1)
        return false;
    Основная структура леммы Гензеля заключается в следующем. (Примечание: непроверенный код; если он не работает, попробуйте t=2 или 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.
    Идея состоит в том, что на каждой итерации вы добавляете один бит в r, "текущий" квадратный корень из x; каждый квадратный корень является точным по модулю все большей и большей степени 2, а именно t/2. В конце r и t/2-r будут квадратными корнями из x по модулю t/2. (Обратите внимание, что если r - это квадратный корень из x, то так же и -r. Это верно даже по модулю чисел, но будьте осторожны, по модулю некоторых чисел вещи могут иметь даже более 2 квадратных корней; в частности, это включает степени 2.) Поскольку наш фактический квадратный корень меньше 2^32, в этот момент мы можем просто проверить, являются ли r или t/2-r действительными квадратными корнями. В моем реальном коде я использую следующий измененный цикл:
    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) );
    Ускорение здесь достигается тремя способами: предварительно вычисленное начальное значение (эквивалентное ~10 итерациям цикла), более ранний выход из цикла и пропуск некоторых значений t. В последней части я смотрю на z = r - x * x и установите t, чтобы быть наибольшей степенью 2, деля z с небольшим уловкой. Это позволяет мне пропустить t значений, которые не повлияли бы на значение r в любом случае. Предварительно вычисленное начальное значение в моем случае выбирает "наименьший положительный" квадратный корень по модулю 8192.

Даже если этот код не работает для вас быстрее, я надеюсь, вам понравятся некоторые идеи, которые он содержит. Далее следует полный проверенный код, включая предварительно вычисленные таблицы.

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

Я довольно поздно на вечеринку, но я надеюсь дать лучший ответ; короче и (при условии, что мой тест верен) также намного быстрее.

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

Первый тест ловит большинство не квадратов быстро. Он использует таблицу из 64 элементов, упакованную в long, поэтому нет затрат на доступ к массиву (проверка косвенности и границ). Для равномерно случайного longесть вероятность окончания здесь 81,25%.

Второй тест ловит все числа, имеющие нечетное число двойок в их факторизации. Метод Long.numberOfTrailingZeros очень быстро, так как он получает JIT-ed в одну инструкцию i86.

После отбрасывания конечных нулей третий тест обрабатывает числа, заканчивающиеся на 011, 101 или 111 в двоичном виде, которые не являются идеальными квадратами. Он также заботится об отрицательных числах, а также обрабатывает 0.

Финальный тест возвращается к double арифметика. Как double имеет только 53 бит мантиссы, преобразование из long в double включает округление для больших значений. Тем не менее, тест верен (если доказательство не верно).

Попытка внедрить идею mod255 не удалась.

Вам нужно будет сделать несколько тестов. Лучший алгоритм будет зависеть от распределения ваших входных данных.

Ваш алгоритм может быть почти оптимальным, но вы можете сделать быструю проверку, чтобы исключить некоторые возможности, прежде чем вызывать подпрограмму с квадратным корнем. Например, посмотрите на последнюю цифру вашего числа в шестнадцатеричном формате, выполнив побитовое "и". Совершенные квадраты могут заканчиваться только 0, 1, 4 или 9 в основании 16, так что для 75% ваших входных данных (при условии, что они распределены равномерно) вы можете избежать вызова квадратного корня в обмен на какое-то очень быстрое переключение битов.

Кип протестировал следующий код, реализующий шестнадцатеричный трюк. При тестировании чисел от 1 до 100 000 000 этот код выполнялся в два раза быстрее оригинала.

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

Когда я тестировал аналогичный код в C++, он на самом деле работал медленнее, чем оригинал. Однако, когда я исключил оператор switch, шестнадцатеричный трюк снова сделал код в два раза быстрее.

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

Исключение оператора switch мало повлияло на код C#.

Я думал об ужасных временах, которые я провел в курсе численного анализа.

И потом я помню, что эта функция кружила по сети из исходного кода 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;
}

Который в основном вычисляет квадратный корень, используя функцию приближения Ньютона (не могу вспомнить точное имя).

Это должно быть удобно и даже быстрее, это из одной из феноменальных игр id!

Он написан на C++, но не должно быть слишком сложно повторно использовать ту же технику в Java, как только вы получите идею:

Первоначально я нашел его по адресу: http://www.codemaestro.com/reviews/9

Метод Ньютона объяснен в Википедии: http://en.wikipedia.org/wiki/Newton%27s_method

Вы можете перейти по ссылке для более подробного объяснения того, как это работает, но если вам все равно, то это примерно то, что я помню из чтения блога и прохождения курса численного анализа:

  • * (long*) &y в основном это функция быстрого преобразования в long, поэтому целые операции могут применяться к необработанным байтам.
  • 0x5f3759df - (i >> 1); линия - это предварительно рассчитанное начальное значение для функции аппроксимации.
  • * (float*) &i преобразует значение обратно в число с плавающей запятой
  • y = y * ( threehalfs - ( x2 * y * y ) ) line снова выполняет итерацию значения по функции.

Функция приближения дает более точные значения, чем больше вы повторяете функцию по результату. В случае с Quake, одна итерация "достаточно хороша", но если бы она была не для вас... тогда вы могли бы добавить столько итераций, сколько вам нужно.

Это должно быть быстрее, потому что это уменьшает количество операций деления, выполняемых в простом квадратном корне, до простого деления на 2 (на самом деле * 0.5F операция умножения) и замените ее несколькими фиксированными числами операций умножения.

Я не уверен, будет ли это быстрее или даже точнее, но вы можете использовать алгоритм магического квадратного корня Джона Кармака, чтобы быстрее решить квадратный корень. Вероятно, вы могли бы легко проверить это для всех возможных 32-битных целых чисел и убедиться, что вы действительно получили правильные результаты, так как это всего лишь приближение. Однако, теперь, когда я думаю об этом, использование двойных чисел также приближенно, так что я не уверен, каким образом это вступит в игру.

Если вы выполните двоичную отбивку, чтобы попытаться найти "правильный" квадратный корень, вы можете довольно легко определить, достаточно ли близко полученное значение, чтобы сказать:

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

Итак, рассчитав n^2, варианты:

  • n^2 = target: сделано, верни истину
  • n^2 + 2n + 1 > target > n^2: вы близки, но не идеальны: верните ложь
  • n^2 - 2n + 1 < target < n^2: то же самое
  • target < n^2 - 2n + 1: бинарная отбивная по нижнему n
  • target > n^2 + 2n + 1: бинарная отбивная на высшем n

(Извините, это использует n как ваше текущее предположение, и target для параметра. Извиняюсь за путаницу!)

Я не знаю, будет ли это быстрее или нет, но стоит попробовать.

РЕДАКТИРОВАТЬ: бинарная отбивная не должна принимать весь диапазон целых чисел, либо (2^x)^2 = 2^(2x)так что, как только вы найдете верхний установленный бит в вашей цели (что может быть сделано с помощью хитрого трюка; я точно забыл, как), вы можете быстро получить диапазон возможных ответов. Имейте в виду, что наивный бинарная отбивная все еще займет всего 31 или 32 итерации.

Я провел собственный анализ нескольких алгоритмов в этой теме и получил некоторые новые результаты. Вы можете увидеть эти старые результаты в истории редактирования этого ответа, но они не точные, так как я допустил ошибку и потратил время на анализ нескольких алгоритмов, которые не являются близкими. Однако, извлекая уроки из нескольких разных ответов, у меня теперь есть два алгоритма, которые сокрушают "победителя" этой темы. Вот основная вещь, которую я делаю иначе, чем все остальные:

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

Однако эта простая строка, которая в большинстве случаев добавляет одну или две очень быстрые инструкции, значительно упрощает switch-case утверждение в одно заявление if. Тем не менее, это может добавить к времени выполнения, если многие из протестированных чисел имеют значительную степень двух факторов.

Алгоритмы ниже следующие:

  • Интернет - опубликованный ответ Кипа
  • Durron - мой модифицированный ответ, использующий однопроходный ответ в качестве основы
  • DurronTwo - мой измененный ответ с использованием двухпроходного ответа (@JohnnyHeggheim) с некоторыми другими незначительными изменениями.

Вот пример времени выполнения, если числа генерируются с использованием 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

А вот пример времени выполнения, если он запускается только для первого миллиона длинных:

 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

Как вы видете, DurronTwo лучше подходит для больших входов, потому что он очень часто использует магический трюк, но затупляется по сравнению с первым алгоритмом и Math.sqrt потому что цифры намного меньше. Между тем, чем проще Durron является огромным победителем, потому что он никогда не должен делиться на 4 много много раз в первом миллионном числе.

Вот 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;
}

А также 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;
}

И мой тестовый жгут: (Требуется 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);
    }
}

ОБНОВЛЕНИЕ: я сделал новый алгоритм, который быстрее в некоторых сценариях, медленнее в других, я получил разные тесты, основанные на разных входах. Если мы посчитаем по модулю 0xFFFFFF = 3 x 3 x 5 x 7 x 13 x 17 x 241, мы можем исключить 97,82% чисел, которые не могут быть квадратами. Это может быть (вроде) сделано в одной строке, с 5 побитовыми операциями:

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

Результирующий индекс либо 1) остаток, 2) остаток + 0xFFFFFFили 3) остаток + 0x1FFFFFE, Конечно, нам нужно иметь таблицу поиска остатков по модулю 0xFFFFFF, что составляет около 3 МБ файла (в этом случае сохраняются как десятичные числа в тексте ascii, не оптимальные, но явно улучшаемые с ByteBuffer и так далее. Но так как это предварительный расчет, это не имеет большого значения. Вы можете найти файл здесь (или создать его самостоятельно):

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

Я загружаю это в boolean массив, как это:

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

Пример выполнения. Это бить Durron (версия первая) в каждом испытании я бежал.

 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

Должно быть намного быстрее использовать метод Ньютона для вычисления корня целочисленного квадрата, затем возвести в квадрат это число и проверить, как вы делаете в своем текущем решении. Метод Ньютона является основой для решения Кармака, упомянутого в некоторых других ответах. Вы должны быть в состоянии получить более быстрый ответ, так как вас интересует только целочисленная часть корня, что позволяет быстрее остановить алгоритм аппроксимации.

Еще одна оптимизация, которую вы можете попробовать: если цифровой корень числа не заканчивается на 1, 4, 7 или 9, число не является идеальным квадратом. Это можно использовать как быстрый способ устранить 60% ваших входных данных перед применением более медленного алгоритма квадратного корня.

Я хочу, чтобы эта функция работала со всеми положительными 64-битными целыми числами со знаком

Math.sqrt() работает с двойными значениями в качестве входных параметров, поэтому вы не получите точных результатов для целых чисел больше 2 ^ 53.

Просто для записи, другой подход заключается в использовании простого разложения. Если каждый фактор разложения четный, то число является идеальным квадратом. Итак, вы хотите увидеть, можно ли разложить число как произведение квадратов простых чисел. Конечно, вам не нужно получать такое разложение, просто чтобы увидеть, существует ли оно.

Сначала создайте таблицу квадратов простых чисел, которые меньше, чем 2^32. Это намного меньше, чем таблица всех целых чисел до этого предела.

Решение тогда будет таким:

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

Я думаю, это немного загадочно. На каждом шаге он проверяет, что квадрат простого числа делит входное число. Если это так, то он делит число на квадрат настолько долго, насколько это возможно, чтобы удалить этот квадрат из простого разложения. Если в результате этого процесса мы пришли к 1, то входное число было разложением квадрата простых чисел. Если квадрат становится больше, чем само число, то этот квадрат или более крупные квадраты никак не могут его разделить, поэтому число не может быть разложением квадратов простых чисел.

Учитывая, что в настоящее время sqrt выполняется аппаратно, и здесь необходимо вычислять простые числа, я думаю, что это решение намного медленнее. Но это должно дать лучшие результаты, чем решение с sqrt, которое не будет работать более 2^54, как говорит mrzl в своем ответе.

Целочисленная задача заслуживает целочисленного решения. таким образом

Выполните бинарный поиск по (неотрицательным) целым числам, чтобы найти наибольшее целое число t, такое что t**2 <= n, Затем проверьте, r**2 = n именно так. Это занимает время O(log n).

Если вы не знаете, как выполнить двоичный поиск натуральных чисел, потому что множество не ограничено, это легко. Вы начинаете с вычисления вашей возрастающей функции F (выше f(t) = t**2 - n) на двоих. Когда вы видите, что это становится положительным, вы нашли верхнюю границу. Затем вы можете сделать стандартный бинарный поиск.

Было отмечено, что последний d цифры идеального квадрата могут принимать только определенные значения. Последний d цифры (в базе b) числа n такой же, как остаток, когда n делится на bd т.е. в нотации С n % pow(b, d),

Это может быть обобщено на любой модуль m т.е. n % m может использоваться для исключения некоторого процента чисел из идеальных квадратов. Модуль, который вы используете в настоящее время, равен 64, что позволяет 12, т.е. 19% остатков, как возможные квадраты. С небольшим кодированием я нашел модуль 110880, который позволяет только 2016, т.е. 1,8% остатков в качестве возможных квадратов. Таким образом, в зависимости от стоимости операции модуля (т. Е. Деления) и поиска в таблице по сравнению с квадратным корнем на вашей машине, использование этого модуля может быть быстрее.

Кстати, если у Java есть способ хранить упакованный массив битов для таблицы поиска, не используйте его. В наши дни 110880 32-битных слов - это не много ОЗУ, и загрузка машинного слова будет быстрее, чем загрузка одного бита.

Следующее упрощение решения maaartinus, похоже, позволяет сократить время выполнения на несколько процентных пунктов, но я недостаточно хорош в тестировании, чтобы произвести тест, которому я могу доверять:

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

Стоит проверить, как пропустить первый тест,

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

повлияет на производительность.

Для производительности вам очень часто приходится идти на некоторые компромиссы. Другие выразили различные методы, однако вы заметили, что хак Кармака был быстрее до определенных значений N. Затем вы должны проверить "n", и если оно меньше, чем число N, используйте хак Кармака, иначе используйте какой-то другой описанный метод в ответах здесь.

Это самая быстрая реализация Java, которую я мог придумать, используя комбинацию методов, предложенных другими в этой теме.

  • Мод-256 тест
  • Неточный тест mod-3465 (избегает целочисленного деления за счет некоторых ложных срабатываний)
  • Квадратный корень с плавающей точкой, округлить и сравнить с входным значением

Я также экспериментировал с этими модификациями, но они не помогли производительности:

  • Дополнительный мод-255 тест
  • Деление входного значения на степени 4
  • Быстрый обратный квадратный корень (для работы при больших значениях N требуется 3 итерации, достаточных для того, чтобы сделать это медленнее, чем аппаратная функция квадратного корня.)

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

}

Вы должны избавиться от 2-степенной части N с самого начала.

2nd Edit Волшебное выражение для м ниже должно быть

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

а не как написано

Конец 2-го редактирования

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

1-е редактирование:

Незначительное улучшение:

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

Конец первого редактирования

Теперь продолжайте как обычно. Таким образом, к тому времени, как вы доберетесь до части с плавающей запятой, вы уже избавились от всех чисел, чья 2-степенная часть нечетна (примерно половина), и тогда вы будете считать только 1/8 того, что осталось. Т.е. вы запускаете часть с плавающей запятой на 6% чисел.

Проект Эйлер упоминается в тегах, и многие из проблем в нем требуют проверки номера >> 2^64. Большинство упомянутых выше оптимизаций не работают легко, когда вы работаете с 80-байтовым буфером.

Я использовал java BigInteger и слегка модифицированную версию метода Ньютона, которая лучше работает с целыми числами. Проблема заключалась в том, что точные квадраты n^2 сходились к (n-1) вместо n, потому что n^2-1 = (n-1)(n+1), и окончательная ошибка была всего на один шаг ниже конечного делителя и алгоритм прекращен. Это было легко исправить, добавив один к исходному аргументу перед вычислением ошибки. (Добавьте два для кубических корней и т. Д.)

Одним из приятных атрибутов этого алгоритма является то, что вы можете сразу сказать, является ли число идеальным квадратом - конечная ошибка (не коррекция) в методе Ньютона будет равна нулю. Простая модификация также позволяет вам быстро вычислить floor(sqrt(x)) вместо ближайшего целого числа. Это удобно с несколькими проблемами Эйлера.

Это доработка от десятичного к двоичному алгоритму старого калькулятора Марчанта (извините, у меня нет ссылки) в Ruby, адаптированном специально для этого вопроса:

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

Вот пример чего-то подобного (пожалуйста, не голосуйте за стиль кодирования / запахи или неуклюжий O/O - это алгоритм, который имеет значение, а C++ не мой родной язык). В этом случае мы ищем остаток == 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;
};

Как уже упоминалось, вызов sqrt не совсем точен, но он интересен и поучителен, так как он не отбрасывает другие ответы с точки зрения скорости. В конце концов, последовательность инструкций на ассемблере для sqrt крошечная. У Intel есть аппаратная инструкция, которая не используется Java, я считаю, потому что она не соответствует IEEE.

Так почему же это медленно? Потому что Java на самом деле вызывает подпрограмму C через JNI, и это на самом деле медленнее, чем вызов подпрограммы Java, которая сама по себе медленнее, чем встроенная. Это очень раздражает, и Java должна была придумать лучшее решение, то есть, при необходимости, создание вызовов библиотеки с плавающей запятой. Ну что ж.

Я подозреваю, что в C++ все сложные альтернативы будут терять скорость, но я не проверял их все. То, что я сделал, и что люди Java найдут полезными, - это простой взлом, расширение тестирования специального случая, предложенного А. Рексом. Используйте одно длинное значение в качестве битового массива, который не проверяется по границам. Таким образом, у вас есть 64-битный логический поиск.

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

Подпрограмма isPerfectSquare5 выполняется примерно на 1/3 времени на моей машине core2 duo. Я подозреваю, что дальнейшие изменения в том же направлении могут в среднем еще больше сократить время, но каждый раз, когда вы проверяете, вы тратите больше тестов на большее устранение, поэтому вы не можете идти слишком далеко по этому пути.

Конечно, вместо того, чтобы иметь отдельный тест для отрицательного значения, вы можете проверить старшие 6 битов таким же образом.

Обратите внимание, что все, что я делаю, это устранение возможных квадратов, но когда у меня есть потенциальный случай, я должен вызвать исходный, встроенный isPerfectSquare.

Процедура init2 вызывается один раз для инициализации статических значений pp1 и pp2. Обратите внимание, что в моей реализации на C++ я использую unsigned long long, поэтому, поскольку вы подписаны, вам придется использовать оператор >>>.

Нет необходимости в проверке массива границ, но оптимизатор Java должен довольно быстро разобраться с этим, поэтому я не виню их за это.

Мне нравится идея использовать почти правильный метод для некоторых входных данных. Вот версия с более высоким "смещением". Код, кажется, работает и проходит мой простой тестовый пример.

Просто замените ваш:

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

код с этим:

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

Учитывая общую длину в битах (хотя здесь я использовал конкретный тип), я попытался разработать упрощенный алгоритм, как показано ниже. Первоначально требуется простая и очевидная проверка на 0,1,2 или<0. Следующее просто в том смысле, что оно не пытается использовать любые существующие математические функции. Большая часть операторов может быть заменена побитовыми операторами. Я не проверял ни с какими контрольными данными все же. Я не являюсь экспертом в области математики или компьютерных алгоритмов, в частности, мне бы очень хотелось, чтобы вы указали на проблему. Я знаю, что есть много шансов на улучшение.

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

Я проверил все возможные результаты, когда наблюдаются последние n бит квадрата. Последовательно исследуя больше битов, можно исключить до 5/6 входных данных. Я на самом деле разработал это для реализации алгоритма факторизации Ферма, и он там очень быстрый.

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

Последний бит псевдокода может использоваться для расширения тестов, чтобы исключить больше значений. Вышеприведенные тесты для k = 0, 1, 2, 3

  • a имеет вид (3 << 2k) - 1
  • b имеет вид (2 << 2k)
  • с имеет вид (2 << 2k + 2) - 1
  • d имеет вид (2 << 2k - 1) * 10

    Сначала он проверяет, имеет ли он квадратную невязку с модулями степени два, затем он проверяет на основе окончательного модуля, а затем использует Math.sqrt для выполнения окончательного теста. Я придумал идею из верхнего поста и попытался ее расширить. Я ценю любые комментарии или предложения.

    Обновление. Используя тест по модулю (modSq) и базе модулей 44352, мой тест выполняется в 96% времени по сравнению с тестом в обновлении OP для чисел до 1 000 000 000.

  • Вот решение "разделяй и властвуй".

    Если квадратный корень из натурального числа (number) натуральное число (solution), вы можете легко определить диапазон для solution на основе количества цифр number:

    • number имеет 1 цифру: solution в диапазоне = 1 - 4
    • number имеет 2 цифры: solution в диапазоне = 3 - 10
    • number имеет 3 цифры: solution в диапазоне = 10 - 40
    • number имеет 4 цифры: solution в диапазоне = 30 - 100
    • number имеет 5 цифр: solution в диапазоне = 100 - 400

    Заметили повторение?

    Вы можете использовать этот диапазон в подходе двоичного поиска, чтобы увидеть, есть ли solution для которого:

    number == solution * solution
    

    Вот код

    Вот мой класс 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;
        }
    
    }
    

    И вот пример того, как его использовать.

    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"
    

    Этот вопрос заставил меня задуматься, поэтому я сделал простой код и представляю его здесь, потому что я думаю, что это интересно, актуально, но я не знаю, насколько полезно. Есть простой алгоритм

          a_n+1 = (a_n + x/a_n)/2
    

    для вычисления квадратных корней, но он предназначен для десятичных дробей. Мне было интересно, что произойдет, если я просто закодирую тот же алгоритм, используя целочисленную математику. Сойдётся ли он вообще к правильному ответу? Я не знал, поэтому написал программу...

          #include <stdio.h>
    #include <stdint.h>
    #include <stdlib.h>
    #include <math.h>
    
    _Bool isperfectsquare(uint64_t x, uint64_t *isqrtx) {
      // NOTE: isqrtx approximate for non-squares. (benchmarked at 162ns 3GHz i5)
      uint32_t i;
      uint64_t ai;
      ai = 1 + ((x & 0xffff000000000000) >> 32) + ((x & 0xffff00000000) >> 24) + ((x & 0xffff0000) >> 16);
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = (ai + x/ai)/2;
      ai = ai & 0xffffffff;
      if (isqrtx != NULL) isqrtx[0] = ai;
      return ai*ai == x;
    }
    
    void main() {
    
      uint64_t x, isqrtx;
      uint64_t i;
      for (i=1; i<0x100000000; i++) {
        if (!isperfectsquare(i*i, &isqrtx)) {
          printf("Failed at %li", i);
          exit(1);
        }
      }
      printf("All OK.\n");
    } 
    

    Итак, получается, что 12 итераций формулы достаточно, чтобы дать правильные результаты для всех 64-битных беззнаковых длинных чисел, которые являются идеальными квадратами, и, конечно же, неквадраты вернут false.

          simon@simon-Inspiron-N5040:~$ time ./isqrt.bin 
    All OK.
    
    real    11m37.096s
    user    11m35.053s
    sys 0m0.272s
    

    Таким образом, 697 с/2^32 составляет примерно 162 нс. Как бы то ни было, функция будет иметь одинаковое время выполнения для всех входных данных. Некоторые из мер, подробно описанных в другом месте в обсуждении, могут ускорить его для неквадратов, проверив последние четыре бита и т. д. Надеюсь, что кто-то найдет это интересным, как и я.

    Квадратный корень числа, учитывая, что число представляет собой полный квадрат.

    Сложность 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;
    }
    

    Вот самый простой и краткий способ, хотя я не знаю, как он сравнивается с точки зрения циклов процессора. Это прекрасно работает, если вы хотите знать, является ли корень целым числом. Если вам действительно важно, является ли оно целым числом, вы также можете понять это. Вот простая (и чистая) функция:

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

    Если вам не нужна микрооптимизация, этот ответ лучше с точки зрения простоты и удобства обслуживания. Если вы будете получать отрицательные числа, возможно, вы захотите использовать Math.abs() для аргумента числа в качестве аргумента Math.sqrt().

    На моем 3,6 ГГц процессоре Intel i7-4790 запуск этого алгоритма на 0–10 000 000 занял в среднем 35–37 наносекунд на расчёт. Я выполнил 10 последовательных прогонов, напечатав среднее время, затрачиваемое на каждый из десяти миллионов расчетов. Каждый полный прогон занимал чуть более 600 мсек.

    Если вы выполняете меньшее количество вычислений, более ранние вычисления занимают немного больше времени.

    Метод Ньютона с целочисленной арифметикой

    Если вы хотите избежать нецелых операций, вы можете использовать метод ниже. Он в основном использует метод Ньютона, модифицированный для целочисленной арифметики.

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

    Эта реализация не может конкурировать с решениями, которые используют Math.sqrt, Однако его производительность может быть улучшена с помощью механизмов фильтрации, описанных в некоторых других публикациях.

    Если скорость вызывает беспокойство, почему бы не выделить из наиболее часто используемых наборов входных данных и их значений таблицу поиска, а затем выполнить любой оптимизированный магический алгоритм, который вы придумали для исключительных случаев?

    Возможно, лучшим алгоритмом для этой проблемы является алгоритм быстрого целочисленного квадратного корня /questions/47657250/ischete-effektivnyij-algoritm-tselochislennogo-kvadratnogo-kornya-dlya-arm-thumb2/47657284#47657284

    Там @Kde утверждает, что трех итераций метода Ньютона было бы достаточно для точности ±1 для 32-битных целых чисел. Конечно, для 64-разрядных целых чисел требуется больше итераций, может быть 6 или 7.

    Расчет квадратных корней по методу Ньютона ужасно быстр... при условии, что начальное значение разумно. Однако разумного начального значения нет, и на практике мы заканчиваем разделением на две части и логарифмическим поведением (2^64).
    Чтобы быть действительно быстрым, нам нужен быстрый способ достичь разумного начального значения, а это значит, что нам нужно погрузиться в машинный язык. Если процессор предоставляет инструкцию типа POPCNT в Pentium, которая подсчитывает начальные нули, мы можем использовать ее, чтобы получить начальное значение с половиной значащих бит. С осторожностью мы можем найти фиксированное количество шагов Ньютона, которое всегда будет достаточно. (Таким образом, отпадает необходимость в цикле и очень быстром исполнении.)

    Второе решение заключается в использовании функции с плавающей запятой, которая может иметь быстрое вычисление sqrt (как, например, сопроцессор i87). Даже экскурсия через exp() и log() может быть быстрее, чем Ньютон, вырождающийся в двоичный поиск. В этом есть один сложный аспект, зависящий от процессора анализ того, что и если впоследствии необходимо усовершенствовать.

    Третье решение решает немного другую проблему, но стоит упомянуть, потому что ситуация описана в этом вопросе. Если вы хотите вычислить большое количество квадратных корней для чисел, которые немного отличаются, вы можете использовать итерацию Ньютона, если вы никогда не инициализируете начальное значение, а просто оставляете его там, где остановились предыдущие вычисления. Я использовал это с успехом по крайней мере в одной проблеме Эйлера.

    Другие вопросы по тегам