Rounding of integer divisions (was Re: linear interpolation / Assembler / ATMega32)

From: Ville Voipio (vvoipio_at_kosh.hut.fi)
Date: 06/13/04 

Date: Sun, 13 Jun 2004 13:52:09 +0300

Paul Keinanen <keinanen@sci.fi> writes:

> If signed multiply and division are used, there should not be a
> problem. If arithmetic shift right is used to handle 2's complement
> negative numbers, the results might not be quite as expected, since
> e.g. -1 >> 4 = -1.

Well, no problem, as:

-17 >> 4 = 11101111b >> 4 = 11111110b = -2
-16 >> 4 = 11110000b >> 4 = 11111111b = -1
...
 -1 >> 4 = 11111111b >> 4 = 11111111b = -1
  0 >> 4 = 00000000b >> 4 = 00000000b = 0
...
 15 >> 4 = 00001111b >> 4 = 00000000b = 0
 16 >> 4 = 00010000b >> 4 = 00000001b = 1

So, everything is rounded down:

-2 <= x < -1 -> -2 (bin average -1.5)
-1 <= x < 0 -> -1 (bin average -0.5)
 0 <= x < 1 -> 0 (bin average 0.5)
 1 <= x < 2 -> 1 (bin average 1.5)

The results are half-step lower than they should be (they should
coincide with the bin averages). This can be corrected by
adding 8 (that is 0.5 * 2^4) to the initial value:

(-9 + 8) >> 4 = 11111111b >> 4 = -1
(-8 + 8) >> 4 = 00000000b >> 4 = 0
...
( 7 + 8) >> 4 = 00001111b >> 4 = 0
( 8 + 8) >> 4 = 00010000b >> 4 = 1

Now the limits are exactly where they shoud be:

-9 / 16 = -0.5625 -> -1
-8 / 16 = -0.5 -> 0
...
 7 / 16 = 0.4375 -> 0
 8 / 16 = 0.5 -> 1

Everything is rounded to the nearest integer, halves are rounded
up.

So, the plain ASR works fine apart from the average -0.5 bias in
the result. The bias can be compensated for in the tabulated data
or by adding the half before shifting.

The worst thing you may do is to unsignedify the number, round
it down, and then put the sign back. Because:

-17 / 16 -> -(17 >> 4) -> -1
-16 / 16 -> -(16 >> 4) -> -1
-15 / 16 -> -(15 >> 4) -> 0
...
  0 / 16 -> ( 0 >> 4) -> 0
...
 15 / 16 -> (15 >> 4) -> 0
 16 / 16 -> (16 >> 4) -> 1

The function becomes dead around zero! The result is zero over
two units, and everything linear becomes non-linear. Cross-over
distortion sets in...

Unfortunately, this is exactly what C does when you write:

  int z, x, y;

  z = x / y;

There is a technical reason for this; division is easier to do,
if you do not need to write different branches for different
signs. C rounds towards zero. Not up, not down, but towards
zero. (This may even be a desired behaviour in some cases,
but with measurements it is a very bad thing.)

One simple way around this is to use biased unsigned divisions:

  // bias by 15:
  z = (x + 15 * y) / y - 15;

If x + 15 * y always >= 0. This is rather easy with constant divisors.
With variable divisors finding a suitable constant is very difficult.
Pick a too low one, and the sum will underflow. Pick a too high one,
and the sum will overflow.
  
The case with the possibility of a negative divisor is even
more complicated. Then the only alternative is to branch the
calculations by the signs before doing anything;

  int x, y, z;
  unsigned int ax, ay, az;
  bool sign;

  // find out the sign of the result and the absolute values of x, y
  if (x < 0)
    {
    ax = -x;
    sign = false;
    }
  else
    {
    ax = x;
    sign = true;
    }

  if (y < 0)
    {
    ay = -y;
    sign = ~sign;
    }
  else
    ay = y;

  // if the sign of the result is negative, round away from zero
  if (sign == false)
    {
    az = (ax + ay - 1) / ay;
    z = -az;
    }
  else
    {
    az = ax / ay;
    z = az;
    }

Yes, there are a lot of extra variables which may be avoided
by using type casts. However, I think it is better avoid the
type casts and let the compiler get rid of the extra variables.

So, twenty-something lines of code to achive a very simple
thing. The good thing is that this does not give a very large
performance hit, as the signed division algorithm does this
anyway.

To make the algorithm round right (towards nearest, halves up)
requires only the following changes:

negative branch:
  az = (ax + ay - 1 - ay/2) / ay

positive branch:
  az = (ax + ay/2) / ay

Even though it is very tempting to say:

  ay - 1 - ay/2 == ay/2 - 1,

don't. Because it isn't. For, e.g., ay = 17:

  17 - 1 - 8 == 8 - 1
           8 == 7

Oops.

I know this is somewhat complicated. The effect of false rounding
is usually quite small. I personally crashed into this well-known
(generally, that is, not well-known by me) problem when I was
looking at a F-transform of some real-world data. There were some
harmonics there shouldn't've been. I scratched through
the skin of my head before finding out the problem.

If off-by-one is bad, off-by-half is not much nicer, either.

- Ville 

-- 
Ville Voipio, Dr.Tech., M.Sc. (EE)

Relevant Pages

  • Re: Division/math bug in perl?
    ... I thought converting from float to int is supposed ... > to give the integer part, which is -3, and not round towards zero, ... The docs for int() specifically say that it truncates towards zero. ...
    (comp.lang.perl.misc)
  • Re: int(2.1) does not work while int(float(2.1)) does
    ... >> Because intmeans round toward zero. ... >> int out of this string. ...
    (comp.lang.python)
  • Re: int(2.1) does not work while int(float(2.1)) does
    ... > Joe Mason wrote: ... > int out of this string. ... So why can't intmean round toward zero? ...
    (comp.lang.python)
  • [PATCH] 2/2 Prezeroing large blocks of pages during allocation Version 4
    ... Have USERZERO and KERNZERO for different types of zero pages to avoid ... This is a patch that makes a step towards merging the modified allocator ... static inline void inc_reserve_count(struct zone *zone, ... +static inline void prep_zero_page(struct page *page, int order, int gfp_flags) ...
    (Linux-Kernel)
  • Re: edit and continue
    ... > Are you using the same sourcesafe I'm using. ... > safe crashing vb6. ... Zero random rollbacks (I can see why this would be a bad ... > int", the professional works on code daily and would much rather type "int ...
    (microsoft.public.vb.general.discussion)