(48G) Miller-Rabin Primality Test

[Thomas Klemm]

Code:

@PRIME? (n -- 1 if n is prime else 0 )
« DUP 1 -
  @TIMES ( a b n -- a * b % n )
  « 1000000 → a b n k
    « a k MOD a OVER -
      b k MOD b OVER -
      → y x v u
      « x u * n NEG MOD
        x v * n MOD
        + n NEG MOD
        y u * n MOD
        + n NEG MOD
        y v *
        + n MOD
      »
    »
  »
  → n m TIMES
  « 0 m
    DO
      SWAP 1 +
      SWAP 2 /
    UNTIL
      DUP 2 MOD
    END → s d
    « { 2 13 23 1662803 }
      IF DUP n POS @ one of the primes in the list
      THEN DROP 1
      ELSE
        IFERR
          WHILE DUP SIZE
          REPEAT
            @POWER ( a -- a ^ d (mod n) )
            1 OVER HEAD d
            WHILE DUP
            REPEAT → r a b
              « IF b 2 MOD
                THEN a r n TIMES EVAL @ r → a * r (mod n)
                ELSE r
                END
                a DUP n TIMES EVAL @ a → a ^ 2 (mod n)
                b 2 / IP
              »
            END
            DROP2
            @ x = a ^ d (mod n)
            IF DUP 1 ==
            THEN DROP 0 @ x == 1 → prime
            ELSE
              IFERR s
                WHILE DUP
                REPEAT
                  SWAP
                  IF DUP m ==
                  THEN 0 DOERR
                  END
                  DUP n TIMES EVAL @ x → x ^ 2 (mod n)
                  SWAP 1 -
                END
              THEN DROP2 0 @ x == n - 1 → prime
              ELSE DROP2 1 @ else → composite
              END
            END
            IF @ composite?
            THEN 0 DOERR
            END
            TAIL
          END
        THEN DROP 0
        ELSE DROP 1
        END
      END
    »
  »
»

Benchmark using HP-48GX
999,999,999,999 → 0 in ~ 8 seconds.
999,999,999,989 → 1 in ~ 31 seconds.
177,777,773 → 1 in ~ 21 seconds.