(42, all flavours) Integer Division - how?

[Werner]

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

[Gerson W. Barbosa]

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.

[Allen]

(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÷ ?

[Werner]

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

[Albert Chan]

(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

[Werner]

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

[Gjermund Skailand]

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

[Werner]

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

[Werner]

Too bad, Gjermund: on Free42,
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

[Gjermund Skailand]

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

[Albert Chan]

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 a
end 


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 + r
end 


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

[Werner]

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

[Werner]

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

[Albert Chan]

(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.

[Albert Chan]

(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-1
end 


[Albert Chan]

(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 - 1
end 


It would be nice if Free42 exposed FMA(a,b,c) to the user ...

[Werner]

(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

[Thomas Okken]

(12-13-2020 06:41 PM)Albert Chan Wrote:  It would be nice if Free42 exposed FMA(a,b,c) to the user ...

Noted...

[Albert Chan]

(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

[Werner]

Ah yes, indeed!
Werner