(42, all flavours) Integer Division - how?

[Albert Chan]

(12-14-2020 10:10 AM)Werner Wrote:  At the risk of being proven wrong again..
as long as b+1 is exact, (and q is not too large), the code for idiv3 works for fractional a and b, too, no?

function idiv3(a,b) -- assumed b > 0
local q, c = floor(a/b), b+1
a, b = a%b - a%c, q%c -- a-b = [2-2c, c-2]
if a==b or a==b-c then return q end
return q-1
end

Yes, fractions still work!

a = (a//b)*b + (a%b) = q*(c-1) + (a%b)

  → a%b - a%c - q%c ≡ 0 (mod c)

If c=b+1 is exact, and q is correct, above identity holds.

And, optimization for (mod c) by 2 comparisons also hold.
Checking range of each terms (note, I use opened end intervals here):

a1 = (a%b) - (a%c) = [0, c-1) + (-c, 0] = (-c, c-1)
b1 = (q%c) = [0, c)

a1 - b1 = (-c,c-1) + (-c,0] = (-2c, c-1)

a1 - b1 ≡ 0 (mod c)     → a1 - b1 = 0, or -c

But, a1 - b1 = -c implied we need room, such that c+ulp(c) ≠ c
So, we switched it around, and test: (a1==b1) or (a1==b1-c)

---

c=b+1 (exactly) is equivalent to test: (b+1-1-b == 0) and (b+1-b-1 == 0)
With b>0, ulp(b) = ulp(c) = some multiples of machine epsilon.

This is lua 5.4 implmentation for floor-mod, from llimits.h

PHP Code:

#define luai_nummod(L,a,b,m)  \
  { (void)L; (m) = l_mathop(fmod)(a,b); \
    if (((m) > 0) ? (b) < 0 : ((m) < 0 && (b) > 0)) (m) += (b); }
#endif 


For our situation, m and b have the same sized ulp.
Even if it needed sign correction, (m) += (b), there no no cancellation errors.

ULP consideration is needed. Unlike truncated fmod, floor-mod might not be exact.
Eli Bendersky's post showed why fmod is exact: Computing Remainders by Doubling

We need to watch out for cancellation errors, when floor-mod fix the sign.

lua> eps = 1e-10
lua> x, y = eps%1, -eps%1
lua> x, y
1e-10      0.9999999999
lua> x+y, 1-y-x                -- x+y rounded to 1, but not exactly equal 1
1            8.274037096265818e-018