我正在寻找确定长值是否为完美平方(即其平方根是另一个整数)的最快方法:

我使用内置的Math.sqrt()以简单的方式完成了这项工作函数,但我想知道是否有一种方法可以通过将自己限制为仅限整数的域。维护查找表是不切实际的(因为平方小于263的231.5个整数)。

下面是我现在做的非常简单明了的方法:

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

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

注意:我在许多Project Euler问题中都使用了这个函数。因此,其他人将永远不必维护此代码。而这种微优化实际上可能会有所不同,因为挑战的一部分是在不到一分钟的时间内完成每一个算法,而在某些问题中,这个函数需要调用数百万次。


我尝试了不同的解决方案:

经过详尽的测试,我发现不需要在Math.sqrt()的结果上加0.5,至少在我的机器上是这样。快速平方根逆运算速度更快,但对于n>=410881,它给出了错误的结果。然而,正如BobbyShaftoe所建议的,我们可以在n<410881时使用FISR黑客。牛顿的方法比Math.sqrt()慢得多。这可能是因为Math.sqr()使用了类似于牛顿方法的东西,但在硬件中实现,所以比Java快得多。此外,牛顿法仍然需要使用双精度。一个经过修改的牛顿方法使用了一些技巧,因此只涉及整数数学,需要一些技巧来避免溢出(我希望这个函数可以处理所有64位有符号的正整数),而且它仍然比math.sqrt()慢。二元斩更慢。这是有意义的,因为二进制斩波平均需要16次才能找到64位数字的平方根。根据John的测试,在C++中使用or语句比使用switch更快,但在Java和C#中,or和switch之间似乎没有区别。我还尝试创建一个查找表(作为64个布尔值的私有静态数组)。然后,我只说if(lookup[(int)(n&0x3F)]){test}else return false;,而不是switch或or语句;。令我惊讶的是,这(只是稍微)慢了一些。这是因为在Java中检查数组边界。


当前回答

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

// 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语句。然而,如果许多被测试的数字具有两个因素的显著幂,则可以增加运行时。

以下算法如下:

互联网-Kip发布的答案Durron-我使用一次通过答案作为基础的修改答案DurronTwo-我使用两遍答案(由@JohnnyHeggheim)进行了修改,并进行了一些其他轻微修改。

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

 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

这里是一个示例运行时,如果它只在前一百万个longs上运行:

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

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

vm: java
trial: 0

正如你所看到的,DurronTwo在大输入方面做得更好,因为它经常使用魔术,但与第一个算法和Math.sqrt相比,它受到了打击,因为数字要小得多。同时,更简单的Durron是一个巨大的赢家,因为在前100万个数字中,它不必多次除以4。

这是Durron:

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

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

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

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

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

还有DurronTwo

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

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

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

还有我的基准线束:(需要谷歌卡尺0.1-rc5)

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

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


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

            return trues;   
        }

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

            return trues;   
        }

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

            return trues;   
        }
    }

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

更新:我做了一个新的算法,在某些情况下更快,在其他情况下更慢,我根据不同的输入获得了不同的基准。如果我们计算模0xFFFFFF=3 x 3 x 5 x 7 x 13 x 17 x 241,我们可以消除97.82%的非平方数。这可以(某种程度上)在一行中完成,有5个按位操作:

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

我将其加载到布尔数组中,如下所示:

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

示例运行时。在我参加的每一次测试中,它都击败了德隆(第一版)。

 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

其他回答

这是我能想到的最快的Java实现,使用了本线程中其他人建议的技术组合。

Mod-256测试不精确的mod-3465测试(避免以某些误报为代价的整数除法)浮点平方根,舍入并与输入值比较

我也尝试了这些修改,但它们对性能没有帮助:

附加mod-255测试将输入值除以4的幂快速逆平方根(要处理高N值,需要3次迭代,足以使其比硬件平方根函数慢。)

public class SquareTester {

    public static boolean isPerfectSquare(long n) {
        if (n < 0) {
            return false;
        } else {
            switch ((byte) n) {
            case -128: case -127: case -124: case -119: case -112:
            case -111: case -103: case  -95: case  -92: case  -87:
            case  -79: case  -71: case  -64: case  -63: case  -60:
            case  -55: case  -47: case  -39: case  -31: case  -28:
            case  -23: case  -15: case   -7: case    0: case    1:
            case    4: case    9: case   16: case   17: case   25:
            case   33: case   36: case   41: case   49: case   57:
            case   64: case   65: case   68: case   73: case   81:
            case   89: case   97: case  100: case  105: case  113:
            case  121:
                long i = (n * INV3465) >>> 52;
                if (! good3465[(int) i]) {
                    return false;
                } else {
                    long r = round(Math.sqrt(n));
                    return r*r == n; 
                }
            default:
                return false;
            }
        }
    }

    private static int round(double x) {
        return (int) Double.doubleToRawLongBits(x + (double) (1L << 52));
    }

    /** 3465<sup>-1</sup> modulo 2<sup>64</sup> */
    private static final long INV3465 = 0x8ffed161732e78b9L;

    private static final boolean[] good3465 =
        new boolean[0x1000];

    static {
        for (int r = 0; r < 3465; ++ r) {
            int i = (int) ((r * r * INV3465) >>> 52);
            good3465[i] = good3465[i+1] = true;
        }
    }

}

我参加聚会已经很晚了,但我希望能提供一个更好的答案;更短,(假设我的基准是正确的)也更快。

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项表,因此没有数组访问成本(间接和边界检查)。对于均匀随机的长,有81.25%的概率在这里结束。

第二个测试捕获因式分解中奇数为2的所有数字。Long.numberOfTrailingZeros方法非常快,因为它被JIT编译成一条i86指令。

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

最后的测试是双倍算术。由于double只有53位尾数,从long到double的转换包括大值的舍入。尽管如此,测试是正确的(除非证明是错误的)。

试图结合mod255的想法并不成功。

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

然后我记得,在Quake源代码中,有一个函数围绕着“网络”旋转:

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

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

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

它基本上使用牛顿近似函数(记不清确切的名字)计算平方根。

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

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

我最初在以下位置找到它:http://www.codemaestro.com/reviews/9

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

您可以通过链接了解更多的工作原理,但如果您不太在意,那么这大概是我在阅读博客和参加数值分析课程时所记得的:

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

在结果上迭代函数的次数越多,逼近函数给出的值就越精确。在Quake的案例中,一次迭代“足够好”,但如果不是为了你。。。然后您可以添加所需的迭代次数。

这应该更快,因为它减少了在简单平方根中执行的除法运算的数量(实际上是一个*0.5F乘法运算),并用一些固定数量的乘法运算代替。

这个问题让我很疑惑,所以我做了一些简单的编码,我在这里介绍它,因为我觉得它很有趣,很相关,但我不知道它有多有用。有一个简单的算法

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

用于计算平方根,但它用于小数。我想知道,如果我只是用整数数学编码相同的算法,会发生什么。它甚至会汇聚到正确的答案上吗?我不知道,所以我写了一个程序。。。

#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <math.h>

_Bool isperfectsquare(uint64_t x, uint64_t *isqrtx) {
  // NOTE: isqrtx approximate for non-squares. (benchmarked at 162ns 3GHz i5)
  uint32_t i;
  uint64_t ai;
  ai = 1 + ((x & 0xffff000000000000) >> 32) + ((x & 0xffff00000000) >> 24) + ((x & 0xffff0000) >> 16);
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = (ai + x/ai)/2;
  ai = ai & 0xffffffff;
  if (isqrtx != NULL) isqrtx[0] = ai;
  return ai*ai == x;
}

void main() {

  uint64_t x, isqrtx;
  uint64_t i;
  for (i=1; i<0x100000000; i++) {
    if (!isperfectsquare(i*i, &isqrtx)) {
      printf("Failed at %li", i);
      exit(1);
    }
  }
  printf("All OK.\n");
} 

因此,事实证明,该公式的12次迭代足以为所有64位无符号长整数(完美平方)提供正确的结果,当然,非平方将返回false。

simon@simon-Inspiron-N5040:~$ time ./isqrt.bin 
All OK.

real    11m37.096s
user    11m35.053s
sys 0m0.272s

因此697s/2^32约为162纳秒。实际上,该函数对于所有输入都具有相同的运行时。讨论中其他地方详细介绍的一些措施可以通过检查最后四位等来加快非正方形的速度。希望有人像我一样觉得这很有趣。

整数问题需要整数解。因此

对(非负)整数进行二进制搜索,以找到最大的整数t,使t**2<=n。然后测试r**2=n是否精确。这需要时间O(log n)。

如果你不知道如何对正整数进行二进制搜索,因为集合是无界的,这很容易。首先计算二次幂的递增函数f(高于f(t)=t**2-n)。当你看到它变为正值时,你已经找到了一个上限。然后可以进行标准的二进制搜索。