+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: General Software Library (/forum-13.html)
+--- Thread: (42S) Quotient and Remainder (/thread-19775.html)
(42S) Quotient and Remainder - Werner - 04-05-2023 09:38 AM
[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: 2e11will 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
RE: (42S) Quotient and Remainder - Albert Chan - 04-05-2023 12:39 PM
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) ?
RE: (42S) Quotient and Remainder - Werner - 04-05-2023 01:11 PM
(04-05-2023 12:39 PM)Albert Chan Wrote: >>> Y, X = 8*10**11+2, 4Hi Albert,
>>> divmod(-Y, X)
(-200000000001L, 2L)
>>> divmod(Y, -X)
(-200000000001L, -2L)
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
RE: (42S) Quotient and Remainder - Albert Chan - 04-05-2023 01:48 PM
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.
RE: (42S) Quotient and Remainder - Werner - 04-05-2023 03:55 PM
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
RE: (42S) Quotient and Remainder - Albert Chan - 04-06-2023 12:03 AM
I added some documentation to pseduo-code, to make algorithm more clear.
Please feel free to correct my errors ...
Quote:t := Y/X;If Y/X not close to integer, floor-divide Q is correct.
Q := FLOOR(t);
R := MOD(Y,X);
if t != Q then return(R,Q);
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;Q is half-way case, rounded-to-even. OP 2 half-way case example:
return Q - MOD(FLOOR(t),2)
>(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)
RE: (42S) Quotient and Remainder - Albert Chan - 04-06-2023 02:55 PM
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 = (qt x) + rt
RHS terms have same sign, thus no cancellation errors.