java - quick - 平方根 導出




整数の平方根が整数かどうかを判断する最速の方法 (20)

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

  1. 私は組み込みのMath.sqrt()関数を使って簡単な方法でやったことがありますが、整数専用のドメインに自分自身を制限することで高速化する方法があるのだろうかと思います。
  2. ルックアップテーブルを維持することは驚くべきことである(なぜなら、平方が2 63未満の2 31.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;
}

注:多くのProject Eulerの問題でこの関数を使用しています。 したがって、誰もこのコードを維持する必要はありません。 そして、この種のマイクロ最適化は実際に差をつけることができます。なぜなら、挑戦の一部は1分以内にすべてのアルゴリズムを実行することであり、この機能はいくつかの問題で何百万回も呼び出される必要があります。

A. Rexによって投稿された新しいソリューションは、さらに高速であることが証明されています。 最初の10億の整数では、元のソリューションが使用した時間の34%しか必要としませんでした。 John Carmackのハックはnの値が小さいほど少し良いですが、このソリューションと比較した場合の利点はかなり小さいです。

ここでA.レックスのソリューションは、Javaに変換されます:

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

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

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

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

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

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

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

私は以下に示すさまざまなソリューションを試しました。

  • 網羅的なテストの結果、私はMath.sqrt()の結果に0.5を加えることは、少なくとも私のマシンでは必要ないことを発見しました。
  • John Carmackのハックは速かったのですが、n = 410881から不正確な結果が出ました。 しかし、 BobbyShaftoe示唆しているBobbyShaftoe 、私たちはCarmackハックをn <410881で使うことができます。
  • ニュートンの方法はMath.sqrt()よりも少し遅いです。 これは、おそらくMath.sqrt()がNewtonのメソッドに似たものを使用しているが、ハードウェアに実装されているためJavaよりもはるかに速いからです。 また、ニュートンの方法は依然としてダブルスの使用を必要とした。
  • 整数の数学だけが含まれるようにいくつかのテクニックを使った修正ニュートンの方法は、オーバーフローを避けるためにいくつかのハックを必要としました(この関数はすべての正の64ビット符号付き整数で動作するようにしたいと思っています)、それでもMath.sqrt()
  • バイナリチョップはさらに遅かった。 バイナリチョップは、平均して64ビットの数値の平方根を求めるために16回のパスを必要とするため、これは理にかなっています。

改善を示した1つの提案はJohn D. Cookによってなされました。 完全な正方形の最後の16進数(すなわち最後の4ビット)は、0,1,4、または9でなければならないことがわかります。これは、75%の数字を可能な四角形として直ちに削除できることを意味します。 このソリューションを実装すると、ランタイムが約50%削減されました。

Johnの提案から、完全な正方形の最後のnビットのプロパティを調べました。 最後の6ビットを分析することによって、最後の6ビットに対して64個の値のうちの12個だけが可能であることが分かった。 これは、数学を使用せずに81%の値を削除できることを意味します。 このソリューションを実装すると、(元のアルゴリズムと比較して)実行時にさらに8%の削減が実現しました。 6ビット以上を分析することは、可能な終了ビットのリストをもたらし、実用的には大きすぎる。

ここで私が使用したコードは元のアルゴリズムで必要とされる時間の42%で実行されます(最初の100万の整数を使って実行されます)。 nが410881未満の場合、元のアルゴリズムで必要とされる時間のわずか29%で実行されます。

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

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

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

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

  default:
    return false;
  }
}

  • Johnのテストによれば、C ++ではswitchを使うよりもusing or statementsが高速ですが、JavaとCでは#とswitch間に違いはありません。
  • 私はルックアップテーブル(64個のブール値のプライベート静的配列)を作成しようとしました。 次に、switch() or statement if(lookup[(int)(n&0x3F)]) { test } else return false;代わりに、 if(lookup[(int)(n&0x3F)]) { test } else return false; 。 私の驚いたことに、これは(わずかに)遅かった。 なぜ私は分からない。 これは、 配列の境界がJavaでチェックされるためです。

私はこの関数をすべての正の64ビット符号付き整数で動作させたい

Math.sqrt()倍精度を入力パラメータとして使用するので、2 ^ 53を超える整数に対して正確な結果は得られません。


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

非整数演算を避けたい場合は、以下の方法を使用できます。これは、基本的に整数計算に変更されたNewton's Methodを使用します。

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


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;

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


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

答えは印象的ですが、私は単純なチェックを見ませんでした:

それがセット(0,1,4,5,6,9)のメンバーであるかどうかを確認してください。そうでなければ、それはおそらく「完璧な正方形」にはなり得ません。

例えば。

4567 - 完璧な正方形にはなれません。


あなたは最初から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;

1回目の編集:

マイナーな改善:

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

最初の編集の終了

今度はいつも通り続ける。このようにして、浮動小数点の部分に到達するまでに、2電源部分が奇数(半分程度)の数をすべて取り除いた後、残った分の1/8しか考慮しません。つまり、数値の6%で浮動小数点部分を実行します。


これは、Rubyで、この質問のために特別に適応された古いMarchant計算アルゴリズムの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;
};

スピードが懸念される場合は、最も一般的に使用される入力とその値をルックアップテーブルに分割し、その例外的なケースに対応した最適化されたマジックアルゴリズムを作成してください。


ベンチマークをしなければならないでしょう。 最良のアルゴリズムは入力の分布に依存します。

あなたのアルゴリズムはほぼ最適ですが、平方根ルーチンを呼び出す前にいくつかの可能性を排除するために素早くチェックしたいかもしれません。 たとえば、数字の最後の桁を16進数で見て、ビット単位で "and"を実行します。 完全な四角形は、基数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#コードにほとんど影響はありませんでした。


整数問題には整数解が必要です。 従って

そのような最大の整数tを見つけるために、(非負の)整数に対してバイナリサーチを行いt**2 <= nます。次に、r**2 = n正確にテストする。これには時間がかかる(log n)。

セットが無制限であるために正の整数をバイナリ検索する方法がわからない場合は、簡単です。あなたf(t) = t**2 - nは2の累乗であなたの増加関数f(上記)を計算することから始めます。それが肯定的になるのを見ると、あなたは上限を見つけました。次に、標準バイナリ検索を実行できます。


私は、Numerical Analysisコースで過ごした恐ろしい時代について考えていました。

そして、私は覚えている、この関数は、ネットの周りに震えソースコードから回っていた:

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://www.codemaestro.com/reviews/9

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

あなたはそれがどのように機能するかについての詳細な説明のためにリンクをたどることができますが、それほど気にしなければ、これはブログを読んで、数値分析コースを取ることからおおよそ覚えています:

  • * (long*) &y基本的にので、整数演算は、生のバイトに適用することができ、高速変換対長の関数です。
  • この0x5f3759df - (i >> 1);線は、近似関数のための予め計算されたシード値である。
  • * (float*) &iバック浮動小数点の値に変換します。
  • そのy = y * ( threehalfs - ( x2 * y * y ) )行は基本的に関数の値を繰り返し反復します。

近似関数を使用すると、より正確な値が得られます。Quakeのケースでは、1回の反復は「十分に良い」ですが、それがあなたのためでない場合は、必要なだけ反復を追加できます。

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


Integer Square Rootを計算するためにhttp://en.wikipedia.org/wiki/Newton%27s_methodを使用する方がはるかに高速でなければなりません。ニュートンの方法は、いくつかの他の答えで言及されているCarmackソリューションの基礎です。近似アルゴリズムを早期に停止できるように、ルートの整数部分だけに興味があるので、より高速な答えを得ることができます。

もう1つの最適化を試すことができます。数字のデジタルルート1,4,7、または9で終わらない場合、数字は完全な正方形ではありません。これは、遅い平方根アルゴリズムを適用する前に、入力の60%を素早く排除するために使用できます。


d完璧な四角形の最後の数字は特定の値をとることしかできないことが指摘されています。d数字の最後の桁(基数b)はn、剰余nを除算したときの剰余と同じb dです。C表記n % pow(b, d)

これは任意のモジュラスmに一般化することができる。n % mいくつかのパーセンテージの数字を完全な正方形で除外するために使うことができます。あなたが現在使用しているモジュラスは64です。余りの19%、可能な正方形として。ちょっとコーディングして、私はモジュラス110880を見つけました。これは2016だけを許します。余りの1.8%が可能な四角である。したがって、モジュラス演算(すなわち、除算)のコストと、マシン上のテーブルルックアップ対平方根に応じて、このモジュラスを使用するほうが高速になる可能性があります。

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


ここでは、CPUサイクルの点でどのように比較するのか分かりませんが、最もシンプルで最も簡潔な方法です。これは、ルートが整数であるかどうかだけを知りたい場合に効果的です。それが整数であるかどうか本当に気にしているなら、それを把握することもできます。ここには単純な(そして純粋な)関数があります:

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

マイクロ最適化が必要ない場合、この答えは単純さと保守性の面で優れています。負の数値を取得する場合は、数値引数にMath.sqrt()引数としてMath.abs()を使用することをお勧めします。

私の3.6Ghzインテルi7-4790 CPUでは、0〜10,000,000でこのアルゴリズムを実行すると、1回の計算で平均35〜37ナノ秒かかりました。私は10回連続して実行し、1,000万sqrtの計算に費やした平均時間を印刷しました。それぞれの合計実行は、完了するのにわずか600msを要した。

少数の計算を実行している場合、以前の計算には少し時間がかかります。


これが以前に言及されているかどうかはわかりません。しかし、私はhere解決策を見つけました:

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

プロジェクトオイラーはタグに記載されており、問題の多くは>> 2 ^ 64という数値をチェックする必要があります。上記のほとんどの最適化は、80バイトのバッファで作業しているときには簡単には機能しません。

私は、Java BigIntegerとNewtonのメソッドを少し修正したものを使用しました。問題は、n ^ 2-1 =(n-1)(n + 1)であり、最終誤差が最終除数の1ステップ下にあるため、正確な正方形n ^ 2がnの代わりに(n-1)アルゴリズムは終了しました。エラーを計算する前に元の引数に1を加えて修正するのは簡単でした。(キューブルーツの場合は2つ追加するなど)

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


一般的なビット長(ここでは特定の型を使用していますが)を考慮して、以下のようにシンプルなアルゴリズムを設計しようとしました。最初は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;  
}  

最後のX桁がNであれば、それをはるかに効率的にパックすることは可能です。私はjavaの32ビット整数を使用し、数字の最後の16ビットをチェックするのに十分なデータを生成します。つまり、2048の16進数のint値です。

...

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でのパフォーマンス低下のために削除されました;改訂履歴を表示するには。)


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

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のべき乗の2乗を持つ正方形の残差があるかどうかをテストし、最終モジュラスに基づいてテストし、最後にMath.sqrtを使用して最終テストを行います。私はトップポストからアイデアを考え出し、それを拡張しようとしました。私は、コメントや提案を感謝します。

    アップデート:モジュラス(modSq)とモジュラスベース44352のテストを使用して、私のテストはOPのアップデートで最大1,000,000,000までの96%の時間で実行されます。


  • 私はこのスレッドでいくつかのアルゴリズムを自分で分析し、いくつかの新しい結果を思いつきました。これらの古い結果はこの回答の編集履歴で見ることができますが、間違いをすると正確ではなく、近くにないいくつかのアルゴリズムを分析する時間を無駄にします。しかし、いくつかの異なる答えからレッスンを引っ張って、私は今このスレッドの "勝者"を押しつぶす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の2乗の重要な要素を持つ場合は、ランタイムを追加する可能性があります。

    以下のアルゴリズムは以下の通りです:

    • インターネット - Kipの回答を掲載
    • Durron - 私の変更された答えは、ベースとしてワンパス答えを使用して
    • DurronTwo - (JonnnyHeggheimによる)2パスの答えを使った私の修正された答え。

    数値が生成されたランタイムの例を次に示します。 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;
    }
    

    と私のベンチマークハーネス:( Googleキャリパー0.1-rc5が必要)

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

    更新:いくつかのシナリオでは速く、他のシナリオでは遅くなるような新しいアルゴリズムを作成しました。異なる入力に基づいて異なるベンチマークを取得しました。モジュロを計算すると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進数として保存されますが、最適ではありませんが、a 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私が走ったすべての試行で勝利しました(バージョン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
    

    私はそれがより速く、正確であるかどうかは分かりませんが、http://www.codemaestro.com/reviews/9アルゴリズムを使って平方根をより速く解くことができます。おそらく、これをすべての可能な32ビット整数で簡単にテストし、実際に正しい結果が得られたかどうかを検証することができます。これは単なるappoximationなのであります。しかし、今私はそれについて考えて、ダブルスを使用して近似しているので、どのように再生になるかわからない。





    perfect-square