檢測完全平方數

黃志斌發表於2017-07-05

尤拉計劃的一些題目需要檢測完全平方數,最簡單的方法是使用以下 C 語言程式:

int isSquare(long n) { long x = (long)sqrt(n); return x * x == n; }

但是,也有更快的方法,即以下的 C 語言程式(exactSquareRoot.c):

#include <math.h>

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

// return: < 0 if x is not a square, r if r*r == x and r >= 0
int exactSquareRoot(long x)
{
  // Quickfail
  if (x < 0 || (x & 2) || (x & 7) == 5 || (x & 11) == 8) return -1;
  if (x == 0) return 0;
  // Check mod 255 = 3 * 5 * 17, for fun
  long y = x;
  y = (y & 4294967295L) + (y >> 32);
  y = (y & 65535) + (y >> 16);
  y = (y & 255) + ((y >> 8) & 255) + (y >> 16);
  if (bad255[y]) return -1;
  // Divide out powers of 4 using binary search
  y = x;
  if ((y & 4294967295L) == 0) y >>= 32;
  if ((y & 65535) == 0) y >>= 16;
  if ((y & 255) == 0) y >>= 8;
  if ((y & 15) == 0) y >>= 4;
  if ((y & 3) == 0) y >>= 2;
  if ((y & 7) != 1) return -1;
  int r = (int)sqrt(x);
  return (r * (long)r == x) ? r : -1;
}

我們來比較一下兩者的速度(a.c):

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

extern int exactSquareRoot(long x);

int isSquare1(long n) { long x = (long)sqrt(n); return x * x == n; }
int isSquare2(long n) { return exactSquareRoot(n) >= 0; }

int count(long a, long n, int (*f)(long))
{
  int z = 0;
  for (long i = 0; i < n; i++) if (f(a + i)) z++;
  return z;
}

int main()
{
  long a = 123456789; long n = 4987654321;
  printf("%ld: %ld\n", a, n); fflush(stdout);
  printf("%d\n", count(a, n, isSquare1)); fflush(stdout);
  printf("%d\n", count(a, n, isSquare2)); fflush(stdout);
}

執行結果:

$ clang -O3 a.c exactSquareRoot.c -lm
$ ./a.out | gnomon -i

  40.4646s | 123456789: 4987654321
  24.9703s | 60381
   0.0069s | 60381
           |
     Total | 65.4442s

這個比較程式統計從 123456789 開始的 4987654321 個整數,計算其中的完全平方數的個數,得到的結果是有 60381 個。前者用時 40.46 秒,後者用時 24.97 秒。

注:

參考資料

  1. Wikipedia: Square number
  2. Stack Overflow: Fastest way to determine if an integer's square root is an integer
  3. 圖靈社群:測量程式的執行時間

相關文章