尽可能快地比较形式(a + sqrt(b））中的两个值

这是一个没有sqrt的版本，尽管我不确定它是否比只有一个sqrt的版本更快(它可能取决于值的分布)。

这是数学(如何删除两个 sqrt)：

ad = a2-a1
bd = b2-b1
a1+sqrt(b1) < a2+sqrt(b2)              // subtract a1
sqrt(b1) < ad+sqrt(b2)              // square it
b1  < ad^2+2*ad*sqrt(b2)+b2    // arrange
ad^2+bd  > -2*ad*sqrt(b2)

在这里，右侧始终为负数。如果左侧为正数，那么我们必须返回 true。

如果左侧为负数，那么我们可以对不等式进行平方：

ad^4+bd^2+2*bd*ad^2 < 4*ad^2*b2

这里要注意的关键是，如果a2>=a1+1000，则is_smaller总是返回true(因为sqrt(b1)的最大值是 1000)。如果a2<=a1+1000，则ad是一个很小的数字，因此ad^4将始终适合 64 位(不需要 128 位算术)。代码如下：

bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
int ad = a2 - a1;
if (ad>1000) {
return true;
}
int bd = b2 - b1;
if (ad*ad+bd>0) {
return true;
}
int ad2 = ad*ad;
return (long long int)ad2*ad2 + (long long int)bd*bd + 2ll*bd*ad2 < 4ll*ad2*b2;
}

编辑：正如Peter Cordes所注意到的，第一个if是不必要的，因为第二个if处理它，所以代码变得更小更快：

bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
int ad = a2 - a1;
int bd = b2 - b1;
if ((long long int)ad*ad+bd>0) {
return true;
}
int ad2 = ad*ad;
return (long long int)ad2*ad2 + (long long int)bd*bd + 2ll*bd*ad2 < 4ll*ad2*b2;
}

我很累，可能犯了一个错误;但我敢肯定，如果我这样做了，有人会指出来。

bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
a_diff = a1-a2;   // May be negative
if(a_diff < 0) {
if(b1 < b2) {
return true;
}
temp = a_diff+sqrt(b1);
if(temp < 0) {
return true;
}
return temp*temp < b2;
} else {
if(b1 >= b2) {
return false;
}
}
//  return a_diff+sqrt(b1) < sqrt(b2);
temp = a_diff+sqrt(b1);
return temp*temp < b2;
}

如果您知道a1 < a2那么它可能会变成：

bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
a_diff = a2-a1;    // Will be positive
if(b1 > b2) {
return false;
}
if(b1 >= a_diff*a_diff) {
return false;
}
temp = a_diff+sqrt(b2);
return b1 < temp*temp;
}

还有用于计算整数 sqrt 的牛顿方法，如下所述 另一种方法是不计算平方根，而是通过二叉搜索搜索 floor(sqrt(n)) ..."只有"1000 个小于 10^6 的全平方数。 这可能性能不佳，但将是一种有趣的方法。我没有测量过这些，但这里有一些例子：

#include <iostream>
#include <array>
#include <algorithm>        // std::lower_bound
#include <cassert>          

bool is_smaller_sqrt(unsigned a1, unsigned b1, unsigned a2, unsigned b2)
{
return a1 + sqrt(b1) < a2 + sqrt(b2);
}
static std::array<int, 1001> squares;
template <typename C>
void squares_init(C& c)
{
for (int i = 0; i < c.size(); ++i)
c[i] = i*i;
}
inline bool greater(const int& l, const int& r)
{
return r < l;
}
inline bool is_smaller_bsearch(unsigned a1, unsigned b1, unsigned a2, unsigned b2)
{
// return a1 + sqrt(b1) < a2 + sqrt(b2)
// find floor(sqrt(b1)) - binary search withing 1000 elems
auto it_b1 = std::lower_bound(crbegin(squares), crend(squares), b1, greater).base();
// find floor(sqrt(b2)) - binary search withing 1000 elems
auto it_b2 = std::lower_bound(crbegin(squares), crend(squares), b2, greater).base();
return (a2 - a1) > (it_b1 - it_b2);
}
unsigned int sqrt32(unsigned long n)
{
unsigned int c = 0x8000;
unsigned int g = 0x8000;
for (;;) {
if (g*g > n) {
g ^= c;
}
c >>= 1;
if (c == 0) {
return g;
}
g |= c;
}
}
bool is_smaller_sqrt32(unsigned a1, unsigned b1, unsigned a2, unsigned b2)
{
return a1 + sqrt32(b1) < a2 + sqrt32(b2);
}
int main()
{
squares_init(squares);
// now can use is_smaller
assert(is_smaller_sqrt(1, 4, 3, 1) == is_smaller_sqrt32(1, 4, 3, 1));
assert(is_smaller_sqrt(1, 2, 3, 3) == is_smaller_sqrt32(1, 2, 3, 3));
assert(is_smaller_sqrt(1000, 4, 1001, 1) == is_smaller_sqrt32(1000, 4, 1001, 1));
assert(is_smaller_sqrt(1, 300, 3, 200) == is_smaller_sqrt32(1, 300, 3, 200));
}

我不确定代数操作是否与整数相结合 算术，必然导致最快的解决方案。您将需要 在这种情况下，许多标量乘法(不是很快)，和/或 分支预测可能会失败，这可能会降低性能。 显然，您必须进行基准测试，以查看在特定情况下哪种解决方案最快。

一种方法制作sqrt更快的方法是将-fno-math-errno选项添加到GCC或CLANG。 在这种情况下，编译器不必检查负输入。 对于 icc，这是默认设置。

通过使用矢量化可以提高性能sqrt指令sqrtpd，而不是标量sqrt指令sqrtsd。 Peter Cordes已经证明，clang能够自动矢量化这段代码， 这样它就会生成这个sqrtpd.

但是，自动矢量化的成功程度在很大程度上取决于正确的编译器设置 以及使用的编译器(CLANG，GCC，ICC等)。对于-march=nehalem或更老的，叮当不会矢量化。

使用以下内部代码可以获得更可靠的矢量化结果，请参见下文。 对于可移植性，我们只假设支持 SSE2，这是 x86-64 基线。

/* gcc -m64 -O3 -fno-math-errno smaller.c                      */
/* Adding e.g. -march=nehalem or -march=skylake might further  */
/* improve the generated code                                  */
/* Note that SSE2 in guaranteed to exist with x86-64           */
#include<immintrin.h>
#include<math.h>
#include<stdio.h>
#include<stdint.h>
int is_smaller_v5(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
uint64_t a64    =  (((uint64_t)a2)<<32) | ((uint64_t)a1); /* Avoid too much port 5 pressure by combining 2 32 bit integers in one 64 bit integer */
uint64_t b64    =  (((uint64_t)b2)<<32) | ((uint64_t)b1); 
__m128i ax      = _mm_cvtsi64_si128(a64);         /* Move integer from gpr to xmm register                  */
__m128i bx      = _mm_cvtsi64_si128(b64);         
__m128d a       = _mm_cvtepi32_pd(ax);            /* Convert 2 integers to double                           */
__m128d b       = _mm_cvtepi32_pd(bx);            /* We don't need _mm_cvtepu32_pd since a,b < 1e6          */
__m128d sqrt_b  = _mm_sqrt_pd(b);                 /* Vectorized sqrt: compute 2 sqrt-s with 1 instruction   */
__m128d sum     = _mm_add_pd(a, sqrt_b);
__m128d sum_lo  = sum;                            /* a1 + sqrt(b1) in the lower 64 bits                     */
__m128d sum_hi  =  _mm_unpackhi_pd(sum, sum);     /* a2 + sqrt(b2) in the lower 64 bits                     */
return _mm_comilt_sd(sum_lo, sum_hi);
}

int is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
return a1+sqrt(b1) < a2+sqrt(b2);
}

int main(){
unsigned a1; unsigned b1; unsigned a2; unsigned b2;
a1 = 11; b1 = 10; a2 = 10; b2 = 10;
printf("smaller?  %i  %i n",is_smaller(a1,b1,a2,b2), is_smaller_v5(a1,b1,a2,b2));
a1 = 10; b1 = 11; a2 = 10; b2 = 10;
printf("smaller?  %i  %i n",is_smaller(a1,b1,a2,b2), is_smaller_v5(a1,b1,a2,b2));
a1 = 10; b1 = 10; a2 = 11; b2 = 10;
printf("smaller?  %i  %i n",is_smaller(a1,b1,a2,b2), is_smaller_v5(a1,b1,a2,b2));
a1 = 10; b1 = 10; a2 = 10; b2 = 11;
printf("smaller?  %i  %i n",is_smaller(a1,b1,a2,b2), is_smaller_v5(a1,b1,a2,b2));
return 0;
}

有关生成的程序集，请参阅此 Godbolt 链接。

在英特尔 Skylake 上进行的简单吞吐量测试中，使用编译器选项gcc -m64 -O3 -fno-math-errno -march=nehalem，我发现了吞吐量is_smaller_v5()比原始is_smaller()好 2.6 倍：6.8 个 CPU 周期与 18 个 CPU 周期，包括环路开销。然而，在一个(太？ 简单的延迟测试，其中输入a1, a2, b1, b2取决于前is_smaller(_v5)的结果，我没有看到任何改进。(39.7 个周期与 39 个周期)。

可能不比其他答案好，但使用了不同的想法(和大量的预分析)。

// Compute approximate integer square root of input in the range [0,10^6].
// Uses a piecewise linear approximation to sqrt() with bounded error in each piece:
//   0 <= x <= 784 : x/28
//   784 < x <= 7056 : 21 + x/112
//   7056 < x <= 28224 : 56 + x/252
//   28224 < x <= 78400 : 105 + x/448
//   78400 < x <= 176400 : 168 + x/700
//   176400 < x <= 345744 : 245 + x/1008
//   345744 < x <= 614656 : 336 + x/1372
//   614656 < x <= 1000000 : (784000+x)/1784
// It is the case that sqrt(x) - 7.9992711366390365897... <= pseudosqrt(x) <= sqrt(x).
unsigned pseudosqrt(unsigned x) {
return 
x <= 78400 ? 
x <= 7056 ?
x <= 764 ? x/28 : 21 + x/112
: x <= 28224 ? 56 + x/252 : 105 + x/448
: x <= 345744 ?
x <= 176400 ? 168 + x/700 : 245 + x/1008
: x <= 614656 ? 336 + x/1372 : (x+784000)/1784 ;
}
// known pre-conditions: a1 < a2, 
//                  0 <= b1 <= 1000000
//                  0 <= b2 <= 1000000
bool is_smaller(unsigned a1, unsigned b1, unsigned a2, unsigned b2) {
// Try three refinements:
// 1: a1 + sqrt(b1) <= a1 + 1000, 
//    so is a1 + 1000 < a2 ?  
//    Convert to a2 - a1 > 1000 .
// 2: a1 + sqrt(b1) <= a1 + pseudosqrt(b1) + 8 and
//    a2 + pseudosqrt(b2) <= a2 + sqrt(b2), 
//    so is  a1 + pseudosqrt(b1) + 8 < a2 + pseudosqrt(b2) ?
//    Convert to a2 - a1 > pseudosqrt(b1) - pseudosqrt(b2) + 8 .
// 3: Actually do the work.
//    Convert to a2 - a1 > sqrt(b1) - sqrt(b2)
// Use short circuit evaluation to stop when resolved.
unsigned ad = a2 - a1;
return (ad > 1000)
|| (ad > pseudosqrt(b1) - pseudosqrt(b2) + 8)
|| ((int) ad > (int)(sqrt(b1) - sqrt(b2)));
}

(我手边没有编译器，所以这可能包含一两个错别字。