Miller-Rabin Primality Test for the HP-48

[Thomas Klemm]

Modular Arithmetic for the HP-48

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

Addition

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 * k²
+ x * v * k
+ y * u * k
+ y * v


Note:
Even though e.g. x * u * k² 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 = 2² 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