Round y = x * x до ближайшего

Question

Округление результата деления на ближайшее целое число равно pretty simple. Но я пытаюсь округлить результат деления, так что последующая операция даст наилучшее приближение. Это лучше всего объяснить простой функции:

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ⁿ?

Редактировать: по многочисленным просьбам Я беру на себя смелость «переводить» вопрос в чистых арифметических терминах (ни одна из этих бит-сдвигающих вещей ...).
Примечание: все деления после этого являются целыми делениями, например 13/3 будут 4.
Проблема: мы не можем вычислить x^2, потому что x большой, поэтому вместо этого мы хотели бы вычислить (x^2)/(2^N).
Для этого вычислим
x_div = X/SQRT (2^N)
, который мы тогда площадь:
у = x_div * x_div

Однако этот результат, как правило, не хватает точного значения (x^2)/(2^N) и ОП предлагает добавить 0,5 * sqrt (2^N) или, может быть, 0,375 * sqrt (2^N), чтобы лучше приблизиться к результату ...

Как ответ Oli Charlesworth предлагает гораздо лучший способ добраться до фактическое значение, считая x^2 as (x_hi + x_lo)^2.

Answer 1

Принимая во внимание, что усечение x_div приведет к ошибке ошибки не более 1, с x_div*x_div ошибка может быть до 1<<(half_digits+2).

Чтобы понять, почему, считают, что мы можем выразить эту квадратуру следующим образом (с использованием long multplication):

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 

Как видим, условия низкого порядка влияют на несколько битов в конечном итоге.

Однако каждое из этих терминов должно вписываться в исходный тип без переполнения. Поэтому при соответствующем переключении и маскировании мы можем получить желаемый результат.

Конечно, компилятор/аппаратное обеспечение делает это все для вас, если вы используете более крупный тип, поэтому вы должны просто сделать это, если у вас есть опция.

Answer 2

Используйте TYPE int для "y". Думаю, это решило бы цель.