整数の平方根が整数であるかどうかを判断する最も速い方法

1497
Kip 2008-11-18 03:43.

long値が完全な正方形であるかどうか(つまり、その平方根が別の整数であるかどうか)を判断するための最速の方法を探しています。

  1. 組み込みMath.sqrt()関数を使用して簡単な方法で実行しましたが、整数のみのドメインに制限することでより高速に実行する方法があるかどうか疑問に思っています。
  2. ルックアップテーブルを維持することは実用的ではありません(平方が2 63未満である約231.5の整数があるため)。

これが私が今やっている非常にシンプルで簡単な方法です:

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

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

注:私はこの関数を多くのプロジェクトオイラーの問題で使用しています。したがって、他の誰もこのコードを維持する必要はありません。そして、この種のマイクロ最適化は実際に違いを生む可能性があります。課題の一部はすべてのアルゴリズムを1分未満で実行することであり、問​​題によってはこの関数を何百万回も呼び出す必要があるためです。


私は問題のさまざまな解決策を試しました:

  • 徹底的なテスト0.5の結果、少なくとも私のマシンでは、Math.sqrt()の結果に追加する必要がないことがわかりました。
  • 平方根逆高速は速かったが、それは、n> = 410881.に対して誤った結果を与えたが、によって提案されたようBobbyShaftoe、我々はN <410881のためFISRハックを使用することができます。
  • ニュートン法は、よりも少し遅かったMath.sqrt()。これはおそらく、Math.sqrt()ニュートン法に似たものを使用しているが、ハードウェアに実装されているため、Javaよりもはるかに高速であるためです。また、ニュートン法では依然としてダブルの使用が必要でした。
  • 整数計算のみが含まれるようにいくつかのトリックを使用した修正ニュートン法は、オーバーフローを回避するためにいくつかのハックを必要とし(この関数をすべての正の64ビット符号付き整数で機能させたい)、それでもMath.sqrt()。よりも低速でした。
  • バイナリチョップはさらに遅くなりました。バイナリチョップは64ビット数の平方根を見つけるために平均16パスを必要とするため、これは理にかなっています。
  • ジョンのテストによると、使用してorステートメントは、より高速なC ++で使用するよりもあるswitchが、JavaやC#での間に違いはないようであるorswitch
  • また、ルックアップテーブルを作成してみました(64ブール値のプライベート静的配列として)。次に、スイッチまたはorステートメントの代わりに、私はただ言うでしょうif(lookup[(int)(n&0x3F)]) { test } else return false;。驚いたことに、これは(ほんの少しだけ)遅くなりました。これは、Javaで配列の境界がチェックされるためです。

30 answers

760
A. Rex 2009-01-09 06:32.

少なくとも私のCPU(x86)とプログラミング言語(C / C ++)では、6ビット+ Carmack + sqrtコードよりも約35%高速に動作する方法を見つけました。特にJavaファクターがどのように機能するかわからないため、結果は異なる場合があります。

私のアプローチは3つあります。

  1. まず、明白な答えを除外します。これには、負の数と最後の4ビットの確認が含まれます。(最後の6つを見ても役に立たなかったことがわかりました。)0についてもyesと答えます。(以下のコードを読むときは、入力が次のようになっていることに注意してください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を法とする正方形であるかどうかを確認します。これは3つの異なる素数の積であるため、mod255の剰余の約1/8のみが正方形です。ただし、私の経験では、モジュロ演算子(%)を呼び出すと、得られるメリットよりもコストがかかるため、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 mod8である必要があります。
    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.
    アイデアは、各反復で、xの「現在の」平方根であるrに1ビットを追加するというものです。各平方根は、2の累乗、つまりt / 2を法として正確です。最後に、rとt / 2-rは、t / 2を法とするxの平方根になります。(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値のスキップの3つの方法で得られます。最後の部分では、を見てz = r - x * x、tをビットトリックでzを除算する2の最大の累乗に設定します。これにより、とにかくrの値に影響を与えなかったt値をスキップできます。私の場合、事前に計算された開始値は、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;
}
398
maaartinus 2013-09-09 07:37.

私はパーティーにかなり遅れていますが、もっと良い答えを提供したいと思っています。より短く、(私のベンチマークが正しいと仮定して)はるかに高速です。

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

最初のテストでは、ほとんどの非正方形をすばやくキャッチします。longにパックされた64項目のテーブルを使用するため、配列アクセスのコスト(間接参照と境界チェック)はありません。均一にランダムなlong場合、ここで終了する確率は81.25%です。

2番目のテストは、因数分解で2の奇数を持つすべての数をキャッチします。このメソッドLong.numberOfTrailingZerosは、単一のi86命令にJITされるため、非常に高速です。

末尾のゼロを削除した後、3番目のテストは、完全な平方ではない2進数の011、101、または111で終わる数値を処理します。また、負の数を考慮し、0も処理します。

最終テストはdouble算術にフォールバックします。double唯一の53ビット仮数を有する、からの変換longにはdouble大きな値の丸め含みます。それにもかかわらず、テストは正しいです(証明が間違っていない限り)。

mod255のアイデアを取り入れようとしても成功しませんでした。

133
John D. Cook 2008-11-18 04:27.

ベンチマークを行う必要があります。最適なアルゴリズムは、入力の分布によって異なります。

アルゴリズムはほぼ最適かもしれませんが、平方根ルーチンを呼び出す前に、いくつかの可能性を除外するために簡単なチェックを行うことをお勧めします。たとえば、ビット単位の「and」を実行して、16進数で数字の最後の桁を確認します。完全な正方形は、基数16で0、1、4、または9でしか終了できないため、入力の75%(均一に分散されていると仮定)では、非常に高速なビットの調整と引き換えに平方根の呼び出しを回避できます。

Kipは、16進トリックを実装する次のコードのベンチマークを行いました。1から100,000,000までの数値をテストする場合、このコードは元のコードの2倍の速度で実行されました。

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ステートメントを削除すると、16進トリックによってコードが2倍高速になります。

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#コードにはほとんど影響しませんでした。

53
chakrit 2008-11-18 03:50.

数値解析のコースで過ごした恐ろしい時間を考えていました。

そして、私は覚えています。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ソフトウェアのゲームの1つです。

これはC ++で書かれていますが、アイデアが浮かんだら、Javaで同じ手法を再利用するのはそれほど難しいことではありません。

私はもともとそれを見つけました: http://www.codemaestro.com/reviews/9

ウィキペディアで説明されているニュートン法: http://en.wikipedia.org/wiki/Newton%27s_method

リンクをたどると、その仕組みの詳細を説明できますが、あまり気にしない場合は、ブログを読んだり、数値解析コースを受講したりして、おおまかに覚えています。

  • これ* (long*) &yは基本的に高速のlong変換関数であるため、整数演算をrawバイトに適用できます。
  • この0x5f3759df - (i >> 1);線は、近似関数の事前に計算されたシード値です。
  • * (float*) &i浮動小数点の値の後ろに変換します。
  • このy = y * ( threehalfs - ( x2 * y * y ) )行は、基本的に関数の値を再度繰り返します。

近似関数は、結果に対して関数を反復するほど、より正確な値を提供します。Quakeの場合、1回の反復で「十分」ですが、それが適切でない場合は、必要な数の反復を追加できます。

これは、単純な平方根で行われる除算演算の数を2による単純な除算(実際には* 0.5F乗算演算)に減らし、代わりにいくつかの固定数の乗算演算に置き換えるため、より高速になるはずです。

38
Kibbee 2008-11-18 03:51.

それがより速いか、あるいは正確であるかどうかはわかりませんが、ジョン・カーマックのマジカル平方根アルゴリズムを使用して平方根をより速く解決することができます。考えられるすべての32ビット整数についてこれを簡単にテストし、実際に正しい結果が得られたことを検証できます。これは単なる概算です。ただ、今考えてみると、ダブルスの使用もおおよそのことなので、どうなるかわかりません。

36
Jon Skeet 2008-11-18 03:50.

バイナリチョップを実行して「正しい」平方根を見つけようとすると、取得した値が次のように十分に近いかどうかをかなり簡単に検出できます。

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

したがって、計算するn^2と、オプションは次のようになります。

  • n^2 = target:完了、trueを返します
  • n^2 + 2n + 1 > target > n^2 :あなたは近くにいますが、完璧ではありません:falseを返します
  • 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回の反復しか必要としません。

24
durron597 2013-06-11 06:39.

このスレッドのいくつかのアルゴリズムについて独自の分析を実行し、いくつかの新しい結果を思いつきました。これらの古い結果はこの回答の編集履歴で確認できますが、私が間違えたため正確ではなく、近くにないいくつかのアルゴリズムの分析に時間を浪費しました。しかし、いくつかの異なる答えから教訓を得て、私は今、このスレッドの「勝者」を粉砕する2つのアルゴリズムを持っています。これが私が他の人とは違うやり方でやる中心的なことです:

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

ただし、ほとんどの場合、1つまたは2つの非常に高速な命令を追加するこの単純な行は、switch-caseステートメントを1つのifステートメントに大幅に単純化します。ただし、テストされた数値の多くに2の累乗の重要な要素がある場合は、ランタイムに追加される可能性があります。

以下のアルゴリズムは次のとおりです。

  • インターネット-Kipの投稿された回答
  • Durron-ワンパス回答をベースとして使用した私の修正された回答
  • DurronTwo -2パス回答(@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

そして、最初の100万ロングのみで実行された場合のサンプルランタイムは次のとおりです。

 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最初の100万の数で何度も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;
}

そして私のベンチマークハーネス:(グーグルキャリパー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つのビット演算を使用して(一種の)1行で実行できます。

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

結果のインデックスは、1)残基、2)残基+ 0xFFFFFF、または3)残基のいずれか+ 0x1FFFFFEです。もちろん、0xFFFFFF約3mbファイルの剰余のルックアップテーブルが必要です(この場合、ASCIIテキストの10進数として保存され、最適ではありませんが、aByteBufferなどで明らかに改善できます。しかし、それは事前計算であるため、そうではありません。 。tは問題でそんなにあなたはここにファイルを見つけることができます(またはそれを自分で生成します):

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私が実行したすべての試行で(バージョン1)を打ち負かしました。

 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
18
Bill the Lizard 2008-11-18 04:58.

現在のソリューションの場合と同様に、ニュートン法を使用して整数平方根を計算し、この数値を2乗して確認する方が、はるかに高速です。ニュートン法は、他のいくつかの回答で言及されているカーマック解の基礎です。ルートの整数部分のみに関心があるため、より迅速な回答を得ることができ、近似アルゴリズムをより早く停止できます。

試すことができるもう1つの最適化:数値の数字根が1、4、7、または9で終わらない場合、その数値は完全な平方ではありません。これは、低速の平方根アルゴリズムを適用する前に、入力の60%を削除する簡単な方法として使用できます。

15
mrzl 2008-11-19 05:57.

この関数をすべての正の64ビット符号付き整数で機能させたい

Math.sqrt()入力パラメータとしてdoubleを使用するため、2 ^ 53より大きい整数では正確な結果が得られません。

13
Cyrille Ka 2008-12-03 00:00.

念のために言うと、別のアプローチは素数分解を使用することです。分解のすべての要素が偶数の場合、その数は完全な二乗です。したがって、必要なのは、ある数が素数の2乗の積として分解できるかどうかを確認することです。もちろん、そのような分解を取得する必要はありません。存在するかどうかを確認するだけです。

まず、2 ^ 32未満の素数の2乗のテーブルを作成します。これは、この制限までのすべての整数のテーブルよりもはるかに小さいです。

その場合の解決策は次のようになります。

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

私はそれが少し不可解だと思います。すべてのステップで、素数の2乗が入力数を除算することを確認します。もしそうなら、それは可能な限り数を正方形で割って、この正方形を素数分解から取り除きます。このプロセスによって1になった場合、入力数は素数の2乗の分解でした。正方形が数自体よりも大きくなると、この正方形またはより大きな正方形がそれを分割する方法がないため、その数を素数の正方形の分解にすることはできません。

今日の平方根はハードウェアで行われ、ここで素数を計算する必要があることを考えると、このソリューションははるかに遅いと思います。しかし、mrzlが彼の答えで述べているように、2 ^ 54を超えて機能しないsqrtを使用したソリューションよりも良い結果が得られるはずです。

11
Hugh Allen 2008-11-29 17:52.

d完全な正方形の最後の桁は特定の値しかとることができないことが指摘されています。最後のd桁(塩基でb数)がn残りのと同じであるnことにより、分割されbd、すなわち、。C表記でn % pow(b, d)

これは、任意の係数mに一般化できます。n % m数の一部を完全な正方形から除外するために使用できます。現在使用しているモジュラスは64で、これにより12が可能になります。余りの19%、可能な限り正方形。少しコーディングするだけで、モジュラス110880が見つかりました。これは、2016年のみを許可します。余りの1.8%を可能な正方形として。したがって、モジュラス演算(つまり除算)のコストと、マシンの平方根に対するテーブルルックアップのコストによっては、このモジュラスを使用する方が速い場合があります。

ちなみに、Javaにルックアップテーブルのビットのパックされた配列を格納する方法がある場合は、それを使用しないでください。110880 32ビットワードは最近RAMが少なく、マシンワードのフェッチはシングルビットのフェッチよりも高速になります。

11
Colonel Panic 2013-10-16 05:48.

整数問題は整数解に値します。したがって、

(非負の)整数で二分探索を実行して、のような最大の整数tを見つけますt**2 <= n。次に、r**2 = n正確にテストします。これには時間O(log n)がかかります。

セットに制限がないために正の整数を二分探索する方法がわからない場合は、簡単です。まずf(t) = t**2 - n、2の累乗で増加する関数f(上記)を計算します。それが正に変わるのを見るとき、あなたは上限を見つけました。次に、標準の二分探索を行うことができます。

10
dfeuer 2014-07-14 00:17.

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;

パフォーマンスに影響します。

9
BobbyShaftoe 2008-12-05 19:36.

パフォーマンスのために、あなたは非常にしばしばいくつかの妥協をしなければなりません。他の人はさまざまな方法を表現していますが、カーマックのハッキングは特定のNの値まで高速であることに気付きました。次に、「n」を確認し、その数N未満の場合は、カーマックのハックを使用するか、説明されている他の方法を使用します。ここの答えで。

8
finnw 2010-05-07 03:29.

これは、このスレッドで他の人が提案した手法の組み合わせを使用して、私が思いつくことができる最速のJava実装です。

  • Mod-256テスト
  • 不正確なmod-3465テスト(誤検出を犠牲にして整数除算を回避します)
  • 浮動小数点平方根、丸めて入力値と比較

私もこれらの変更を試しましたが、パフォーマンスには役立ちませんでした。

  • 追加のmod-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;
        }
    }

}
7
David Lehavi 2009-01-02 12:12.

最初からNの2乗部分を取り除く必要があります。

2回目の編集以下のmの魔法の表現は次のようになります

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;

最初の編集:

マイナーな改善:

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

最初の編集の終わり

いつものように続けます。このように、浮動小数点部分に到達するまでに、2乗部分が奇数(約半分)であるすべての数値をすでに取り除き、残りの1/8のみを考慮します。つまり、浮動小数点部分を数値の6%で実行します。

7
bgiles 2009-03-26 05:23.

プロジェクトオイラーはタグで言及されており、その問題の多くは番号を確認する必要があります>> 2^64。上記の最適化のほとんどは、80バイトのバッファーを使用している場合は簡単に機能しません。

Java BigIntegerと、整数でより適切に機能するニュートン法のわずかに変更されたバージョンを使用しました。問題は、正確な二乗が理由ではなくにn^2収束し、最終的なエラーが最終的な除数の1ステップ下にあり、アルゴリズムが終了したことでした。エラーを計算する前に元の引数に1を追加することで、簡単に修正できました。(立方根などに2つ追加します)(n-1)nn^2-1 = (n-1)(n+1)

このアルゴリズムの優れた属性の1つは、数値が完全な平方であるかどうかをすぐに判断できることです。ニュートン法の最終誤差(補正ではない)はゼロになります。簡単な変更floor(sqrt(x))により、最も近い整数の代わりにすばやく計算することもできます。これは、いくつかのオイラー問題で便利です。

6
Brent.Longborough 2009-01-02 12:27.

これは、Rubyで、この質問に特に適合した、古いマーチャント計算機アルゴリズム(申し訳ありませんが、参照はありません)の10進数から2進数への作り直しです。

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;
};
6
hydrodog 2009-03-12 03:25.

すでに述べたように、sqrt呼び出しは完全に正確ではありませんが、速度の点で他の答えを吹き飛ばさないことは興味深く、有益です。結局のところ、sqrtのアセンブリ言語命令のシーケンスはごくわずかです。Intelにはハードウェア命令がありますが、IEEEに準拠していないため、Javaでは使用されていないと思います。

では、なぜ遅いのですか?Javaは実際にはJNIを介してCルーチンを呼び出しているため、Javaサブルーチンを呼び出すよりも実際には遅く、インラインで呼び出すよりも遅くなります。これは非常に煩わしいので、Javaはより良い解決策を考え出す必要があります。つまり、必要に応じて浮動小数点ライブラリ呼び出しを組み込む必要があります。しかたがない。

C ++では、すべての複雑な選択肢が速度を落とすと思いますが、すべてをチェックしていません。私がしたこと、そしてJavaの人々が役立つと思うことは、単純なハックであり、A。Rexによって提案された特別なケースのテストの拡張です。境界チェックされていないビット配列として単一のlong値を使用します。そうすれば、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は、私のcore2duoマシンで約1/3の時間で実行されます。同じ線に沿ってさらに微調整することで、平均して時間をさらに短縮できると思いますが、チェックするたびに、より多くのテストと引き換えに、より多くの除去を行うため、その道をあまり遠くに行くことはできません。

確かに、ネガティブを個別にテストするのではなく、同じ方法で上位6ビットをチェックできます。

私がしているのは可能な正方形を排除することだけですが、潜在的なケースがある場合は、元のインライン化されたisPerfectSquareを呼び出す必要があることに注意してください。

init2ルーチンは、pp1とpp2の静的な値を初期化するために1回呼び出されます。C ++での実装では、unsigned long longを使用しているため、署名されているため、>>>演算子を使用する必要があることに注意してください。

配列を境界チェックする本質的な必要はありませんが、Javaのオプティマイザーはこれを非常に迅速に把握する必要があるため、私はそれらを非難しません。

6
Jonny Heggheim 2009-05-25 15:22.

一部の入力でほぼ正しい方法を使用するというアイデアが好きです。これは、「オフセット」が高いバージョンです。コードは機能しているようで、私の単純なテストケースに合格しています。

交換するだけです:

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);
}
6
nabam serbang 2010-12-04 15:22.

一般的なビット長を考慮して(ここでは特定のタイプを使用しましたが)、以下のように単純なアルゴリズムを設計しようとしました。最初は、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;  
}  
5
Fractaly 2011-12-04 18:30.

正方形の最後のnビットが観察されたときに、考えられるすべての結果を確認しました。より多くのビットを連続して調べることにより、最大5/6の入力を排除できます。私は実際にこれをFermatの因数分解アルゴリズムを実装するように設計しましたが、そこでは非常に高速です。

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)
  • cの形式は(2 << 2k + 2)-1
  • dの形式は(2 << 2k-1)* 10

    最初に、2の累乗のモジュラスを持つ平方残差があるかどうかをテストし、次に最終モジュラスに基づいてテストし、次にMath.sqrtを使用して最終テストを実行します。私はトップポストからアイデアを思いつき、それを拡張しようとしました。コメントや提案をいただければ幸いです。

    更新:モジュラス(modSq)と44352のモジュラスベースによるテストを使用すると、私のテストは、最大1,000,000,000の数のOPの更新の時間の96%で実行されます。

  • 2
    aventurin 2018-10-18 12:49.

    整数演算を使用したニュートン法

    整数以外の操作を避けたい場合は、以下の方法を使用できます。基本的に、整数演算用に修正されたニュートン法を使用します。

    /**
     * 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。ただし、他のいくつかの投稿で説明されているフィルタリングメカニズムを使用することで、パフォーマンスを向上させることができます。

    2
    MWB 2018-12-28 07:13.

    これが分割統治法の解決策です。

    自然数(number)の平方根が自然数()のsolution場合solution、次の桁数に基づいて範囲を簡単に決定できますnumber

    • number1桁あり:solution範囲内= 1〜4
    • number2桁あります:solution範囲内= 3-10
    • number3桁の数字:solution範囲内= 10〜40
    • number4桁の数字:solution範囲内= 30-100
    • number5桁の数字: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"
    
    2
    Sajjad Ali Vayani 2019-11-01 21:53.

    数が完全な平方であるとすると、数の平方根。

    複雑さは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;
    }
    
    1
    Elijah 2008-11-18 13:29.

    速度が懸念される場合は、最も一般的に使用される入力のセットとその値をルックアップテーブルに分割してから、例外的なケースに対して思いついた最適化された魔法のアルゴリズムを実行してみませんか?

    1
    paulmurray 2009-02-22 20:22.

    最後のX桁がNの場合、それよりもはるかに効率的に「完全な正方形にすることはできません」をパックできるはずです。Java 32ビット整数を使用し、数値の最後の16ビットをチェックするのに十分なデータを生成します。これは2048個の16進整数値です。

    ..。

    OK。少し超えた数論に遭遇したか、コードにバグがあります。いずれにせよ、ここにコードがあります:

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

    結果は次のとおりです。

    (ed:prettify.jsのパフォーマンスが低いために削除されました。リビジョン履歴を表示して確認してください。)

    1
    Albert van der Horst 2018-12-27 16:42.

    ニュートン法による平方根の計算は、開始値が妥当であれば、非常に高速です。ただし、妥当な開始値はなく、実際には、二分法とlog(2 ^ 64)の動作で終わります。
    本当に速くするには、妥当な開始値に到達するための速い方法が必要です。つまり、機械語に降りる必要があります。プロセッサがPentiumでPOPCNTのような命令を提供する場合、先行ゼロをカウントし、それを使用して、有効ビットの半分の開始値を持つことができます。注意して、常に十分な固定数のニュートンステップを見つけることができます。(したがって、ループして非常に高速に実行する必要があります。)

    2番目の解決策は、浮動小数点機能を使用することです。浮動小数点機能は、(i87コプロセッサーのように)高速な平方根計算を行う場合があります。exp()およびlog()を介したエクスカーションでさえ、ニュートンが二分探索に縮退するよりも高速な場合があります。これには注意が必要な側面があります。プロセッサに依存して、後で改善が必要かどうかを分析します。

    3番目の解決策は、わずかに異なる問題を解決しますが、状況が質問で説明されているため、言及する価値があります。わずかに異なる数の非常に多くの平方根を計算する場合、開始値を再初期化せずに、前の計算を中断したままにしておくと、ニュートン反復を使用できます。私はこれを使用して、少なくとも1つのオイラー問題で成功しました。

    0
    Celestial M Weasel 2008-11-18 03:48.

    整数のサイズが有限であることを考えると、速度が必要な場合は、(a)パラメーターをサイズごとに(たとえば、最大ビットセットごとにカテゴリに)分割し、完全な正方形の配列に対して値をチェックするのが最も簡単な方法だと思います。その範囲内。

    0
    Ben 2009-03-12 03:18.

    カーマックスの方法に関しては、もう一度繰り返すだけで非常に簡単なようです。これにより、精度の桁数が2倍になるはずです。結局のところ、これは非常に切り捨てられた反復法であり、ニュートンの方法であり、最初の推測は非常に優れています。

    あなたの現在のベストに関して、私は2つのマイクロ最適化を見ます:

    • mod255を使用して、チェック後にチェックと0を移動します。
    • 通常(75%)の場合のすべてのチェックをスキップするために、4の分割累乗を再配置します。

    すなわち:

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

    シンプルな方がいいかもしれません

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

    明らかに、各チェックポイントでいくつの数がカリングされるかを知ることは興味深いでしょう-私はむしろ、チェックが本当に独立していることを疑っています。

    Related questions

    MORE COOL STUFF

    ケイト・ブランシェットは3日間一緒に夫と一緒に寝て、25年経ってもまだ夫と結婚しています

    ケイト・ブランシェットは3日間一緒に夫と一緒に寝て、25年経ってもまだ夫と結婚しています

    ケイト・ブランシェットは、夫に会ったとき、典型的な交際のアドバイスに逆らいました。

    マイケルシーンが非営利の俳優である理由

    マイケルシーンが非営利の俳優である理由

    マイケルシーンは非営利の俳優ですが、それは正確にはどういう意味ですか?

    ホールマークスターのコリンエッグレスフィールドがRomaDramaLiveでスリル満点のファンと出会う![エクスクルーシブ]

    ホールマークスターのコリンエッグレスフィールドがRomaDramaLiveでスリル満点のファンと出会う![エクスクルーシブ]

    特徴的なスターのコリン・エッグレスフィールドは、RomaDrama Liveでのスリル満点のファンとの出会いについて料理しました!加えて、大会での彼のINSPIREプログラム。

    「たどりつけば」をオンラインでストリーミングできない理由

    「たどりつけば」をオンラインでストリーミングできない理由

    ノーザンエクスポージャーが90年代の最も人気のある番組の1つになった理由を確認するには、Blu-rayまたはDVDプレーヤーをほこりで払う必要があります。

    バイオニック読書はあなたをより速く読むことができますか?

    バイオニック読書はあなたをより速く読むことができますか?

    BionicReadingアプリの人気が爆発的に高まっています。しかし、それは本当にあなたを速読術にすることができますか?

    ドミニカのボイリング湖:アクセスは簡単ではありませんが、ハイキングする価値があります

    ドミニカのボイリング湖:アクセスは簡単ではありませんが、ハイキングする価値があります

    ドミニカのボイリング湖は、世界で2番目に大きいボイリング湖です。そこにたどり着くまでのトレッキングは大変で長いですが、努力する価値は十分にあります。

    私たちの水をきれいに保つのを助けるためにあなたの髪を寄付してください

    私たちの水をきれいに保つのを助けるためにあなたの髪を寄付してください

    サロンからのヘアトリミングや個人的な寄付は、油流出を吸収して環境を保護するのに役立つマットとして再利用できます。

    ホワイトハウスの最も記憶に残る結婚式を見てください

    ホワイトハウスの最も記憶に残る結婚式を見てください

    過去200年以上の間にホワイトハウスで結婚したのはほんの数人です。彼らは誰でしたか、そしてそこで結婚式を獲得するために何が必要ですか?

    グッドジョブ、ESPN

    グッドジョブ、ESPN

    今日のホワイトハウスの記者会見で、報道官のサラ・ハッカビー・サンダースは、スポーツセンターのホストであるイェメル・ヒルのドナルド・トランプに関する最近のツイートについてコメントするよう求められ、大統領と彼の政策を白人至上主義者として説明した。ヒルは彼女のつぶやきのためにESPNによって公に叱責されました。

    アマゾンからのこのレミントンツールセールで髪を整える、スタイルを整える、乾かす、または取り除く

    アマゾンからのこのレミントンツールセールで髪を整える、スタイルを整える、乾かす、または取り除く

    基本的に頭のあらゆる部分から髪の毛を整えたり、ブロードライしたり、まっすぐにしたり、脱毛したりする必要がある場合は、このレミントンゴールドボックスが最適です。今日だけ、Amazonは、すでに人気のあるShortcut Pro Self-HaircutKitやPearlPro Ceramic Flat Ironのように、グルーミングをはるかに簡単にするヘアツールをマークダウンしています。

    カナダの元桂冠詩人が、史上最高の文化の盗用でトゥパックとマヤアンジェロウから盗んだ

    カナダの元桂冠詩人が、史上最高の文化の盗用でトゥパックとマヤアンジェロウから盗んだ

    トゥパックシャクール(ティムモーゼンフェルダー/ゲッティイメージズ); マヤアンジェロウ(マーティンゴッドウィン/ゲッティイメージズ)移動、テイラースウィフト。ケンドールとカイリーは、必要な数の座席を持っています。

    テスラは、ハリケーンイルマによる避難を容易にするために、フロリダでの車両の範囲を拡大しています

    テスラは、ハリケーンイルマによる避難を容易にするために、フロリダでの車両の範囲を拡大しています

    写真:Tesla Motorsフロリダに住んでいて、Tesla Model S、Model X 60、またはModel 60Dを使用している場合、車の自律性は50 km(約30マイル)高くなります。これは失敗ではなく、ハリケーンイルマによる避難作業を容易にするために会社自身が命じた自治権の一時的な延長です。

    Zendaya Wishes Boyfriend Tom Holland Happy Birthday with Cuddly Photo: He 'Makes Me the Happiest'

    Zendaya Wishes Boyfriend Tom Holland Happy Birthday with Cuddly Photo: He 'Makes Me the Happiest'

    Zendaya shared a sweet photo in honor of boyfriend Tom Holland's 26th birthday Wednesday

    小さな女性:脳卒中を患った後に病院から解放されたアトランタのジューシーな赤ちゃん:「まだ癒し」

    小さな女性:脳卒中を患った後に病院から解放されたアトランタのジューシーな赤ちゃん:「まだ癒し」

    シーレン「Ms.JuicyBaby」ピアソンは、先月脳卒中で入院した後、「もう一度たくさんのことをする方法を学ばなければならない」ため、言語療法を受けていることを明らかにしました。

    エマストーンは彼女のクリフサイドマリブビーチハウスを420万ドルでリストアップしています—中を見てください!

    エマストーンは彼女のクリフサイドマリブビーチハウスを420万ドルでリストアップしています—中を見てください!

    オスカー受賞者の世紀半ばの家には、3つのベッドルーム、2つのバス、オーシャンフロントの景色があります。

    ジーニー・メイ・ジェンキンスは、母乳育児の経験の中で、彼女は「本当に、本当に落ち込んでいる」と言います

    ジーニー・メイ・ジェンキンスは、母乳育児の経験の中で、彼女は「本当に、本当に落ち込んでいる」と言います

    ジーニー・メイ・ジェンキンスは、生後4か月の娘、モナコに母乳育児をしていると語った。

    投資ノート:Bioscout AU$300万シード

    投資ノート:Bioscout AU$300万シード

    Bioscoutは、農家を運転席に置くという使命を負っています。Artesian(GrainInnovate)やUniseedと並んで、最新のシードラウンドでチームを支援できることをうれしく思います。問題真菌症による重大な作物の損失は、農民にとって試練であることが証明されています。

    リトルマーケットリサーチ1| 2022年のクイックグリンプス遠隔医療市場

    リトルマーケットリサーチ1| 2022年のクイックグリンプス遠隔医療市場

    遠隔医療は、パンデミック後の時代では新しいものではなく、時代遅れの分野でもありません。しかし、業界を詳しく見ると、需要と供給の強力な持続可能性と、米国で絶え間ない革命となる強力な潜在的成長曲線を示しています。

    スタートアップ資金調達環境:タイのスタートアップエコシステムの次は何ですか?

    スタートアップ資金調達環境:タイのスタートアップエコシステムの次は何ですか?

    2021年は、世界的なベンチャーキャピタル(VC)の資金調達にとって記録的な年でした。DealStreetAsiaによると、東南アジアも例外ではなく、この地域では年間で記録的な25の新しいユニコーンが採掘されました。

    ムーアの法則を超えて

    ムーアの法則を超えて

    計算に対する私たちの欲求とムーアの法則が提供できるものとの間には、指数関数的に増大するギャップがあります。私たちの文明は計算に基づいています—建築と想像力の現在の限界を超える技術を見つけなければなりません。

    Language