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 sk(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 = sk(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 | sk(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 )
s1(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
s7(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
s0(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
N1/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.
-s0(n) is the number of divisors of n.
-s1(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 sk(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 =
35 74 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 = 24 35 72 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:
sk(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 = 24 35 72 17
d) Pollard rho factorization
-This method is a probabilistic algorithm to find one factor of an integer n:
• Let x0 = y0 = a random integer between 0 and n-2 and f(x) = x2 + c ( c # 0 and # -2 , c = 1 in the routine below )
• Compute xk = f(xk-1)
, yk = f [f(yk-1)]
and d = gcd ( | xk - yk |
, 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 = xk R03 = yk
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
< 1010 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---- |
0 | 1 | 2 | 3 | 4 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3 |
After executing this subroutine, CPU register C contains:
MS| --------------- Mantissa -------------- | --S&X---- |
0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 2 | 3 | 4 | 0 | 0 | 0 |
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 = 1010 or 1+1010
...
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 2r where q is odd
• Let a > 1 , n is called "strong probable prime to base a" or "a-strong probable prime"
if
aq = 1 mod n
or aq.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
= 999912 )
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 < 1012 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 < 1010 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 sk(n) = Sumd | n d k
Formula: if n = p1r1
....... p1mrm , sk(n) = Producti
[ pik(ri+1) - 1 ] / ( pik -
1 ) = Producti ( 1 + pik + pi2k
+ ...... + pik.ri )
Data Registers: R00: temp
R01 = k R02 = sk(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 | sk(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"
>>>> s1(n) = 11577384
--- Execution time = 6s ---
7.010 ENTER^
2
R/S >>>>
s2(n) = 1.610126288 1013
Note:
sk(1) = 1 for all k
b) Program#2
-If n is small, we can compute sk(n) directly:
Data Registers: R00: temp
R01 = n R02 = k R03 = sk(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 | sk(n) |
Examples:
25 ENTER^
4 XEQ "SKN2"
>>>> s4(25) = 391251
( in 3 seconds )
1000 ENTER^
0
R/S >>>>
s0(1000) = 16
( in 13 seconds )
999999 ENTER^
7
R/S >>>> s7(999999)
~ 1.000451735 E42 ( in 4mn22s )
c) Program#3 (
multiprecision , n < 100000 )
-"SKN3" computes sk(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 sk(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 s7(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 )