hp41programs

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)  2 Focal Programs
   b)  2 M-Code Routines
   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

 
Latest Update:    §2a) & 2b)
 

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

           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)  2 Focal Programs
 

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

Notes:

-If n is a prime, execution time is about  sqrt(n) / 25 seconds

-After testing n mod 2 , n mod 3 & n mod 5, this program skips the multiples of 2 , 3 & 5.
-We can also skip the multiples of 7 after testing n mod 7
-The routine becomes slightly faster, but much longer:


Data Registers:   R00 unused   R01 = 2  R02 = 4  R03 = 6  R04 = 8  R05 = 10
Flags: /
Subroutines: /
 
-Line 325 is a three-byte GTO 01
-All the GTO 02 before line 244 should be three-byte GTOs, but it's perhaps not really useful...

01 LBL "PR?"  02 STO Y  03 4  04 STO 02  05 6  06 STO 03  07 +  08 STO 05  09 8  10 STO 04  11 -  12 STO 01  13 MOD  14 X=0?  15 GTO 02  16 CLX  17 3  18 MOD  19 X=0?  20 GTO 02  21 CLX  22 5  23 MOD  24 X=0?  25 GTO 02  26 CLX  27 7  28 MOD  29 X=0?  30 GTO 02  31 SIGN  32 SIGN  33 LBL 01  34 X<> L  35 RCL 05            36 +  37 MOD  38 X=0?  39 GTO 02  40 X<> L  41 RCL 01  42 +  43 MOD  44 X=0?  45 GTO 02  46 X<> L  47 RCL 02  48 +

49 MOD  50 X=0?  51 GTO 02  52 X<> L  53 RCL 01  54 +  55 MOD  56 X=0?  57 GTO 02  58 X<> L  59 RCL 02  60 +  61 MOD  62 X=0?  63 GTO 02  64 X<> L  65 RCL 03  66 +  67 MOD  68 X=0?  69 GTO 02  70 X<> L  71 RCL 01  72 +  73 MOD  74 X=0?  75 GTO 02  76 X<> L  77 RCL 03  78 +  79 MOD  80 X=0?  81 GTO 02  82 X<> L  83 RCL 02            84 +  85 MOD  86 X=0?  87 GTO 02  88 X<> L  89 RCL 01  90 +  91 MOD  92 X=0?  93 GTO 02  94 X<> L  95 RCL 02  96 +

97 MOD  98 X=0?  99 GTO 02 100 X<> L 101 RCL 03 102 + 103 MOD 104 X=0? 105 GTO 02 106 X<> L 107 RCL 03 108 + 109 MOD 110 X=0? 111 GTO 02 112 X<> L 113 RCL 01 114 + 115 MOD 116 X=0? 117 GTO 02 118 X<> L 119 RCL 03 120 + 121 MOD 122 X=0? 123 GTO 02 124 X<> L 125 RCL 02 126 + 127 MOD 128 X=0? 129 GTO 02 130 X<> L 131 RCL 01           132 + 133 MOD 134 X=0? 135 GTO 02 136 X<> L 137 RCL 03 138 + 139 MOD 140 X=0? 141 GTO 02 142 X<> L 143 RCL 02 144 +

145 MOD 146 X=0? 147 GTO 02 148 X<> L 149 RCL 03 150 + 151 MOD 152 X=0? 153 GTO 02 154 X<> L 155 RCL 04 156 + 157 MOD 158 X=0? 159 GTO 02 160 X<> L 161 RCL 02 162 + 163 MOD 164 X=0? 165 GTO 02 166 X<> L 167 RCL 01 168 + 169 MOD 170 X=0? 171 GTO 02 172 X<> L 173 RCL 02 174 + 175 MOD 176 X=0? 177 GTO 02 178 X<> L 179 RCL 01           180 + 181 MOD 182 X=0? 183 GTO 02 184 X<> L 185 RCL 02 186 + 187 MOD 188 X=0? 189 GTO 02 190 X<> L 191 RCL 04 192 +

193 MOD 194 X=0? 195 GTO 02 196 X<> L 197 RCL 03 198 + 199 MOD 200 X=0? 201 GTO 02 202 X<> L 203 RCL 02 204 + 205 MOD 206 X=0? 207 GTO 02 208 X<> L 209 RCL 03 210 + 211 MOD 212 X=0? 213 GTO 02 214 X<> L 215 RCL 01 216 + 217 MOD 218 X=0? 219 GTO 02 220 X<> L 221 RCL 02 222 + 223 MOD 224 X=0? 225 GTO 02 226 X<> L 227 RCL 03           228 + 229 MOD 230 X=0? 231 GTO 02 232 X<> L 233 RCL 01 234 + 235 MOD 236 X=0? 237 GTO 02 238 X<> L 239 RCL 03 240 +

241 MOD 242 X=0? 243 GTO 02 244 X<> L 245 RCL 03 246 + 247 MOD 248 X=0? 249 GTO 02 250 X<> L 251 RCL 02 252 + 253 MOD 254 X=0? 255 GTO 02 256 X<> L 257 RCL 01 258 + 259 MOD 260 X=0? 261 GTO 02 262 X<> L 263 RCL 02 264 + 265 MOD 266 X=0? 267 GTO 02 268 X<> L 269 RCL 03 270 + 271 MOD 272 X=0? 273 GTO 02 274 X<> L 275 RCL 01           276 + 277 MOD 278 X=0? 279 GTO 02 280 X<> L 281 RCL 03 282 + 283 MOD 284 X=0? 285 GTO 02 286 X<> L 287 RCL 02 288 +

289 MOD 290 X=0? 291 GTO 02 292 X<> L 293 RCL 01 294 + 295 MOD 296 X=0? 297 GTO 02 298 X<> L 299 RCL 02 300 + 301 MOD 302 X=0? 303 GTO 02 304 X<> L 305 RCL 01 306 + 307 MOD 308 X=0? 309 GTO 02 310 X<> L 311 RCL 05 312 + 313 MOD 314 X=0? 315 GTO 02 316 X<> L 317 RCL 01           318 + 319 MOD 320 X=0? 321 GTO 02 322 X<> L 323 X^2 324 X<Y? 325 GTO 01 326 RDN 327 SIGN 328 LBL 02 329 CLX 330 LASTX 331 X<Y? 332 CLX 333 X=Y? 334 SIGN 335 END

         ( 447 bytes / SIZE 006 )

 

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

Example:

   9999999967   XEQ "PR?"   >>>>   1    thus   9999999967  is a prime                ---Execution time = 52 mn---


Note:

-If n is a prime, execution time is about  sqrt(n) / 32 seconds

     b) 2 M-Code Routines
 

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

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  FA84  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 = n   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. 

-Here is another M-Code routine - a little slower - but independant from the address in your ROM:

0BF   "?"
012   "R"
010   "P"
0F8   C=X
361   ?NCXQ          checks for alpha data
050   14D8              and execute SETDEC
05E   C= | C |
088   C
0ED  =
064   INT(C)
0E8  X=C
2F9  C=
060  SQRT(C)
088   C
0ED  =
064   INT(C)
1A8  N=C
244   CLRF 9
0F8  C=X
10E  A=C ALL
04E  C
35C  =
090  2
128  L=C
34E  ?A=0
3A0  ?NC RTN
044   C
070   =
171  
064  AmodC
2EE  ?C#0
0C3  JNC+24d
0F8  C=X
10E  A=C ALL
04E  C
35C  =
0D0  3
128  L=C
044   C
070   =
171  
064  AmodC
2EE  ?C#0
063  JNC+12d
0F8  C=X
10E  A=C ALL
04E  C
35C  =
150  5
128  L=C
044  C
070   =
171  
064  AmodC
2EE  ?C#0
13B  JNC+39d
226  C=C+1 S&X
0EE  B<>C  ALL
1B8  C=N
070  N=C ALL
0F8  C=X
006  A=0 S&X
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
0F0  C<>N ALL
066  A<>B S&X
1A6  A=A-1 S&X
01F  JC+03
066  A<>B S&X
39B  JNC-13d
0E8  X=C
0EE  B<>C ALL
04E  C=0
01C  PT=3
050  LD@PT-1
128  L=C
0DC  PT=10
090   C
190
110
090
110  =
090
110
190
130
007  26424246 E7
228  P=C
013  JNC+02
1EB  JNC+61d
238   C=P                 LOOP1
168   M=C                LOOP2
00E  A=0 ALL
01C  PT=3
0A2  A<>C @PT
138  C=L
21A  C=A+C M
128  L=C
10E  A=C ALL
0FC  RCR 10
0AE  A>C ALL
3EE  LSHFA  ALL
32E  ?A<B ALL
3F7  JC-02
38E  RSHFA  ALL
07A  A<>B M
19A  A=A-B M
3FB  JNC-01
13A  A=A+B M
3BA  RSHFB  M
31A  ?A<C M
3DB  JNC-05
0F8  C=X
0EE  B<>C ALL
35A  ?A#0 M
05B  JNC+11d
178  C=M
3DA  RSHFC M
266  C=C-1  S&X
323  JNC-28d              GTO LOOP2
138  C=L
10E  A=C ALL
0B0  C=N ALL
31A  ?A<C M
2F7  JNC-34d             GTO LOOP1
248  SETF 9
130  LDI S&X
001  001
0EE  B<>C ALL
070  N=C ALL
138  C=L
10E  A=C ALL
130  LDI S&X
009  009
35C  PT=12
342  ?A#0 @PT
027  JC+04
266  C=C-1 S&X
3FA  LSHFA  M
3E3  JNC-04
0BA  A<>C M
0F0  C<>N ALL
066  A<>B S&X
1A6  A=A-1 S&X
01F  JC+03
066 A<>B S&X
38B  JNC-15d
128  L=C
0B0  C<>N ALL
0E8  X=C
24C  ?FSET 9
04F  JC+09
1B8  C=N
10E  A=C ALL
138  C=L
2BE  C=-C
01D  C=
061  A+C
2FE  ?C<0
03B  JNC+07
0F8  C=X
128  L=C
04E  C
35C  =
050  1
013  JNC+02
04E  C=0 ALL
0E8  X=C
345  ?NCGO
042  CLA

( 175 words )

 

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

  X-output = 1   if n is a prime or if n = 1       L-output = n   in this case
  X-output = 0   otherwise                             L-output = divisor d

Examples:

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

Notes:

-PR? first takes the absolute value and the integer part of X-input.
-Alpha "register" is cleared  ( synthetic registers  M N O P )
-Registers Y Z T are 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	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"  >>>>  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	sk(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 )