(42, all flavours) Integer Division - how?
12-19-2020, 10:03 AM
Post: #41
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-19-2020, 02:55 PM (This post was last modified: 12-19-2020 02:56 PM by Albert Chan.)
Post: #42
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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.
12-19-2020, 04:42 PM (This post was last modified: 12-19-2020 08:14 PM by Werner.)
Post: #43
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-19-2020, 09:21 PM
Post: #44
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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.
12-19-2020, 09:56 PM (This post was last modified: 12-21-2020 10:04 PM by Albert Chan.)
Post: #45
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
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 qend
12-20-2020, 02:14 AM
Post: #46
 robve Senior Member Posts: 458 Joined: Sep 2020
RE: (42, all flavours) Integer Division - how?
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

"I count on old friends to remain rational"
12-20-2020, 03:22 AM (This post was last modified: 01-08-2021 03:05 PM by Albert Chan.)
Post: #47
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
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
12-20-2020, 08:46 AM
Post: #48
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
(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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-20-2020, 10:09 AM
Post: #49
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
(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 qend

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

Cheers, Werner

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-20-2020, 11:06 AM
Post: #50
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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.
12-20-2020, 11:50 AM
Post: #51
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
I know but HP calcs have decimal floating point ;-)
Werner

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-20-2020, 08:37 PM (This post was last modified: 12-20-2020 08:45 PM by Albert Chan.)
Post: #52
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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.
12-21-2020, 07:51 AM
Post: #53
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
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

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-21-2020, 03:20 PM
Post: #54
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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
12-21-2020, 07:35 PM
Post: #55
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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
12-23-2020, 02:59 PM (This post was last modified: 12-23-2020 09:34 PM by Albert Chan.)
Post: #56
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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 Qqend

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.
12-24-2020, 11:12 AM
Post: #57
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
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 caseend
12-24-2020, 12:34 PM
Post: #58
 Werner Senior Member Posts: 887 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Yes, this is the first one that is simpler than idiv3, in RPN, and a lot!
Thanks, Albert!
Werner

41CV†,42S,48GX,49G,DM42,DM41X,17BII,15CE,DM15L,12C,16CE
12-25-2020, 01:07 PM
Post: #59
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
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
01-08-2021, 05:24 PM
Post: #60
 Albert Chan Senior Member Posts: 2,682 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(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|)
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)