+- HP Forums (https://www.hpmuseum.org/forum)
+-- Forum: HP Software Libraries (/forum-10.html)
+--- Forum: HP-41C Software Library (/forum-11.html)
+--- Thread: (41) Exponentiation & Residue Reduction (/thread-14798.html)
(41) Exponentiation & Residue Reduction - SlideRule - 04-06-2020 01:31 PM
Pages 341 & 344 from Number Theory in Science and Communication, Second Enlarged Edition, © Sprinler·Veriag Berlin Heidelberg 1984 and 1986, ISBN 978-3-662-22246-1 (eBoook)
"A. A Calculator Program for Exponentiation and Residue Reduction
We list here a program for calculating the least nonnegative remainder of aⁿ modulo m on a Hewlett-Packard 41C or 41CV pocket calculator. In other words, in our acute bracket notation, the calculator calculates
The 3 variables, the base a, the exponent n and the modulus m are entered into the calculator in that order.
To call the program, which is labeled "AN", from storage, one presses
GTO "AN" .
To calculate, say,
{2³⁴⁰}₃₄₁ ,
one proceeds by pressing the following buttons:
2
ENTER
340
ENTER
341
To start the program, one presses
R/S
After 2 seconds the HP 41 display shows in rapid succession
8.
6.
4.
2.
This is the binary decomposition of the exponent 340.
Check: 28 + 26 + 24 + 22 = 340. Check!
Next, the computer will start the necessary repeated squarings and reductions modulo 341, which will take about 7 seconds. Then the display will show the end result, a 1 (without a decimal point!).
Thus,
{2³⁴⁰)₃₄₁ = 1 .
In other words, 341 is a pseudoprime to the base 2.
Similarly, after pressing 2, ENTER, 170, ENTER, 341, R/S, the display shows in succession:
7. 5. 3. 1. 1
where that last display (the 1 without the decimal point) tells us that
{2¹⁷⁰}₃₄₁ = 1 .
Proceeding in the same manner, we find
{2⁸⁵}₃₄₁ = 32 .
i.e., 341 is not a strong pseudoprime to the base 2. Also, by using 3 as a base we find
{3³⁴⁰}₃₄₁ = 56 .
Thus, 341 is certainly not an absolute pseudoprime (Carmichael number).
However, for the modulus 2821, we find
{2²⁸²⁰}₂₈₂₁ = {3²⁸²⁰}₂₈₂₁ = 1 .
two of the many steps necessary to show that 2821 is an absolute pseudoprime or a Carmichael number
Note that our little calculator with a limited accuracy and a 10-digit display, in calculating (3²⁸²⁰}₂₈₂₁, has coped with a number having 1346 decimal digits! This has been made possible by the frequent intermediate modulo reduction that the program employs.
Listing for "AN"
Comment Step Code
call program 01 LBL "AN"
decimal point 02 SF 29
store modulus 03 STO 18
get exponent 04 RDN
store exponent 05 STO 17
get base 06 RDN
store base 07 STO 16
constant 08 2
store 2 09 STO 01
take logarithm 10 LN
store log 11 STO 15
display no fractions 12 FIX 0
constant 13 1
store 1 14 STO 14
constant 15 0
store 0 16 STO 02
subroutine for calculating 17 LBL 10
binary representation of n 18 RCL 17
19 LN
20 RCL 15
21 /
add small constant to 22 0.000000001
avoid inaccurate rounding 23 +
24 INT
25 ENTER
26 ST- IND 01
display the binary 27 PSE
exponents of n 28 1
29 ST+ 01
30 RDN
31 STO IND 01
32 2
33 x <> y
34 y^x
35 ST- 17
36 RCL 17
binary representation 37 x = 0?
of n completed? 38 GTO 11
39 GTO 10
subroutine for executing 40 LBL 11
repeated squaring 41 RCL IND 01
42 x = 0?
43 GTO 12
44 RCL16
45 RCL 18
46 2
47 /
48 x < > y
49 x > y?
50 XEQ 16
51 x²
52 RCL 18
take remainder 53 MOD
modulo n 54 STO 16
55 1
56 ST- IND 01
57 GTO 11
subroutine for cal- 58 LBL 12
culating intermediate 59 RCL 16
squaring results 60 ST* 14
61 RCL 14
62 RCL 18
63 MOD
64 STO 14
65 1
66 ST- 01
67 RCL 01
68 2
69 x = y?
70 GTO 13
71 GTO 11
72 LBL 16
73 RCL 18
74 -
75 RTN
Subroutine for 76 LBL 13
recalling and 77 RCL 14
displaying residue 78 CF 29
display residue 79 STOP
ready to start over 80 GTO "AN"
81 RTN
82 END"
BEST!
SlideRule
RE: (41) Exponentiation & Residue Reduction - Werner - 05-08-2020 01:58 PM
Stack-only, without synthetics, 50 bytes
Calculate R := B^X mod M
Code:
L X Y Z T
01>LBL"Z^YMOD" M X B
02 ABS
03 1
04 RDN
05 RDN
06>LBL 02 M X B R
08 2
09 ST/ Y
10 X<> L M X/2 B R
11 X<>Y
12 INT
13 ST- L FP(X) IP(X) M B R
14 RDN
15 X<> L
16 X=0? M FP(X) B R X
17 GTO 00
18 RDN @ R := MOD(R*B,M);
19 STx Y
20 X<>Y M R.B B X
21 LASTX
22 MOD
23 X<> Z
24>LBL 00 M - B R X
25 RDN @ B := MOD(B^2,M);
26 STx X M B^2 R X
27 LASTX
28 MOD M B R X X
29 RUP
30 X!=0? M X B R
31 GTO 02
32 RCL Z
33 ENDCheers, Werner
RE: (41) Exponentiation & Residue Reduction - Albert Chan - 05-11-2020 02:55 PM
(05-08-2020 01:58 PM)Werner Wrote: Stack-only, without synthetics, 50 bytes
Calculate R := B^X mod M ...
Amazing !
I touch up a bit so that it can copy/paste directly into Free42 (52 bytes)
Code:
01▸LBL "Z↑YMOD"
02 ABS
03 1
04 R↓
05 R↓
06▸LBL 02
07 2
08 STO÷ ST Y
09 X<> ST L
10 X<>Y
11 IP
12 STO- ST L
13 R↓
14 X<> ST L
15 X=0?
16 GTO 00
17 R↓
18 STO× ST Y
19 X<>Y
20 LASTX
21 MOD
22 X<> ST Z
23▸LBL 00
24 R↓
25 STO× ST X
26 LASTX
27 MOD
28 R↑
29 X≠0?
30 GTO 02
31 RCL ST Z
32 ENDExample 123^456 mod 789:
123 [Enter] 456 [Enter] 789 [XEQ] "Z↑YMOD" → 699
This version does binary ladder expoentiation from right to left, like this:
Code:
(define (pow-mod b x m r)
(if (fxzero? x) r
(pow-mod [mod (* b b) m] [fx/ x 2] m
[if (even? x) r (mod (* b r) m)] )))scheme> (trace pow-mod)
(pow-mod)
scheme> (pow-mod 123 456 789 1)
|(pow-mod 123 456 789 1)
|(pow-mod 138 228 789 1)
|(pow-mod 108 114 789 1)
|(pow-mod 618 57 789 1)
|(pow-mod 48 28 789 618)
|(pow-mod 726 14 789 618)
|(pow-mod 24 7 789 618)
|(pow-mod 576 3 789 630)
|(pow-mod 396 1 789 729)
|(pow-mod 594 0 789 699)
|699
699
RE: (41) Exponentiation & Residue Reduction - Thomas Klemm - 04-26-2022 03:17 AM
Let us start with the following Python program:
Code:
def pow(b, x, m):
r = 1
while x:
if x % 2:
r = b * r % m
b = b ** 2 % m
x = IP(x / 2)
return rWith the help of the Python to RPN - source code converter it is translated to:
Code:
LBL "Z^Y%X" // Label. Identify programs and routines for execution and branching. Parameter: local or global label.
RTN
LBL A // def pow
XEQ 66 // reorder 3 params for storage
STO 00 // param: b
RDN
STO 01 // param: x
RDN
STO 02 // param: m
RDN
1
STO 03 // r
LBL 00 // while
RCL 01 // x
X≠0? // while true?
GTO 01 // while body
GTO 02 // resume
LBL 01 // while body
RCL 01 // x
2
MOD
X≠0? // if true?
GTO 03 // if body
GTO 04 // resume
LBL 03 // if body
RCL 00 // b
RCL 03 // r
*
RCL 02 // m
MOD
STO 03 // r
LBL 04 // resume
RCL 00 // b
X↑2
RCL 02 // m
MOD
STO 00 // b
RCL 01 // x
2
/
IP // Returns the integer part of x
STO 01 // x
GTO 00 // while
LBL 02 // resume
RCL 03 // r
RTN // return
LBL 50 // PyRPN Support Library of
"-Utility Funcs-"
RTN // ---------------------------
LBL 66 // Reverse params. (a,b,c) -> (c,b,a)
X<>Y
RDN
RDN
X<>Y
RDN
RTNLet's clear it up a bit:
- remove local label A
- remove reordering of params for storage
- remove empty support Library
- invert logic in condition
Instead of:
Code:
X≠0? // while true?
GTO 01 // while body
GTO 02 // resume
LBL 01 // while bodyCode:
X=0? // while not?
GTO 02 // resumeWith these changes we are already down at 51 bytes and could call it a day:
Code:
00 { 51-Byte Prgm }
01▸LBL "Z↑Y%X"
02 STO 02
03 R↓
04 STO 01
05 R↓
06 STO 00
07 1
08 STO 03
09▸LBL 00
10 RCL 01
11 X=0?
12 GTO 02
13 RCL 01
14 2
15 MOD
16 X=0?
17 GTO 04
18 RCL 00
19 RCL 03
20 ×
21 RCL 02
22 MOD
23 STO 03
24▸LBL 04
25 RCL 00
26 X↑2
27 RCL 02
28 MOD
29 STO 00
30 RCL 01
31 2
32 ÷
33 IP
34 STO 01
35 GTO 00
36▸LBL 02
37 RCL 03
38 ENDBut can we do better?
You may have noticed that line 13 could be removed since x is still present from line 10.
Furthermore, division by 2 and modulo 2 of x can be combined.
Code:
00 { 47-Byte Prgm }
01▸LBL "Z↑Y%X"
02 STO 02
03 R↓
04 STO 01
05 R↓
06 STO 00
07 1
08 STO 03
09▸LBL 00
10 RCL 01
11 X=0?
12 GTO 02
13 2
14 ÷
15 ENTER
16 IP
17 STO 01
18 X=Y?
19 GTO 04
20 RCL 00
21 RCL 03
22 ×
23 RCL 02
24 MOD
25 STO 03
26▸LBL 04
27 RCL 00
28 X↑2
29 RCL 02
30 MOD
31 STO 00
32 GTO 00
33▸LBL 02
34 RCL 03
35 ENDSo far we haven't bothered with registers or keeping results on the stack.
But if we leave x on the stack, we save 2 STO and 1 RCL instructions at the cost of 2 RDN instructions.
With this we not only save a register, but also another byte:
Code:
00 { 46-Byte Prgm }
01▸LBL "Z↑Y%X" ; b x m -- b^x (mod m)
02 STO 00 ; m
03 1
04 STO 01 ; r
05 R↑
06 STO 02 ; b
07 R↑ ; x
08▸LBL 00 ; while
09 X=0? ; x == 0
10 GTO 02 ; end while
11 2
12 ÷ ; x / 2
13 ENTER
14 IP ; x // 2
15 X=Y? ; x is even
16 GTO 01 ; end if
17 RCL 01 ; r
18 RCL 02 ; b
19 ×
20 RCL 00 ; m
21 MOD ; r × b % m
22 STO 01 ; r
23 R↓ ; x
24▸LBL 01 ; end if
25 RCL 02 ; b
26 X↑2
27 RCL 00 ; m
28 MOD ; b^2 % m
29 STO 02 ; b
30 R↓ ; x
31 GTO 00 ; while
32▸LBL 02 ; end while
33 RCL 01 ; r
34 ENDFeel free to replace the 3 remaining registers with the synthetic registers M, N and O.
Errata: (41) Exponentiation & Residue Reduction - Thomas Klemm - 04-26-2022 10:59 AM
Errata
342 Appendix
This is the binary decomposition of the exponent 340.
Check: \(2^8 + 2^6 + 2^4 + 2^2 = 340\). Check!
Next, the computer will start the necessary repeated squarings and
reductions modulo 341, which will take about 7 seconds. Then the display
will show the end result, a 1 (without a decimal point!).
Thus,
\(
\left< 2^{340} \right>_{341} = 1
\).
In other words, 341 is a pseudoprime to the base 2.
Similarly, after pressing 2, ENTER, 170, ENTER, 341, R/S, the display
shows in succession:
7. 5. 3. 1. 1
where that last display (the 1 without the decimal point) tells us that
\(
\left< 2^{170} \right>_{341} = 1
\).
Proceeding in the same manner, we find
\(
\left< 2^{85} \right>_{341} = 32
\),
i.e., 341 is not a strong pseudoprime to the base 2. Also, by using 3 as a base
we find
\(
\left< 3^{340} \right>_{341} = 56
\).
Thus, 341 is certainly not an absolute pseudoprime (Carmichael number).
However, for the modulus 2821, we find
\(
\left< 2^{2820} \right>_{2821} = \left< 3^{2820} \right>_{2821} = 1
\),
two of the many steps necessary to show that 2821 is an absolute pseudoprime
or a Carmichael number.
Note that our little calculator with a limited accuracy and a 10-digit
display, in calculating \(\left< 2^{2820} \right>_{2821}\) has coped with a number having 1346
decimal digits! This has been made possible by the frequent intermediate
modulo reduction that the program employs.
Appendix 343
Original program for the HP-41C:
Code:
01 LBL "AN" ; call program
02 SF 29 ; decimal point
03 STO 18 ; store modulus
04 RDN ; get exponent
05 STO 17 ; store exponent
06 RDN ; get base
07 STO 16 ; store base
08 2 ; constant
09 STO 01 ; store 2
10 LN ; take logarithm
11 STO 15 ; store log
12 FIX 0 ; display no fractions
13 1 ; constant
14 STO 14 ; store 1
15 0 ; constant
16 STO 02 ; store 0
17 LBL 10 ; subroutine for calculating
18 RCL 17 ; binary representation of ?
19 LN
20 RCL 15
21 /
22 0.000000001 ; add small constant to
23 + ; avoid inaccurate rounding
24 INT
25 ENTER
26 ST- IND 01
27 PSE ; display the binary
28 1 ; exponents of ?
29 ST+ 01
30 RDN
31 STO IND 01
32 2
33 X<>Y
34 Y^X
35 ST- 17
36 RCL 17
37 X=0? ; binary representation
38 GTO 11 ; of ? completed?
39 GTO 10
40 LBL 11 ; subroutine for executing
41 RCL IND 01 ; repeated squaring
42 X=0?
43 GTO 12
44 RCL 16
45 RCL 18
46 2
47 /
48 X<>Y
49 X>Y?
50 XEQ 16
51 X^2
52 RCL 18
53 MOD ; take remainder
54 STO 16 ; modulo ?
55 1
56 ST- IND 01
57 GTO 11
58 LBL 12 ; subroutine for cal-
59 RCL 16 ; culating intermediate
60 ST* 14 ; squaring results
61 RCL 14
62 RCL 18
63 MOD
64 STO 14
65 1
66 ST- 01
67 RCL 01
68 2
69 X=Y?
70 GTO 13
71 GTO 11
72 LBL 16
73 RCL 18
74 -
75 RTN
76 LBL 13 ; Subroutine for
77 RCL 14 ; recalling and
78 CF 29 ; displaying residue
79 STOP ; display residue
80 GTO "AN" ; ready to start over
81 RTN
82 ENDTranslated program for the HP-42S:
Code:
00 { 141-Byte Prgm }
01▸LBL "AN"
02 SF 29
03 STO 18
04 R↓
05 STO 17
06 R↓
07 STO 16
08 2
09 STO 01
10 LN
11 STO 15
12 FIX 00
13 1
14 STO 14
15 0
16 STO 02
17▸LBL 10
18 RCL 17
19 LN
20 RCL 15
21 ÷
22 0.000000001
23 +
24 IP
25 ENTER
26 STO- IND 01
27 PSE
28 1
29 STO+ 01
30 R↓
31 STO IND 01
32 2
33 X<>Y
34 Y↑X
35 STO- 17
36 RCL 17
37 X=0?
38 GTO 11
39 GTO 10
40▸LBL 11
41 RCL IND 01
42 X=0?
43 GTO 12
44 RCL 16
45 RCL 18
46 2
47 ÷
48 X<>Y
49 X>Y?
50 XEQ 16
51 X↑2
52 RCL 18
53 MOD
54 STO 16
55 1
56 STO- IND 01
57 GTO 11
58▸LBL 12
59 RCL 16
60 STO× 14
61 RCL 14
62 RCL 18
63 MOD
64 STO 14
65 1
66 STO- 01
67 RCL 01
68 2
69 X=Y?
70 GTO 13
71 GTO 11
72▸LBL 16
73 RCL 18
74 -
75 RTN
76▸LBL 13
77 RCL 14
78 CF 29
79 STOP
80 GTO "AN"
81 RTN
82 ENDResources
Errata: (41) Exponentiation & Residue Reduction - Thomas Klemm - 04-26-2022 04:05 PM
(04-06-2020 01:31 PM)SlideRule Wrote: Proceeding in the same manner, we findThat sentence didn't make any sense to me, so I looked it up in the original:
{2⁸⁵}₃₄₁ = 1 .
i.e., 341 is not a strong pseudoprime to the base 2.
(04-26-2022 10:59 AM)Thomas Klemm Wrote: Proceeding in the same manner, we findThis makes a lot more sense.
\(
\left< 2^{85} \right>_{341} = 32
\),
i.e., 341 is not a strong pseudoprime to the base 2.
Most of the original PDF can be copied, but among other things, the formula from above is missing and so had to be pasted in manually by the OP.
![[Image: duty_calls.png]](https://imgs.xkcd.com/comics/duty_calls.png)
I hope these fixes will make the article easier to read and also allow you to run the original program.
For those interested in the topic:
- Miller-Rabin Primality Test for the HP-48
- (48G) Miller-Rabin Primality Test
- [HP35s] Program for prime number (Wheel Sieve and Miller-Rabin)
- Primes & Number Theoretic Functions for the HP-41
- Miller-Rabin Primality Test for the HP-42S