(42, all flavours) Integer Division - how?
|
12-19-2020, 10:03 AM
Post: #41
|
|||
|
|||
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
|
|||
|
|||
RE: (42, all flavours) Integer Division - how?
(12-19-2020 10:03 AM)Werner Wrote: But I rewrite it as 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
|
|||
|
|||
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 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
|
|||
|
|||
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
|
|||
|
|||
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 } 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. |
|||
12-20-2020, 02:14 AM
Post: #46
|
|||
|
|||
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
|
|||
|
|||
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 } 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
|
|||
|
|||
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: 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
|
|||
|
|||
RE: (42, all flavours) Integer Division - how?
(12-19-2020 03:14 AM)Albert Chan Wrote: 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
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
RE: (42, all flavours) Integer Division - how?
(12-21-2020 07:51 AM)Werner Wrote: Hi Albert, 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
|
|||
|
|||
RE: (42, all flavours) Integer Division - how?
(12-21-2020 03:20 PM)Albert Chan Wrote: a = q*b + s*h 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
|
|||
|
|||
RE: (42, all flavours) Integer Division - how?
(12-17-2020 08:53 AM)Werner Wrote: I wanted to use integer division to correctly split 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) 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
|
|||
|
|||
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) |
|||
12-24-2020, 12:34 PM
Post: #58
|
|||
|
|||
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
|
|||
|
|||
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
|
|||
|
|||
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 } 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: