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




平方根 導出 (24)

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

私は、 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でチェックされるためです。

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

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

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

  • Mod-256テスト
  • 不正確なmod-3465テスト(いくつかの偽陽性を犠牲にして整数除算を避ける)
  • 浮動小数点平方根、丸め、入力値との比較

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

  • 追加のmod-255テスト
  • 入力値を4の累乗で割る
  • Fast Inverse Square Root(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;
        }
    }

}

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


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つのビットをフェッチするよりも高速になります。


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

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


Carmacメソッドに関しては、もう一度繰り返すだけで、正確さの桁数を2倍にすることは非常に簡単なようです。結局のところ、極端に切り詰められた反復的な方法である - ニュートンは、最初の非常に良い推測である。

あなたの現在のベストに関しては、2つのマイクロ最適化があります。

  • mod255を使用してチェックの後にチェックを移動します
  • 通常の(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;

明らかに、各チェックポイントでどれだけの数が選ばれるかを知ることは面白いでしょう。チェックが本当に独立していると疑いがあります。


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

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


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

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

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


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

私のアプローチは三倍です:

  1. まず、明白な答えを除外します。 これには、負の数と最後の4ビットを見ます。 (私は最後の6つを見ても助けにはならなかったことにint64 xました)。私も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の2乗であるかどうかをチェックします。これは3つの異なる素数の積であるため、残差mod 255の約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. 最後に、 Henselの補題に似た方法を使用して平方根を計算しようとします。 (私はそれが直接当てはまるとは思わないが、いくつかの変更を加えても動作する)。これを行う前に、2のすべてのべき乗を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.
    考え方は、各繰り返しで、xの "現在の"平方根であるrに1ビットを加えることです。 各平方根は、2の大きな累乗、すなわちt / 2をモジュロで正確に表します。 最後に、rとt / 2-rはx / 2を法とする平方根になります。 (rがxの平方根ならば、-rもそうであることに注意してください。これはモジュロ数でも当てはまりますが、いくつかのモジュロ数に注意してください.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をビットのトリックで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;
}

私はいくつかの入力に対してほぼ正しい方法を使うという考えが好きです。ここにはより高い「オフセット」を持つバージョンがあります。コードはうまくいくと思われ、私の単純なテストケースをパスします。

あなたの代わりに:

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

私は、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乗算演算)に至るまで、ナイーブな平方根で行われる除算演算の数を減らし、その代わりにいくつかの固定数の乗算演算で置き換えるので、より高速にすべきです。


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

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

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

例えば。

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


私は、正方形の最後の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%の時間で実行されます。


  • ここでは、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を要した。

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


    私はパーティーにかなり遅れていますが、より良い答えを提供したいと思います。 (私のbenchmarkが正しいと仮定して)短く、またはるかにfaster

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

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

    2番目の検定は、分解の2倍の奇数を持つすべての数をキャッチします。 Long.numberOfTrailingZerosメソッドはJIT-edを単一のi86命令に変換するので非常に高速です。

    末尾の0を削除した後、3番目のテストでは、完全な正方形ではない2進数で011,101、または111で終わる番号を処理します。 また、負の数を扱い、0も扱います。

    最後のテストはdouble算術に戻ります。 doubleは53ビットの仮数しかないので、 longからdoubleへの変換には大きな値の丸めが含まれます。 それにもかかわらず、テストは正しい( proofが間違っていない限り)。

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


    "正しい"平方根を見つけようとバイナリチョップを行うと、あなたが得た値が近くにあるかどうかを簡単に検出できます:

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

    計算したn^2ので、オプションは次のとおりです。

    • n^2 = target :done、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回の反復しか必要としません。


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

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

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


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

    非整数演算を避けたい場合は、以下の方法を使用できます。これは、基本的に整数計算に変更された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。ただし、他の投稿のいくつかに記載されているフィルタリングメカニズムを使用することで、そのパフォーマンスを向上させることができます。


    あなたの整数が有限のサイズであることを考えれば、私は、最も簡単な方法は、(a)大きさでパラメータを分割すること(例えば、最も大きなビットセットでカテゴリに分割する)、そして完全な四角形の配列その範囲内にある。


    プロジェクトオイラーはタグに記載されており、問題の多くは>> 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))を素早く計算することもできます。これはいくつかのオイラーの問題に便利です。


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






    perfect-square