(42S) Quotient and Remainder

[Werner]

[Updated to use floor-divide. X and/or Y may be negative, and Y = Q*X+R holds when possible]
The 41C version was easy, but here we have to solve the halfway problem, eg:
Y: 8e11+2
X: 4
must return
Q: 2E11
R: 2

but

Y: 8e11+6
X: 4
must return
Q: 2E11+1
R: 2

In the first case, the division rounds down, in the second it rounds up.

Thanks to Albert Chan for the halfway correcting algorithm.
(which I adapted so I could keep on using IP instead of FLOOR)

split Y = Q*X + R

t := X/Y;
Q := INT(t);
R := MOD(Y,X);
if t#Q then return(R,Q);
t := X - R;
if t<R then Q:= Q - 1;
if t=R then Q:= Q - MOD(FLOOR(MOD(Y,X*10)/X),2); /* equals IP(MOD(MOD(Y,X*10)/X,2)) */
return(R,Q);

Code:

@ given X and Y, determine Quotient DIV(Y,X) and Remainder MOD(Y,X)
@               L       X       Y       Z       T
@ In:                   X       Y
@ Out:          X       R       Q
00 { 69-Byte Prgm }
01▸LBL "RQ" @           X       Y
02 RCL ST Y
03 RCL÷ ST Y @         Y/X      X       Y
04 ENTER
05 IP @                 Q      Q..      X       Y
06 X=Y?
07 GTO 00
08 X>Y? @ FLOOR
09 DSE ST X @ always skips
10 X≠Y? @ NOP
11 R↓
12 R↓ @                 X       Y       Q
13 MOD @        X       R       Q
14 RTN
15▸LBL 00 @             Q       Q       X       Y
16 R^ @                 Y       Q       Q       X
17 STO ST Z
18 R^ @                 X       Y       Q       Y
19 MOD @        X       R       Q       Y       Y
20 LASTX
21 RCL- ST Y @  X      X-R      R       Q       Y
22 X=Y?
23 GTO 02
24 X<Y?
25 +/-
26 X<0?
27 DSE ST Z @ will skip when Q is negative
28▸LBL 00
29 R↓
30 RTN
31▸LBL 02 @             R       R       Q       Y
32 R↓
33 X<> ST Z @           Y       Q       R       R
34 20
35 RCL× ST T @         X.10     Y       Q       R
36 MOD
37 R^
38 STO+ ST X
39 ÷
40 2
41 MOD
42 IP
43 -
44 R↓
45 +
46 -
47 +/-
48 END

examples:

[code]Y: 4e11+39       Q: 1e10
X: 40            R: 39

Y: 2e12          Q: 666 666 666 666
X: 3             R: 2

Y: 5e23          Q: 5e11
X: 1e12-1        R: 5e11

Y: 2e23+3e12     Q: 5e11+7
X: 4e11          R: 2e11

will work for X and Y integers, abs(X)<1e12 and abs(Y)<=X*1e12
I confess I haven't tried it with fractional inputs, but that should also work. Counterexamples welcome, I think ;-)
As with the 41 code, this routine will work when Q can be represented with 12 digits eg Y=1e40 and X=6 will not work.

Cheers, Werner

[Albert Chan]

Hi, Werner

It had been 2 years. I provided 2 versions for floor-divide

#1: we assumed "normal" division was round-to-nearest, half-way-to-even.
The trick is to reverse engineer how it was done, and reverse the "damage".

(12-20-2020 03:22 AM)Albert Chan Wrote:  To handle signs of a, b, correction test: remainder(a,b)/b < 0 and q-1 or q

#2: more direct way, based from FMA(-q, b, a) = -q*b + a
Fixing q (for floor-divide) is similar: FMA(-q, b, a)/b < 0 and q-1 or q

The downside is that 42S does not have FMA (almost no calculator does ...)
Free42/Plus42 only had it recently (Free42 version ≥ 2.5.23, all Plus42 versions)

(01-08-2021 05:24 PM)Albert Chan Wrote:  FMA based DIV behave the same as my REM based DIV.

Click on the green arrow for floor-divide code.

---

Note that Free42/Plus42 (42S too?) MOD is really floor-mod. Q should be floor-quotient to match.
Example, Python use floor-divide/floor-mod: divmod(y,x) == (y//x, y%x)

>>> Y, X = 8*10**11+2, 4
>>> divmod(-Y, X)
(-200000000001L, 2L)
>>> divmod(Y, -X)
(-200000000001L, -2L)

OP code, it seems Q is truncated quotient, R is floor-mod.
Or, if you want truncated quotient, R should match that.
Either way, invariant Y = Q*X+R should hold.

(-Y) (X) "RQ" --> (-2e11) (+2)      ?
(Y) (-X) "RQ" --> (-2e11) (-2)       ?

[Werner]

(04-05-2023 12:39 PM)Albert Chan Wrote:  >>> Y, X = 8*10**11+2, 4
>>> divmod(-Y, X)
(-200000000001L, 2L)
>>> divmod(Y, -X)
(-200000000001L, -2L)

Hi Albert,
My RQ routine matches the results above? remark that I do use
Q:= Q - MOD(FLOOR(MOD(Y,X*10)/X),2);
but I switched FLOOR and MOD around, then I could use IP:
Q := Q - IP(MOD(MOD(Y,X*10)/X,2)) ;

the results of MOD(Y,X*10)/X and the corresponding adjustment of Q (0 or 1) are
-3.5 - 0
-2.5 - 1
-1.5 - 0
-0.5 - 1
0.5 - 0
1.5 - 1
2.5 - 0

whether you now do MOD(FLOOR(X),2) or IP(MOD(X,2)) is the same.
ah but the initial IP should also be FLOOR then I see.
(don't get me wrong, but your explanations are sometimes a bit short, and I don't understand, or understand wrongly ;-)

I wanted to:
- use only the stack
- return X in LASTX
- be fast for simple cases, exceptions last.

That ruled out many of your suggestion algorithms, eg the one with remainder, or FMA.
I still have to do the Free42/DM42 one where I will be able to use FMA;

I will correct the routines, thanks

Cheers, Werner

[Albert Chan]

Hi, Werner

Sorry if I was a bit unclear.

Python code only tried to show invariant, Y = Q*X + R, always hold.
It does not ask for QR code to match it.

Truncated quotient is fine, but should match with truncated remainder.

Quote:I still have to do the Free42/DM42 one where I will be able to use FMA;

I like the FMA approach. It is simple and direct.

(Y-Q*X)/X = Y/X - Q

The expression is simply the fractional part that was rounded-off.
For floor-divide, simply make sure fractional part is non-negative.

[Werner]

So, to work for negative X and/or Y, the algorithm should become:
(changes in red)

t := X/Y;
Q := FLOOR(t);
R := MOD(Y,X);
if t#Q then return(R,Q);
t := X - R;
if ABS(t)<ABS(R) then Q:= Q - 1;
if t=R then Q:= Q - MOD(FLOOR(MOD(Y,X*10)/X),2); /* equals IP(MOD(MOD(Y,X*10)/X,2)) */
return(R,Q);

Sigh. Back to the drawing board.[update: done]
Werner

[Albert Chan]

I added some documentation to pseduo-code, to make algorithm more clear.
Please feel free to correct my errors ...

Quote:t := Y/X;
Q := FLOOR(t);
R := MOD(Y,X);
if t != Q then return(R,Q);

If Y/X not close to integer, floor-divide Q is correct.
R is floor-mod, always correct, just come for the ride.

Quote:t := X - R;

Y = Q*X + R = (Q+1)*X + (R-X) = (Q+1)*X + (-t)

The rest of the code is to decide remainder(Y,X) = R or -t
If remainder is (-t), we had over-estimated Q by 1

Quote:if ABS(t)<ABS(R) then Q := Q - 1;

With floor-mod, sign(X) = sign(R) = sign(t)
If Z = remainder(Y,X)/X = (-t)/X < 0, then Q over-estimated

Example, with a 2 significant digits calculator
80/7 = 11.428... ≈ 11. remainder(80,7) = +3 --> floor-divide Q = 11
81/7 = 11.571... ≈ 12. remainder(81,7) = −3 --> floor-divide Q = 12 - 1 = 11

Quote:if R != t then return (R,Q)

If |Z| < 1/2, not half-way case, we are done.

Quote:t := MOD(Y,X*10) / X;
return Q - MOD(FLOOR(t),2)

Q is half-way case, rounded-to-even. OP 2 half-way case example:

>(8e11+2) / 4
 200000000000
>(8e11+6) / 4
 200000000002

Z = ±1/2. We zoom-in to decide Z sign.

We could scale by 2, but scale with system base is error-free, simply an exponent shift.
This made t numerator calculation exact. (MOD part is exact)

t numerator divide by denominator may generate rounding errors, but it does not matter.
With half-way case, we expected t = .5, 1.5, 2.5, 3.5, ...

Since we scaled by 10, Z sign from t = +- +- +- +- +-
In other words, Z = -1/2 if FLOOR(t) =  1,  3,  5,  7,  9

--> Z = 1/2 * (-1)^FLOOR(t)

[Albert Chan]

Floor-divide may be more useful, in the sense that it is trivial to convert to other kinds.

Code:

def sign(x): return (x>0) - (x<0)
def qr_floor(y,x, s=1): return divmod(s*y,x)
def qr_trunc(y,x, s=1): q,r = divmod(s*abs(y),abs(x)); return sign(x*y)*q, sign(y)*r
def qr_ceil(y,x): q,r = qr_floor(y,x,-1); return -q,-r
def qr_away(y,x): q,r = qr_trunc(y,x,-1); return -q,-r

>>> y, x = 10, 7
>>> funcs = (qr_floor, qr_trunc, qr_ceil, qr_away)
>>> args = ((y,x), (-y,x), (y,-x), (-y,-x))
>>> for f in funcs: print f.__name__, [f(*a) for a in args]
...
qr_floor [(1, 3), (-2, 4), (-2, -4), (1, -3)]
qr_trunc [(1, 3), (-1, -3), (-1, 3), (1, -3)]
qr_ceil  [(2, -4), (-1, -3), (-1, 3), (2, 4)]
qr_away  [(2, -4), (-2, 4), (-2, -4), (2, 4)]

---

Truncated quotient, because of symmetry, is easier to implement.
It also has the accuracy advantage: y = (q_t x) + r_t

RHS terms have same sign, thus no cancellation errors.