hellerve / silly

A silly implementation of fixed point arithmetic (mirror because I linked here on my blog)
1 stars 1 forks source link

silly_div : I can't understand why b is 33. I'm silly #1

Open gnahznux opened 3 years ago

gnahznux commented 3 years ago

I wonder why b is 33 instead of 32.

silly silly_div(silly x, silly y) {
  // disregard the sign
  uint64_t sign_bit = 1UL<<63;
  uint64_t rem = x & ~sign_bit;
  uint64_t div = y & ~sign_bit;

  uint64_t quo = 0UL;
  int b = 33; // 64 / 2 + 1

  while (rem && b >= 0) {
      uint64_t d = rem / div;
      rem %= div;
      quo += d << b;

      rem <<= 1;
      --b;
  }

  uint64_t res = quo >> 1;
  // add the sign
  res |= (x & sign_bit) ^ (y & sign_bit);

  return res;
}
hellerve commented 3 years ago

b is the number of bits we go through. The upper bound for this (unless rem shortcuts us) is <width> / 2 + 1, where width is the width of the number. Since we know that that width is always 64 bits, we can just hardcode it in our function.

I encourage you to read through my blog post about this to get a better sense as to why.

gnahznux commented 3 years ago

thanks for your reply !!! I have read your blog before asking the question. It's very impressive. That the upper bound is 64/2 +1 means we go through 33 bit at the most? Can you explain specific? And after while loop , quo >> 1 is weird . If b is 32 , we don't need >> 1 after while loop, am I right?

hellerve commented 3 years ago

I haven’t thought about this code in quite some time, but you seem to be correct. IIRC the whole thing has something to do with rounding (also the increment of quo that I commented out in the code). I can’t replicate the behavior we’re fixing here, however, so the only answer that I have now is: I forgot why I did what I did there (I realize that this is a disappointment, but such is life with code that you wrote a few years ago).

gnahznux commented 3 years ago

Don't feel upset about it (doge). I will tell you the answer when I figure it out.