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  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 )