(42, all flavours) Integer Division - how?
12-11-2020, 10:48 AM
Post: #1
 Werner Senior Member Posts: 768 Joined: Dec 2013
(42, all flavours) Integer Division - how?
Anybody know how to implement integer division in a fast and reliable way?
And, no, / IP is not the answer..

Illustrating my point with a hypothetical 2-digit calculator, then:
79 DIV 40 = 1
200 DIV 3 = 66
890 DIV 99 = 8
2300 DIV 40 = 57

the goal is to be able to split Y=Q*X+R when possible, i.e. up to Cc00 = Qq*Xx + Rr with Cc and Rr < Xx, and Cc00 DIV Xx should give Qq (all single letters are half-length integers, or single digits in the case of the 2-digit calculator)

42S equivalents
(4e11+39) DIV 40 = 1e10 (/ IP = 1e10+1)
2e12 DIV 3 = 666666666666 (666666666667)
5e23 DIV (1e12-1) = 5e11 (5e11+1)
(2e23+3e12) DIV 4e11 = 5e11+7 (5e11+8) (round-to-even, and a particularly difficult one to get right)

Free42 equivalents
(4e33+39) DIV 40 = 1e32
2e34 DIV 3 = 666...666

I CAN do it, but it's not remotely pretty.
It would be a welcome addition to Free42 ;-)

Cheers, Werner
12-11-2020, 03:12 PM
Post: #2
 Gerson W. Barbosa Senior Member Posts: 1,498 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Code:
 00 { 17-Byte Prgm } 01▸LBL "INT÷" 02 RCL ST Y 03 X<>Y 04 MOD 05 STO- ST Y 06 X<> ST L 07 ÷ 08 END

This gets only your first example right on the HP-42S. So does INT÷ on the hp 33s.
12-11-2020, 04:47 PM
Post: #3
 Allen Member Posts: 196 Joined: Aug 2014
RE: (42, all flavours) Integer Division - how?
(12-11-2020 10:48 AM)Werner Wrote:  Anybody know how to implement integer division in a fast and reliable way?
And, no, / IP is not the answer..

Illustrating my point with a hypothetical 2-digit calculator, then:
79 DIV 40 = 1
200 DIV 3 = 66
890 DIV 99 = 8
2300 DIV 40 = 57

Have you tried BASE÷ ?

17bii | 32s | 32sii | 41c | 41cv | 41cx | 42s | 48g | 48g+ | 48gx | 50g | 30b

12-11-2020, 04:56 PM
Post: #4
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Base/ doesn’t work for 7.2 DIV 5 for instance, but also not for the Free42 examples, as it is limited to 36 bits in the 42S and 64 in Free42.
Werner
12-11-2020, 07:32 PM
Post: #5
 Albert Chan Senior Member Posts: 2,142 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(12-11-2020 10:48 AM)Werner Wrote:  Anybody know how to implement integer division in a fast and reliable way?

I don't know if this is fast, but FMA should work.

10 DEF FND(A,B)
20 Q=FLOOR(A/B)
30 B1=B*1000001 @ B1=B+B1-B1 @ B=B-B1
40 Q1=Q*1000001 @ Q1=Q+Q1-Q1 @ Q=Q-Q1
50 FND=FLOOR((A-Q1*B1-Q*B1-B*Q1-B*Q)/(B+B1))+Q+Q1
60 END DEF

>RUN
>FND(4E11+39,40)
10000000000
>FND(2E12,3)
666666666666
>FND(5E23,1E12-1)
500000000000
>FND(2E23+3E12,4E11)
500000000007
12-11-2020, 08:04 PM
Post: #6
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Yes that looks like Dekker’s double precision routines. That works, but in RPN it doesn’t look quite so elegant. I was hoping there would be a shorter, simpler way.
Werner
12-11-2020, 09:28 PM
Post: #7
 Gjermund Skailand Member Posts: 52 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Hi.
Here is one version fulfilling the test examples on the Free42.
It uses the DOT function to force calculation with all digits prior to rounding.
I am not sure how reliable it actually is.
It uses one extra stack level

00 {28-Byte Prgm }
01 LBL "IDIV2"
02 X<>Y
03 -0.5
04 RCLX ST Z
05 COMPLEX
06 X<>Y
07 1/X
08 ENTER
09 COMPLEX
10 DOT
11 IP
12 END

Best regards
Gjermund
12-11-2020, 09:46 PM
Post: #8
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Never too old to learn! I didn’t know DOT worked on complex numbers, too. I will have to take a look at this, but one thing’s for sure: DOT does not use extended precision, as it does in the 42S.
Werner
12-12-2020, 10:36 AM
Post: #9
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
4e33+6 DIV 4 should be 1e33+1, and your routine returns 1e33+2
But it does work on a 42S, where DOT uses 15 intermediate digits:
4e11+6 DIV 4 = 1e11 + 1
well, - for this particular example. But since it basically does Y/X - 1/2,
it fails for eg 7 DIV 3
Nevertheless, I learned something and it will be put to good use!
Werner
12-12-2020, 09:31 PM
Post: #10
 Gjermund Skailand Member Posts: 52 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Yeah, it turned out to be a bad idea. I also found out that the HP50g and the 42S behaves differently for the DOT and CROSS when using complex numbers for 2D. HP50g will not allow it.
best regards
Gjermund
12-13-2020, 01:24 AM
Post: #11
 Albert Chan Senior Member Posts: 2,142 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
There is the "cheating" way, by temporarily setting the rounding mode
I coded mathx.setround() for this purpose

PHP Code:
require 'mathx'function idiv1(a,b)    mathx.setround(-1)  -- round downwards    a = mathx.rint(a/b)    mathx.setround(0)   -- restore round-to-nearest    return aend

Another way is to correct the quotient of a and b (assumed both integers)
Here, we assumed q, Q may have errors of ±1

a = q*b + r = q*c - q + r ﻿ ﻿ ﻿ ﻿ ﻿ , where c = b+1
a = Q*c + R

0 = (q-Q)*c - q + r - R﻿ ﻿ ﻿ ﻿ ﻿ ﻿ → -q + r - R ≡ 0 (mod c)

Since r and R can be calculated with MOD, we can correct for q
Assuming |q| < 2^53, this is the code:
PHP Code:
function idiv2(a,b) -- assumed b > 0    local q, r, c = floor(a/b), a%b, b+1    r = (r+1 - q%c - a%c) % c - 1    return q + rend

lua> a = 0x1p72
lua> for b=1e6+1, 1e6+9 do
: ﻿ ﻿ ﻿ ﻿ ﻿ q = idiv1(a,b) -- reference
: ﻿ ﻿ ﻿ ﻿ ﻿ print(b, q - floor(a/b), q - idiv2(a,b))
: ﻿ ﻿ ﻿ ﻿ ﻿ end

1000001 ﻿ ﻿ ﻿ ﻿ -1 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000002 ﻿ ﻿ ﻿ ﻿ -1 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000003 ﻿ ﻿ ﻿ ﻿ -1 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000004 ﻿ ﻿ ﻿ ﻿ ﻿ 0 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000005 ﻿ ﻿ ﻿ ﻿ ﻿ 0 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000006 ﻿ ﻿ ﻿ ﻿ ﻿ 0 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000007 ﻿ ﻿ ﻿ ﻿ ﻿ 0 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000008 ﻿ ﻿ ﻿ ﻿ ﻿ 0 ﻿ ﻿ ﻿ ﻿ ﻿ 0
1000009 ﻿ ﻿ ﻿ ﻿ -1 ﻿ ﻿ ﻿ ﻿ ﻿ 0
12-13-2020, 08:12 AM
Post: #12
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Yes. YES!
Thanks a million, Albert, this is what I was looking for!
(the first part, cheating, doesn't apply to 41,42,Free42 of course)
Werner
12-13-2020, 12:09 PM
Post: #13
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
All that is left is to turn it in a routine..First try:
Code:
 { 38-Byte Prgm } @     X       Y       Z       T  LBL "DIV" @            b       a  RCL ST Y  RCL ST Y  /  IP @                   q       b       a       a  Rv  MOD @          b       r       a       q       q  X<>Y  1  STO+ ST Z  RCL+ ST L @            c       a       r+1     q  MOD  STO- ST Y @    c       R       r+1-R   q       q  X<> ST T  LASTX  MOD @          c       q%c     r+1-R   q       q  STO- ST Y  X<> ST L @             c       r'      q       q  MOD  +  DSE ST X  END

I can recover b too, if needed, but that's it.

Werner
12-13-2020, 01:18 PM
Post: #14
 Albert Chan Senior Member Posts: 2,142 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(12-13-2020 01:24 AM)Albert Chan Wrote:  return q + r

We had assumed q never overflow, which results in |r| ≤ 1

However, if q already overflowed, calculated r is basically garbage.

We should return q + sign(r), to limit the damage.
In other words, if q overflow calculator precision, don't correct the quotient.
12-13-2020, 04:21 PM (This post was last modified: 12-13-2020 04:22 PM by Albert Chan.)
Post: #15
 Albert Chan Senior Member Posts: 2,142 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(12-13-2020 01:18 PM)Albert Chan Wrote:  We had assumed q never overflow, which results in |r| ≤ 1

Assuming q is correctly rounded (mode round-to-nearest), r = 0 or -1

The code is simplified to correct, or not correct.

PHP Code:
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-1end
12-13-2020, 06:41 PM
Post: #16
 Albert Chan Senior Member Posts: 2,142 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(12-13-2020 04:21 PM)Albert Chan Wrote:  Assuming q is correctly rounded (mode round-to-nearest), r = 0 or -1

This allowed optimization of my FMA (Fused-Multiply-Add) code.
Without using MOD, this is the fastest of all.

Bonus: it removed the integer arguments requirement.
PHP Code:
function idiv4(a,b)     -- asuumed b > 0    local q = floor(a/b)    local b1 = b * (0x1p27+1)    local q1 = q * (0x1p27+1)    b1 = b + b1 - b1    -- hi bits    q1 = q + q1 - q1    local b0 = b - b1   -- lo bits    local q0 = q - q1    if a - q1*b1 - q1*b0 - q0*b1 >= q0*b0 then return q end    return q - 1end

It would be nice if Free42 exposed FMA(a,b,c) to the user ...
12-13-2020, 07:19 PM
Post: #17
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
(12-13-2020 04:21 PM)Albert Chan Wrote:  Assuming q is correctly rounded (mode round-to-nearest), r = 0 or -1

Only for a and b positive?
For a negative, it doesn’t work, eg
(12 digits calc)
a=-4e11-6
b=4
DIV returns -1e11-2

Werner
12-13-2020, 07:43 PM
Post: #18
 Thomas Okken Senior Member Posts: 1,828 Joined: Feb 2014
RE: (42, all flavours) Integer Division - how?
(12-13-2020 06:41 PM)Albert Chan Wrote:  It would be nice if Free42 exposed FMA(a,b,c) to the user ...

Noted...
12-13-2020, 08:51 PM
Post: #19
 Albert Chan Senior Member Posts: 2,142 Joined: Jul 2018
RE: (42, all flavours) Integer Division - how?
(12-13-2020 07:19 PM)Werner Wrote:
(12-13-2020 04:21 PM)Albert Chan Wrote:  Assuming q is correctly rounded (mode round-to-nearest), r = 0 or -1
Only for a and b positive?
For a negative, it doesn’t work, eg
(12 digits calc)
a=-4e11-6
b=4
DIV returns -1e11-2

I defined IDIV matching MOD behavior: a = b * IDIV(a,b) + MOD(a,b)

Free42, binary and decimal, uses floor-mod: (a MOD b) has sign of b
To match it, IDIV(a,b) = floor(a/b)

So, above DIV is correct: floor((-4e11 - 6)/4) = floor(-1e11 - 1.5) = -1e11 - 2

Python also define it this way, see Why Python Integer Division Floors

>>> a, b = -4*10**11-6, 4
>>> q, r = a//b, a%b
>>> print q, r, q*b+r
-100000000002 ﻿ ﻿ ﻿ ﻿ ﻿ 2 ﻿ ﻿ ﻿ ﻿ ﻿ -400000000006
12-13-2020, 08:59 PM
Post: #20
 Werner Senior Member Posts: 768 Joined: Dec 2013
RE: (42, all flavours) Integer Division - how?
Ah yes, indeed!
Werner
 « Next Oldest | Next Newest »