题 确定整数平方根是否为整数的最快方法


我正在寻找确定是否最快的方法 long value是一个完美的正方形(即它的平方根是另一个整数):

  1. 我通过使用内置的Math.sqrt()以简单的方式完成了它 功能,但我想知道是否有办法更快地做到这一点 将自己限制为仅限整数的域。
  2. 维护查找表是不切实际的(因为有关于 231.5 平方小于2的整数63)。

这是我现在正在做的非常简单直接的方式:

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

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

注意:我在许多人中使用此功能 项目欧拉 问题。所以没有人会不得不维护这段代码。这种微优化实际上可以产生影响,因为部分挑战是在不到一分钟内完成每个算法,并且在某些问题中需要将此函数调用数百万次。


发布的新解决方案 A.雷克斯 已被证明更快。在前10亿个整数的运行中,该解决方案仅需要原始解决方案使用的34%的时间。虽然约翰卡马克的黑客攻击对于小的价值来说好一点 ñ,与此解决方案相比的好处非常小。

以下是转换为Java的A. Rex解决方案:

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

我尝试了下面介绍的不同解决方案。

  • 经过详尽的测试,我发现添加了 0.5 不需要Math.sqrt()的结果,至少不在我的机器上。
  • 约翰卡马克黑客 更快,但它从n = 410881开始给出了错误的结果。但是,正如所建议的那样 BobbyShaftoe,我们可以使用Carmack hack for n <410881。
  • 牛顿的方法慢了一点 Math.sqrt()。这可能是因为 Math.sqrt() 使用与牛顿方法类似的东西,但在硬件中实现,因此它比在Java中快得多。此外,牛顿的方法仍然需要使用双打。
  • 一个改进的牛顿方法,它使用了一些技巧,只涉及整数数学,需要一些黑客来避免溢出(我希望这个函数适用于所有正64位有符号整数),它仍然比 Math.sqrt()
  • 二进制斩甚至更慢。这是有道理的,因为二进制斩波平均需要16遍才能找到64位数的平方根。

确实改进的一个建议是由 约翰D.库克。您可以观察到完美正方形的最后一个十六进制数字(即最后4位)必须是0,1,4或9.这意味着75%的数字可以立即作为可能的正方形消除。实施此解决方案后,运行时间缩短了约50%。

根据John的建议,我调查了最后一个的属性 ñ 一个完美的正方形。通过分析最后6位,我发现最后6位中64个值中只有12个是可能的。这意味着可以在不使用任何数学的情况下消除81%的值。实施此解决方案使运行时间减少了8%(与我原来的算法相比)。分析超过6位导致可能的结束位列表太大而不实用。

这是我使用的代码,它在原始算法所需的42%的时间内运行(基于前1亿个整数的运行)。对于的值 ñ 小于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的测试,使用 or 语句在C ++中比使用a更快 switch,但在Java和C#之间似乎没有区别 or 和 switch
  • 我还尝试制作一个查找表(作为64个布尔值的私有静态数组)。然后而不是开关或 or 声明,我只想说 if(lookup[(int)(n&0x3F)]) { test } else return false;。令我惊讶的是,这只是(稍微)慢了。 我不知道为什么。  这是因为 在Java中检查数组边界

1210


起源


这是Java代码,其中int == 32位,long == 64位,两者都是有符号的。 - Kip
@Shreevasta:我已经对大值(大于2 ^ 53)进行了一些测试,并且你的方法给出了一些误报。遇到的第一个是n = 9007199326062755,这不是一个完美的正方形,而是作为一个返回。 - Kip
请不要称它为“John Carmack hack”。他没想出来。 - user9282
@mamama - 也许,但这归功于他。亨利福特并没有发明汽车,赖特兄弟没有发明这架飞机,并且Galleleo并不是第一个发现地球围绕太阳旋转的人......这个世界是由偷来的发明组成的(和爱)。 - Robert Fraser
通过使用类似的东西,你可能会在'quickfail'中获得极小的速度提升 ((1<<(n&15))|65004) != 0而不是三个单独的检查。 - Nabb


答案:


我发现一种方法比你的6bits + Carmack + sqrt代码快〜35%,至少用我的CPU(x86)和编程语言(C / C ++)。您的结果可能会有所不同,特别是因为我不知道Java因素将如何发挥作用。

我的方法有三个方面:

  1. 首先,筛选出明显的答案。这包括负数并查看最后4位。 (我发现看到最后六个没有帮助。)我也回答是0.(在阅读下面的代码时,请注意我的输入是 int64 x。)
    if( x < 0 || (x&2) || ((x & 7) == 5) || ((x & 11) == 8) )
        return false;
    if( x == 0 )
        return true;
  2. 接下来,检查它是否为方形模255 = 3 * 5 * 17.因为这是三个不同质数的乘积,所以只有大约1/8的mod 255残差是正方形。但是,根据我的经验,调用模运算符(%)的成本高于获得的效益,因此我使用涉及255 = 2 ^ 8-1的位技巧来计算残差。 (无论好坏,我没有使用从单词中读取单个字节的技巧,只是按位 - 并且移位。)
    int64 y = x;
    y = (y & 4294967295LL) + (y >> 32); 
    y = (y & 65535) + (y >> 16);
    y = (y & 255) + ((y >> 8) & 255) + (y >> 16);
    // At this point, y is between 0 and 511.  More code can reduce it farther.
    

602



哇!我将尝试将其转换为Java并进行比较,并对结果进行准确性检查。我会让你知道我发现了什么。 - Kip
哇,这很美。我之前看到Hensel解除了(计算多项式的根模数为素数)但是我甚至没有意识到引理可以一直小心地降低以计算数字的平方根;这是......令人振奋:) - ShreevatsaR
@nightcracker没有。 9 < 0 => false, 9&2 => 0, 9&7 == 5 => false, 9&11 == 8 => false。 - primo
Maartinus发布了 快2倍的解决方案 (并且更短)在下面,稍后,似乎没有得到太多的爱。 - Jason C
通过滤除明显的正方形,似乎可以获得不同解决方案中的许多速度优势。有没有人通过Maartinus的解决方案对过滤掉的情况进行基准测试,然后只使用sqrt函数作为内置函数? - user1914292


我参加晚会很晚,但我希望能提供更好的答案;更短和(假设我的 基准 是正确的)也很多 更快

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%。

第二个测试捕获了因子分解中具有奇数个二进制数的所有数字。方法 Long.numberOfTrailingZeros 因为它被JIT编辑成单个i86指令,所以速度非常快。

在删除尾随零之后,第三个测试以二进制形式处理以011,101或111结尾的数字,这些数字不是完美的正方形。它也关心负数,也处理0。

最后的测试可以追溯到 double 算术。如 double 只有53位尾数, 转换自 long 至 double 包括舍入大价值。尽管如此,测试是正确的(除非 证明 是错的)。

试图纳入mod255的想法是不成功的。


262



隐含的移位值掩盖有点......邪恶。你知道为什么那个Java规范吗? - dfeuer
特别是关于您的代码的问题:为什么需要检查奇数是否结束 001?是不是由处理 goodMask 测试? - dfeuer
@dfeuer我猜有两个原因:1。更多的转移毫无意义。 2.就像HW工作一样,任何使用按位运算的人都对性能感兴趣,所以做其他事情都是错误的。 -  该 goodMask 测试做到了,但它做到了 之前 正确的转变。所以你必须重复它,但这样它更简单,AFAIK更快,同样好。 - maaartinus
@Sebastian一个可能更好的测试: if ((x & (7 | Integer.MIN_VALUE)) != 1) return x == 0;。 - maaartinus
“因为双尾只有56位尾数” - >我会说它更有可能有一个 53位 一。 也 - chux


你必须做一些基准测试。最佳算法取决于输入的分布。

您的算法可能几乎是最优的,但您可能希望在调用平方根例程之前快速检查以排除某些可能性。例如,通过逐位“和”来查看十六进制数字的最后一位数字。完美的正方形只能在基数16中以0,1,4或9结束。因此,对于75%的输入(假设它们是均匀分布的),您可以避免调用平方根来换取一些非常快速的比特。

Kip对实现hex技巧的以下代码进行了基准测试。测试1到100,000,000时,此代码的运行速度是原始代码的两倍。

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

    switch((int)(n & 0xF))
    {
    case 0: case 1: case 4: case 9:
        long tst = (long)Math.sqrt(n);
        return tst*tst == n;

    default:
        return false;
    }
}

当我在C ++中测试类似代码时,它实际上比原始代码运行得慢。但是,当我删除switch语句时,十六进制技巧再次使代码快两倍。

int isPerfectSquare(int n)
{
    int h = n & 0xF;  // h is the last hex "digit"
    if (h > 9)
        return 0;
    // Use lazy evaluation to jump out of the if statement as soon as possible
    if (h != 2 && h != 3 && h != 5 && h != 6 && h != 7 && h != 8)
    {
        int t = (int) floor( sqrt((double) n) + 0.5 );
        return t*t == n;
    }
    return 0;
}

消除switch语句对C#代码几乎没有影响。


118



这很聪明......不会想到这一点 - warren
关于尾随位的好点。我会尝试将该测试与其他一些评论结合起来。 - PeterAllenWebb
我对前1亿个整数进行了基准测试。这大约减少了所需的时间。 - Kip
精湛的解决方案想知道你是怎么想出来的吗?是一个相当成熟的原则还是你想出来的东西? :d - Jeel Shah
@LarsH没有必要添加0.5,请参阅我的解决方案以获取证据的链接。 - maaartinus


我在思考我在数值分析课程中度过的可怕时光。

然后我记得,这个函数在Quake源代码的网络周围盘旋:

float Q_rsqrt( float number )
{
  long i;
  float x2, y;
  const float threehalfs = 1.5F;

  x2 = number * 0.5F;
  y  = number;
  i  = * ( long * ) &y;  // evil floating point bit level hacking
  i  = 0x5f3759df - ( i >> 1 ); // wtf?
  y  = * ( float * ) &i;
  y  = y * ( threehalfs - ( x2 * y * y ) ); // 1st iteration
  // y  = y * ( threehalfs - ( x2 * y * y ) ); // 2nd iteration, this can be removed

  #ifndef Q3_VM
  #ifdef __linux__
    assert( !isnan(y) ); // bk010122 - FPE?
  #endif
  #endif
  return y;
}

这基本上使用牛顿的近似函数计算平方根(不能记住确切的名称)。

它应该是可用的,甚至可能更快,它来自一个非凡的id软件的游戏!

它是用C ++编写的,但是一旦你得到这个想法,在Java中重用相同的技术应该不会太难:

我最初发现它: http://www.codemaestro.com/reviews/9

牛顿的方法在维基百科上解释: http://en.wikipedia.org/wiki/Newton%27s_method

您可以按照链接获取有关其工作原理的更多说明,但如果您不在乎,那么这大致是我在阅读博客和参加数值分析课程时所记得的:

  • * (long*) &y 基本上是一个快速转换为长函数,因此可以对原始字节应用整数运算。
  • 0x5f3759df - (i >> 1); line是近似函数的预先计算的种子值。
  • * (float*) &i 将值转换回浮点。
  • y = y * ( threehalfs - ( x2 * y * y ) ) line bascially再次迭代函数的值。

在函数迭代结果时,近似函数给出的值越精确。在Quake的情况下,一次迭代“足够好”,但如果它不适合你......那么你可以根据需要添加尽可能多的迭代。

这应该更快,因为它减少了天真平方根生成的除法运算的数量,直到简单除以2(实际上是a * 0.5F 乘法运算)并用几个固定数量的乘法运算代替它。


41



应该注意,这返回1 / sqrt(数字),而不是sqrt(数字)。我做了一些测试,这从n = 410881开始失败:当实际平方根为641时,John Carmack魔法公式返回642.00104。 - Kip
你可以看一下关于快速平方根的Chris Lomonts论文: lomont.org/Math/Papers/2003/InvSqrt.pdf 它使用与此处相同的技术,但具有不同的幻数。该论文解释了为什么选择了神奇的数字。 - Dan
也, beyond3d.com/content/articles/8 和 beyond3d.com/content/articles/15 阐明了这种方法的起源。它通常归功于John Carmack,但似乎原始代码(可能)由Gary Tarolli,Greg Walsh和其他人编写。 - Dan
你也不能用Java打字浮动和整数。 - Antimony
@Antimony谁说? FloatToIntBits 和 IntToFloatBits 自java 1.0.2以来一直存在。 - corsiKa


我不确定它是否会更快,甚至更准确,但你可以使用 约翰卡马克的魔法平方根,算法更快地解决平方根。您可以轻松地对所有可能的32位整数进行测试,并验证您实际上得到了正确的结果,因为它只是一个适应性。然而,现在我考虑一下,使用双打也是近似的,所以我不确定它会如何发挥作用。


33



适用于常规Carmack参考的+1! - Mitch Wheat
我相信卡马克的伎俩近来毫无意义。内置的sqrt指令比以前快得多,所以你最好只执行常规的平方根并测试结果是否为int。一如既往,以此为基准。 - jalf
这打破了n = 410881,John Carmack魔法公式返回642.00104,当实际平方根为641时。 - Kip
我最近在Java游戏中使用了Carmack的技巧,它非常有效,速度提高了大约40%,所以它仍然很有用,至少在Java中是这样。 - finnw
@Robert Fraser是整体帧率+ 40%。游戏中有一个粒子物理系统,它占据了几乎所有可用的CPU周期,由平方根函数和舍入到最接近整数函数(我还使用类似的比特杂乱的黑客进行了优化)。 - finnw


如果你做一个二进制印章试图找到“正确”的平方根,你可以很容易地检测出你得到的值是否足够接近:

(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 :你很亲密,但并不完美:回归虚假
  • 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次迭代。


30



我的钱就是采用这种方法。避免调用sqrt(),因为它正在计算一个完整的平方根,你只需要前几个数字。 - PeterAllenWebb
另一方面,如果浮点是在专用的FP单元中完成的,它可能正在使用各种有趣的技巧。我不想在没有基准测试的情况下下注:)(我今晚可能会尝试使用C#,只是为了看...) - Jon Skeet
硬件sqrts现在实际上非常快。 - Adam Rosenfield


我对这个线程中的几个算法进行了自己的分析,并得出了一些新的结果。你可以在这个答案的编辑历史中看到那些旧的结果,但它们不准确,因为我犯了一个错误,浪费时间分析几个不接近的算法。然而,从几个不同的答案中吸取教训,我现在有两种算法可以粉碎这个线程的“赢家”。这是我做的核心事情与其他人不同:

// This is faster because a number is divisible by 2^4 or more only 6% of the time
// and more than that a vanishingly small percentage.
while((x & 0x3) == 0) x >>= 2;
// This is effectively the same as the switch-case statement used in the original
// answer. 
if((x & 0x7) != 1) return false;

然而,这条简单的线路,大部分时间都添加了一个或两个非常快速的指令,大大简化了 switch-case 声明成一个if语句。但是,如果许多测试数字具有显着的两个幂因子,它可以添加到运行时。

以下算法如下:

  • 互联网  - 基普的答案
  • Durron  - 我使用一遍答案作为基础修改了答案
  • DurronTwo  - 我使用双通答案(@JohnnyHeggheim)修改了答案,并进行了一些其他的修改。

如果使用生成数字,这是一个示例运行时 Math.abs(java.util.Random.nextLong())

 0% Scenario{vm=java, trial=0, benchmark=Internet} 39673.40 ns; ?=378.78 ns @ 3 trials
33% Scenario{vm=java, trial=0, benchmark=Durron} 37785.75 ns; ?=478.86 ns @ 10 trials
67% Scenario{vm=java, trial=0, benchmark=DurronTwo} 35978.10 ns; ?=734.10 ns @ 10 trials

benchmark   us linear runtime
 Internet 39.7 ==============================
   Durron 37.8 ============================
DurronTwo 36.0 ===========================

vm: java
trial: 0

这是一个示例运行时,如果它仅在第一百万个长度上运行:

 0% Scenario{vm=java, trial=0, benchmark=Internet} 2933380.84 ns; ?=56939.84 ns @ 10 trials
33% Scenario{vm=java, trial=0, benchmark=Durron} 2243266.81 ns; ?=50537.62 ns @ 10 trials
67% Scenario{vm=java, trial=0, benchmark=DurronTwo} 3159227.68 ns; ?=10766.22 ns @ 3 trials

benchmark   ms linear runtime
 Internet 2.93 ===========================
   Durron 2.24 =====================
DurronTwo 3.16 ==============================

vm: java
trial: 0

如你看到的, DurronTwo 对于大型输入来说效果更好,因为它经常使用魔术技巧,但与第一种算法相比遭到破坏 Math.sqrt 因为这些数字要小得多。同时,更简单 Durron 是一个巨大的赢家,因为它永远不会在前一百万个数字中多次划分4次。

这里的 Durron

public final static boolean isPerfectSquareDurron(long n) {
    if(n < 0) return false;
    if(n == 0) return true;

    long x = n;
    // This is faster because a number is divisible by 16 only 6% of the time
    // and more than that a vanishingly small percentage.
    while((x & 0x3) == 0) x >>= 2;
    // This is effectively the same as the switch-case statement used in the original
    // answer. 
    if((x & 0x7) == 1) {

        long sqrt;
        if(x < 410881L)
        {
            int i;
            float x2, y;

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

            sqrt = (long)(1.0F/y);
        } else {
            sqrt = (long) Math.sqrt(x);
        }
        return sqrt*sqrt == x;
    }
    return false;
}

DurronTwo

public final static boolean isPerfectSquareDurronTwo(long n) {
    if(n < 0) return false;
    // Needed to prevent infinite loop
    if(n == 0) return true;

    long x = n;
    while((x & 0x3) == 0) x >>= 2;
    if((x & 0x7) == 1) {
        long sqrt;
        if (x < 41529141369L) {
            int i;
            float x2, y;

            x2 = x * 0.5F;
            y = x;
            i = Float.floatToRawIntBits(y);
            //using the magic number from 
            //http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
            //since it more accurate
            i = 0x5f375a86 - (i >> 1);
            y = Float.intBitsToFloat(i);
            y = y * (1.5F - (x2 * y * y));
            y = y * (1.5F - (x2 * y * y)); //Newton iteration, more accurate
            sqrt = (long) ((1.0F/y) + 0.2);
        } else {
            //Carmack hack gives incorrect answer for n >= 41529141369.
            sqrt = (long) Math.sqrt(x);
        }
        return sqrt*sqrt == x;
    }
    return false;
}

我的基准测试工具:(需要谷歌卡尺0.1-rc5)

public class SquareRootBenchmark {
    public static class Benchmark1 extends SimpleBenchmark {
        private static final int ARRAY_SIZE = 10000;
        long[] trials = new long[ARRAY_SIZE];

        @Override
        protected void setUp() throws Exception {
            Random r = new Random();
            for (int i = 0; i < ARRAY_SIZE; i++) {
                trials[i] = Math.abs(r.nextLong());
            }
        }


        public int timeInternet(int reps) {
            int trues = 0;
            for(int i = 0; i < reps; i++) {
                for(int j = 0; j < ARRAY_SIZE; j++) {
                    if(SquareRootAlgs.isPerfectSquareInternet(trials[j])) trues++;
                }
            }

            return trues;   
        }

        public int timeDurron(int reps) {
            int trues = 0;
            for(int i = 0; i < reps; i++) {
                for(int j = 0; j < ARRAY_SIZE; j++) {
                    if(SquareRootAlgs.isPerfectSquareDurron(trials[j])) trues++;
                }
            }

            return trues;   
        }

        public int timeDurronTwo(int reps) {
            int trues = 0;
            for(int i = 0; i < reps; i++) {
                for(int j = 0; j < ARRAY_SIZE; j++) {
                    if(SquareRootAlgs.isPerfectSquareDurronTwo(trials[j])) trues++;
                }
            }

            return trues;   
        }
    }

    public static void main(String... args) {
        Runner.main(Benchmark1.class, args);
    }
}

更新: 我已经制作了一种新算法,在某些情况下更快,在其他情况下更慢,我根据不同的输入得到了不同的基准。如果我们计算模数 0xFFFFFF = 3 x 3 x 5 x 7 x 13 x 17 x 241,我们可以消除不能为正方形的97.82%的数字。这可以(一种)在一行中完成,具有5个按位操作:

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

得到的指数是1)残留物,2)残留物 + 0xFFFFFF,或3)残留物 + 0x1FFFFFE。当然,我们需要有一个模数残差的查找表 0xFFFFFF,这是一个3mb文件(在这种情况下存储为ascii文本十进制数字,不是最佳的,但明确可以改进 ByteBuffer 等等。但由于这是预先计算,所以这并不重要。 你可以在这里找到这个文件 (或自己生成):

public final static boolean isPerfectSquareDurronThree(long n) {
    if(n < 0) return false;
    if(n == 0) return true;

    long x = n;
    while((x & 0x3) == 0) x >>= 2;
    if((x & 0x7) == 1) {
        if (!goodLookupSquares[(int) ((n & 0xFFFFFFl) + ((n >> 24) & 0xFFFFFFl) + (n >> 48))]) return false;
        long sqrt;
        if(x < 410881L)
        {
            int i;
            float x2, y;

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

            sqrt = (long)(1.0F/y);
        } else {
            sqrt = (long) Math.sqrt(x);
        }
        return sqrt*sqrt == x;
    }
    return false;
}

我加载到一个 boolean 像这样的数组:

private static boolean[] goodLookupSquares = null;

public static void initGoodLookupSquares() throws Exception {
    Scanner s = new Scanner(new File("24residues_squares.txt"));

    goodLookupSquares = new boolean[0x1FFFFFE];

    while(s.hasNextLine()) {
        int residue = Integer.valueOf(s.nextLine());
        goodLookupSquares[residue] = true;
        goodLookupSquares[residue + 0xFFFFFF] = true;
        goodLookupSquares[residue + 0x1FFFFFE] = true;
    }

    s.close();
}

示例运行时。它击败了 Durron (第一版)我跑的每一次试验。

 0% Scenario{vm=java, trial=0, benchmark=Internet} 40665.77 ns; ?=566.71 ns @ 10 trials
33% Scenario{vm=java, trial=0, benchmark=Durron} 38397.60 ns; ?=784.30 ns @ 10 trials
67% Scenario{vm=java, trial=0, benchmark=DurronThree} 36171.46 ns; ?=693.02 ns @ 10 trials

  benchmark   us linear runtime
   Internet 40.7 ==============================
     Durron 38.4 ============================
DurronThree 36.2 ==========================

vm: java
trial: 0

21



巨大的查找表似乎不是一个好主意。缓存未命中比x86硬件sqrt指令(~20个周期)慢(~100到150个周期)。在吞吐量方面,您可以维持许多未完成的缓存未命中,但您仍在驱逐其他有用的数据。如果它比任何其他选项快得多,那么一个巨大的查找表将是值得的,并且此功能是整个程序性能的主要因素。 - Peter Cordes


它的使用速度应该快得多 牛顿的方法 计算 整数平方根,然后将此数字平方并检查,就像在当前解决方案中一样。牛顿方法是其他一些答案中提到的Carmack解决方案的基础。您应该能够获得更快的答案,因为您只对根的整数部分感兴趣,允许您更快地停止近似算法。

您可以尝试的另一个优化:如果 数字根 一个数字不会结束 数字是1,4,7或9  一个完美的广场。在应用较慢的平方根算法之前,这可以用作消除60%输入的快速方法。


14



数字根在计算上严格等同于模数,因此应该与其他模数方法一起考虑,例如mod 16和mod 255。 - Christian Oudard
你确定数字根相当于模数吗?如链接所解释的那样,它似乎完全不同。请注意,列表是1,4,7,9而不是1,4,5,9。 - Fractaly
十进制系统中的数字根相当于使用模9(well dr(n)= 1 +((n-1)mod 9);因此也有轻微的移位)。数字0,1,4,5,9用于模16,0,1,4,7用于模9 - 对应于数字根的1,4,7,9。 - Hans Olsson