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

【问题标题】：Comparing sqrt(n) with the rational p/q将 sqrt(n) 与有理 p/q 进行比较
【发布时间】：2018-07-13 19:00:40
【问题描述】：

给你一个整数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 由于溢出而变大（通常是 p 或 q em> >= 2^32 当使用 64 位整数时）。

解决方案 4：将解决方案 3 与 bignum 库/在具有未绑定整数的编程语言中结合使用。

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

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

【问题讨论】：

浮点不一定比整数慢。见stackoverflow.com/questions/2550281/…
5 的反例：n = 4; p = 2; q = 1; 预期结果：sqrt(n) == p/q; 实际结果：(1 / 2) * 4 整数除法为4，而2/1 为2，产生4 <= 2这是错误的。
您认为一种解决方案比另一种解决方案更好的标准是什么？或者标准是什么？正如您的问题所暗示的那样，您是否只关心可以 32 位或 64 位表示的整数？
@HighPerformanceMark 标准将是效率和优雅。关于整数：问题只出现在我的设置中的 64 位。如果我有 32 位或更少的整数，我可以使用 64 位来进行计算并且我没有溢出问题。如果我有超过 64 位的整数，那么溢出也不会成为问题（因为我将使用 bignum 库）。但问题仍然可以转化为任何大小：当您的输入有 n 位，但您的计算不能使用超过 2*n-1 位时。跨度>
在不知道 n,p,q 上界的情况下，答案是 4，使用 bignum

标签： performance math floating-point

【解决方案1】：

正如我所评论的，一个简单而优雅的解决方案是使用 bignum，特别是如果是内置的，或者在所选语言中很容易获得。它可以不受 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) 和f s = 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.

因此，我们得到了差值d 的近似值，一阶：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 */
}

如您所见，它相对较短，应该是有效的，但很难破译，所以至少从这个 POV 来看，我不会认为这个解决方案比琐碎的 bignum 更好。

【讨论】：

【解决方案2】：

您可以考虑使用 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);

【讨论】：