(41) Exponentiation & Residue Reduction - Printable Version +- 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 END Cheers, 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 END Example 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 r With 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 RTN Let's clear it up a bit: remove local label A remove reordering of params for storage remove empty support Library invert logic in condition In the last step, the condition is inverted. Instead of: Code: X≠0?            // while true? GTO 01          // while body GTO 02          // resume LBL 01          // while body We use: Code: X=0?            // while not? GTO 02          // resume This saves us both a LBL and GTO instruction. With 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 END We also saved a byte by choosing a shorter name. But 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 END So 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 END Feel 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 END Translated 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 END Resources 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 find     {2⁸⁵}₃₄₁ = 1 . i.e., 341 is not a strong pseudoprime to the base 2. That sentence didn't make any sense to me, so I looked it up in the original: (04-26-2022 10:59 AM)Thomas Klemm Wrote:  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. This makes a lot more sense. 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. 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: