Rounded division by power of 2

I'm implementing a quantization algorithm from a textbook. I'm at a point where things pretty much work, except I get off-by-one errors when round开发者_开发百科ing. This is what the textbook has to say about that:

Rounded division by 2^p may be carried out by adding an offset and right-shifting by p bit positions

Now, I get the bit about the right shift, but what offset are they talking about?

Here's my sample code:

def scale(x, power2=16):
    if x < 0:
        return -((-x) >> power2)
    else:
        return x >> power2
def main():
    inp = [ 12595827, -330706, 196605, -387168, -274244, 377496, -241980, 
            -545272,  -196605, 24198,   196605,  193584, 104858, 424683,
            -40330,     41944 ]
    expect = [ 192, -5, 3, -6, -4, 5, -3, -8, -3, 0, 3, 3, 1, 6, 0, 0 ]
    actual = map(scale, inp)
    for i in range(len(expect)):
        if actual[i] == expect[i]:
            continue
        print 'inp: % 8d expected: % 3d actual: % 3d err: %d' % (inp[i], 
                expect[i], actual[i], expect[i] - actual[i])
if __name__ == '__main__':
    main()

I'm checking for negative input as bit shifting a negative integer appears to be implementation-dependent.

My output:

inp:   196605 expected:   3 actual:   2 err: 1
inp:  -387168 expected:  -6 actual:  -5 err: -1
inp:  -196605 expected:  -3 actual:  -2 err: -1
inp:   196605 expected:   3 actual:   2 err: 1
inp:   193584 expected:   3 actual:   2 err: 1

What is the offset that is mentioned in the textbook, and how can I use it to get rid of this error?

The shift will truncate. The shift is a binary operator operating. I'm using square brackets to denote the base here:

196605[10] = 101111111111111111[2]
101111111111111111[2] >> 16[10] = 10[2] = 2[10]

To perform correct rounding you need to add half of your divisor before doing the shift.

101111111111111111[2] + 1000000000000000[2] >> 16[10] = 110111111111111111[2] >> 16[10] = 11[2] = 3[10]

Shifting by p gives division by 2^p rounded down (truncated).

If you want to divide by 2^p but round to closest integer, do:

shift-right by (p-1)
add 1
shift-right 1

In your code:

def scale(x, power2=16):
    if x < 0:
        return -((((-x) >> (power2-1)) + 1) >> 1)
    else:
        return ((x >> (power2-1)) + 1) >> 1

As suspected, your algorithm is not actually rounded, but truncating division. What's more is that there are errors in your textbook as well. So even if you fix your algorithm, you wouldn't get the expected results.

In order to confirm that the results are in fact flawed, you can try running your code with a correct, float-based rounded division function:

def scale(x, power2=16):
    divider = float(1<<power2)
    result = round(x/divider)
    return result

Yet we get the following errors:

inp:   377496 expected:   5 actual:   6 err: -1
inp:  -241980 expected:  -3 actual:  -4 err: 1
inp:   104858 expected:   1 actual:   2 err: -1
inp:   -40330 expected:   0 actual:  -1 err: 1
inp:    41944 expected:   0 actual:   1 err: -1

By calculating correct results for rounding division, we can confirm that these expectations, in fact, are faulty:

377496 / 65536 = 5,7601 -> should round to 6
104858 / 65536 = 1,600 -> should round to 2
-241980 / 65536 = -3,692 -> should round to -4
104858 / 65536 = 1,600 -> should round to 2
-40330 / 65536 = -0,6154 -> should round to -1
41994 / 65536 = 0,641 -> should round to 1

Thus, if rounded division is what you really want, your list of expected values should be:

expect = [ 192, -5, 3, -6, -4, 6, -4, -8, -3, 0, 3, 3, 2, 6, -1, 1 ]

The "expected" answers are just not consistent with one of the possible rounding methods (down, nearest, up) and that's obvious from the positive dividends, before we consider the complications introduced by negative dividends.

dividend  exp  float div
   24198    0   0.3692322 DN
   41944    0   0.6400146 D
  104858    1   1.6000061 D
  193584    3   2.9538574  NU
  196605    3   2.9999542  NU
  377496    5   5.7601318 D
  424683    6   6.4801483 DN
12595827  192 192.1970673 DN

So down gets 6 out of 8, nearest gets 5, and up gets only 2.

What textbook? It's "Name'n'Shame" time!

Update after further experimentation:

If you add 8192 in just before you do the truncating division by 65536, you get the "expected" results. No other power of 2 in (512, ..., 32768) has the same effect.

Describing this as adding an offset to bias rounding downwards is somewhat unfortunate.

Rewrite: The object is to round to the NEAREST integer, but introduce a bias towards zero (smaller absolute integers). Rounding to the nearest would be done by adding in 32768 before the truncating division. Using a smaller "offset" than 32768 gives the desired bias effect. If the offset is a power of 2 e.g 2**k, it can be done by: shift k bits, add 1, shift 16-k bits.