hp41programs

primes Primes & Number Theoretic Functions for the HP-41

Overview

1°) Factorization

   a) Primes, Divisor Functions, Euler Function
   b) Primes, Euler Function, Moebius Function - program#1
   c) Primes, Euler Function, Moebius Function - program#2
   d) Pollard rho factorization

2°) Primality Tests

   a) Program#1
   b) M-Code Routine
   c) a-Strong Probable Primes
   d) Other Primality tests
   e) Next Prime

3°) Divisor Functions

   a) Program#1
   b) Program#2
   c) Program#3 ( multiprecision , n < 100000 )
   d) List of the Divisors

4°) Moebius Function

a) Program#1
b) Program#2

5°) Euler Totient Function

1°) Factorization

a) Primes, Divisor Functions, Euler Function

-The following program displays the factorization of any integer n ( 1 < n < E10 )
and returns s_k(n) = the sum of the k-th powers of the divisors of n is in X-register and in R06.
and phi(n) = the Euler function in Y-register and in R04.

-phi(n) is the number of integers not exceeding and relatively prime to n.

Registers: R00 , R07 , R08: temp

                           R01 = 2 ; R02 = 4 ; R03 = 6     ( these numbers are stored in data registers because digit entry is very slow on the HP-41 )
                           R04 = phi(N)
                           R05 = k ; R06 = s_k(n)
Flag: F29
Subroutines: /

-The append character is denoted ~

01 LBL "PRF"
02 SIGN
03 STO 06
04 X<>Y
05 STO 05
06 LASTX
07 ENTER^
08 STO 04
09 6
10 STO 03
11 4
12 STO 02
13 -
14 STO 01
15 FIX 0
16 CF 29
17 CLA
18 ARCL Y
19 "~="
20 MOD
21 X=0?
22 XEQ 02
23 CLX
24 3
25 MOD
26 X=0?
27 XEQ 02
28 CLX
29 5
30 MOD
31 X=0?

32 XEQ 02
33 SIGN
34 SIGN
35 LBL 01
36 X<> L
37 RCL 03
38 +
39 MOD
40 X=0?
41 XEQ 02
42 X<> L
43 RCL 02
44 +
45 MOD
46 X=0?
47 XEQ 02
48 X<> L
49 RCL 01
50 +
51 MOD
52 X=0?
53 XEQ 02
54 X<> L
55 RCL 02
56 +
57 MOD
58 X=0?
59 XEQ 02
60 X<> L
61 RCL 01
62 +

63 MOD
64 X=0?
65 XEQ 02
66 X<> L
67 RCL 02
68 +
69 MOD
70 X=0?
71 XEQ 02
72 X<> L
73 RCL 03
74 +
75 MOD
76 X=0?
77 XEQ 02
78 X<> L
79 RCL 01
80 +
81 MOD
82 X=0?
83 XEQ 02
84 X<> L
85 X^2
86 X<Y?
87 GTO 01
88 SIGN
89 X=Y?
90 GTO 05
91 X<>Y
92 ARCL X
93 "~^1"

94 ST/ 04
95 DSE X
96 ST* 04
97 R^
98 RCL 05
99 Y^X
100 1
101 +
102 ST* 06
103 AVIEW
104 GTO 05
105 LBL 02
106 STO 00
107 LBL 03
108 ISG 00
108 CLX
110 X<> L
111 ST/ Y
112 ST/ Z
113 ST/ T
114 MOD
115 X=0?
116 GTO 03
117 ARCL L
118 "~^"
119 ARCL 00
120 X<> L
121 ST/ 04
122 SIGN
123 X#Y?
124 "~*"

125 CLX
126 LASTX
127 DSE X
128 ST* 04
129 X<> L
130 STO 07
131 RCL 05
132 Y^X
133 SIGN
134 STO 08
135 LBL 04
136 LASTX
137 *
138 ST+ 08
139 DSE 00
140 GTO 04
141 X<> 08
142 ST* 06
143 X<> 07
144 SIGN
145 AVIEW
146 RTN
147 LBL 05
148 FIX 4
149 SF 29
150 RCL 04
151 RCL 06
152 END

( 232 bytes / SIZE 009 )

STACK	INPUTS	OUTPUTS
Y	k	phi(n)
X	n	s_k(n)

Example1:

1 ENTER^
3238704 XEQ "PRF" displays ( successively )

           3238704=2^4*
           3238704=2^4*3^5*
           3238704=2^4*3^5*7*2*
           3238704=2^4*3*5*7^2*17^1        ( if there are more than 24 characters, the left part of the alpha string will be gradually truncated )

s₁(3238704) = 11577384 in X-register and in R06
and phi(3238704) = 870912 in Y-register and in R04

Example2:

7 ENTER^
24 R/S produces 24^1=2^3*3*1

s₇(24) = 4624699020 in registers X and R06
and phi(24) = 8 in registers Y and R04

Example3:

0 ENTER^
999983 R/S yields 999983=999983^1 ( in 42 seconds ) thus, 999983 is prime

s₀(999983) = 2 in registers X and R04
phi(999983) = 999982 in registers Y and R06

Notes:

-If N is a prime, execution time is approximately N^1/2 / 25 seconds.
-The greatest prime < E10 is 9999999967 ( 1h06m with an HP-41CX )
-If you set flag F21, the program will stop at each AVIEW.
-s₀(n) is the number of divisors of n.
-s₁(n) is the sum of the divisors of n.

b) Primes, Euler Function, Moebius Function - program#1

-"PRF" has some drawbacks: the factors are lost, and we must execute the routine several times if we have to calculate s_k(n) for several k-values.
-The following variant stores the factors into contiguous registers from register R07 and it computes Moebius function µ(n).

µ(0) = 0 µ(1) = 1

otherwise,

µ(n) = 0 if n is a multiple of the square of a prime
µ(n) = (-1)^k if n is the product of k distinct primes.

-The factors themselves are stored as follows: if, for instance n = 3⁵ 7⁴ R07 = 3.005 R08 = 7.004
-In other words, the exponents are divided by 1000 and added to the corresponding prime factor p.
-If p > 9999999 , the fractional part will be zero but should be read 0.001

Registers: R00 to R03: temp
R04 = phi(N) R05 = µ(n) R06 = 7.eee R07 , R08 , ... the factors
Flags: /
Subroutines: /

01 LBL "PMF"
02 ENTER^
03 STO 04
04 1
05 X=Y?
06 RTN
07 STO 05
08 CLX
09 X=Y?
10 RTN
11 CLX
12 6
13 STO 03
14 STO 06
15 4
16 STO 02
17 -
18 STO 01
19 MOD
20 X=0?
21 XEQ 02
22 CLX
23 3
24 MOD
25 X=0?
26 XEQ 02
27 CLX
28 5
29 MOD
30 X=0?

31 XEQ 02
32 SIGN
33 SIGN
34 LBL 01
35 X<> L
36 RCL 03
37 +
38 MOD
39 X=0?
40 XEQ 02
41 X<> L
42 RCL 02
43 +
44 MOD
45 X=0?
46 XEQ 02
47 X<> L
48 RCL 01
49 +
50 MOD
51 X=0?
52 XEQ 02
53 X<> L
54 RCL 02
55 +
56 MOD
57 X=0?
58 XEQ 02
59 X<> L
60 RCL 01

61 +
62 MOD
63 X=0?
64 XEQ 02
65 X<> L
66 RCL 02
67 +
68 MOD
69 X=0?
70 XEQ 02
71 X<> L
72 RCL 03
73 +
74 MOD
75 X=0?
76 XEQ 02
77 X<> L
78 RCL 01
79 +
80 MOD
81 X=0?
82 XEQ 02
83 X<> L
84 X^2
85 X<Y?
86 GTO 01
87 SIGN
88 X=Y?
89 GTO 04
90 ST+ 06

91 RDN
92 ST/ 04
93 DSE X
94 ST* 04
95 CLX
96 E-3
97 +
98 STO IND 06
99 RCL 05
100 CHS
101 STO 05
102 GTO 04
103 LBL 02
104 STO 00
105 ISG 06
106 LBL 03
107 ISG 00
108 CLX
109 X<> L
110 ST/ Y
111 ST/ Z
112 ST/ T
113 MOD
114 X=0?
115 GTO 03
116 CLX
117 LASTX
118 ST/ 04
119 DSE X
120 ST* 04

121 X<> L
122 RCL 00
123 E3
124 /
125 +
126 STO IND 06
127 INT
128 RCL 05
129 CHS
130 DSE 00
131 CLX
132 STO 05
133 CLX
134 +
135 SIGN
136 RTN
137 LBL 04
138 RCL 05
139 RCL 04
140 RCL 06
141 E3
142 /
143 7
144 +
145 STO 06
146 END

( 210 bytes )

STACK	INPUTS	OUTPUTS
Z	/	µ(n)
Y	/	phi(n)
X	n	7.eee

Exception: if n = 0 or 1 , X-output = 0 or 1 respectively

Example1:

n = 3238704 XEQ "PMF" >>>>       7.010                                    --- Execution time = 9s ---
                                               RDN    phi(n) = 870912
                                               RDN     µ(n) =   0

{ R07 R08 R09 R10 } = { 2.004 3.005 7.002 17.001 } whence 3238704 = 2⁴ 3⁵ 7² 17

Example2:

n =   999983      R/S           >>>>        7.007                                    --- Execution time = 42s ---
                                             RDN    phi(n) = 999982
                                             RDN     µ(n) = -1

R07 = 999983.001 thus 999983 is a prime

Notes:

s_k(n) may then be computed by one of the programs listed in paragraph 3
SIZE 016 is always enough

The initial number may be recomputed from the control number of the factorization by this short routine:

01 LBL "F-N"
02 1
03 LBL 01
04 RCL IND Y
05 FRC
06 E3
07 *
08 RCL IND Z
09 INT
10 X<>Y
11 X=0?
12 ISG X
13 INT
14 Y^X
15 *
16 ISG Y
17 GTO 01
18 END

( 32 bytes )

STACK	INPUT	OUTPUT
X	ddd.eee	n

-With the example above,

7.010 XEQ "F-N" >>>> n = 3238704

c) Primes, Euler Function, Moebius Function - program#2

-A faster program may of course be written if you use an M-code routine ( cf §2b )

Registers: R00: temp

R01 = phi(N) R02 = µ(n) R03 = 4.eee R04 , R05 , ... the factors
Flags: /
Subroutines: /

01 LBL "PMF2"
02 ENTER^
03 ENTER^
04 STO 01
05 3
06 STO 03
07 SIGN
08 STO 02
09 X=Y?
10 RTN
11 *
12 X=0?
13 RTN
14 GTO 03
15 LBL 01
16 CLX

17 STO 00
18 ISG 03
19 LBL 02
20 ISG 00
21 CLX
22 X<> L
23 ST/ Y
24 ST/ Z
25 ST/ T
26 MOD
27 X=0?
28 GTO 02
29 CLX
30 LASTX
31 ST/ 01
32 DSE X

33 ST* 01
34 X<> L
35 RCL 00
36 E3
37 /
38 +
39 STO IND 03
40 INT
41 RCL 02
42 CHS
43 DSE 00
44 CLX
45 STO 02
46 +
47 CLX
48 +

49 LBL 03
50 PR?
51 X<=0?
52 GTO 01
53 X=Y?
54 GTO 04
55 ST+ 03
56 RDN
57 ST/ 01
58 DSE X
59 ST* 01
60 CLX
61 E-3
62 +
63 STO IND 03
64 RCL 02

65 CHS
66 STO 02
67 LBL 04
68 RCL 02
69 RCL 01
70 RCL 03
71 E3
72 /
73 4
74 +
75 STO 03
76 END

( 111 bytes )

STACK	INPUTS	OUTPUTS
Z	/	µ(n)
Y	/	phi(n)
X	n	4.eee

Exception: if n = 0 or 1 , X-output = 0 or 1 respectively

Example:

n = 3238704 XEQ "PMF2" >>>>   4.007                                    --- Execution time = 8s ---
                                                RDN    phi(n) = 870912
                                                RDN     µ(n) =   0

{ R04 R05 R06 R07 } = { 2.004 3.005 7.002 17.001 } whence 3238704 = 2⁴ 3⁵ 7² 17

d) Pollard rho factorization

-This method is a probabilistic algorithm to find one factor of an integer n:

• Let x₀ = y₀ = a random integer between 0 and n-2 and f(x) = x² + c ( c # 0 and # -2 , c = 1 in the routine below )

• Compute x_k = f(x_k-1) , y_k = f [f(y_k-1)] and d = gcd ( | x_k - y_k | , n )
Until d > 1

• If d < n we have found a non-trivial divisor
• If d = n failure: execute "POLL" again!

Registers: • R00 = random number ( Register R00 is to be initialized before executing "POLL" )

R01 = n R02 = x_k R03 = y_k
Flags: /
Subroutine: "SQM" ( cf "Modular Arithmetic for the HP-41" )

01 LBL "POLL"
02 STO 01
03 DSE X
04 RCL 00
05 R-D
06 FRC
07 STO 00
08 *
09 INT

10 STO 02
11 STO 03
12 LBL 01
13 RCL 02
14 RCL 01
15 XEQ "SQM"
16 1
17 +
18 STO 02

19 RCL 03
20 RCL 01
21 XEQ "SQM"
22 1
23 +
24 RCL 01
25 XEQ "SQM"
26 1
27 +

28 STO 03
29 RCL 02
30 -
31 ABS
32 RCL 01
33 LBL 02
34 MOD
35 LASTX
36 X<>Y

37 X#0?
38 GTO 02
39 SIGN
40 X<>Y
41 X=Y?
42 GTO 01
43 END

( 67 bytes / SIZE 004 )

STACK	INPUT	OUTPUT
X	n	d

if d < n success!
if d = n failure!

Example: Find a factor of n = 9976979

If you store 2 into R00

2 STO 00
9976979 XEQ "POLL" >>>> d = 997 --- Execution time = 2mn06s ---

-Thus, we've found a non-trivial factor... but not very quickly!
-According to the random seed in R00, the execution time seems to fall between 73s and 3m44s with n = 9976979
-A large factor is sometimes found more quickly than with "PRF" but not always.
-This routine is essentially given for completeness.

2°) Primality Tests

a) Program#1

-"PR?" takes a positive integer n in X-register,
and returns 1 if n is a prime or if n = 1 and 0 otherwise
-The smallest divisor d is in L-register.

-"PR?" and "PMF" use the same method.

Data Registers: R00 = 2 R01 = 4 R02 = 6
Flags: /
Subroutines: /

01 LBL "PR?"
02 STO Y
03 6
04 STO 02
05 4
06 STO 01
07 -
08 STO 00
09 MOD
10 X=0?
11 GTO 02
12 CLX
13 3
14 MOD
15 X=0?
16 GTO 02
17 CLX
18 5

19 MOD
20 X=0?
21 GTO 02
22 SIGN
23 SIGN
24 LBL 01
25 X<> L
26 RCL 02
27 +
28 MOD
29 X=0?
30 GTO 02
31 X<> L
32 RCL 01
33 +
34 MOD
35 X=0?
36 GTO 02

37 X<> L
38 RCL 00
39 +
40 MOD
41 X=0?
42 GTO 02
43 X<> L
44 RCL 01
45 +
46 MOD
47 X=0?
48 GTO 02
49 X<> L
50 RCL 00
51 +
52 MOD
53 X=0?
54 GTO 02

55 X<> L
56 RCL 01
57 +
58 MOD
59 X=0?
60 GTO 02
61 X<> L
62 RCL 02
63 +
64 MOD
65 X=0?
66 GTO 02
67 X<> L
68 RCL 00
69 +
70 MOD
71 X=0?
72 GTO 02

73 X<> L
74 X^2
75 X<Y?
76 GTO 01
77 RDN
78 SIGN
79 LBL 02
80 CLX
81 LASTX
82 X<Y?
83 CLX
84 X=Y?
85 SIGN
86 END

( 116 bytes / SIZE 003 )

STACK	INPUTS	OUTPUTS
Y	/	n
X	n	1 or 0
L	/	d

Example:

n =     999983       XEQ "PR?" >>>>   1      so 999983 is a prime                    --- Execution time =   39s   ---
n = 9999865009        R/S         >>>>   0      so n is not a prime                       --- Execution time = 6mn47s ---
                              and the smallest non-trivial divisor = L = 10007

b) M-Code Routine

-"PR?" takes an integer n in X-register and returns:

1 in X and n in L if n is a prime or if n = 1
0 in X and the smallest positive non-trivial divisor d in L otherwise

-However, if n > 9999999999 X-output = -1 and L-output = 2

and if n = 0 , X-output = L-output = 0

-Two other M-code routines are used to format the positive integers < 10¹⁰ as needed in CPU-register C
and to perform the reverse operation at the end
-Both assume that the CPU is in decimal mode.

Format

006    A=0 S&X             at the address F723 in my ROM
39C   PT=0
1A2   A=A-1 @PT
366   ?A#C S&X
01B   JNC +03
3DA RSHFC M
3E3   JNC -04
046   C=0 S&X
3E0   RTN

For instance, if C = 1234 = 1.234 E3 , this number is coded as follows

MS| --------------- Mantissa -------------- | --S&X----

After executing this subroutine, CPU register C contains:

MS| --------------- Mantissa -------------- | --S&X----

Reverse Operation

2FA ?C#0 M                at the address F753 in my ROM
3A0 ?NC RTN
130   LDI S&X
009   009
35C PT=12
11A A=C M
342   ?A#0 @PT
027   JC +04
266   C=C-1 S&X
3FA LSHFA M
3E3   JNC -04
0BA A<>C M
3E0   RTN

M-Code Routine "PR?"

>>>> All the words written in red are to be changed according to the addresses of the subroutines sub1 sub2 sub3 sub4
and the 2 subroutines above in your own ROM

0BF   "?"
012   "R"
010   "P"
0F8   READ3 (X)     the first executable word is at the address FA10 in my ROM
361   ?NCXQ          checks for alpha data
050   14D8              and execute SETDEC
05E   C=0 MS          C= | C |
088   SETF5                C
0ED ?NCXQ              =
064   193B               INT(C)
0E8   WRIT3 (X)
128   WRIT4 (L)
2EE   ?C#0 ALL
3A0   ?NC RTN
10E   A=C ALL
130   LDI S&X
010   010d
306   ?A<C S&X
057   JC +10d
04E   C=0 ALL           C
35C   PT=12               =
090   LD@PT- 2         2
128   WRIT4 (L)
2BE   C=-C-1 MS      C
35C   PT=12               =
050    LD@PT- 1       -1
0E8   WRIT3 (X)
3E0   RTN
0F8   READ3 (X)
00E   A=0 ALL           A
35C   PT=12               =
162   A=A+1 @PT     1
36E   ?A#C ALL
3A0   ?NC RTN
2F9   ?NCXQ             C=
060    18BE              sqrt(C)
088   SETF5                C
0ED ?NCXQ              =
064   193B               INT(C)
08D ?NCXQ           ?ncxq
3DC F723                format
33C   RCR 1
158   M=C ALL
0F8   READ3 (X)
08D ?NCXQ           ?ncxq
3DC F723                format
0E8   WRIT3 (X)
2DC PT=13
3D4   PT=PT-1
2E2   ?C#0 @PT
3F3   JNC -02
04E   C=0 ALL
090   LD@PT- 2
10E   A=C ALL
128   WRIT4 (L)
3DC PT=PT+1
3DC PT=PT+1
21D   ?NCXQ          ?ncxq
3E8   FA87                sub4
3D4   PT=PT-1
00E   A=0 ALL
162   A=A+1 @PT
3DC PT=PT+1
138   READ4 (L)
211   ?NCXQ           ?ncxq
3E8   FA84                sub3
06E   A<>B ALL
138   READ4 (L)
211   ?NCXQ           ?ncxq
3E8   FA84                sub3
3D4   PT=PT-1
00E   A=0 ALL
162   A=A+1 @PT
3DC PT=PT+1
013   JNC +02
06E   A<>B ALL      ----------loop-----------
138   READ4 (L)
209   ?NCXQ           ?ncxq
3E8   FA82                sub1
06E   A<>B ALL
138   READ4 (L)
20D   ?NCXQ           ?ncxq
3E8   FA83                sub2
06E   A<>B ALL
138   READ4 (L)
211   ?NCXQ           ?ncxq
3E8   FA84                sub3
06E   A<>B ALL
138   READ4 (L)
20D   ?NCXQ           ?ncxq
3E8   FA83                sub2
06E   A<>B ALL
138   READ4 (L)
211   ?NCXQ           ?ncxq
3E8   FA84                sub3
06E   A<>B ALL
138   READ4 (L)
20D   ?NCXQ           ?ncxq
3E8   FA83                sub2
06E   A<>B ALL
138   READ4 (L)
209   ?NCXQ           ?ncxq
3E8   FA82                sub1
06E   A<>B ALL
138   READ4 (L)
211   ?NCXQ           ?ncxq
3E8   FA84                sub3
10E   A=C ALL
198   C=M ALL
30E   ?A<C ALL          if potential divisor < sqrt(n)
2EF   JC -35d=-23h     goto loop
0F8   READ3 (X)         otherwise
0BB JNC +23d=17h    n is a prime
38E   RSHFA ALL
3CE RSHFC ALL
128   WRIT4 (L)
033   JNC +06
14E   A=A+C ALL     sub1              at the address FA82 in my ROM
14E   A=A+C ALL     sub2              at the address FA83 in my ROM
14E   A=A+C ALL     sub3              at the address FA83 in my ROM
342   ?A#0 @PT
3C7   JC -08
08E   B=A ALL          sub4              at the address FA87 in my ROM
0F8   READ3 (X)
0AE A<>C ALL
1CE A=A-C ALL      these lines calculate A mod C
3FB JNC -01
14E A=A+C ALL
3CE RSHFC ALL
2E6   ?C#0 S&X
3DB JNC -05
34E   ?A#0 ALL
360   ?C RTN             if A mod C # 0 we've not found a divisor
020   XQ->GO
2FC RCR 13
10E   A=C ALL
0F8   READ3 (X)
0EE   C<>B ALL
04E   C=0 ALL
32E   ?A<B ALL
01F   JC +03
35C PT=12
050   LD@PT- 1
0E8   WRIT3 (X)
0AE A<>C ALL
14D ?NCXQ               the divisor is restored
3DC F753                   in the standard format
128   WRIT4 (L)
3E0   RTN

STACK	INPUTS	OUTPUTS
X	n	1 or 0 or -1
L	/	n or d or 2

X-output = 1   if n is a prime or if n = 1       L-output = 1   in this case
X-output = 0   otherwise if n < E10            L-output = divisor
X-output = -1 if n > 9999999999              L-output = 2

Examples:

99999999967 XEQ "PR?" gives 1 in X and 9999999967 in L --- Execution time = 7mn09s ---
100140049 XEQ "PR?" gives 0 in X and 10007 in L --- Execution time = 39s ---

Remarks:

PR? first takes the absolute value and the integer part of X-input.
It returns -1 if n > 9999999999 because if n is the result of a calculation, we don't really know if n = 10¹⁰ or 1+10¹⁰ ...
So, you can decide how to take into account d = 2 in L-register according to your preference.

The stack registers Y Z T and the synthetic registers M N O P Q are unchanged.
CPU register N is also unused.

c) a-Strong Probable Primes

-"SPP?" checks if an odd integer n > 2 satisfies Miller-Rabin primality test: n is written n = 1 + q 2^r where q is odd

• Let a > 1 , n is called "strong probable prime to base a" or "a-strong probable prime"

if a^q = 1 mod n
or a^{q.2^i} = -1 mod n for an integer i ( 0 <= i < r )

Data Registers: R00 thru R02 are used by "M^" R03 = a R04 = n R05 = r , r-1 , ..........
Flags: /
Subroutines: "SQM" "M^" ( cf "Modular Arithmetic for the HP-41" )

01 LBL "SPP?"
02 STO 03
03 X<>Y
04 STO 01
05 STO 04
06 1
07 X=Y?
08 RTN
09 ST+ X
10 X=Y?
11 GTO 03
12 MOD

13 X=0?
14 RTN
15 ST- 01
16 CLX
17 STO 05
18 GTO 00
19 LBL 01
20 LASTX
21 ST/ 01
22 ISG 05
23 LBL 00
24 RCL 01

25 LASTX
26 MOD
27 X=0?
28 GTO 01
29 RCL 03
30 RCL 01
31 RCL 04
32 XEQ "M^"
33 1
34 X=Y?
35 RTN
36 X<>Y

37 GTO 00
38 LBL 02
39 RCL 04
40 XEQ "SQM"
41 LBL 00
42 X=0?
43 RTN
44 LASTX
45 DSE X
46 X=Y?
47 GTO 03
48 SIGN

49 X=Y?
50 GTO 00
51 X<>Y
52 DSE 05
53 GTO 02
54 LBL 00
55 CLX
56 RTN
57 LBL 03
58 SIGN
59 END

( 88 bytes / SIZE 006 )

STACK	INPUTS	OUTPUTS
Y	n	/
X	a > 1	1 or 0

n must be non-negative

Examples:

9999999967 ENTER^
2 XEQ "SPP?" >>>> 1 --- Execution time = 1m20s ---

9998200081 ENTER^
3 XEQ "SPP?" >>>> 0 --- Execution time = 1m18s ---

-So, 9999999967 is 2-strong probable prime but we are not sure that it's a prime!
9998200081 is not 3-strong probable prime whence 9998200081 is certainly composite ( in fact 9998200081 = 99991² )

Remarks:

"SPP?" also returns 1 if n = 1 or 2
and 0 if n is even # 2

d) Other Primality Tests

-We can use Miller-Rabin test for several a-values to be sure that n is a prime, for instance:

• If n < 10¹² and if n is strong probable prime to bases 2 , 13 , 23 and 1662803 then n is a prime

-Thus, this criterion could also be applied with a calculator working with 12 digits.

Data Registers: R00 thru R05: temp ( R04 = n )
Flags: /
Subroutine: "SPP?"

01 LBL "PR2?"
02 2
03 XEQ "SPP?"
04 X=0?
05 RTN
06 RCL 04

07 13
08 X<Y?
09 X=0?
10 GTO 00
11 XEQ "SPP?"
12 X=0?

13 RTN
14 RCL 04
15 23
16 X<Y?
17 X=0?
18 GTO 00

19 XEQ "SPP?"
20 X=0?
21 RTN
22 RCL 04
23 1662803
24 X#Y?

25 XEQ "SPP?"
26 LBL 00
27 X#0?
28 SIGN
29 END

( 68 bytes / SIZE 006 )

STACK	INPUT	OUTPUT
X	n >= 0	1 or 0

X-output = 1 if n is a prime or if n = 1
X-output = 0 otherwise

Examples:

9999999967 XEQ "PR2?" >>>> 1 --- Execution time = 5m28s ---

9998200081 XEQ "PR2?" >>>> 0 --- Execution time = 1m18s ---

1234567899 XEQ "PR2?" >>>> 0 --- Execution time = 1m07s ---

-This program is even faster than the M-code routine listed in §2b for "large" primes or when the smallest non-trivial divisor is large.
-On the other hand, it may be much slower if n has a "small" divisor.

-Other similar tests can be used:

     • If n < 1373653 is 2- and 3-strong probable prime, then n is prime.
     • If n < 25326001 is 2-, 3- and 5-strong probable prime, then n is prime.
     • If n < 118670087467 is 2-, 3-, 5- and 7-strong probable prime, and if n # 3215031751 then n is prime

-In fact, only 8 integers < 10¹⁰ are simultaneously 2-, 3- and 5-strong pseudoprimes ( cf http://www.research.att.com/~njas/sequences/A056915 )

namely: 25326001 161304001 960946321 1157839381 3215031751 3697278427 5764643587 6770862367

-So we can avoid calling "SPP?" for a = 7 if we check that n is different from these 8 composite numbers.

Data Registers: R00 thru R05: temp ( R04 = n )
Flags: /
Subroutine: "SPP?"

01 LBL "PR3?"
02 STO 04
03 25326001
04 -
05 X=0?
06 RTN
07 RCL 04
08 161304001
09 -
10 X=0?
11 RTN
12 RCL 04
13 960946321

14 -
15 X=0?
16 RTN
17 RCL 04
18 1157839381
19 -
20 X=0?
21 RTN
22 RCL 04
23 3215031751
24 -
25 X=0?
26 RTN

27 RCL 04
28 3697278427
29 -
30 X=0?
31 RTN
32 RCL 04
33 5764643587
34 -
35 X=0?
36 RTN
37 RCL 04
38 6770862367
39 -

40 X=0?
41 RTN
42 RCL 04
43 2
44 XEQ "SPP?"
45 X=0?
46 RTN
47 RCL 04
48 3
49 X<Y?
50 X=0?
51 GTO 00
52 XEQ "SPP?"

53 X=0?
54 RTN
55 1373653
56 RCL 04
57 X<Y?
58 GTO 00
59 5
60 XEQ "SPP?"
61 LBL 00
62 X#0?
63 SIGN
64 END

( 164 bytes / SIZE 006 )

STACK	INPUT	OUTPUT
X	n >= 0	1 or 0

X-output = 1 if n is a prime or if n = 1
X-output = 0 otherwise

Example:

9999999967 XEQ "PR3?" >>>> 1 --- Execution time = 4m09s ---

e) Next Prime

-Given an integer n < 9999999967 , "NPRM" returns the smallest prime p > n

Data Registers: R00 to R02: temp R03 = n , .... , p
Flags: /
Subroutine: "PR?" or "PR2?" or "PR3?"

-Replace R03 by R06 if you want to use "PR2?" or "PR3?" instead of "PR?"

01 LBL "NPRM"
02 STO 03
03 2
04 X>Y?
05 RTN
06 MOD
07 X=0?
08 DSE 03
09 LBL 01
10 RCL 03
11 2
12 +
13 STO 03
14 XEQ "PR?"
15 X=0?
16 GTO 01
17 RCL 03
18 END

( 33 bytes / SIZE 004 )

STACK	INPUT	OUTPUT
X	n	p

Examples:

12346 XEQ "NPRM" >>>> p = 12347
12347 R/S >>>> p = 12373

-If you have an M-code routine like "PR?" ( listed in §2-b) ),

   replace line 02 by STO Y
   replace line 08 by DSE Y
   replace line 10 by CLX
   replace line 13 by ENTER^
   replace line 14 by PR?
   replace line 17 by X<>Y

( 31 bytes / SIZE 000 )

3°) Divisor Functions

a) Program#1

-After factorizing n > 1 with "PMF" or "PMF2", "SKN" can be used to compute s_k(n) = Sum_{d | n} d^k

Formula: if n = p₁^r1 ....... p_1m^rm , s_k(n) = Product_i [ p_i^k(ri+1) - 1 ] / ( p_i^k - 1 ) = Product_i ( 1 + p_i^k + p_i^2k + ...... + p_i^k.ri )

Data Registers: R00: temp R01 = k R02 = s_k(n)
Flags: /
Subroutines: /

01 LBL "SKN"
02 STO 01
03 X<>Y
04 STO 00
05 1
06 STO 02

07 LBL 01
08 RCL IND 00
09 FRC
10 E3
11 *
12 RCL IND 00

13 INT
14 RCL 01
15 Y^X
16 SIGN
17 ENTER^
18 LBL 02

19 LASTX
20 *
21 ST+ Y
22 DSE Z
23 GTO 02
24 X<>Y

25 ST* 02
26 ISG 00
27 GTO 01
28 RCL 02
29 END

( 46 bytes )

STACK	INPUTS	OUTPUTS
Y	bbb.eee	/
X	k	s_k(n)

where bbb.eee is the control number of the factorization bbb > 002

Example:

n = 3238704 "PMF" gives 7.010

7.010 ENTER^
1 XEQ "SKN" >>>> s₁(n) = 11577384 --- Execution time = 6s ---

7.010 ENTER^
2 R/S >>>> s₂(n) = 1.610126288 10¹³

Note:

s_k(1) = 1 for all k

b) Program#2

-If n is small, we can compute s_k(n) directly:

Data Registers: R00: temp R01 = n R02 = k R03 = s_k(n)
Flags: /
Subroutines: /

01 LBL "SKN2"
02 STO 02
03 CLX
04 STO03
05 X<>Y
06 STO 01
07 SQRT

08 INT
09 STO 00
10 LBL 01
11 RCL 01
12 RCL 00
13 MOD
14 X#0?

15 GTO 02
16 RCL 01
17 LASTX
18 /
19 STO Y
20 RCL 02
21 Y^X

22 ST+ 03
23 CLX
24 RCL 00
25 X=Y?
26 GTO 02
27 LASTX
28 Y^X

29 ST+ 03
30 LBL 02
31 DSE 00
32 GTO 01
33 RCL 03
34 END

( 50 bytes / SIZE 004 )

STACK	INPUTS	OUTPUTS
Y	n > 0	/
X	k	s_k(n)

Examples:

25 ENTER^
4 XEQ "SKN2" >>>> s₄(25) = 391251 ( in 3 seconds )

1000 ENTER^
0 R/S >>>> s₀(1000) = 16 ( in 13 seconds )

999999 ENTER^
7 R/S >>>> s₇(999999) ~ 1.000451735 E42 ( in 4mn22s )

c) Program#3 ( multiprecision , n < 100000 )

-"SKN3" computes s_k(n) exactly and stores the result ( by groups of 5 digits ) in registers R09 R10 ...... from the right to the left.
-Unlike the program listed in paragraph 1, k must be an integer greater than 1 and n cannot exceed 100,000

Data Registers: R00 = bbb.eee = control number of the result , R01 = n , R02 = k , R03 thru R08: temp

Rbb to Ree = the digits of s_k(n) Ree+1 to Ree+1+ee-bb: temp
Flags: /
Subroutines: /

01 LBL "SKN3"
02 STO 01
03 X<>Y
04 STO 02
05 Y^X
06 LOG
07 E5
08 STO 08
09 LOG
10 /
11 INT
12 .1
13 %
14 9.009
15 +
16 STO 00
17 CLRGX
18 X<>Y
19 1
20 +
21 1.001

22 *
23 +
24 STO 03
25 RCL 01
26 SQRT
27 INT
28 STO 04
29 LBL 01
30 RCL 01
31 RCL 04
32 /
33 FRC
34 X#0?
35 GTO 04
36 LASTX
37 STO 07
38 XEQ 10
39 RCL 07
40 RCL 04
41 STO 07
42 X#Y?

43 XEQ 10
44 LBL 04
45 DSE 04
46 GTO 01
47 RCL 00
48 RTN
49 LBL 10
50 RCL 02
51 STO 05
52 RCL 03
53 CLRGX
54 RCL 07
55 STO IND Y
56 DSE 05
57 LBL 02
58 ST* IND Y
59 ISG Y
60 GTO 02
61 RCL 03
62 0
63 LBL 09

64 RCL IND Y
65 +
66 RCL X
67 RCL 08
68 MOD
69 STO IND Z
70 -
71 RCL 08
72 /
73 ISG Y
74 GTO 09
75 RCL 03
76 RCL 07
77 DSE 05
78 GTO 02
79 RCL 00
80 STO 05
81 RCL 03
82 STO 06
83 CLX
84 LBL 03

85 RCL IND 05
86 +
87 RCL IND 06
88 +
89 STO Y
90 RCL 08
91 MOD
92 STO IND 05
93 -
94 RCL 08
95 /
96 ISG 05
97 ""
98 ISG 06
99 GTO 03
100 LASTX
101 *
102 DSE 05
103 ""
104 ST+ IND 05
105 END

( 154 bytes )

STACK	INPUTS	OUTPUTS
Y	k	/
X	n	bbb.eee

where bbb.eee = control number of the result ( in registers Rbb to Ree )

Example:

7 ENTER^
99999 XEQ "SKN3" >>>> 9.015 ( execution time = 7mn03s )

we find R09 = 61408 , R10 = 13064 , R11 = 95537 , R12 = 84451 , R13 = 30096 , R14 = 74265 , R15 = 100038

whence s₇(99999) = 100038,74265,30096,84451,95537,13064,61408

Notes:

-The last group of digits ( in register Ree ) may have more than 5 digits.
-If you want to reverse the content of registers Rbb thru Ree, in order to get the groups of digits from the left to the right,
add XEQ "RVL" RCL 00 after line 47 where "RVL" is listed in "Miscellaneous Short Routines for the HP-41"

d) List of the Divisors

Data Registers: R00 = eee , R01 = 1st divisor , R02 = 2nd divisor , ........... , Ree = 1
Flags: /
Subroutines: /

01 LBL "LDIV"
02 ENTER^
03 ENTER^
04 STO 01
05 2
06 STO 00
07 /
08 INT
09 LBL 01
10 MOD
11 X#0?
12 GTO 01
13 CLX
14 LASTX
15 STO IND 00
16 ISG 00
17 LBL 01
18 X<> L
19 DSE X
20 GTO 01
21 SIGN
22 ST- 00
23 RCL 00
24 E3
25 /
26 +
27 END

( 44 bytes )

STACK	INPUTS	OUTPUTS
Y	/	n
X	n > 1	1.eee

where 1.eee = control number of the list of all the divisors ( in registers R01 to Ree )

Example:

117 XEQ "LDIV" >>>> 1.006 --- Execution time = 11s ---

and { R01 R02 R03 R04 R05 R06 } = { 117 39 13 9 3 1 }

-Add STO 00 XEQ "RVL" RCL 00 after line 26 if you want to get the divisors in increasing order.
( "RVL" is listed in "Miscellaneous Short Routines for the HP-41" )

4°) Moebius Function

a) Program#1

"PMF" calculates Moebius function but when we have to compute µ(n) only, the program may be stopped as soon as µ(n) = 0 ( if µ(n) = 0 , line 30 )

Data Registers: R00 - R01: temp R02 = µ(n)
Flags: /
Subroutines: /

01 LBL "MOEB"
02 STO 01
03 1
04 X=Y?
05 RTN
06 STO 02
07 X<>Y
08 X=0?
09 RTN
10 2
11 MOD
12 X#0?

13 GTO 00
14 LBL 01
15 STO 00
16 LBL 02
17 ISG 00
18 CLX
19 RCL 01
20 LASTX
21 /
22 STO 01
23 LASTX
24 MOD

25 X=0?
26 GTO 02
27 RCL 02
28 CHS
29 DSE 00
30 CLX
31 STO 02
32 X=0?
33 RTN
34 LBL 00
35 2
36 LASTX

37 X=Y?
38 DSE L
39 LBL 03
40 RCL 01
41 RCL 01
42 LASTX
43 2
44 +
45 MOD
46 X=0?
47 GTO 01
48 X<> L

49 X^2
50 X<Y?
51 GTO 03
52 SIGN
53 X=Y?
54 GTO 00
55 RCL 02
56 CHS
57 STO 02
58 LBL 00
59 RCL 02
60 END

( 78 bytes / SIZE 003 )

STACK	INPUT	OUTPUT
X	n	µ(n)

Examples:

12 XEQ "MOEB" >>>> µ(12) = 0
105 XEQ "MOEB" >>>> µ(105) = -1

-Of course, lines 39 to 51 are not the fastest way to seek divisors
-The following variant uses the M-code routine "PR?" listed in paragraph 2-b

b) Program#2

01 LBL "MOEB2"
02 STO 01
03 1
04 X=Y?
05 RTN
06 STO 02
07 X<>Y
08 X=0?
09 RTN

10 GTO 03
11 LBL 01
12 STO 00
13 LBL 02
14 ISG 00
15 CLX
16 RCL 01
17 LASTX
18 /

19 STO 01
20 LASTX
21 MOD
22 X=0?
23 GTO 02
24 RCL 02
25 CHS
26 DSE 00
27 CLX

28 STO 02
29 X=0?
30 RTN
31 LBL 03
32 RCL 01
33 PR?
34 X=0?
35 GTO 01
36 RCL 01

37 X=Y?
38 GTO 00
39 RCL 02
40 CHS
41 STO 02
42 LBL 00
43 RCL 02
44 END

( 61 bytes / SIZE 003 )

STACK	INPUT	OUTPUT
X	n	µ(n)

Examples:

12 XEQ "MOEB2" >>>> µ(12) = 0
105 XEQ "MOEB2" >>>> µ(105) = -1

-"MOEB" ( or "MOEB2" ) is called as a subroutine by "CYCLO2" to calculate cyclotomic polynomials

5°) Euler Totient Function

-This short routine can be used if n is very small

Data Registers: R00 = phi(n) R01 = n R02 = n-1 , n-2 , .... , 0
Flags: /
Subroutines: /

01 LBL "PHI"
02 STO 01
03 STO 02
04 CLX
05 STO 00
06 DSE 02
07 LBL 01
08 RCL 01
09 RCL 02
10 LBL 02
11 MOD
12 LASTX
13 X<>Y
14 X#0?
15 GTO 02
16 SIGN
17 X=Y?
18 ST+ 00
19 DSE 02
20 GTO 01
21 RCL 00
22 END

( 35 bytes / SIZE 003 )

STACK	INPUT	OUTPUT
X	n	phi(n)

Examples:

1 XEQ "PHI" >>>> phi(1) = 1
360 R/S >>>> phi(360) = 96 --- Execution time = 3mn27s ---

-The programs listed in paragraph 1 are of course much faster.

References:

[1] Abramowitz and Stegun - "Handbook of Mathematical Functions" - Dover Publications - ISBN 0-486-61272-4
[2] Chris Caldwell - http://primes.utm.edu/
[3] Jean-Paul Delahaye - "Merveilleux Nombres Premiers" - Belin - ISBN 2-84245-017-5 ( in French )