将sqrt(n)与有理数p/q进行比较

Question

您将获得一个整数n和一个有理数p/q（其中p和q是整数）。

如何比较sqrt(n)和p/q？

解决方案1：sqrt(n) <= (double)p / q。这应该可以工作，但调用sqrt比仅使用乘法/除法慢。

解决方案2：(double)n * q * q <= p * p。更好，但我认为因为我们使用浮点数，如果p/q非常接近sqrt(n)，我们可能会得到错误的答案。此外，它需要将整数转换为浮点数，这比仅使用整数稍微慢一些。

解决方案3：n*q*q <= p*p。更好，但如果p和q变得很大，就会遇到问题，因为会发生溢出（通常，如果在使用64位整数时p或q >= 2^32）。

解决方案4：使用具有无限制整数的bignum库/编程语言的解决方案3。

解决方案5：(q/p)*n <= p/q。成功避免了任何溢出问题，但我不确定在所有情况下都正确，因为涉及整数除法...

所以...我愿意采用解决方案2或4，但我想知道是否有聪明的技巧来解决这个问题，或者是否有证明（或反例）表明解决方案5有效（或无效）。

- R2B2

---

浮点数并不一定比整数慢。请参见https://dev59.com/PnE85IYBdhLWcg3w6oB0 - llllllllll

反例为5：n = 4; p = 2; q = 1; 预期结果：sqrt(n) == p/q; 实际结果：(1 / 2) * 4 使用整数除法得到的结果是 4，而 2/1 的结果是 2，因此 4 <= 2 是错误的。 - Niet the Dark Absol

你认为一个解决方案比另一个解决方案更好的标准是什么？或者说有哪些标准？正如你的问题所暗示的，你是否只关注于32位或64位可表示的整数？ - High Performance Mark

@HighPerformanceMark 标准应该是效率和优雅。关于整数：在我的设置中，只有64位才真正出现问题。如果我有32位或更少位的整数，我可以使用64位来进行计算，而且不会有溢出问题。如果我有超过64位的整数，再次溢出也不会成为问题（因为我将使用bignum库）。但问题仍然可以转换为任何大小：当您的输入具有_n_位时，但您不能在计算中使用超过_2*n-1_位。 - R2B2

没有关于n、p、q的边界的先前知识，答案是4，使用bignum。 - aka.nice

Answer 1

如我所评论的，一个简单而优雅的解决方案是使用大数，特别是如果内置或在选择的语言中容易获得。它将在n，p，q上不受限制的工作。

当以下条件成立时，我将开发基于IEEE浮点的备用解决方案：

n，p，q都可以用给定的浮点精度准确表示（例如，在单个或双精度IEEE 754中为24或53位）

可用融合乘加。

我将注明f为浮点类型，f（x）为值x转换为f，可能已舍入到最接近的浮点数，偶数绑定。

fsqrt(x)将表示精确平方根的浮点近似值。

让f x = fsqrt（f（n）），f y = f(p) / f(q)。

根据IEEE 754属性，x和y都是精确结果的最近浮点数，并且从我们的初步条件中n = f(n)，p = f(p)，q = f(q) 。

因此，如果x < y，则问题得到解决sqrt（n）< p / q。

如果x > y，则也解决了这个问题sqrt（n）> p / q。

否则，如果x == y，则我们不能立即判断......

让我们注释残留f r = fma（x，x，-n）和fs = fma（y，q，-p）。

我们有r = x * x - n和s = y * q - p。因此，s / q = y-p / q（精确操作，而不是浮点数）。

现在，我们可以比较残余误差。 (p / q) ^ 2 = y ^ 2-2 * y * s / q +（s / q）^ 2。它与n = x ^ 2-r相比如何？

n-（p / q）^ 2 = 2 * y * s / q-r-（s / q）^ 2。

因此，我们有一个一阶的差值近似f d = 2 * y * s / f(q) - r。所以这里是一个类似于C的原型：

int sqrt_compare(i n,i p,i q)
/* answer -1 if sqrt(n)<p/q, 0 if sqrt(n)==p/q, +1 if sqrt(n)>p/q */
/* n,p,q are presumed representable in f exactly */
{
  f x=sqrt((f) n);
  f y=(f) p / (f) q;
  if(x<y) return -1;
  if(x>y) return +1;
  f r=fma(x,x,-(f) n);
  f s=fma(y,(f) q,-(f) p);
  f d=y*s/(f) q - r;
  if(d<0) return -1;
  if(d>0) return +1;
  if(r==0 && s==0) return 0; /* both exact and equal */
  return -1; /* due to 2nd order */
}

如您所见，它相对较短，应该是高效的，但很难解读，因此至少从这个角度来看，我不会认为这个解决方案比普通bignum更好。

- aka.nice

Answer 2

你可以考虑使用两倍大小的整数，选择解决方案3。

n * uint2n_t{q} * q <= uint2n_t{p} * p

如果 n * q * q 溢出，那么此处会溢出，但在这种情况下，您无论如何都会返回 false。

uint2n_t nqq;
bool overflow = __builtin_mul_overflow(uint2n_t{n} * q, q, &nqq);
(!overflow) && (uint2n_t{n} * q * q <= uint2n_t{p} * p);

- Veedrac