(42, all flavours) Integer Division - how?

[Werner]

Since b+b may overflow, I tend to stick to idiv3:

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

But I rewrite it as

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

Turning it into RPN, the latter can be done using only the stack as I need 1 less variable.
(it's part of a Quotient-and-remainder routine a=q*b+r, and I return b in L, q in X and r in Y)

Werner

[Albert Chan]

(12-19-2020 10:03 AM)Werner Wrote:  But I rewrite it as

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

If a ≥ b > 0, the rewrite work.

If a < 0, q < 0, q-c might overflow
If |a| < b, it may hit with floor-mod cancellation errors.

For tiny a = ε < b, rewritten code had:
a1 = a%b = ε
b1 = (q-c+a%c)%c = (-c+ε) % c = c - (c-ε), might not get back ε

lua> a, b, c, q = 1e-10, 2, 3, 0
lua> a1 = a%b
lua> b1 = (q - c + a%c)%c
lua> a1, b1 --> inequality force correction, returning wrong q-1 = -1
1e-010      1.000000082740371e-010

Original version work for all a, if q = ceil(a/b), or IP(a/b)

For tiny a = ε < b, original code had:
a1 = a%b - a%c = ε - ε = 0
b1 = q%c = 0, or 1

Both variables are exactly calculated. Thus correction code will work.

[Werner]

I always use IP.
41/42 doesn't have native FLOOR or CEIL.

So, this would be the 42 version:
In: X Y
Out: Q R, X in LASTX

Code:

00 { 70-Byte Prgm } @     X       Y       Z       T
01▸LBL "RQ" @             X       Y
02 RCL ST Y @             Y       X       Y
03 RCL÷ ST Y @           Y/X      X       Y
04 ENTER
05 IP @                  [Q]      Q       X       Y
06 X=Y?
07 GTO 00
08 R↓
09 R↓ @                   X       Y      [Q]
10 MOD
11 X<>Y @         X      [Q]      R
12 RTN
13▸LBL 00 @               Q       Q       X       Y
14 CLX
15 + @                    Q       X       Y       Y
16 R↓ @                   X       Y       Y       Q
17 MOD @          X       R       Y       Q       Q
18 LASTX @                X       R       Y       Q
19 RCL- ST Y @    X      X-R      R       Y       Q
20 X>Y?
21 GTO 00
22 X=Y?
23 GTO 01
24 DSE ST T
25▸LBL 00
26 X<> ST T @     X       Q       R
27 RTN
28▸LBL 01 @       X       R       R       Y       Q
29 R↓
30 X<>Y @         X       Y       R       Q
31 1
32 RCL+ ST L @            C       Y       R       Q
33 MOD @          C      Y%C      R       Q       Q
34 LASTX
35 - @            C     Y%C-C     R       Q       Q
36 R^
37 STO+ ST Y
38 X<> ST L @             C       T       R       Q
39 MOD @          C      T%C      R       Q       Q
40 X≠Y?
41 DSE ST T
42▸LBL 00
43 CLX
44 1
45 STO- ST L
46 X<> ST T @     X       Q       R
47 END

A version for the 41 would be at the same time simpler and harder: the 41 doesn't do round-to even (rounds up the halfway case), but doesn't have recall arithmetic.

Werner

[Albert Chan]

(12-19-2020 10:03 AM)Werner Wrote:  Since b+b may overflow, I tend to stick to idiv3:

We can generalized, to r = a % (b*base), to avoid rounding errors (base must be even)

Halfway case, r = k*b, possible k = 0.5, 1.5, 2.5, 3.5, 4.5, ... (base-0.5)
Since we can subtract multiples of 2b, still keeping even quotient, we have:

k = 0.5, 2.5, 4.5 ... ⇒ remainder(a,b) = +b/2
k = 1.5, 3.5, 5.5 ... ⇒ remainder(a,b) = −b/2

→ if floor(k) is odd then q := q - 1

lua> a, b = 460, 40
lua> q = rint(a/b) -- rint(11.5) = 12 (round-to-even)
lua> for base = 2, 12, 2 do
:      k = a % (b*base) / b
:      print(base, k, q - floor(k)%2)
:      end

2   1.5   11
4   3.5   11
6   5.5   11
8   3.5   11
10  1.5   11
12 11.5   11


Of course, it is best if Free42 had REM function.
Both IEEE double and decimal128 already have them built-in.

[Albert Chan]

In the mean time, this is my implementation for remainder(a,b), for Free42 Decimal.
Since a = q*b + r = (-q)*(-b) + r, remainder(a,b) = remainder(a,|b|)

That's why I set b = |b|, before the code even started.
This allowed REM work with signed arguments, both a and b.

Code:

00 { 44-Byte Prgm }
01▸LBL "REM"    @  L     X     Y     Z     T
02 ABS          @  All code below, b >= 0
03 RCL ST Y
04 RCL ST Y     @        b     a     b     a
05 MOD          @  b   a%b     b     a     a
06 STO- ST L
07 LASTX        @    b-a%b   a%b     b     a
08 X<>Y         @      a%b b-a%b     b     a
09 X<Y?
10 RTN          @ positive remainder
11 X=Y?
12 GTO 00
13 X<>Y         @ negative remainder
14 +/-
15 RTN
16▸LBL 00       @ halfway case
17 R↓
18 R↓
19 10
20 ×            @   b*base     a   b/2   b/2
21 MOD
22 X<>Y
23 STO+ ST X
24 ÷            @        k   b/2   b/2   b/2
25 IP
26 -1
27 X<>Y         @    IP(k)    -1   b/2   b/2
28 Y↑X
29 ×
30 END

Example: (for halfway case, quotient must be even)

 459 40 XEQ "REM"   →  19 =  459 - 11*40
 460 40 XEQ "REM"   → -20 =  460 - 12*40
 461 40 XEQ "REM"   → -19 =  461 - 12*40
-459 40 XEQ "REM"   → -19 = -459 + 11*40
-460 40 XEQ "REM"   →  20 = -460 + 12*40
-461 40 XEQ "REM"   →  19 = -461 + 12*40

9 PI XEQ "MOD" →  2.716814692820413523074713233440994
9 PI XEQ "REM" → -0.424777960769379715387930149838509

Since REM based from MOD, it is numerically better if we pull out a's sign.
This way, we avoided floor-mod subtraction cancellation errors, when it "fix" sign.

remainder(a,b) = sign(a) * remainder(|a|, b)

We can use REM to implement DIV, like this

(12-18-2020 12:08 AM)Albert Chan Wrote:  remainder() based idiv is very elegant.

PHP Code:

function idiv6(a, b)    -- assumed b > 0
    local r = a/b
    local q = floor(r)
    if q == r and remainder(a, b) < 0 then q=q-1 end
    return q
end 


[robve]

Using the int function (or floor function), the quotient Q and remainder R of Y/X is simply:

Q = sign(X)*int(Y/abs(X))
R = Y-Q*X

This satisfies Y=Q*X+R, 0<=R<|X|, for any Y and X!=0.

I've always used this formulation, rather than fiddling with MOD.

Is there a downside to this method, I mean in general, not specific to the 42?

- Rob

[Albert Chan]

DIV code, using REM. It work with signed arguments, both a and b.

Code:

00 { 42-Byte Prgm }
01▸LBL "DIV"    @  L     X     Y     Z     T
02 RCL ST Y
03 RCL ST Y     @        b     a     b     a
04 ÷
05 ENTER
06 IP           @    IP(q)     q     b     a
07 X>Y?
08 DSE ST X     @ floor(q) < 0, always skip
09 X<Y?
10 RTN
11 LSTO "q"
12 R↓
13 R↓
14 LSTO "b"
15 XEQ "REM"
16 RCL÷ "b"     @ remainder(a,b)/b
17 ENTER
18 SIGN         @ note: SIGN(0)=1
19 X>0?
20 CLX
21 RCL+ "q"     @ remainder(a,b)/b < 0 and q-1 or q
22 END

To handle signs of a, b, correction test: remainder(a,b)/b < 0 and q-1 or q

idiv(a, b) = idiv(sign(b)*a, |b|)      -- previously, code required positive divisor

remainder(a,b)/b
= remainder(a,|b|) / (sign(b)*|b|)
= sign(b) * remainder(a, |b|) / |b|
= remainder(sign(b)*a, |b|) / |b|

The code gives the same answer, as if we manually made b positive.

Note: remainder(a,b)/b = fractional part of quotient, if integer part is closest to a/b, not floored.

Example:

4e33 39 + STO "u" 40 STO "v" XEQ "DIV"      → X = 1e32, Y = -0.025

With Y < 0, this meant closest integer for u/v is X+1
(-u)/(-v) gives exactly the same result of u/v, as expected.
u/v = X+1 - Y = 1e32 + 0.975

(-u)/v gives X = -1e32 - 1, Y = 0.025

With Y >= 0, closest integer for (-u)/v = X
u/(-v) gives exactly the same result of (-u)/v, as expected.
(-u)/v = X + Y = -1e32 - 0.975

Note: if u/v != floor(u/v), the code skipped REM, and return floor(u/v), u/v

100 16 XEQ "DIV"      → X = 6, Y = 6.25

Update: removed the unneeded NOP after DSE

[Werner]

(12-20-2020 02:14 AM)robve Wrote:  Using the int function (or floor function), the quotient Q and remainder R of Y/X is simply:

Q = sign(X)*int(Y/abs(X))
R = Y-Q*X

This satisfies Y=Q*X+R, 0<=R<|X|, for any Y and X!=0.

In unlimited arithmetic, yes.
But it fails when Q is rounded up, or when Q*X carries more significant digits than the maximum (12 for a 42S or 34 for Free42).

For instance:
42S: Y: 1e12 X: 6
Y: 9.8696044e21 X: 31415926536

You might consider these examples farfetched and contrived, but when you try to implement a multiprecision division (as I mentioned in one of the previous posts here), it is what you need to get right.

Cheers, Werner

[Werner]

(12-19-2020 03:14 AM)Albert Chan Wrote:  

PHP Code:

function idiv6c(a, b)   -- assumed b > 0
    local r = a/b
    local q = floor(r)
    if q ~= r then return q end
    r = (a%b) * 2
    if r == 0 then r = a % (b+b) end
    return r > b and q-1 or q
end 


2 remarks:
- shouldn't that be if r== b instead of if r==0 ?
- and, a%b*2 may overflow again?

Cheers, Werner

[Albert Chan]

(12-20-2020 10:09 AM)Werner Wrote:  - shouldn't that be if r== b instead of if r==0 ?

Yes. I made a typo, and will correct it. Thanks.
For testing halfway, a%b == b/2, or (a%b)*2 == b

Quote:- and, a%b*2 may overflow again?

I assume "overflow" meant lost of precision ?

This is one nice thing about binary math, *2 or /2, we do not lose precision.
Scaling just update the exponent, similar to *10 or /10 for decimal system.

[Werner]

I know but HP calcs have decimal floating point ;-)
Werner

[Albert Chan]

(12-11-2020 10:48 AM)Werner Wrote:  (2e23+3e12) DIV 4e11 = 5e11+7 (5e11+8) (round-to-even, and a particularly difficult one to get right)

Round-to-even behavior actually make halfway case easy to get floor-divide.

a, b = 2e23 + 3e12, 4e11
a ≡ 0 + 2e11 ≡ b/2, (mod b), indeed a halfway case.
Problem is reduced to finding s = sign of remainder(a,b)

Let q = rounded-to-nearest of a/b, which we know is correct
Let h = halfway = b/2

a = q*b + s*(b/2) = (2*q+s)*h

While 3|h, reduce both a, h by 3. This made mod 3 test (below), work all the time.
Without this step, we get into modulo fractions equivalent of divide-by-zero. (*)

a ≡ (2*q + s)*h ≡ (s - q)*h (mod 3)
s ≡ a/h + q (mod 3)

Mod 3 is same as Mod 9, then Mod 3. Just add the digits, ignoring powers-of-tens.
For above example, h = 2e11 ≡ 2 ≡ -1 (mod 3), not divisible by 3. So, we proceed:

a = 2e23 + 3e12 ≡ 2+3 ≡ -1 (mod 3)
q = 5e11+8 ≡ 5 + 8 ≡ 13 ≡ 1 (mod 3)

→ s ≡ -1/-1 + 1 ≡ 2 ≡ -1 (mod 3)
→ floor-divide(a/b) = q - (s<0) = q - 1 = 5e11 + 7


(*) modulo arithmetics with fractions actually meant multiply by its inverse.
x ≡ a/b (mod m)      ⇔      x ≡ a*b^-1 (mod m)      ⇔      b*x ≡ a (mod m)

Example: 5/4 ≡ -1/1 ≡ -1 (mod 3)
By keeping the mod 3 range of -1,0,1, mod 3 fractions is same as normal fractions.

[Werner]

Hi Albert,
can that then not be simplified to verifying whether

a = (1-q)*h (mod 3)

and if not equal q := q-1

works for

2300 = 57*40 + 20 ( 2 ^= (1-58)%3 * 20%3, so q-1)
and
-2300 = -58*40 + 20 (1 = ((1+58)%3 *20%3 )%3; so q=-58 is correct)

Werner

[Albert Chan]

(12-21-2020 07:51 AM)Werner Wrote:  Hi Albert,
can that then not be simplified to verifying whether

a = (1-q)*h (mod 3)

and if not equal q := q-1

Yes ! For halfway case, s = ±1. If s≠1 then s=-1
I "divided" both side by h, only to show mod 3 test cannot have h with factor of 3.

With assumed s=1, we can also do mod 9 test.
mod 9 test allowed b having 1 factor of 3, instead of none at all.

a = q*b + s*h
a-h ≡ q*b (mod 9), if s=1

We start with sum-of-digits, to expand for non-integer inputs. (e.g. 1.2e-9 → 1+2 → 3)
Redoing previous post examples, but scale up both a, b by 3

lua> a, b = 2300*3, 40*3
lua> q = rint(a/b)
lua> a, b, q, a%b
6900      120      58      60

This is halfway case, with b%9 = 3 ≠ 0. We do mod 9 test.

a-h ≡ (6+9) - 6 ≡ 6 - 6 ≡ 0
q*b ≡ (5+8)*3 ≡ 4*3 ≡ 12 ≡ 3
→ s≠1, return q-1 = 57

lua> a, q = -a, -q -- 2nd example
lua> a, b, q, a%b
-6900      120      -58      60

a-h ≡ -(6+9) - 6 ≡ -6 - 6 ≡ -12 ≡ 6
q*b ≡ -(5+8)*3 ≡ -4*3 ≡ -12 ≡ 6
→ s=1, return q = -58

---

Even if initial q was wrongly rounded, we still get the right correction.

1st example, but q wrongly rounded 57.5 to 57:
q*b ≡ (5+7)*3 ≡ 3*3 ≡ 0      → = a-h (mod 9), return q = 57

2nd example, but q wrongly rounded to -57.5 to -57:
q*b = -(5+7)*3 ≡ 3*3 ≡ 0     → ≠ a-h (mod 9), return q-1 = -58

[Albert Chan]

(12-21-2020 03:20 PM)Albert Chan Wrote:  a = q*b + s*h
a-h ≡ q*b (mod 9), if s=1

I was being stupid, why mod 9 test, when we could do mod 10 ?
Again, b cannot have factors of 10 (nonzero last digit)

Redoing previous post examples, but with mod-10 test

a/b = 6900 / 120 = 690 / 12 = 57.5
q = 58
h = b/2 = 6

a-h ≡ 0 - 6 ≡ 4
q·b ≡ 8 * 2 ≡ 6
→ s≠1, return q-1 = 57

2nd example a = -a; q = -q

a-h ≡ (-0) - 6 ≡ 4
q·b ≡ (-8) * 2 ≡ 4
→ s=1, return q = -58

[Albert Chan]

(12-17-2020 08:53 AM)Werner Wrote:  I wanted to use integer division to correctly split
Cc00 = Qq*Xx + Rr, where each letter is a half-length integer and Cc<Xx

Revisiting DIV(Cc00, Xx), using mod 10 test for halfway case.

Note: the code assumed half-length = 1.
This can be easily generalize, by adjust code for last nonzero digit of Rr.

PHP Code:

function idiv7(Cc00, Xx)
    local Qq, Rr = d2(Cc00/Xx), Cc00%Xx -- simulate 2 digits calculator
    local c = floor(Qq)
    if Qq ~= c then return c end        -- Fast-path
    if Rr < Xx-Rr then return c end     -- Cc00/Xx rounded-down to c
    if Rr > Xx-Rr then return c-1 end   -- Cc00/Xx rounded-up to c
    c = Rr % 10                         -- halfway: Cc00 = (2*Qq +/- 1) * Rr
    if c == 0 then c = Rr/10 end        -- last nonzero digit of Rr
    if 0 ~= ((Qq%10)*2-9)*c % 10 then Qq = Qq-1 end
    return Qq
end 


Since Lua is not really a 2-digits calculator, I added d2():
d2 = function(x) return string.format('%.2g',x) + 0 end

The code confirmed correct for all possible 2-digits combinations of (Cc, Xx), Cc<Xx
Proof that above half-way test work generally, for half-length integers.

Cc00 / Xx = (Cc00/100) / (Xx/100)
But, Xx < 100, so factoring out 10's always have numerator ends in 0.

For halfway case, taking mod 10 both side, we have:

Cc00 = (2*Qq ± 1) * Rr
0 ≡ (2*Qq ± 1) * c (mod 10), where c = last non-zero decimal digit of Rr

Update:

If Xx had all factors of 10 removed, Xx must still be even.
Rr = Xx/2 ends in 1,2,3,4, 6,7,8,9. In other words, c ≠ 0, or 5

If (2*Qq + 1)*c ≡ 0, then (2*Qq - 1)*c ≡ -2*c ≠ 0 (mod 10)

Passing +1 test implied failing -1 test, and vice versa.

[Albert Chan]

Since c ≠ 0, or 5 (mod 10), we can rewrite the equality as:

0 ≡ (2*Qq ± 1) * (2c) (mod 10), where c = last non-zero decimal digit of Rr
0 ≡ (2*Qq ± 1) * c (mod 5)

+1 case: Qq ≡ 2 (mod 5)
−1 case: Qq ≡ 3 (mod 5)

But, -1 case should return Qq-1.
This meant half-way case floor-divide quotient always in 5n+2 form

PHP Code:

function idiv7b(Cc00, Xx)
    local Qq, Rr = d2(Cc00/Xx), Cc00%Xx
    local c = floor(Qq)
    if Qq ~= c then return c end
    if Rr < Xx-Rr then return c end     -- Cc00/Xx rounded-down to c
    if Rr > Xx-Rr then return c-1 end   -- Cc00/Xx rounded-up to c
    return c - (c-2)%5                  -- halfway case
end 


[Werner]

Yes, this is the first one that is simpler than idiv3, in RPN, and a lot!
Thanks, Albert!
Werner

[Albert Chan]

Previously, we had tried mod 3, mod 9, mod 10 to get halfway floor-divide quotient.

a = (2*q+1) * h, where h = b/2, q is floor-divide quotient of a/b

Mod 2 seems impossible. With (2*q+1) odd, a, h have same parity.

Say, b = 2^k * b', where b' is odd. Reduced both side by 2^(k-1), and take mod 4, we have

a' ≡ (2*q+1) * h' (mod 4)

With a', h' both odd, a' ≡ ± 1 (mod 4), and so is h'

If a' ≡ +h' (mod 4), then 2*q+1 ≡ 1 (mod 4)      → q ≡ 0 (mod 2)
If a' ≡ -h' (mod 4), then 2*q+1 ≡ -1 (mod 4)      → q ≡ 1 (mod 2)

So, we can use mod 2 to get floor-divided quotient after all

\(\qquad\qquad\large q ≡ {a'-h' \over 2} ≡ {a-h \over 2^k} \normalsize \pmod 2 \)

Redo previous post example, using mod 2:

>>> from gmpy2 import *
>>> a, b = 6900, 120 # a/b = 57.5, a halfway case
>>> k = bit_scan1(b) # k = 3, see Find First Set
>>> q, h = 58, 60
>>> q - ((q&1) ^ (bit_test(a,k) ^ bit_test(h,k)))
57

[Albert Chan]

(12-20-2020 03:22 AM)Albert Chan Wrote:  DIV code, using REM. It work with signed arguments, both a and b.

The latest Free42 Windows version (2.5.23, 1/6/2021) now has FMA.
Thank you. Thomas Okken

PI PI PI X^2 +/- XEQ "FMA"      → -1.37075656833106266883964580072991e-34

FMA based DIV behave the same as my REM based DIV.

Code:

00 { 36-Byte Prgm }
01▸LBL "DIV"    @  L     X     Y     Z     T
02 RCL ST Y
03 RCL ST Y     @        b     a     b     a
04 ÷
05 ENTER
06 IP           @    IP(q)     q     b     a
07 X>Y?
08 DSE ST X     @ floor(q) < 0, always skip
09 X<Y?
10 RTN          @       q      q     b     a
11 X<> ST Z     @       b      q     q     a
12 LSTO "b"
13 +/-
14 R↑           @       a     -b     q     q
15 FMA          @   a-b*q      q     q     q
16 RCL÷ "b"
17 X<0?
18 DSE ST Y
19 NOP
20 X<>Y
21 END

Note: (a-b*q)/b = (sign(b)*a - |b|*q) / |b|

In other words, DIV(a,b) = DIV(sign(b)*a, |b|)