Miller-Rabin Primality Test for the HP-48
07-01-2018, 08:19 PM
Post: #1
 Thomas Klemm Senior Member Posts: 1,447 Joined: Dec 2013
Miller-Rabin Primality Test for the HP-48
Modular Arithmetic for the HP-48

We assume that n is the modulus: 0 < n < 1,000,000,000,000

Code:
PLUS ( a b n -- (a + b) % n )

When we add two integers a and b their sum a + b may overflow when it's bigger than 999,999,999,999.
However if we can assume that both values are between 0 and a modulus n we can calculate a - n + b instead.

This leads to the following program:
Code:
« → a b n   « a b +     IF DUP n <     THEN     ELSE DROP a n - b +     END   » »

Example

This allows to calculate:

500,000,000,000 + 500,000,000,001 = 2 (mod 999,999,999,999)

Otherwise the last digit of the sum is cancelled and the result was 1.

Note:
The result is between 0 and n.

Multiplication

Code:
TIMES ( a b n -- a * b % n )

To multiply two numbers we can translate them to the base k = 1,000,000.

a = x * k + y
b = u * k + v

We can then calculate their product:

a * b
= (x * k + y) * (u * k + v)
= x * u * k2
+ x * v * k
+ y * u * k
+ y * v

Note:
Even though e.g. x * u * k2 is bigger than 1'000'000'000'000 the HP-48 calculates MOD n correctly.
We just have to make sure that the product x * u doesn't overflow.

Code:
« 1000000 → a b n k   « a k / IP a k MOD     b k / IP b k MOD     → x y u v     « x u * k SQ * n MOD       x v * k * n MOD n PLUS       y u * k * n MOD n PLUS       y v * n MOD n PLUS     »   » »

Power

Code:
POWER ( a b n -- a ^ b % n )

We will use the binary method for modular exponentiation.

Python
Code:
def power(a, b, n):     r = 1     while b > 0:         if b % 2 == 1:             r = r * a % n         a = a * a % n         b /= 2     return r

Note:
Python provides the function pow(a, b, n) out of the box.

UserRPL
Code:
« → n   « 1 ROT ROT     WHILE DUP     REPEAT → r a b       «         IF b 2 MOD         THEN a r n TIMES         ELSE r         END         a DUP n TIMES         b 2 / IP       »     END     DROP2   » »

Miller-Rabin Primality Test

If you only want to check if a number n is prime and don't need to find the factors there are more efficient algorithms to do so.
One of them is the Miller-Rabin Primality Test.

Here we use the deterministic variant.

Composite

Check with witness a if n is composite.
If the result is 1 we can be sure that n is composite.
If the result is 0 the number may still be composite.
In this case a is called a strong liar for n.

Example:

n = 221
n - 1 = 220 = 22 55

Code:
174 2 55 221 COMPOSITE? → 0

Thus 174 is a strong liar for 221.

But:

Code:
137 2 55 221 COMPOSITE? → 1

Thus we now can be sure that 221 = 13 * 17 is indeed composite.

Code:
COMPOSITE? ( a s d n -- true if n is composite )

Python
Code:
def composite(a, s, d, n):     x = pow(a, d, n)     if x == 1:         return False     for i in range(s):         if x == n - 1:             return False         x = x * x % n     return True

UserRPL
Code:
« DUP 1 - → a s d n m   « a d n POWER     IF DUP 1 ==     THEN DROP 0     ELSE       IFERR s         WHILE DUP         REPEAT           SWAP           IF DUP m ==           THEN 0 DOERR           END           DUP n TIMES           SWAP 1 -         END       THEN DROP2 0       ELSE DROP2 1       END     END   » »

Prime?

Code:
PRIME? ( n -- true if n is prime )

According to Deterministic variants:

Quote:if n < 1,122,004,669,633, it is enough to test a = 2, 13, 23, and 1662803;

Python
Code:
def prime(n):     d, s = n - 1, 0     while not d % 2:         d, s = d >> 1, s + 1     return not any(composite(a, s, d, n) for a in (2, 13, 23, 1662803))

UserRPL
Code:
« → n    « 0 n 1 -     WHILE DUP 2 MOD NOT     REPEAT       SWAP 1 +       SWAP 2 /     END → s d     « { 2 13 23 1662803 }       IFERR         WHILE DUP SIZE         REPEAT           IF DUP HEAD s d n COMPOSITE?           THEN 0 DOERR           END           TAIL         END       THEN DROP 0       ELSE DROP 1       END     »   » »

Examples:

Code:
999999999997 PRIM? 0

In fact: 999999999997 = 5507×181587071

Code:
999999999989 PRIM? 1

Code:
177777773 PRIM? 1

It would be interesting to see how long this takes on a real HP-28.
We can use the division by 0 trick to raise an error since DOERR is missing.

Cheers
Thomas
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)