Re: [PATCH 2/5] dynlink.c: compute modulus via magic multiplication

[Alexander Monakov <amonakov@ispras.ru>]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 1719 bytes --]

> > How can I do the last step, 'x-v*div' without it?
> 
> Ah yes. Would it be preferable to have that in struct udiv then,
> though? Then the caller never has to care about the divisor, and in
> the case where umod doesn't get inlined, the call will be more
> efficient (fewer args).

I doubt that; the caller needs that value ('nbuckets') anyway.

[...]
> Using the post-mul add rather than the saturated increment would make
> this work for 0xffffffff too.

post-mul add is problematic on 32-bit
 
> Another obvious solution is not using the +32 offset so that the right
> shift can just be 31, but that pessimizes the code on 32-bit archs
> quite a bit.

I agree that a special-case for a power-of-two divisor will fire too rarely,
so figuring a replacement out would be nice.

If you (or anyone) want to play with ideas, I'm attaching my test driver for
magic division that tests 2^32-1 inputs for a divisor given on the command
line (the main loop looks odd, but my goal was to have gcc vectorize it).

> > p->s1 check is for this reason: shift are relatively costly, and s1 is rarely
> > non-zero, so try to skip the shift if possible; in the rare case the pre-shift
> > is non-zero, the check allows to skip the saturating increment operation.
> 
> Shifts are costly? Are we talking about P4-era junk? ;-)

I primarily had modern Intel cores in mind where as far as I can see shifts by
%ecx have latency 2.  But even ignoring that, there's the second part of my
argument.

Anyway it's not trivial to measure an impact on the same "modern Intel cores",
so if anybody can say how that looks on a different platform, I'd appreciate
that.  Removing my if-else there is obviously beneficial for code size.

Alexander

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: Type: TEXT/x-c; name=udiv.c, Size: 1923 bytes --]

#include <stdlib.h>
static void precalc(unsigned div, unsigned *mulp, int *prep, int *postp, int *incp)
{
  int bits = 0;
again:
  *prep = bits;
  unsigned pow = 1u<<31, quo = pow / div, rem = pow % div;

  int log=0, exp;
  unsigned tmp=div;
  do log++; while (tmp >>= 1);

  unsigned down_mul=0;
  for (exp=0, pow=1u<<bits; ; exp++, pow<<=1) {
    int adj = rem >= div - rem;
    quo *= 2; rem *= 2;
    if (adj) {quo++; rem -= div; }

    if (exp >= log-bits || div-rem <= pow)
      break;

    if (!down_mul && rem <= pow) {
      down_mul = quo; *postp = exp;
    }
  }
  if (exp < log) {
    *mulp = quo+1;
    *postp = exp;
    *incp = 0;
  } else if (div & 1) {
    *mulp = down_mul;
    *incp = 1;
  } else {
    do bits++; while (!((div >>= 1) & 1));
    goto again;
  }
}
int main(int argc, char *argv[])
{
  unsigned div = atoi(argv[1]);
  if (!(div&(div-1)))
    return 0;

  unsigned mul;
  int pre, post, inc;
  precalc(div, &mul, &pre, &post, &inc);
  unsigned s32 = 32;
  if (sizeof(long) == 8) {
    s32 += post;
    post = 0;
  }
  unsigned rem, quo, val;
  char bad=0;
#define loop(pre, inc, post) ({    \
  for (val=~0u, quo=val/div, rem=val%div+1; quo+1; quo--, rem=div)      \
    for (; rem; rem--, val--) {  \
      unsigned v = val;            \
      v >>= pre;                   \
      if (v+inc) v+=inc;           \
      v = (1ull * v * mul) >> s32; \
      v >>= post;                  \
      bad |= (v != quo);      \
    } })
  if (post)
    if (inc)
      if (pre)
        return 2;
      else
        loop(0, inc, post);
    else
      if (pre)
        loop(pre, 0, post);
      else
        loop(0, 0, post);
  else
    if (inc)
      if (pre)
        return 2;
      else
        loop(0, inc, 0);
    else
      if (pre)
        loop(pre, 0, 0);
      else
        loop(0, 0, 0);
  return bad;
}