舍入y = x*x到最近

Question

将除法结果舍入到最接近的整数非常简单.但我试图围绕除法的结果,以便后续操作将给出最佳近似值.最好的解释是一个简单的函数:

const int halfbits = std::numeric_limits<unsigned int>::digits / 2;
unsigned x = foo(); // Likely big.
unsigned x_div = x >> halfbits; // Rounded down
unsigned y = x_div * x_div; // Will fit due to division.

我可以x_div通过添加来舍入到最近1<<(halfbits-1).但由于x²不是线性函数,因此y通常不能正确舍入.有(x*x) >> (halfbits*2)没有使用更大类型的简单而更准确的计算方法？

我认为添加3<<(halfbits-3)到x_div会改进舍入,但无法证明这是最佳解决方案.此外,这可以推广为xⁿ？

编辑:按照大众的要求,我冒昧地用纯算术术语"翻译"这个问题(这些都不是C位移位的东西......).
注意:此后的所有除法都是整数除数,例如13/3就是4.
问题:我们无法计算x ^ 2,因为x很大,所以我们想要计算(x ^ 2)/(2 ^ N).
为此,我们计算
x_div = x/sqrt(2 ^ N)
,然后我们
求平方:y = x_div*x_div

然而,这个结果通常短于确切的值,(x^2)/(2^N)并且OP建议添加0.5*sqrt(2 ^ N)或者可能是0.375*sqrt(2 ^ N)以更好地近似结果...

正如Oli Charlesworth答案所暗示的那样,通过考虑x^2为(x_hi + x_lo)^ 2 ,可以更好地获得实际值.

Answer 1

而x_div截断将导致误差幅度最多为1,x_div*x_div误差可能达到1<<(half_digits+2).

要了解原因,请考虑我们可以如下表达这种平方(使用长多重):

x * x = (x_lo + x_hi) * (x_lo + x_hi)
      = x_lo^2 + x_hi^2 + 2*x_lo*x_hi

其中x_lo和x_hi分别是下半部分和上半部分x.有了一些不错的ASCII艺术,我们可以考虑这些如何排列:

 MSB     :      :      :     LSB
  +------+------+      :      :
  |    x_hi^2   |      :      :
  +-----++-----++-----+:      :
  :     | 2*x_lo*x_hi |:      :
  :     +------++-----++------+
  :      :      |    x_lo^2   |
  :      :      +------+------+
  :<----------->:      :      :
    Our result

我们可以看到,低阶项会影​​响最终结果中的多个位.

但是,这些术语中的每一个都应该适合原始类型而不会溢出.因此,通过适当的移位和屏蔽,我们可以获得所需的结果.

当然,如果您使用更大的类型,编译器/硬件会为您完成所有操作,因此如果您有选项,您应该这样做.