Little help in understanding expmed.c::choose_multiplier()

[Denis Vlasenko]

I noticed that gcc's div-by-constant optimization is a bit
suboptimal and want to improve it. I submitted it to
bugzilla:
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=28417

but I think "big guys" have no time for such a low-impact thing.
I want to do it myself.

However, I need a little helt in understanding
what "precision" parameter is for in this function.
The comments in the source are unfortunately are too terse for me.
My current understanding of the functions internals are below.
Please read on or jump to [*] marker below if you are impatient.

So the routine tries to find such m so that:

x div d == x * m >> POST_SHIFT

and m is at max N+1 bit.

As far as I can see, at first routine picks a largest possible
POST_SHIFT:

lgup = ceil(log2(d)); POST_SHIFT = N + lgup;

Then it finds mlow and mhigh so that:

x * mlow >> POST_SHIFT
    <= x * m >> POST_SHIFT
	<= x * mhigh >> POST_SHIFT

Using v >> POST_SHIFT == v div (1<<(N+lgup)) we can rewrite it as
["div" denotes an integer division, unlike real division "/"):

x * mlow div (1<<(N+lgup))
    <= x * m div (1<<(N+lgup))
	<= x * mhigh div (1<<(N+lgup))

"Precise" m is equal to (1<<(N+lgup))/d (a real value, not integer).
Substitute it into middle expr. Then obviously
mlow = floor(m) = (1<<(N+lgup)) div d.

For mhigh: for arbitrary unsigned a,b: (a+b-1) div b >= a/b.
So we want to add a value bigger than d-1 to the dividend.
Let's add 1<<N:

x * ((1<<(N+lgup)) div d) div (1<<(N+lgup))
    <=  x * (1<<(N+lgup))/d div (1<<(N+lgup))
	<= x * (((1<<(N+lgup)) + 1<<N) div d) div (1<<(N+lgup))

Then routine checks whether POST_SHIFT can be lowered.

The inefficiency here is that routine does not take into account
_at which value of x_  x * m >> POST_SHIFT will fail.
It will always fail at some x = n*d - 1, and for large d
these values are rather scarce -> we may be lucky and not hit
such a value!

An example:

choose_multiplier(d=1577682821,n=32,precision=32) returns
        *post_shift_ptr=31,multiplier=5846151023
whereas optimal one is
        *post_shift_ptr=27,multiplier=365384439

and it is correct wrt algorithm: mlow=5846151022 < m < mhigh=5846151023.

But the catch is that 5846151024 _too_ will "by chance of d being large"
work for any 32-bit x, and 5846151024/16 = 365384439 will also work
(because of all those zeroes in low-order bits)!

IOW: "mlow < m < mhigh" algorithm is not optimal.
It misses potentially better values.

[*]
I have a better alrorithm. See attachment if you are curious.

I'd like to put it instead of current one, but I don't
understand the role of the "precision" parameter in the current code:
Currently comment says:

   Choose a minimal N + 1 bit approximation to 1/D that can be used to
   replace division by D, and put the least significant N bits of the result
   in *MULTIPLIER_PTR and return the most significant bit.

Ok, this one I understand.

   The width of operations is N (should be <= HOST_BITS_PER_WIDE_INT),

Width of _which_ ops? x * m is N * N = N bits? or N * N = 2N? Or what??

   the needed precision is in PRECISION (should be <= N).

The "needed" precision? Does this mean that that x*m>>shift is allowed
to deviate by +/- 1<<(N-PRECISION) - 1 from true result x/d? Or what??

I'd prefer much more verbose comments here.... HEEEEEELP :)
Maybe an example of the call with N and PRECISION which
are not equal to eaqch other and not equal to 32 will be helpful.

Gory details of choose_multiplier(d=1577682821,n=32,precision=32)
are below:

  lgup = ceil_log2 (d); //// 31
  pow = n + lgup; //// 63
  pow2 = n + lgup - precision; //// 31

  /* mlow = 2^(N + lgup)/d */
 if (pow >= HOST_BITS_PER_WIDE_INT) //// yes
    {
      nh = (HOST_WIDE_INT) 1 << (pow - HOST_BITS_PER_WIDE_INT); //// 1<<31
      nl = 0;
    }
  else
    {
      nh = 0;
      nl = (unsigned HOST_WIDE_INT) 1 << pow;
    }
//// 1<<63 / d: mlow=5846151022
  div_and_round_double (TRUNC_DIV_EXPR, 1, nl, nh, d, (HOST_WIDE_INT) 0,
                        &mlow_lo, &mlow_hi, &dummy1, &dummy2);

  /* mhigh = (2^(N + lgup) + 2^N + lgup - precision)/d */
  if (pow2 >= HOST_BITS_PER_WIDE_INT) //// no
    nh |= (HOST_WIDE_INT) 1 << (pow2 - HOST_BITS_PER_WIDE_INT);
  else
    nl |= (unsigned HOST_WIDE_INT) 1 << pow2; //// 1<<31
//// (1<<63 + 1<<31) / d: mhigh=5846151023 (5846151023.661466 in fp)
  div_and_round_double (TRUNC_DIV_EXPR, 1, nl, nh, d, (HOST_WIDE_INT) 0,
                        &mhigh_lo, &mhigh_hi, &dummy1, &dummy2);

--
vda

-------------- next part --------------
A non-text attachment was scrubbed...
Name: find_fast_div.c
Type: text/x-csrc
Size: 3830 bytes
Desc: not available
URL: <https://gcc.gnu.org/pipermail/gcc/attachments/20060725/688723b9/attachment.bin>