[ prog / sol / mona ]

prog


What are you working on?

165 2020-12-31 15:10

When trying to perform integer division of some number n by some divisor d, we may instead wish to approximate this with multiplication by some magic factor followed by integer division with a power of the base of the number system of our hardware. The last step can be performed using a right shift operation, and the multiply-shift combination is likely to be faster than performing integer division. For example our base might be 2 and our divisor 10, as when trying to produce the decimal string representing our number n. The correctness condition we impose is that we must obtain the same result as if we performed the integer division. Since the right shift operation rounds downward on positive integers, our approximation must be no less than the result of the integer division, and may only be greater by an amount that does not overflow into the result for the successor integer. As an example, we cannot allow the result of the integer division 20//10 to be approximated by anything less since downward rounding will follow, and we cannot allow the approximation for the integer division 29//10 to reach as high as 30//10. The correctness condition can be expressed as:

n    n * magic   n + 1
- <= --------- < -----
d    basepower     d  

Taking the cross product of the first inequality yields

basepower <= d * magic

The n has simplified away so this inequality can be made to hold regardless of n. We are left with a d-multiple that is no less than basepower. As an example we might take basepower to be 2^16=65536 and magic to be 6554 which gives

n // 10 <= n * 6554 // 65536

The second inequality yields

n * magic * d < basepower * (n + 1)

Since magic * d is a d-multiple that is no less than basepower, we make this explicit by writing

magic * d = basepower + over

Substitution and simplification yield

n * over < basepower

As an example we might have

6554 * 10 = 65536 + 4

which yields

n * 4 < 65536

Several conclusions can be drawn. One is that minimizing 'over' helps satisfy this inequality. This is equivalent to choosing magic * d to be the smallest d-multiple no less than basepower. Another conclusion is that if 'over' is not 0, it will be at least 1 and thus basepower will be greater than n. The main conclusion, however, is that n has not simplified away and therefore in the general case the inequality will not hold for arbitrary n. For any choice of magic factor and basepower where 'over' is not 0, we can always find a large enough n to break this inequality. It is for this reason that the process of picking an appropriate magic factor and basepower does not depend merely on the base and the divisor, it also depends on the maximum value of n upto which we expect our method to work. To go over this limit we need to find a higher basepower and magic factor to satisfy both inequaliies of the correctness condition. This strategy of iterative search can be implemented as follows:

def diim_basebits (base, target, strict):
   shift = 1
   exp   = base

   while exp < target:
      exp   *= base
      shift += 1

   if strict and (exp == target):
      exp   *= base
      shift += 1

   return (shift, exp)

def work_diim (base, div, limit, *, debug = False):
   msg = print

   if debug:
      msg ("base {:,d} div {:,d} limit {:,d}".format (
         base, div, limit))

   shift, exp = diim_basebits (base, div, False)
   more       = True
   tries      = 0

   while more:
      tries   += 1
      exptodiv = divmod (exp, div)

      if exptodiv [1]:
         factor = exptodiv [0] + 1
         over   = div - exptodiv [1]
      else:
         factor = exptodiv [0]
         over   = 0

      product = over * limit

      if debug:
         msg ("try {:,d} shift {:,d} exp {:,d}".format (
            tries, shift, exp))
         msg ("  factor {:,d} over {:,d} product {:,d}".format (
            factor, over, product))

      if product < exp:
         more = False
      else:
         exp   *= base
         shift += 1

   product     = limit * factor
   pbits, pexp = diim_basebits (base, product, True)

   if debug:
      msg ("x div {:,d} -> (x * {:,d}) >>({:,d}) {:,d}".format (
         div, factor, base, shift))
      msg ("  basebits {:,d} max {:,d}".format (
         pbits, product))

   return {
      "factor":   factor,
      "shift":    shift,
      "basebits": pbits,
   }
199


VIP:

do not edit these